branchandbound_tools.py 18.7 KB
Newer Older
Flavio Baccari's avatar
Flavio Baccari committed
1 2 3 4 5 6 7 8 9 10
#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Created on Wed Nov 22 23:50:46 2017

@author: Flavio
"""
from __future__ import print_function, division
import numpy as np
import copy as cp
Peter Wittek's avatar
Peter Wittek committed
11
from numpy import dot
Flavio Baccari's avatar
Flavio Baccari committed
12
from operator import itemgetter
Peter Wittek's avatar
Peter Wittek committed
13 14 15
from time import time
from ncpol2sdpa import flatten, SdpRelaxation, get_monomials
from ncpol2sdpa.nc_utils import get_support
Flavio Baccari's avatar
Flavio Baccari committed
16 17


18
def get_triang_ineqs(spins):
Peter Wittek's avatar
Peter Wittek committed
19 20 21
    """Generate the list of triangle inequalities involving the given spin
    variables

22 23 24 25
    :param spins: list of spin variables
    :type spins: list of sympy.core.symbol.Symbol
    :returns: list of sympy.core.add.Add
    """
Peter Wittek's avatar
Peter Wittek committed
26

Flavio Baccari's avatar
Flavio Baccari committed
27
    ineq = []
Peter Wittek's avatar
Peter Wittek committed
28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45

    # Run over all the triples
    for i in range(len(spins) - 2):
        for j in range(i + 1, len(spins) - 1):
            for k in range(j + 1, len(spins)):

                # Add all the different sign choices
                ineq.append(spins[i] * spins[j] + spins[i] * spins[k] +
                            spins[j] * spins[k] + 1)
                ineq.append(spins[i] * spins[j] - spins[i] * spins[k] -
                            spins[j] * spins[k] + 1)
                ineq.append(-spins[i] * spins[j] - spins[i] * spins[k] +
                            spins[j] * spins[k] + 1)
                ineq.append(-spins[i] * spins[j] + spins[i] * spins[k] -
                            spins[j] * spins[k] + 1)
    return ineq


46 47 48 49 50 51 52 53 54 55 56 57 58 59 60
def get_det_guess(s_variables, cliques, Mom, one, two):
    """Obtains the closest spin configuration from the moment matrix resulting from 
    the lower bound SDP
    
    :param s_variables: list of spin variables
    :type s_variables: list of sympy.core.symbol.Symbol
    :param cliques: list of the cliques of the chordal extension
    :type cliques: list of lists of sympy.core.symbol.Symbol
    :param Mom: the block moment matrix at level 1
    :type Mom: list of numpy.ndarray
    :param one: coefficients of the one-body terms in the hamiltonian 
    :type one: list of float
    :param two: coefficients of the two-body terms in the hamiltonian as a matrix
    :type one: list of list of float
    
Flavio Baccari's avatar
Flavio Baccari committed
61

62 63
    :returns: the deterministic configuration final_config as list of floats
    """
Peter Wittek's avatar
Peter Wittek committed
64

65
    #taking a Choleski decomposition of each block in the moment matrix
Flavio Baccari's avatar
Flavio Baccari committed
66 67 68
    Red = []
    configs = []

Peter Wittek's avatar
Peter Wittek committed
69 70 71
    for c, M in enumerate(Mom):

        B = np.linalg.cholesky(M + (10**(-6)) * np.identity(len(M)))
Flavio Baccari's avatar
Flavio Baccari committed
72
        BT = B.T.conj()
Peter Wittek's avatar
Peter Wittek committed
73 74 75
        # Get the smallest size of the vectors

        diff = np.linalg.norm(M - np.dot(B, BT))
76
        l = len(flatten(s_variables))
Flavio Baccari's avatar
Flavio Baccari committed
77 78

        while diff < 10**(-5):
Peter Wittek's avatar
Peter Wittek committed
79

Flavio Baccari's avatar
Flavio Baccari committed
80
            l -= 1
Peter Wittek's avatar
Peter Wittek committed
81 82 83 84 85
            tmp = B[::, 0:l]
            diff = np.linalg.norm(M - np.dot(tmp, tmp.T.conj()))

        d = l + 1
        Red.append(B[::, 0:d])
Flavio Baccari's avatar
Flavio Baccari committed
86

Peter Wittek's avatar
Peter Wittek committed
87
        configs.append([np.sign(dot(B[0, 0:d], el)) for el in B[1::, 0:d]])
Flavio Baccari's avatar
Flavio Baccari committed
88

Peter Wittek's avatar
Peter Wittek committed
89 90 91 92 93
    # Extract the overall configuration
    sub_conf = {}

    for b, block in enumerate(cliques):
        for e, el in enumerate(block):
Flavio Baccari's avatar
Flavio Baccari committed
94 95 96 97 98

            sub_conf[el] = configs[b][e]

    config = []

99
    for el in flatten(s_variables):
Peter Wittek's avatar
Peter Wittek committed
100

Flavio Baccari's avatar
Flavio Baccari committed
101
        config.append(sub_conf[el])
Peter Wittek's avatar
Peter Wittek committed
102 103

    bound = get_bound_from_det(flatten(config), one, two)
Flavio Baccari's avatar
Flavio Baccari committed
104 105

    final_config = cp.copy(config)
Peter Wittek's avatar
Peter Wittek committed
106
    # Check if the bound can be improved by flipping a single spin
Flavio Baccari's avatar
Flavio Baccari committed
107 108

    for i in range(len(config)):
Peter Wittek's avatar
Peter Wittek committed
109

Flavio Baccari's avatar
Flavio Baccari committed
110 111 112
        flip_config = cp.deepcopy(config)
        flip_config[i] = cp.deepcopy(-flip_config[i])

Peter Wittek's avatar
Peter Wittek committed
113
        bound_tmp = get_bound_from_det(flatten(flip_config), one, two)
Flavio Baccari's avatar
Flavio Baccari committed
114 115

        if bound_tmp < bound:
Peter Wittek's avatar
Peter Wittek committed
116

Flavio Baccari's avatar
Flavio Baccari committed
117
            bound = cp.deepcopy(bound_tmp)
Peter Wittek's avatar
Peter Wittek committed
118 119
            final_config = cp.deepcopy(flip_config)

120
    return final_config
Peter Wittek's avatar
Peter Wittek committed
121

Flavio Baccari's avatar
Flavio Baccari committed
122

Peter Wittek's avatar
Peter Wittek committed
123
def get_bound_from_det(config, one, two):
124 125 126 127 128 129 130 131 132 133 134 135
    """Gets the energy of the corresponding configuration
    
    :param config: the spin configuration
    :type config: list of float
    :param one: coefficients of the one-body terms in the hamiltonian 
    :type one: list of float
    :param two: coefficients of the two-body terms in the hamiltonian as a matrix
    :type one: list of list of float
    

    :returns: the upper bound to the ground state energy as float
    """
Peter Wittek's avatar
Peter Wittek committed
136
    # Config has to be flattened
Flavio Baccari's avatar
Flavio Baccari committed
137 138

    bound = 0
Peter Wittek's avatar
Peter Wittek committed
139 140 141 142 143 144 145 146 147 148 149 150 151

    # One-body
    for i, el in enumerate(config):

        bound += one[i] * el

    # Two-body

    for i, el1 in enumerate(config):
        for j, el2 in enumerate(config[i + 1::], i + 1):

            bound += two[i][j] * el1 * el2

Flavio Baccari's avatar
Flavio Baccari committed
152 153 154
    return bound


155
def get_strengthen_low(s_variables, substitutions, hamiltonian, cliques,
Peter Wittek's avatar
Peter Wittek committed
156
                       threshold, solverparameters, verbose):
157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178
    """Runs the relaxation for the first time and determines how many triangle inequalities
    to add to get a reasonably good lower bound. Outputs the corresponding lower bound 
    as well

    :param s_variables: list of spin variables
    :type s_variables: list of sympy.core.symbol.Symbol
    :param substitutions: substitution to be applied in the generation of the moment matrix at the level of operators
    :type substitutions: dict of items of polynomials of sympy.core.symbol.Symbol
    :param hamiltonian: hamiltonian of the problem as a symbolic polynomial of the spins
    :type hamiltonian: sympy.core.add.Add
    :param cliques: list of the cliques of the chordal extension
    :type cliques: list of lists of sympy.core.symbol.Symbol
    :param threshold: threshold clique size n_t under which generats blocks at level 2
    :type threshold: int
    :param solverparameters: parameters for the SDP solver
    :type solverparameters: dict
    :param verbose: verbosity level
    :type verbose: int

    :returns: sdp as ncpol2sdpa sdp class, list of triangle inequalities to add as 
    list of sympy.core.add.Add
    """
179
    N = len(flatten(s_variables))
Flavio Baccari's avatar
Flavio Baccari committed
180 181 182
    n_ineqs = N
    new_ineqs = []
    ineqs = []
Peter Wittek's avatar
Peter Wittek committed
183 184
    # Generate the mixed level with blocks al level 2 only if the are below
    # the threshold size
Flavio Baccari's avatar
Flavio Baccari committed
185
    monomial_blocks = []
Peter Wittek's avatar
Peter Wittek committed
186

Flavio Baccari's avatar
Flavio Baccari committed
187
    for el in cliques:
Peter Wittek's avatar
Peter Wittek committed
188

Flavio Baccari's avatar
Flavio Baccari committed
189
        if len(el) <= threshold:
Peter Wittek's avatar
Peter Wittek committed
190 191

            monomial_blocks.append(get_monomials(el, 2))
Flavio Baccari's avatar
Flavio Baccari committed
192
        else:
Peter Wittek's avatar
Peter Wittek committed
193 194 195 196
            monomial_blocks.append(get_monomials(el, 1))

    # Solve the first SDP

197
    sdp = SdpRelaxation(flatten(s_variables), verbose=verbose)
Peter Wittek's avatar
Peter Wittek committed
198 199 200 201 202 203 204
    sdp.get_relaxation(
        -1,
        substitutions=substitutions,
        objective=hamiltonian,
        momentinequalities=flatten(ineqs),
        extramonomials=monomial_blocks)

Flavio Baccari's avatar
Flavio Baccari committed
205
    sdp.solve(solver="mosek", solverparameters=solverparameters)
Peter Wittek's avatar
Peter Wittek committed
206

Flavio Baccari's avatar
Flavio Baccari committed
207 208
    bound_init = cp.deepcopy(sdp.primal)
    bound_new = cp.deepcopy(sdp.primal)
209

Peter Wittek's avatar
Peter Wittek committed
210 211 212 213
    # Get the moment matrix without the identity and computing the triangular
    # inequalities

    lengths = [len(el) for el in cliques]
Flavio Baccari's avatar
Flavio Baccari committed
214
    improv = 1
Peter Wittek's avatar
Peter Wittek committed
215 216 217

    while (improv > 0.005):

Flavio Baccari's avatar
Flavio Baccari committed
218 219
        bound_old = cp.deepcopy(sdp.primal)
        ineqs_values = []
Peter Wittek's avatar
Peter Wittek committed
220

Flavio Baccari's avatar
Flavio Baccari committed
221
        for c in range(len(cliques)):
Peter Wittek's avatar
Peter Wittek committed
222 223 224 225 226 227

            M = cp.deepcopy(sdp.x_mat[c][1:lengths[c] + 1, 1:lengths[c] + 1])

            for i in range(lengths[c] - 2):
                for j in range(i + 1, lengths[c] - 1):
                    for k in range(j + 1, lengths[c]):
Flavio Baccari's avatar
Flavio Baccari committed
228
                        tmp = []
Peter Wittek's avatar
Peter Wittek committed
229 230 231 232 233 234 235 236

                        tmp.append(M[i, j] + M[j, k] + M[i, k] + 1)
                        tmp.append(M[i, j] - M[j, k] - M[i, k] + 1)
                        tmp.append(-M[i, j] + M[j, k] - M[i, k] + 1)
                        tmp.append(-M[i, j] - M[j, k] + M[i, k] + 1)

                        ineqs_values.append([[c, i, j, k], min(tmp)])

237
        # Select the first N most violated ineqs and add the corresponding 4
Peter Wittek's avatar
Peter Wittek committed
238 239
        # constraints to the problem

Flavio Baccari's avatar
Flavio Baccari committed
240
        for _ in range(n_ineqs):
Peter Wittek's avatar
Peter Wittek committed
241

Flavio Baccari's avatar
Flavio Baccari committed
242 243
            b = [el[1] for el in ineqs_values]
            index = min(enumerate(b), key=itemgetter(1))[0]
Peter Wittek's avatar
Peter Wittek committed
244
            [c, i, j, k] = ineqs_values[index][0]
Flavio Baccari's avatar
Flavio Baccari committed
245
            ineqs_values.remove(ineqs_values[index])
Peter Wittek's avatar
Peter Wittek committed
246 247 248 249 250 251

            variables = [
                flatten(cliques[c])[i],
                flatten(cliques[c])[j],
                flatten(cliques[c])[k]
            ]
Flavio Baccari's avatar
Flavio Baccari committed
252
            new_ineqs.append(get_triang_ineqs(variables))
Peter Wittek's avatar
Peter Wittek committed
253 254

        # Solve the new sdp
255

Flavio Baccari's avatar
Flavio Baccari committed
256
        time0 = time()
Peter Wittek's avatar
Peter Wittek committed
257

258
        sdp = SdpRelaxation(flatten(s_variables), verbose=verbose)
Peter Wittek's avatar
Peter Wittek committed
259 260 261 262 263 264 265
        sdp.get_relaxation(
            -1,
            substitutions=substitutions,
            objective=hamiltonian,
            momentinequalities=flatten(new_ineqs),
            extramonomials=monomial_blocks)

Flavio Baccari's avatar
Flavio Baccari committed
266
        sdp.solve(solver="mosek", solverparameters=solverparameters)
Peter Wittek's avatar
Peter Wittek committed
267

Flavio Baccari's avatar
Flavio Baccari committed
268
        bound_new = cp.deepcopy(sdp.primal)
Peter Wittek's avatar
Peter Wittek committed
269
        improv = abs(bound_new - bound_old) / abs(bound_init)
Peter Wittek's avatar
Peter Wittek committed
270

271
        #if the relative improvement is good enough add the new inequalities
Peter Wittek's avatar
Peter Wittek committed
272

Flavio Baccari's avatar
Flavio Baccari committed
273
        if (improv > 0.005):
Peter Wittek's avatar
Peter Wittek committed
274

Flavio Baccari's avatar
Flavio Baccari committed
275
            ineqs = cp.deepcopy(new_ineqs)
Peter Wittek's avatar
Peter Wittek committed
276 277 278 279 280

    return sdp, ineqs


def get_Mom(sdp, cliques):
281 282 283 284 285 286 287 288 289 290 291
    """Extract the level 1 moment matrix from the solved sdp problem

    :param spd: SDP problem generated by ncpol2sdpa
    :type sdp: ncpol2sdpa.sdp_relaxation.SdpRelaxation
    :param cliques: list of the cliques of the chordal extension
    :type cliques: list of lists of sympy.core.symbol.Symbol 


    :returns: the level 1 the block moment matrix Mom as list of numpy.ndarray, 
    the average spin values as list of numpy.ndarray 
    """
Flavio Baccari's avatar
Flavio Baccari committed
292

Peter Wittek's avatar
Peter Wittek committed
293
    # Extract the moment matrix for the cliques
Flavio Baccari's avatar
Flavio Baccari committed
294

Peter Wittek's avatar
Peter Wittek committed
295
    lengths = [len(el) for el in cliques]
Flavio Baccari's avatar
Flavio Baccari committed
296 297
    Mom = []
    spins = []
Peter Wittek's avatar
Peter Wittek committed
298

Flavio Baccari's avatar
Flavio Baccari committed
299
    for c in range(len(cliques)):
Peter Wittek's avatar
Peter Wittek committed
300 301 302 303 304 305 306 307 308

        M = cp.copy(sdp.x_mat[c][0:lengths[c] + 1, 0:lengths[c] + 1])
        Mom.append(sdp.x_mat[c][0:lengths[c] + 1, 0:lengths[c] + 1])

        spins.append(M[1:, 0])

    return Mom, spins


309
def get_up_and_low(s_variables, substitutions, ineqs, hamiltonian, one, two,
Peter Wittek's avatar
Peter Wittek committed
310
                   cliques, threshold, solverparameters, verbose):
311 312 313 314 315 316 317 318 319 320 321 322 323 324 325 326 327 328 329 330 331 332
    """Gets upper and lower bound for the given instance of branching

    :param s_variables: list of spin variables
    :type s_variables: list of sympy.core.symbol.Symbol
    :param substitutions: substitution to be applied in the generation of the moment matrix at the level of operators
    :type substitutions: dict of items of polynomials of sympy.core.symbol.Symbol
    :param ineqs: triangle inequalities to be added to improve the lower bound
    :type ineqs: list of sympy.core.add.Add
    :param hamiltonian: hamiltonian of the problem as a symbolic polynomial of the spins
    :type hamiltonian: sympy.core.add.Add
    :param one: coefficients for the one-body terms in the hamiltonian
    :type one: list of floats
    :param two: coefficients for the two-body terms in the hamiltonian
    :type two: list of floats
    :param cliques: list of the cliques of the chordal extension
    :type cliques: list of lists of sympy.core.symbol.Symbol
    :param threshold: threshold clique size n_t under which generats blocks at level 2
    :type threshold: int
    :param solverparameters: parameters for the SDP solver
    :type solverparameters: dict
    :param verbose: verbosity level
    :type verbose: int
Peter Wittek's avatar
Peter Wittek committed
333

334 335 336
    :returns: lower bound as float, upper bound as float, spin expectation values
    extracted from the moment matrix for each clique as list of array of floats,
    the deterministic configuration for each spin as list of floats
Peter Wittek's avatar
Peter Wittek committed
337 338
    """

Flavio Baccari's avatar
Flavio Baccari committed
339 340 341 342 343
    time0 = time()

    monomial_blocks = []

    for el in cliques:
Peter Wittek's avatar
Peter Wittek committed
344

Flavio Baccari's avatar
Flavio Baccari committed
345
        if len(el) <= threshold:
Peter Wittek's avatar
Peter Wittek committed
346 347

            monomial_blocks.append(get_monomials(el, 2))
Flavio Baccari's avatar
Flavio Baccari committed
348
        else:
Peter Wittek's avatar
Peter Wittek committed
349 350
            monomial_blocks.append(get_monomials(el, 1))

Flavio Baccari's avatar
Flavio Baccari committed
351
    #solving the first sdp
Peter Wittek's avatar
Peter Wittek committed
352

353
    sdp = SdpRelaxation(flatten(s_variables), verbose=verbose)
Peter Wittek's avatar
Peter Wittek committed
354 355 356 357 358 359 360
    sdp.get_relaxation(
        -1,
        substitutions=substitutions,
        objective=hamiltonian,
        momentinequalities=flatten(ineqs),
        extramonomials=monomial_blocks)

Flavio Baccari's avatar
Flavio Baccari committed
361
    sdp.solve(solver="mosek", solverparameters=solverparameters)
Peter Wittek's avatar
Peter Wittek committed
362

Flavio Baccari's avatar
Flavio Baccari committed
363
    print("The lower bound is " + str(sdp.primal))
Peter Wittek's avatar
Peter Wittek committed
364 365
    print("Time taken for lower is " + str(time() - time0))

Flavio Baccari's avatar
Flavio Baccari committed
366
    time0 = time()
Peter Wittek's avatar
Peter Wittek committed
367 368
    [Mom, spins] = get_Mom(sdp, cliques)

Flavio Baccari's avatar
Flavio Baccari committed
369
    config = get_det_guess(s_variables, cliques, Mom, one, two)
Peter Wittek's avatar
Peter Wittek committed
370 371 372 373 374 375
    bound = get_bound_from_det(flatten(config), one, two)
    print("Time taken for upper is " + str(time() - time0))

    return sdp.primal, bound, spins, config


Peter Wittek's avatar
Peter Wittek committed
376
def get_branching_node(spins, s_variables, cliques, eqs):
377 378 379 380 381 382 383 384 385 386 387 388 389 390
    """Gets the spin on which to perform the next branching. Works with an "easy-first"
    choice

    :param spins: the spin average values extracted from the moment matrix
    :type spins: list of floats
    :param s_variables: list of spin variables
    :type s_variables: list of sympy.core.symbol.Symbol
    :param cliques: list of the cliques of the chordal extension
    :type cliques: list of lists of sympy.core.symbol.Symbol
    :param eqs: list of already implemented branching choices
    :type eqs: list of sympy.core.add.Add

    :returns: the spin variable of the branching node as sympy.core.symbol.Symbol
    """
Peter Wittek's avatar
Peter Wittek committed
391

Flavio Baccari's avatar
Flavio Baccari committed
392
    det = 1
Peter Wittek's avatar
Peter Wittek committed
393 394 395 396 397
    sup = get_sup(s_variables, eqs)

    for c, block in enumerate(spins):
        for i, el in enumerate(block):

Flavio Baccari's avatar
Flavio Baccari committed
398
            if el not in sup:
Peter Wittek's avatar
Peter Wittek committed
399

Flavio Baccari's avatar
Flavio Baccari committed
400
                if abs(el) < det:
Peter Wittek's avatar
Peter Wittek committed
401

Flavio Baccari's avatar
Flavio Baccari committed
402
                    det = cp.deepcopy(abs(el))
Flavio Baccari's avatar
Flavio Baccari committed
403
                    branch_node = cp.deepcopy(cliques[c][i])
Flavio Baccari's avatar
Flavio Baccari committed
404

Flavio Baccari's avatar
Flavio Baccari committed
405
    return branch_node
Peter Wittek's avatar
Peter Wittek committed
406 407


408 409 410 411 412 413 414 415 416 417 418
def get_as_coeffs(hamiltonian, s_variables):
    """Extracts the one and two-body coefficients from the symbolic hamiltonian
    
    :param hamiltonian: hamiltonian of the problem as a symbolic polynomial of the spins
    :type hamiltonian: sympy.core.add.Add
    :param s_variables: list of spin variables
    :type s_variables: list of sympy.core.symbol.Symbol
    
    :returns: the one-body coeffs as list of float, the two-body coeffs as list 
    of list of float
    """
Peter Wittek's avatar
Peter Wittek committed
419

Flavio Baccari's avatar
Flavio Baccari committed
420
    one = []
Peter Wittek's avatar
Peter Wittek committed
421
    two = []
Flavio Baccari's avatar
Flavio Baccari committed
422
    pol = hamiltonian.as_poly()
423
    mon = flatten(s_variables)
Peter Wittek's avatar
Peter Wittek committed
424 425 426

    for i, el in enumerate(mon):

Flavio Baccari's avatar
Flavio Baccari committed
427 428
        one.append(pol.coeff_monomial(el))
        tmp = []
Peter Wittek's avatar
Peter Wittek committed
429 430
        for j, el2 in enumerate(mon):

Flavio Baccari's avatar
Flavio Baccari committed
431
            if j < i:
Peter Wittek's avatar
Peter Wittek committed
432

Flavio Baccari's avatar
Flavio Baccari committed
433 434
                tmp.append(0)
            else:
Peter Wittek's avatar
Peter Wittek committed
435 436

                tmp.append(pol.coeff_monomial(el * el2))
Flavio Baccari's avatar
Flavio Baccari committed
437 438
        two.append(tmp)

Peter Wittek's avatar
Peter Wittek committed
439 440 441
    return one, two


442 443 444 445 446 447 448 449 450 451
def get_sup(s_variables, eqs):
    """Gets the spin variables on which the branching has alreadu been performed
    
    :param s_variables: list of spin variables
    :type s_variables: list of sympy.core.symbol.Symbol
    :param eqs: list of already implemented branching choices
    :type eqs: list of sympy.core.add.Add
    
    :returns: support of the branching equations as list of sympy.core.symbol.Symbol
    """
Peter Wittek's avatar
Peter Wittek committed
452

Flavio Baccari's avatar
Flavio Baccari committed
453
    sup = []
Peter Wittek's avatar
Peter Wittek committed
454

Flavio Baccari's avatar
Flavio Baccari committed
455
    for el in eqs:
Peter Wittek's avatar
Peter Wittek committed
456

457 458
        s = get_support(flatten(s_variables), el)
        sup.append(dot(s[1], flatten(s_variables)))
Peter Wittek's avatar
Peter Wittek committed
459

Flavio Baccari's avatar
Flavio Baccari committed
460 461 462
    return sup


Peter Wittek's avatar
Peter Wittek committed
463 464 465 466 467 468 469
def get_groundBandB(s_variables,
                    substitutions,
                    hamiltonian,
                    cliques,
                    threshold,
                    solverparameters,
                    verbose=0):
470 471
    """Runs the whole branch&bound algorithm for the given hamiltonian

Flavio Baccari's avatar
Flavio Baccari committed
472 473
    :param s_variables: list of spin variables
    :type s_variables: list of sympy.core.symbol.Symbol
474
    :param substitutions: substitution to be applied in the generation of the moment matrix at the level of operators
475 476 477 478 479 480 481 482 483
    :type substitutions: dict of items of polynomials of sympy.core.symbol.Symbol
    :param hamiltonian: hamiltonian of the problem as a symbolic polynomial of the spins
    :type hamiltonian: sympy.core.add.Add
    :param cliques: list of the cliques of the chordal extension
    :type cliques: list of lists of sympy.core.symbol.Symbol
    :param threshold: threshold clique size n_t under which generats blocks at level 2
    :type threshold: int
    :param solverparameters: parameters for the SDP solver
    :type solverparameters: dict
Flavio Baccari's avatar
Flavio Baccari committed
484 485
    :param verbose: verbosity level
    :type verbose: int
Peter Wittek's avatar
Peter Wittek committed
486

487 488 489
    :returns z_low: final lower bound as float
    :returns z_up: final upper bound as float
    :returns final_config: final ground state configuration as list of floats
490
    """
Peter Wittek's avatar
Peter Wittek committed
491

Flavio Baccari's avatar
Flavio Baccari committed
492
    [one, two] = get_as_coeffs(hamiltonian, s_variables)
Peter Wittek's avatar
Peter Wittek committed
493 494
    # First root

Flavio Baccari's avatar
Flavio Baccari committed
495 496 497
    eqs = []
    tree = []
    print("Starting initial optimisation")
Peter Wittek's avatar
Peter Wittek committed
498 499 500

    # Get the first lower bound and choosing the list of triangle ineqs to add

Peter Wittek's avatar
Peter Wittek committed
501
    [sdp, ineqs] = get_strengthen_low(
Flavio Baccari's avatar
Flavio Baccari committed
502
        s_variables,
Peter Wittek's avatar
Peter Wittek committed
503 504 505 506 507 508
        substitutions,
        hamiltonian,
        cliques,
        threshold,
        solverparameters,
        verbose=0)
Peter Wittek's avatar
Peter Wittek committed
509 510 511

    # Extract the corresponding upper bound and spin configuration

Peter Wittek's avatar
Peter Wittek committed
512
    [Mom, spins] = get_Mom(sdp, cliques)
Flavio Baccari's avatar
Flavio Baccari committed
513
    config = get_det_guess(s_variables, cliques, Mom, one, two)
Peter Wittek's avatar
Peter Wittek committed
514
    final_config = cp.deepcopy(config)
Flavio Baccari's avatar
Flavio Baccari committed
515

Peter Wittek's avatar
Peter Wittek committed
516
    z_up = get_bound_from_det(flatten(config), one, two)
Flavio Baccari's avatar
Flavio Baccari committed
517
    z_low = sdp.primal
Peter Wittek's avatar
Peter Wittek committed
518

Flavio Baccari's avatar
Flavio Baccari committed
519 520
    print("The lower bound is " + str(z_low))
    print("The upper bound is " + str(z_up))
Peter Wittek's avatar
Peter Wittek committed
521

Peter Wittek's avatar
Peter Wittek committed
522 523 524 525
    # Go for the branching with easy first

    if ((z_low - z_up) / abs(z_up) < -0.00001):

Peter Wittek's avatar
Peter Wittek committed
526
        node = get_branching_node(spins, s_variables, cliques, eqs)
Peter Wittek's avatar
Peter Wittek committed
527 528 529 530
        tree = [[z_low, eqs, substitutions, node]]

    while (len(tree) > 0):

Flavio Baccari's avatar
Flavio Baccari committed
531 532
        lis = [el[0] for el in tree]
        index = min(enumerate(lis), key=itemgetter(1))[0]
Peter Wittek's avatar
Peter Wittek committed
533
        [z_low, eqs, subs, node] = tree[index]
Flavio Baccari's avatar
Flavio Baccari committed
534
        tree.remove(tree[index])
Peter Wittek's avatar
Peter Wittek committed
535 536 537

        # Branch +1

Flavio Baccari's avatar
Flavio Baccari committed
538 539 540 541
        eqs_new = cp.deepcopy(eqs)
        eqs_new.append(node - 1)
        subs_new = cp.deepcopy(subs)
        subs_new[node] = 1
Peter Wittek's avatar
Peter Wittek committed
542 543

        [z_low, tmp, spins, config] = get_up_and_low(
Flavio Baccari's avatar
Flavio Baccari committed
544
            s_variables,
Peter Wittek's avatar
Peter Wittek committed
545 546 547 548 549 550 551 552 553 554 555
            subs_new,
            ineqs,
            hamiltonian,
            one,
            two,
            cliques,
            threshold,
            solverparameters,
            verbose=0)

        # Substitute the upper bound if better
Flavio Baccari's avatar
Flavio Baccari committed
556
        if tmp < z_up:
Peter Wittek's avatar
Peter Wittek committed
557

Flavio Baccari's avatar
Flavio Baccari committed
558 559
            z_up = cp.deepcopy(tmp)
            final_config = cp.deepcopy(config)
Peter Wittek's avatar
Peter Wittek committed
560 561 562 563 564

        # Include this node in the branching if the lower bound is meaningful

        if ((z_low - z_up) / abs(z_up) < -0.00001):

Peter Wittek's avatar
Peter Wittek committed
565
            node_new = get_branching_node(spins, s_variables, cliques, eqs_new)
Peter Wittek's avatar
Peter Wittek committed
566 567 568 569
            tree.append([z_low, eqs_new, subs_new, node_new])

        # branch -1

Flavio Baccari's avatar
Flavio Baccari committed
570 571 572 573
        eqs_new = cp.deepcopy(eqs)
        eqs_new.append(node + 1)
        subs_new = cp.deepcopy(subs)
        subs_new[node] = -1
Peter Wittek's avatar
Peter Wittek committed
574 575

        [z_low, tmp, spins, config] = get_up_and_low(
Flavio Baccari's avatar
Flavio Baccari committed
576
            s_variables,
Peter Wittek's avatar
Peter Wittek committed
577 578 579 580 581 582 583 584 585 586 587
            subs_new,
            ineqs,
            hamiltonian,
            one,
            two,
            cliques,
            threshold,
            solverparameters,
            verbose=0)

        # Substitute the upper bound if better
Flavio Baccari's avatar
Flavio Baccari committed
588
        if tmp < z_up:
Peter Wittek's avatar
Peter Wittek committed
589

Flavio Baccari's avatar
Flavio Baccari committed
590 591
            z_up = cp.deepcopy(tmp)
            final_config = cp.deepcopy(config)
Peter Wittek's avatar
Peter Wittek committed
592 593

        # Include this node in the branching if the lower bound is meaningful
Flavio Baccari's avatar
Flavio Baccari committed
594
        print("The upper bound is " + str(z_up))
Peter Wittek's avatar
Peter Wittek committed
595 596
        if ((z_low - z_up) / abs(z_up) < -0.00001):

Peter Wittek's avatar
Peter Wittek committed
597
            node_new = get_branching_node(spins, s_variables, cliques, eqs_new)
Peter Wittek's avatar
Peter Wittek committed
598 599 600 601
            tree.append([z_low, eqs_new, subs_new, node_new])

        # Check all the previous nodes to see if we can eliminate some with
        # the new z_up
Flavio Baccari's avatar
Flavio Baccari committed
602 603 604

        tree_old = cp.deepcopy(tree)
        tree = []
Peter Wittek's avatar
Peter Wittek committed
605

Flavio Baccari's avatar
Flavio Baccari committed
606 607
        for el in tree_old:

Peter Wittek's avatar
Peter Wittek committed
608
            if ((el[0] - z_up) / abs(z_up) < -0.00001):
Flavio Baccari's avatar
Flavio Baccari committed
609

Peter Wittek's avatar
Peter Wittek committed
610
                tree.append(el)
Flavio Baccari's avatar
Flavio Baccari committed
611

Peter Wittek's avatar
Peter Wittek committed
612
    return z_low, z_up, final_config