branchandbound_tools.py 18.7 KB
 Flavio Baccari committed Jun 28, 2018 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 committed Jun 29, 2018 11 ``````from numpy import dot `````` Flavio Baccari committed Jun 28, 2018 12 ``````from operator import itemgetter `````` Peter Wittek committed Jun 29, 2018 13 14 15 ``````from time import time from ncpol2sdpa import flatten, SdpRelaxation, get_monomials from ncpol2sdpa.nc_utils import get_support `````` Flavio Baccari committed Jun 28, 2018 16 17 `````` `````` Flavio Baccari committed Jun 28, 2018 18 ``````def get_triang_ineqs(spins): `````` Peter Wittek committed Jun 29, 2018 19 20 21 `````` """Generate the list of triangle inequalities involving the given spin variables `````` Flavio Baccari committed Jun 28, 2018 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 committed Jun 29, 2018 26 `````` `````` Flavio Baccari committed Jun 28, 2018 27 `````` ineq = [] `````` Peter Wittek committed Jun 29, 2018 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 `````` Flavio Baccari committed Jul 15, 2018 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 committed Jun 28, 2018 61 `````` `````` Flavio Baccari committed Jul 15, 2018 62 63 `````` :returns: the deterministic configuration final_config as list of floats """ `````` Peter Wittek committed Aug 03, 2018 64 `````` `````` Flavio Baccari committed Jul 15, 2018 65 `````` #taking a Choleski decomposition of each block in the moment matrix `````` Flavio Baccari committed Jun 28, 2018 66 67 68 `````` Red = [] configs = [] `````` Peter Wittek committed Jun 29, 2018 69 70 71 `````` for c, M in enumerate(Mom): B = np.linalg.cholesky(M + (10**(-6)) * np.identity(len(M))) `````` Flavio Baccari committed Jun 28, 2018 72 `````` BT = B.T.conj() `````` Peter Wittek committed Jun 29, 2018 73 74 75 `````` # Get the smallest size of the vectors diff = np.linalg.norm(M - np.dot(B, BT)) `````` Flavio Baccari committed Jul 15, 2018 76 `````` l = len(flatten(s_variables)) `````` Flavio Baccari committed Jun 28, 2018 77 78 `````` while diff < 10**(-5): `````` Peter Wittek committed Jun 29, 2018 79 `````` `````` Flavio Baccari committed Jun 28, 2018 80 `````` l -= 1 `````` Peter Wittek committed Jun 29, 2018 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 committed Jun 28, 2018 86 `````` `````` Peter Wittek committed Jun 29, 2018 87 `````` configs.append([np.sign(dot(B[0, 0:d], el)) for el in B[1::, 0:d]]) `````` Flavio Baccari committed Jun 28, 2018 88 `````` `````` Peter Wittek committed Jun 29, 2018 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 committed Jun 28, 2018 94 95 96 97 98 `````` sub_conf[el] = configs[b][e] config = [] `````` Flavio Baccari committed Jul 15, 2018 99 `````` for el in flatten(s_variables): `````` Peter Wittek committed Jun 29, 2018 100 `````` `````` Flavio Baccari committed Jun 28, 2018 101 `````` config.append(sub_conf[el]) `````` Peter Wittek committed Jun 29, 2018 102 103 `````` bound = get_bound_from_det(flatten(config), one, two) `````` Flavio Baccari committed Jun 28, 2018 104 105 `````` final_config = cp.copy(config) `````` Peter Wittek committed Jun 29, 2018 106 `````` # Check if the bound can be improved by flipping a single spin `````` Flavio Baccari committed Jun 28, 2018 107 108 `````` for i in range(len(config)): `````` Peter Wittek committed Jun 29, 2018 109 `````` `````` Flavio Baccari committed Jun 28, 2018 110 111 112 `````` flip_config = cp.deepcopy(config) flip_config[i] = cp.deepcopy(-flip_config[i]) `````` Peter Wittek committed Jun 29, 2018 113 `````` bound_tmp = get_bound_from_det(flatten(flip_config), one, two) `````` Flavio Baccari committed Jun 28, 2018 114 115 `````` if bound_tmp < bound: `````` Peter Wittek committed Jun 29, 2018 116 `````` `````` Flavio Baccari committed Jun 28, 2018 117 `````` bound = cp.deepcopy(bound_tmp) `````` Peter Wittek committed Jun 29, 2018 118 119 `````` final_config = cp.deepcopy(flip_config) `````` Flavio Baccari committed Jul 15, 2018 120 `````` return final_config `````` Peter Wittek committed Jun 29, 2018 121 `````` `````` Flavio Baccari committed Jun 28, 2018 122 `````` `````` Peter Wittek committed Jun 29, 2018 123 ``````def get_bound_from_det(config, one, two): `````` Flavio Baccari committed Jul 15, 2018 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 committed Jun 29, 2018 136 `````` # Config has to be flattened `````` Flavio Baccari committed Jun 28, 2018 137 138 `````` bound = 0 `````` Peter Wittek committed Jun 29, 2018 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 committed Jun 28, 2018 152 153 154 `````` return bound `````` Flavio Baccari committed Jun 29, 2018 155 ``````def get_strengthen_low(s_variables, substitutions, hamiltonian, cliques, `````` Peter Wittek committed Jun 29, 2018 156 `````` threshold, solverparameters, verbose): `````` Flavio Baccari committed Jun 29, 2018 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 """ `````` Flavio Baccari committed Jun 29, 2018 179 `````` N = len(flatten(s_variables)) `````` Flavio Baccari committed Jun 28, 2018 180 181 182 `````` n_ineqs = N new_ineqs = [] ineqs = [] `````` Peter Wittek committed Jun 29, 2018 183 184 `````` # Generate the mixed level with blocks al level 2 only if the are below # the threshold size `````` Flavio Baccari committed Jun 28, 2018 185 `````` monomial_blocks = [] `````` Peter Wittek committed Jun 29, 2018 186 `````` `````` Flavio Baccari committed Jun 28, 2018 187 `````` for el in cliques: `````` Peter Wittek committed Jun 29, 2018 188 `````` `````` Flavio Baccari committed Jun 28, 2018 189 `````` if len(el) <= threshold: `````` Peter Wittek committed Jun 29, 2018 190 191 `````` monomial_blocks.append(get_monomials(el, 2)) `````` Flavio Baccari committed Jun 28, 2018 192 `````` else: `````` Peter Wittek committed Jun 29, 2018 193 194 195 196 `````` monomial_blocks.append(get_monomials(el, 1)) # Solve the first SDP `````` Flavio Baccari committed Jun 29, 2018 197 `````` sdp = SdpRelaxation(flatten(s_variables), verbose=verbose) `````` Peter Wittek committed Jun 29, 2018 198 199 200 201 202 203 204 `````` sdp.get_relaxation( -1, substitutions=substitutions, objective=hamiltonian, momentinequalities=flatten(ineqs), extramonomials=monomial_blocks) `````` Flavio Baccari committed Jun 28, 2018 205 `````` sdp.solve(solver="mosek", solverparameters=solverparameters) `````` Peter Wittek committed Jun 29, 2018 206 `````` `````` Flavio Baccari committed Jun 28, 2018 207 208 `````` bound_init = cp.deepcopy(sdp.primal) bound_new = cp.deepcopy(sdp.primal) `````` Flavio Baccari committed Jul 15, 2018 209 `````` `````` Peter Wittek committed Jun 29, 2018 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 committed Jun 28, 2018 214 `````` improv = 1 `````` Peter Wittek committed Jun 29, 2018 215 216 217 `````` while (improv > 0.005): `````` Flavio Baccari committed Jun 28, 2018 218 219 `````` bound_old = cp.deepcopy(sdp.primal) ineqs_values = [] `````` Peter Wittek committed Jun 29, 2018 220 `````` `````` Flavio Baccari committed Jun 28, 2018 221 `````` for c in range(len(cliques)): `````` Peter Wittek committed Jun 29, 2018 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 committed Jun 28, 2018 228 `````` tmp = [] `````` Peter Wittek committed Jun 29, 2018 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)]) `````` Flavio Baccari committed Jun 29, 2018 237 `````` # Select the first N most violated ineqs and add the corresponding 4 `````` Peter Wittek committed Jun 29, 2018 238 239 `````` # constraints to the problem `````` Flavio Baccari committed Jun 28, 2018 240 `````` for _ in range(n_ineqs): `````` Peter Wittek committed Jun 29, 2018 241 `````` `````` Flavio Baccari committed Jun 28, 2018 242 243 `````` b = [el[1] for el in ineqs_values] index = min(enumerate(b), key=itemgetter(1))[0] `````` Peter Wittek committed Jun 29, 2018 244 `````` [c, i, j, k] = ineqs_values[index][0] `````` Flavio Baccari committed Jun 28, 2018 245 `````` ineqs_values.remove(ineqs_values[index]) `````` Peter Wittek committed Jun 29, 2018 246 247 248 249 250 251 `````` variables = [ flatten(cliques[c])[i], flatten(cliques[c])[j], flatten(cliques[c])[k] ] `````` Flavio Baccari committed Jun 28, 2018 252 `````` new_ineqs.append(get_triang_ineqs(variables)) `````` Peter Wittek committed Jun 29, 2018 253 254 `````` # Solve the new sdp `````` Flavio Baccari committed Jul 15, 2018 255 `````` `````` Flavio Baccari committed Jun 28, 2018 256 `````` time0 = time() `````` Peter Wittek committed Jun 29, 2018 257 `````` `````` Flavio Baccari committed Jun 29, 2018 258 `````` sdp = SdpRelaxation(flatten(s_variables), verbose=verbose) `````` Peter Wittek committed Jun 29, 2018 259 260 261 262 263 264 265 `````` sdp.get_relaxation( -1, substitutions=substitutions, objective=hamiltonian, momentinequalities=flatten(new_ineqs), extramonomials=monomial_blocks) `````` Flavio Baccari committed Jun 28, 2018 266 `````` sdp.solve(solver="mosek", solverparameters=solverparameters) `````` Peter Wittek committed Jun 29, 2018 267 `````` `````` Flavio Baccari committed Jun 28, 2018 268 `````` bound_new = cp.deepcopy(sdp.primal) `````` Peter Wittek committed Jun 29, 2018 269 `````` improv = abs(bound_new - bound_old) / abs(bound_init) `````` Peter Wittek committed Aug 03, 2018 270 `````` `````` Flavio Baccari committed Jun 29, 2018 271 `````` #if the relative improvement is good enough add the new inequalities `````` Peter Wittek committed Aug 03, 2018 272 `````` `````` Flavio Baccari committed Jun 28, 2018 273 `````` if (improv > 0.005): `````` Peter Wittek committed Jun 29, 2018 274 `````` `````` Flavio Baccari committed Jun 28, 2018 275 `````` ineqs = cp.deepcopy(new_ineqs) `````` Peter Wittek committed Jun 29, 2018 276 277 278 279 280 `````` return sdp, ineqs def get_Mom(sdp, cliques): `````` Flavio Baccari committed Jul 15, 2018 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 committed Jun 28, 2018 292 `````` `````` Peter Wittek committed Jun 29, 2018 293 `````` # Extract the moment matrix for the cliques `````` Flavio Baccari committed Jun 28, 2018 294 `````` `````` Peter Wittek committed Jun 29, 2018 295 `````` lengths = [len(el) for el in cliques] `````` Flavio Baccari committed Jun 28, 2018 296 297 `````` Mom = [] spins = [] `````` Peter Wittek committed Jun 29, 2018 298 `````` `````` Flavio Baccari committed Jun 28, 2018 299 `````` for c in range(len(cliques)): `````` Peter Wittek committed Jun 29, 2018 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 `````` Flavio Baccari committed Jun 29, 2018 309 ``````def get_up_and_low(s_variables, substitutions, ineqs, hamiltonian, one, two, `````` Peter Wittek committed Jun 29, 2018 310 `````` cliques, threshold, solverparameters, verbose): `````` Flavio Baccari committed Jun 29, 2018 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 committed Jun 29, 2018 333 `````` `````` Flavio Baccari committed Jun 29, 2018 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 committed Jun 29, 2018 337 338 `````` """ `````` Flavio Baccari committed Jun 28, 2018 339 340 341 342 343 `````` time0 = time() monomial_blocks = [] for el in cliques: `````` Peter Wittek committed Jun 29, 2018 344 `````` `````` Flavio Baccari committed Jun 28, 2018 345 `````` if len(el) <= threshold: `````` Peter Wittek committed Jun 29, 2018 346 347 `````` monomial_blocks.append(get_monomials(el, 2)) `````` Flavio Baccari committed Jun 28, 2018 348 `````` else: `````` Peter Wittek committed Jun 29, 2018 349 350 `````` monomial_blocks.append(get_monomials(el, 1)) `````` Flavio Baccari committed Jun 28, 2018 351 `````` #solving the first sdp `````` Peter Wittek committed Jun 29, 2018 352 `````` `````` Flavio Baccari committed Jun 29, 2018 353 `````` sdp = SdpRelaxation(flatten(s_variables), verbose=verbose) `````` Peter Wittek committed Jun 29, 2018 354 355 356 357 358 359 360 `````` sdp.get_relaxation( -1, substitutions=substitutions, objective=hamiltonian, momentinequalities=flatten(ineqs), extramonomials=monomial_blocks) `````` Flavio Baccari committed Jun 28, 2018 361 `````` sdp.solve(solver="mosek", solverparameters=solverparameters) `````` Peter Wittek committed Jun 29, 2018 362 `````` `````` Flavio Baccari committed Jun 28, 2018 363 `````` print("The lower bound is " + str(sdp.primal)) `````` Peter Wittek committed Jun 29, 2018 364 365 `````` print("Time taken for lower is " + str(time() - time0)) `````` Flavio Baccari committed Jun 28, 2018 366 `````` time0 = time() `````` Peter Wittek committed Jun 29, 2018 367 368 `````` [Mom, spins] = get_Mom(sdp, cliques) `````` Flavio Baccari committed Jul 16, 2018 369 `````` config = get_det_guess(s_variables, cliques, Mom, one, two) `````` Peter Wittek committed Jun 29, 2018 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 committed Aug 03, 2018 376 ``````def get_branching_node(spins, s_variables, cliques, eqs): `````` Flavio Baccari committed Jul 15, 2018 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 committed Jun 29, 2018 391 `````` `````` Flavio Baccari committed Jun 28, 2018 392 `````` det = 1 `````` Peter Wittek committed Aug 03, 2018 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 committed Jun 28, 2018 398 `````` if el not in sup: `````` Peter Wittek committed Aug 03, 2018 399 `````` `````` Flavio Baccari committed Jul 16, 2018 400 `````` if abs(el) < det: `````` Peter Wittek committed Aug 03, 2018 401 `````` `````` Flavio Baccari committed Jul 16, 2018 402 `````` det = cp.deepcopy(abs(el)) `````` Flavio Baccari committed Jun 28, 2018 403 `````` branch_node = cp.deepcopy(cliques[c][i]) `````` Flavio Baccari committed Jul 16, 2018 404 `````` `````` Flavio Baccari committed Jun 28, 2018 405 `````` return branch_node `````` Peter Wittek committed Jun 29, 2018 406 407 `````` `````` Flavio Baccari committed Jul 15, 2018 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 committed Aug 03, 2018 419 `````` `````` Flavio Baccari committed Jun 28, 2018 420 `````` one = [] `````` Peter Wittek committed Jun 29, 2018 421 `````` two = [] `````` Flavio Baccari committed Jun 28, 2018 422 `````` pol = hamiltonian.as_poly() `````` Flavio Baccari committed Jul 15, 2018 423 `````` mon = flatten(s_variables) `````` Peter Wittek committed Jun 29, 2018 424 425 426 `````` for i, el in enumerate(mon): `````` Flavio Baccari committed Jun 28, 2018 427 428 `````` one.append(pol.coeff_monomial(el)) tmp = [] `````` Peter Wittek committed Jun 29, 2018 429 430 `````` for j, el2 in enumerate(mon): `````` Flavio Baccari committed Jun 28, 2018 431 `````` if j < i: `````` Peter Wittek committed Jun 29, 2018 432 `````` `````` Flavio Baccari committed Jun 28, 2018 433 434 `````` tmp.append(0) else: `````` Peter Wittek committed Jun 29, 2018 435 436 `````` tmp.append(pol.coeff_monomial(el * el2)) `````` Flavio Baccari committed Jun 28, 2018 437 438 `````` two.append(tmp) `````` Peter Wittek committed Jun 29, 2018 439 440 441 `````` return one, two `````` Flavio Baccari committed Jul 15, 2018 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 committed Aug 03, 2018 452 `````` `````` Flavio Baccari committed Jun 28, 2018 453 `````` sup = [] `````` Peter Wittek committed Jun 29, 2018 454 `````` `````` Flavio Baccari committed Jun 28, 2018 455 `````` for el in eqs: `````` Peter Wittek committed Jun 29, 2018 456 `````` `````` Flavio Baccari committed Jul 15, 2018 457 458 `````` s = get_support(flatten(s_variables), el) sup.append(dot(s[1], flatten(s_variables))) `````` Peter Wittek committed Jun 29, 2018 459 `````` `````` Flavio Baccari committed Jun 28, 2018 460 461 462 `````` return sup `````` Peter Wittek committed Aug 03, 2018 463 464 465 466 467 468 469 ``````def get_groundBandB(s_variables, substitutions, hamiltonian, cliques, threshold, solverparameters, verbose=0): `````` Flavio Baccari committed Jun 29, 2018 470 471 `````` """Runs the whole branch&bound algorithm for the given hamiltonian `````` Flavio Baccari committed Jun 29, 2018 472 473 `````` :param s_variables: list of spin variables :type s_variables: list of sympy.core.symbol.Symbol `````` Flavio Baccari committed Jun 29, 2018 474 `````` :param substitutions: substitution to be applied in the generation of the moment matrix at the level of operators `````` Flavio Baccari committed Jun 29, 2018 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 committed Jun 29, 2018 484 485 `````` :param verbose: verbosity level :type verbose: int `````` Peter Wittek committed Jun 29, 2018 486 `````` `````` Flavio Baccari committed Jun 29, 2018 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 `````` Flavio Baccari committed Jun 29, 2018 490 `````` """ `````` Peter Wittek committed Jun 29, 2018 491 `````` `````` Flavio Baccari committed Jun 29, 2018 492 `````` [one, two] = get_as_coeffs(hamiltonian, s_variables) `````` Peter Wittek committed Jun 29, 2018 493 494 `````` # First root `````` Flavio Baccari committed Jun 28, 2018 495 496 497 `````` eqs = [] tree = [] print("Starting initial optimisation") `````` Peter Wittek committed Jun 29, 2018 498 499 500 `````` # Get the first lower bound and choosing the list of triangle ineqs to add `````` Peter Wittek committed Jun 29, 2018 501 `````` [sdp, ineqs] = get_strengthen_low( `````` Flavio Baccari committed Jun 29, 2018 502 `````` s_variables, `````` Peter Wittek committed Jun 29, 2018 503 504 505 506 507 508 `````` substitutions, hamiltonian, cliques, threshold, solverparameters, verbose=0) `````` Peter Wittek committed Jun 29, 2018 509 510 511 `````` # Extract the corresponding upper bound and spin configuration `````` Peter Wittek committed Jun 29, 2018 512 `````` [Mom, spins] = get_Mom(sdp, cliques) `````` Flavio Baccari committed Jul 16, 2018 513 `````` config = get_det_guess(s_variables, cliques, Mom, one, two) `````` Peter Wittek committed Jun 29, 2018 514 `````` final_config = cp.deepcopy(config) `````` Flavio Baccari committed Jun 28, 2018 515 `````` `````` Peter Wittek committed Jun 29, 2018 516 `````` z_up = get_bound_from_det(flatten(config), one, two) `````` Flavio Baccari committed Jun 28, 2018 517 `````` z_low = sdp.primal `````` Peter Wittek committed Jun 29, 2018 518 `````` `````` Flavio Baccari committed Jun 28, 2018 519 520 `````` print("The lower bound is " + str(z_low)) print("The upper bound is " + str(z_up)) `````` Peter Wittek committed Jun 29, 2018 521 `````` `````` Peter Wittek committed Jun 29, 2018 522 523 524 525 `````` # Go for the branching with easy first if ((z_low - z_up) / abs(z_up) < -0.00001): `````` Peter Wittek committed Aug 03, 2018 526 `````` node = get_branching_node(spins, s_variables, cliques, eqs) `````` Peter Wittek committed Jun 29, 2018 527 528 529 530 `````` tree = [[z_low, eqs, substitutions, node]] while (len(tree) > 0): `````` Flavio Baccari committed Jun 28, 2018 531 532 `````` lis = [el[0] for el in tree] index = min(enumerate(lis), key=itemgetter(1))[0] `````` Peter Wittek committed Jun 29, 2018 533 `````` [z_low, eqs, subs, node] = tree[index] `````` Flavio Baccari committed Jun 28, 2018 534 `````` tree.remove(tree[index]) `````` Peter Wittek committed Jun 29, 2018 535 536 537 `````` # Branch +1 `````` Flavio Baccari committed Jun 28, 2018 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 committed Jun 29, 2018 542 543 `````` [z_low, tmp, spins, config] = get_up_and_low( `````` Flavio Baccari committed Jun 29, 2018 544 `````` s_variables, `````` Peter Wittek committed Jun 29, 2018 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 committed Jun 28, 2018 556 `````` if tmp < z_up: `````` Peter Wittek committed Jun 29, 2018 557 `````` `````` Flavio Baccari committed Jun 28, 2018 558 559 `````` z_up = cp.deepcopy(tmp) final_config = cp.deepcopy(config) `````` Peter Wittek committed Jun 29, 2018 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 committed Aug 03, 2018 565 `````` node_new = get_branching_node(spins, s_variables, cliques, eqs_new) `````` Peter Wittek committed Jun 29, 2018 566 567 568 569 `````` tree.append([z_low, eqs_new, subs_new, node_new]) # branch -1 `````` Flavio Baccari committed Jun 28, 2018 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 committed Jun 29, 2018 574 575 `````` [z_low, tmp, spins, config] = get_up_and_low( `````` Flavio Baccari committed Jun 29, 2018 576 `````` s_variables, `````` Peter Wittek committed Jun 29, 2018 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 committed Jun 28, 2018 588 `````` if tmp < z_up: `````` Peter Wittek committed Jun 29, 2018 589 `````` `````` Flavio Baccari committed Jun 28, 2018 590 591 `````` z_up = cp.deepcopy(tmp) final_config = cp.deepcopy(config) `````` Peter Wittek committed Jun 29, 2018 592 593 `````` # Include this node in the branching if the lower bound is meaningful `````` Flavio Baccari committed Jun 28, 2018 594 `````` print("The upper bound is " + str(z_up)) `````` Peter Wittek committed Jun 29, 2018 595 596 `````` if ((z_low - z_up) / abs(z_up) < -0.00001): `````` Peter Wittek committed Aug 03, 2018 597 `````` node_new = get_branching_node(spins, s_variables, cliques, eqs_new) `````` Peter Wittek committed Jun 29, 2018 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 committed Jun 28, 2018 602 603 604 `````` tree_old = cp.deepcopy(tree) tree = [] `````` Peter Wittek committed Jun 29, 2018 605 `````` `````` Flavio Baccari committed Jun 28, 2018 606 607 `````` for el in tree_old: `````` Peter Wittek committed Jun 29, 2018 608 `````` if ((el[0] - z_up) / abs(z_up) < -0.00001): `````` Flavio Baccari committed Jun 28, 2018 609 `````` `````` Peter Wittek committed Jun 29, 2018 610 `````` tree.append(el) `````` Flavio Baccari committed Jun 28, 2018 611 `````` `````` Peter Wittek committed Jun 29, 2018 612 `` return z_low, z_up, final_config``