specific_analyses.model_free.optimisation

1 ############################################################################### 2 # # 3 # Copyright (C) 2003-2014 Edward d'Auvergne # 4 # Copyright (C) 2007 Gary S Thompson (https://gna.org/users/varioustoxins) # 5 # # 6 # This file is part of the program relax (http://www.nmr-relax.com). # 7 # # 8 # This program is free software: you can redistribute it and/or modify # 9 # it under the terms of the GNU General Public License as published by # 10 # the Free Software Foundation, either version 3 of the License, or # 11 # (at your option) any later version. # 12 # # 13 # This program is distributed in the hope that it will be useful, # 14 # but WITHOUT ANY WARRANTY; without even the implied warranty of # 15 # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the # 16 # GNU General Public License for more details. # 17 # # 18 # You should have received a copy of the GNU General Public License # 19 # along with this program. If not, see <http://www.gnu.org/licenses/>. # 20 # # 21 ############################################################################### 22 23 # Module docstring. 24 """The model-free analysis optimisation functions.""" 25 26 # Python module imports. 27 from math import pi 28 from minfx.generic import generic_minimise 29 from minfx.grid import grid, grid_point_array 30 from numpy import array, dot, float64 31 from re import match 32 33 # relax module imports. 34 import lib.arg_check 35 from lib.errors import RelaxError, RelaxInfError, RelaxMultiVectorError, RelaxNaNError 36 from lib.float import isNaN, isInf 37 from lib.optimisation import test_grid_ops 38 from lib.physical_constants import return_gyromagnetic_ratio 39 from multi import Memo, Result_command, Slave_command 40 from pipe_control import pipes 41 from pipe_control.interatomic import return_interatom_list 42 from pipe_control.mol_res_spin import return_spin, return_spin_from_index, spin_loop 43 from specific_analyses.model_free.model import determine_model_type 44 from specific_analyses.model_free.parameters import assemble_param_vector, disassemble_param_vector 45 from target_functions.mf import Mf 46 47

48 -def disassemble_result(param_vector=None, func=None, iter=None, fc=None, gc=None, hc=None, warning=None, spin=None, sim_index=None, model_type=None, scaling=None, scaling_matrix=None):

49 """Disassemble the optimisation results. 50 51 @keyword param_vector: The model-free parameter vector. 52 @type param_vector: numpy array 53 @keyword func: The optimised chi-squared value. 54 @type func: float 55 @keyword iter: The number of optimisation steps required to find the minimum. 56 @type iter: int 57 @keyword fc: The function count. 58 @type fc: int 59 @keyword gc: The gradient count. 60 @type gc: int 61 @keyword hc: The Hessian count. 62 @type hc: int 63 @keyword warning: Any optimisation warnings. 64 @type warning: str or None 65 @keyword spin: The spin container. 66 @type spin: SpinContainer instance or None 67 @keyword sim_index: The Monte Carlo simulation index. 68 @type sim_index: int or None 69 @keyword model_type: The model-free model type, one of 'mf', 'local_tm', 'diff', or 70 'all'. 71 @type model_type: str 72 @keyword scaling: If True, diagonal scaling is enabled during optimisation to 73 allow the problem to be better conditioned. 74 @type scaling: bool 75 @keyword scaling_matrix: The diagonal, square scaling matrix. 76 @type scaling_matrix: numpy diagonal matrix 77 """ 78 79 # No result. 80 if param_vector == None: 81 return 82 83 # Alias the current data pipe. 84 cdp = pipes.get_pipe() 85 86 # Catch infinite chi-squared values. 87 if isInf(func): 88 raise RelaxInfError('chi-squared') 89 90 # Catch chi-squared values of NaN. 91 if isNaN(func): 92 raise RelaxNaNError('chi-squared') 93 94 # Scaling. 95 if scaling: 96 param_vector = dot(scaling_matrix, param_vector) 97 98 # Check if the chi-squared value is lower. This allows for a parallelised grid search! 99 if sim_index == None: 100 # Get the correct value. 101 chi2 = None 102 if (model_type == 'mf' or model_type == 'local_tm') and hasattr(cdp, 'chi2'): 103 chi2 = spin.chi2 104 if (model_type == 'diff' or model_type == 'all') and hasattr(cdp, 'chi2'): 105 chi2 = cdp.chi2 106 107 # Spin text. 108 spin_text = '' 109 if spin != None and hasattr(spin, '_spin_ids') and len(spin._spin_ids): 110 spin_text = " for the spin '%s'" % spin._spin_ids[0] 111 112 # No improvement. 113 if chi2 != None and func >= chi2: 114 print("Discarding the optimisation results%s, the optimised chi-squared value is higher than the current value (%s >= %s)." % (spin_text, func, chi2)) 115 116 # Exit! 117 return 118 119 # New minimum. 120 else: 121 print("Storing the optimisation results%s, the optimised chi-squared value is lower than the current value (%s < %s)." % (spin_text, func, chi2)) 122 123 # Disassemble the parameter vector. 124 disassemble_param_vector(model_type, param_vector=param_vector, spin=spin, sim_index=sim_index) 125 126 # Monte Carlo minimisation statistics. 127 if sim_index != None: 128 # Sequence specific minimisation statistics. 129 if model_type == 'mf' or model_type == 'local_tm': 130 131 # Chi-squared statistic. 132 spin.chi2_sim[sim_index] = func 133 134 # Iterations. 135 spin.iter_sim[sim_index] = iter 136 137 # Function evaluations. 138 spin.f_count_sim[sim_index] = fc 139 140 # Gradient evaluations. 141 spin.g_count_sim[sim_index] = gc 142 143 # Hessian evaluations. 144 spin.h_count_sim[sim_index] = hc 145 146 # Warning. 147 spin.warning_sim[sim_index] = warning 148 149 # Global minimisation statistics. 150 elif model_type == 'diff' or model_type == 'all': 151 # Chi-squared statistic. 152 cdp.chi2_sim[sim_index] = func 153 154 # Iterations. 155 cdp.iter_sim[sim_index] = iter 156 157 # Function evaluations. 158 cdp.f_count_sim[sim_index] = fc 159 160 # Gradient evaluations. 161 cdp.g_count_sim[sim_index] = gc 162 163 # Hessian evaluations. 164 cdp.h_count_sim[sim_index] = hc 165 166 # Warning. 167 cdp.warning_sim[sim_index] = warning 168 169 # Normal statistics. 170 else: 171 # Sequence specific minimisation statistics. 172 if model_type == 'mf' or model_type == 'local_tm': 173 # Chi-squared statistic. 174 spin.chi2 = func 175 176 # Iterations. 177 spin.iter = iter 178 179 # Function evaluations. 180 spin.f_count = fc 181 182 # Gradient evaluations. 183 spin.g_count = gc 184 185 # Hessian evaluations. 186 spin.h_count = hc 187 188 # Warning. 189 spin.warning = warning 190 191 # Global minimisation statistics. 192 elif model_type == 'diff' or model_type == 'all': 193 # Chi-squared statistic. 194 cdp.chi2 = func 195 196 # Iterations. 197 cdp.iter = iter 198 199 # Function evaluations. 200 cdp.f_count = fc 201 202 # Gradient evaluations. 203 cdp.g_count = gc 204 205 # Hessian evaluations. 206 cdp.h_count = hc 207 208 # Warning. 209 cdp.warning = warning

210 211

212 -def grid_search_config(num_params, spin=None, spin_id=None, lower=None, upper=None, inc=None, scaling_matrix=None, verbosity=1):

213 """Configure the grid search. 214 215 @param num_params: The number of parameters in the model. 216 @type num_params: int 217 @keyword spin: The spin data container. 218 @type spin: SpinContainer instance 219 @keyword spin_id: The spin identification string. 220 @type spin_id: str 221 @keyword lower: The lower bounds of the grid search which must be equal to the 222 number of parameters in the model. 223 @type lower: array of numbers 224 @keyword upper: The upper bounds of the grid search which must be equal to the 225 number of parameters in the model. 226 @type upper: array of numbers 227 @keyword inc: The increments for each dimension of the space for the grid 228 search. The number of elements in the array must equal to the 229 number of parameters in the model. 230 @type inc: array of int 231 @keyword scaling_matrix: The diagonal, square scaling matrix. 232 @type scaling_matrix: numpy diagonal matrix 233 @keyword verbosity: A flag specifying the amount of information to print. The 234 higher the value, the greater the verbosity. 235 @type verbosity: int 236 """ 237 238 # Test the grid search options. 239 test_grid_ops(lower=lower, upper=upper, inc=inc, n=num_params) 240 241 # If inc is a single int, convert it into an array of that value. 242 if isinstance(inc, int): 243 inc = [inc]*num_params 244 245 # Set up the default bounds. 246 if not lower: 247 # Init. 248 lower = [] 249 upper = [] 250 251 # Determine the model type. 252 model_type = determine_model_type() 253 254 # Minimisation options for diffusion tensor parameters. 255 if model_type == 'diff' or model_type == 'all': 256 # Get the diffusion tensor specific configuration. 257 grid_search_diff_bounds(lower, upper) 258 259 # Model-free parameters (residue specific parameters). 260 if model_type != 'diff': 261 # The loop. 262 if spin: 263 loop = [spin] 264 else: 265 loop = spin_loop(spin_id) 266 267 # Loop over the spins. 268 for spin in loop: 269 # Skip deselected residues. 270 if not spin.select: 271 continue 272 273 # Get the spin specific configuration. 274 grid_search_spin_bounds(spin, lower, upper) 275 276 # Diagonal scaling of minimisation options. 277 lower_new = [] 278 upper_new = [] 279 for i in range(num_params): 280 lower_new.append(lower[i] / scaling_matrix[i, i]) 281 upper_new.append(upper[i] / scaling_matrix[i, i]) 282 283 # Return the minimisation options. 284 return inc, lower_new, upper_new

285 286

287 -def grid_search_diff_bounds(lower, upper):

288 """Set up the default grid search bounds the diffusion tensor. 289 290 This method appends the default bounds to the lower and upper lists. 291 292 @param lower: The lower bound list to append to. 293 @type lower: list 294 @param upper: The upper bound list to append to. 295 @type upper: list 296 """ 297 298 # Spherical diffusion {tm}. 299 if cdp.diff_tensor.type == 'sphere': 300 lower.append(1.0 * 1e-9) 301 upper.append(12.0 * 1e-9) 302 303 # Spheroidal diffusion {tm, Da, theta, phi}. 304 if cdp.diff_tensor.type == 'spheroid': 305 # tm. 306 lower.append(1.0 * 1e-9) 307 upper.append(12.0 * 1e-9) 308 309 # Da. 310 if cdp.diff_tensor.spheroid_type == 'prolate': 311 lower.append(0.0) 312 upper.append(1e7) 313 elif cdp.diff_tensor.spheroid_type == 'oblate': 314 lower.append(-1e7) 315 upper.append(0.0) 316 else: 317 lower.append(-1e7) 318 upper.append(1e7) 319 320 # theta. 321 lower.append(0.0) 322 upper.append(pi) 323 324 # phi. 325 lower.append(0.0) 326 upper.append(pi) 327 328 # Ellipsoidal diffusion {tm, Da, Dr, alpha, beta, gamma}. 329 elif cdp.diff_tensor.type == 'ellipsoid': 330 # tm. 331 lower.append(1.0 * 1e-9) 332 upper.append(12.0 * 1e-9) 333 334 # Da. 335 lower.append(0.0) 336 upper.append(1e7) 337 338 # Dr. 339 lower.append(0.0) 340 upper.append(1.0) 341 342 # alpha. 343 lower.append(0.0) 344 upper.append(pi) 345 346 # beta. 347 lower.append(0.0) 348 upper.append(pi) 349 350 # gamma. 351 lower.append(0.0) 352 upper.append(pi)

353 354

355 -def grid_search_spin_bounds(spin, lower, upper):

356 """Set up the default grid search bounds for a single spin. 357 358 This method appends the default bounds to the lower and upper lists. The ordering of the 359 lists in min_options matches that of the params list in the spin container. 360 361 @param spin: A SpinContainer object. 362 @type spin: class instance 363 @param lower: The lower bound list to append to. 364 @type lower: list 365 @param upper: The upper bound list to append to. 366 @type upper: list 367 """ 368 369 # Loop over the model-free parameters. 370 for i in range(len(spin.params)): 371 # Local tm. 372 if spin.params[i] == 'local_tm': 373 lower.append(1.0 * 1e-9) 374 upper.append(12.0 * 1e-9) 375 376 # {S2, S2f, S2s}. 377 elif match('s2', spin.params[i]): 378 lower.append(0.0) 379 upper.append(1.0) 380 381 # {te, tf, ts}. 382 elif match('t', spin.params[i]): 383 lower.append(0.0) 384 upper.append(500.0 * 1e-12) 385 386 # Rex. 387 elif spin.params[i] == 'rex': 388 lower.append(0.0) 389 upper.append(5.0 / (2.0 * pi * cdp.spectrometer_frq[cdp.ri_ids[0]])**2) 390 391 # Bond length. 392 elif spin.params[i] == 'r': 393 lower.append(1.0 * 1e-10) 394 upper.append(1.05 * 1e-10) 395 396 # CSA. 397 elif spin.params[i] == 'csa': 398 lower.append(-120 * 1e-6) 399 upper.append(-200 * 1e-6) 400 401 # Unknown option. 402 else: 403 raise RelaxError("Unknown model-free parameter.")

404 405

406 -def minimise_data_setup(data_store, min_algor, num_data_sets, min_options, spin=None, sim_index=None):

407 """Set up all the data required for minimisation. 408 409 @param data_store: A data storage container. 410 @type data_store: class instance 411 @param min_algor: The minimisation algorithm to use. 412 @type min_algor: str 413 @param num_data_sets: The number of data sets. 414 @type num_data_sets: int 415 @param min_options: The minimisation options array. 416 @type min_options: list 417 @keyword spin: The spin data container. 418 @type spin: SpinContainer instance 419 @keyword sim_index: The optional MC simulation index. 420 @type sim_index: int 421 @return: An insane tuple. The full tuple is (ri_data, ri_data_err, equations, param_types, param_values, r, csa, num_frq, frq, num_ri, remap_table, noe_r1_table, ri_types, num_params, xh_unit_vectors, diff_type, diff_params) 422 @rtype: tuple 423 """ 424 425 # Initialise the data structures for the model-free function. 426 data_store.ri_data = [] 427 data_store.ri_data_err = [] 428 data_store.equations = [] 429 data_store.param_types = [] 430 data_store.param_values = None 431 data_store.r = [] 432 data_store.csa = [] 433 data_store.num_frq = [] 434 data_store.frq = [] 435 data_store.num_ri = [] 436 data_store.remap_table = [] 437 data_store.noe_r1_table = [] 438 data_store.ri_types = [] 439 data_store.gx = [] 440 data_store.gh = [] 441 data_store.num_params = [] 442 data_store.xh_unit_vectors = [] 443 if data_store.model_type == 'local_tm': 444 data_store.mf_params = [] 445 elif data_store.model_type == 'diff': 446 data_store.param_values = [] 447 448 # Set up the data for the back_calc function. 449 if min_algor == 'back_calc': 450 # The spin data. 451 data_store.ri_data = [0.0] 452 data_store.ri_data_err = [0.000001] 453 data_store.equations = [spin.equation] 454 data_store.param_types = [spin.params] 455 data_store.csa = [spin.csa] 456 data_store.num_frq = [1] 457 data_store.frq = [[min_options[3]]] 458 data_store.num_ri = [1] 459 data_store.remap_table = [[0]] 460 data_store.noe_r1_table = [[None]] 461 data_store.ri_types = [[min_options[2]]] 462 data_store.gx = [return_gyromagnetic_ratio(spin.isotope)] 463 464 # The interatomic data. 465 interatoms = return_interatom_list(data_store.spin_id) 466 for i in range(len(interatoms)): 467 # No relaxation mechanism. 468 if not interatoms[i].dipole_pair: 469 continue 470 471 # The surrounding spins. 472 if data_store.spin_id != interatoms[i].spin_id1: 473 spin_id2 = interatoms[i].spin_id1 474 else: 475 spin_id2 = interatoms[i].spin_id2 476 spin2 = return_spin(spin_id2) 477 478 # The data. 479 data_store.r = [interatoms[i].r] 480 data_store.gh = [return_gyromagnetic_ratio(spin2.isotope)] 481 if data_store.model_type != 'local_tm' and cdp.diff_tensor.type != 'sphere': 482 data_store.xh_unit_vectors = [interatoms[i].vector] 483 else: 484 data_store.xh_unit_vectors = [None] 485 486 # Count the number of model-free parameters for the spin index. 487 data_store.num_params = [len(spin.params)] 488 489 # Loop over the number of data sets. 490 for j in range(num_data_sets): 491 # Set the spin index and get the spin, if not already set. 492 if data_store.model_type == 'diff' or data_store.model_type == 'all': 493 spin_index = j 494 spin, data_store.spin_id = return_spin_from_index(global_index=spin_index, return_spin_id=True) 495 496 # Skip deselected spins. 497 if not spin.select: 498 continue 499 500 # Skip spins where there is no data or errors. 501 if not hasattr(spin, 'ri_data') or not hasattr(spin, 'ri_data_err'): 502 continue 503 504 # Make sure that the errors are strictly positive numbers. 505 for ri_id in cdp.ri_ids: 506 # Skip missing data. 507 if not ri_id in spin.ri_data_err: 508 continue 509 510 # Alias. 511 err = spin.ri_data_err[ri_id] 512 513 # Checks. 514 if err != None and err == 0.0: 515 raise RelaxError("Zero error for spin '%s' for the relaxation data ID '%s', minimisation not possible." % (data_store.spin_id, ri_id)) 516 elif err != None and err < 0.0: 517 raise RelaxError("Negative error of %s for spin '%s' for the relaxation data ID '%s', minimisation not possible." % (err, data_store.spin_id, ri_id)) 518 519 # The relaxation data optimisation structures. 520 data = relax_data_opt_structs(spin, sim_index=sim_index) 521 522 # Append the data. 523 data_store.ri_data.append(data[0]) 524 data_store.ri_data_err.append(data[1]) 525 data_store.num_frq.append(data[2]) 526 data_store.num_ri.append(data[3]) 527 data_store.ri_types.append(data[4]) 528 data_store.frq.append(data[5]) 529 data_store.remap_table.append(data[6]) 530 data_store.noe_r1_table.append(data[7]) 531 if sim_index == None or data_store.model_type == 'diff': 532 data_store.csa.append(spin.csa) 533 else: 534 data_store.csa.append(spin.csa_sim[sim_index]) 535 536 # Repackage the spin data. 537 data_store.equations.append(spin.equation) 538 data_store.param_types.append(spin.params) 539 data_store.gx.append(return_gyromagnetic_ratio(spin.isotope)) 540 541 # Repackage the interatomic data. 542 interatoms = return_interatom_list(data_store.spin_id) 543 for i in range(len(interatoms)): 544 # No relaxation mechanism. 545 if not interatoms[i].dipole_pair: 546 continue 547 548 # The surrounding spins. 549 if data_store.spin_id != interatoms[i].spin_id1: 550 spin_id2 = interatoms[i].spin_id1 551 else: 552 spin_id2 = interatoms[i].spin_id2 553 spin2 = return_spin(spin_id2) 554 555 # The data. 556 data_store.gh.append(return_gyromagnetic_ratio(spin2.isotope)) 557 if sim_index == None or data_store.model_type == 'diff' or not hasattr(interatoms[i], 'r_sim'): 558 data_store.r.append(interatoms[i].r) 559 else: 560 data_store.r.append(interatoms[i].r_sim[sim_index]) 561 562 # Vectors. 563 if data_store.model_type != 'local_tm' and cdp.diff_tensor.type != 'sphere': 564 # Check that this is a single vector! 565 if lib.arg_check.is_num_list(interatoms[i].vector[0], raise_error=False): 566 raise RelaxMultiVectorError(data_store.spin_id) 567 568 # Store the vector. 569 data_store.xh_unit_vectors.append(interatoms[i].vector) 570 571 # No vector. 572 else: 573 data_store.xh_unit_vectors.append(None) 574 575 # Stop - only one mechanism is current supported. 576 break 577 578 # Model-free parameter values. 579 if data_store.model_type == 'local_tm': 580 pass 581 582 # Count the number of model-free parameters for the spin index. 583 data_store.num_params.append(len(spin.params)) 584 585 # Repackage the parameter values for minimising just the diffusion tensor parameters. 586 if data_store.model_type == 'diff': 587 data_store.param_values.append(assemble_param_vector(model_type='mf')) 588 589 # Convert to numpy arrays. 590 for k in range(len(data_store.ri_data)): 591 data_store.ri_data[k] = array(data_store.ri_data[k], float64) 592 data_store.ri_data_err[k] = array(data_store.ri_data_err[k], float64) 593 594 # Diffusion tensor type. 595 if data_store.model_type == 'local_tm': 596 data_store.diff_type = 'sphere' 597 else: 598 data_store.diff_type = cdp.diff_tensor.type 599 600 # Package the diffusion tensor parameters. 601 data_store.diff_params = None 602 if data_store.model_type == 'mf': 603 # Spherical diffusion. 604 if data_store.diff_type == 'sphere': 605 data_store.diff_params = [cdp.diff_tensor.tm] 606 607 # Spheroidal diffusion. 608 elif data_store.diff_type == 'spheroid': 609 data_store.diff_params = [cdp.diff_tensor.tm, cdp.diff_tensor.Da, cdp.diff_tensor.theta, cdp.diff_tensor.phi] 610 611 # Ellipsoidal diffusion. 612 elif data_store.diff_type == 'ellipsoid': 613 data_store.diff_params = [cdp.diff_tensor.tm, cdp.diff_tensor.Da, cdp.diff_tensor.Dr, cdp.diff_tensor.alpha, cdp.diff_tensor.beta, cdp.diff_tensor.gamma] 614 elif min_algor == 'back_calc' and data_store.model_type == 'local_tm': 615 # Spherical diffusion. 616 data_store.diff_params = [spin.local_tm]

617 618

619 -def relax_data_opt_structs(spin, sim_index=None):

620 """Package the relaxation data into the data structures used for optimisation. 621 622 @param spin: The spin container to extract the data from. 623 @type spin: SpinContainer instance 624 @keyword sim_index: The optional MC simulation index. 625 @type sim_index: int 626 @return: The structures ri_data, ri_data_err, num_frq, num_ri, ri_ids, frq, remap_table, noe_r1_table. 627 @rtype: tuple 628 """ 629 630 # Initialise the data. 631 ri_data = [] 632 ri_data_err = [] 633 ri_labels = [] 634 frq = [] 635 remap_table = [] 636 noe_r1_table = [] 637 638 # Loop over the relaxation data. 639 for ri_id in cdp.ri_ids: 640 # Skip missing data. 641 if ri_id not in spin.ri_data: 642 continue 643 644 # The Rx data. 645 if sim_index == None: 646 data = spin.ri_data[ri_id] 647 else: 648 data = spin.ri_data_sim[ri_id][sim_index] 649 650 # The errors. 651 err = spin.ri_data_err[ri_id] 652 653 # Missing data, so don't add it. 654 if data == None or err == None: 655 continue 656 657 # Append the data and error. 658 ri_data.append(data) 659 ri_data_err.append(err) 660 661 # The labels. 662 ri_labels.append(cdp.ri_type[ri_id]) 663 664 # The frequencies. 665 if cdp.spectrometer_frq[ri_id] not in frq: 666 frq.append(cdp.spectrometer_frq[ri_id]) 667 668 # The remap table. 669 remap_table.append(frq.index(cdp.spectrometer_frq[ri_id])) 670 671 # The NOE to R1 mapping table. 672 noe_r1_table.append(None) 673 674 # The number of data sets. 675 num_ri = len(ri_data) 676 677 # Fill the NOE to R1 mapping table. 678 for i in range(num_ri): 679 # If the data corresponds to 'NOE', try to find if the corresponding R1 data. 680 if cdp.ri_type[cdp.ri_ids[i]] == 'NOE': 681 for j in range(num_ri): 682 if cdp.ri_type[cdp.ri_ids[j]] == 'R1' and cdp.spectrometer_frq[cdp.ri_ids[i]] == cdp.spectrometer_frq[cdp.ri_ids[j]]: 683 noe_r1_table[i] = j 684 685 # Return the structures. 686 return ri_data, ri_data_err, len(frq), num_ri, ri_labels, frq, remap_table, noe_r1_table

687 688

689 -def reset_min_stats():

690 """Reset all the minimisation statistics. 691 692 All global and spin specific values will be set to None. 693 """ 694 695 # Global stats. 696 if hasattr(cdp, 'chi2'): 697 cdp.chi2 = None 698 cdp.iter = None 699 cdp.f_count = None 700 cdp.g_count = None 701 cdp.h_count = None 702 cdp.warning = None 703 704 # Spin specific stats. 705 for spin in spin_loop(): 706 if hasattr(spin, 'chi2'): 707 spin.chi2 = None 708 spin.iter = None 709 spin.f_count = None 710 spin.g_count = None 711 spin.h_count = None 712 spin.warning = None

713 714

715 -def spin_print(spin_id, verbosity):

716 """Print out some header text for the spin. 717 718 @param spin_id: The spin ID string. 719 @type spin_id: str 720 @param verbosity: The amount of information to print. The higher the value, the greater the verbosity. 721 @type verbosity: int 722 """ 723 724 # Some extra spacing for verbose printouts. 725 if verbosity >= 2: 726 print("\n\n") 727 728 # The header. 729 string = "Fitting to spin " + repr(spin_id) 730 print("\n\n" + string) 731 print(len(string) * '~')

732 733 734

735 -class MF_memo(Memo):

736 """The model-free memo class. 737 738 Not quite a momento so a memo. 739 """ 740

741 - def __init__(self, model_free=None, model_type=None, spin=None, sim_index=None, scaling=None, scaling_matrix=None):

742 """Initialise the model-free memo class. 743 744 This memo stores the model-free class instance so that the disassemble_result() method can be called to store the optimisation results. The other args are those required by this method but not generated through optimisation. 745 746 @keyword model_free: The model-free class instance. 747 @type model_free: specific_analyses.model_free.Model_free instance 748 @keyword spin: The spin data container. If this argument is supplied, then the spin_id argument will be ignored. 749 @type spin: SpinContainer instance 750 @keyword sim_index: The optional MC simulation index. 751 @type sim_index: int 752 @keyword scaling: If True, diagonal scaling is enabled. 753 @type scaling: bool 754 @keyword scaling_matrix: The diagonal, square scaling matrix. 755 @type scaling_matrix: numpy diagonal matrix 756 """ 757 758 # Execute the base class __init__() method. 759 super(MF_memo, self).__init__() 760 761 # Store the arguments. 762 self.model_free = model_free 763 self.model_type = model_type 764 self.spin = spin 765 self.sim_index = sim_index 766 self.scaling = scaling 767 self.scaling_matrix = scaling_matrix

768 769 770

771 -class MF_minimise_command(Slave_command):

772 """Command class for standard model-free minimisation.""" 773

774 - def __init__(self):

775 """Initialise the base class.""" 776 777 # Execute the base class __init__() method. 778 super(MF_minimise_command, self).__init__()

779 780

781 - def optimise(self):

782 """Model-free optimisation. 783 784 @return: The optimisation results consisting of the parameter vector, function value, iteration count, function count, gradient count, Hessian count, and warnings. 785 @rtype: tuple of numpy array, float, int, int, int, int, str 786 """ 787 788 # Minimisation. 789 results = generic_minimise(func=self.mf.func, dfunc=self.mf.dfunc, d2func=self.mf.d2func, args=(), x0=self.opt_params.param_vector, min_algor=self.opt_params.min_algor, min_options=self.opt_params.min_options, func_tol=self.opt_params.func_tol, grad_tol=self.opt_params.grad_tol, maxiter=self.opt_params.max_iterations, A=self.opt_params.A, b=self.opt_params.b, full_output=True, print_flag=self.opt_params.verbosity) 790 791 # Return the minfx results unmodified. 792 return results

793 794

795 - def run(self, processor, completed):

796 """Setup and perform the model-free optimisation.""" 797 798 # Initialise the function to minimise. 799 self.mf = Mf(init_params=self.opt_params.param_vector, model_type=self.data.model_type, diff_type=self.data.diff_type, diff_params=self.data.diff_params, scaling_matrix=self.data.scaling_matrix, num_spins=self.data.num_spins, equations=self.data.equations, param_types=self.data.param_types, param_values=self.data.param_values, relax_data=self.data.ri_data, errors=self.data.ri_data_err, bond_length=self.data.r, csa=self.data.csa, num_frq=self.data.num_frq, frq=self.data.frq, num_ri=self.data.num_ri, remap_table=self.data.remap_table, noe_r1_table=self.data.noe_r1_table, ri_labels=self.data.ri_types, gx=self.data.gx, gh=self.data.gh, h_bar=self.data.h_bar, mu0=self.data.mu0, num_params=self.data.num_params, vectors=self.data.xh_unit_vectors) 800 801 # Print out. 802 if self.opt_params.verbosity >= 1 and (self.data.model_type == 'mf' or self.data.model_type == 'local_tm'): 803 spin_print(self.data.spin_id, self.opt_params.verbosity) 804 805 # Preform optimisation. 806 results = self.optimise() 807 808 # Disassemble the results list. 809 param_vector, func, iter, fc, gc, hc, warning = results 810 811 processor.return_object(MF_result_command(processor, self.memo_id, param_vector, func, iter, fc, gc, hc, warning, completed=False))

812 813

814 - def store_data(self, data, opt_params):

815 """Store all the data required for model-free optimisation. 816 817 @param data: The data used to initialise the model-free target function class. 818 @type data: class instance 819 @param opt_params: The parameters and data required for optimisation using minfx. 820 @type opt_params: class instance 821 """ 822 823 # Store the data. 824 self.data = data 825 self.opt_params = opt_params

826 827 828

829 -class MF_grid_command(MF_minimise_command):

830 """Command class for the model-free grid search.""" 831

832 - def __init__(self):

833 """Initialise all the data.""" 834 835 # Execute the base class __init__() method. 836 super(MF_grid_command, self).__init__()

837 838

839 - def optimise(self):

840 """Model-free grid search. 841 842 @return: The optimisation results consisting of the parameter vector, function value, iteration count, function count, gradient count, Hessian count, and warnings. 843 @rtype: tuple of numpy array, float, int, int, int, int, str 844 """ 845 846 # Normal grid search. 847 if not hasattr(self.opt_params, 'subdivision'): 848 results = grid(func=self.mf.func, args=(), num_incs=self.opt_params.inc, lower=self.opt_params.lower, upper=self.opt_params.upper, A=self.opt_params.A, b=self.opt_params.b, verbosity=self.opt_params.verbosity) 849 850 # Subdivided grid. 851 else: 852 results = grid_point_array(func=self.mf.func, args=(), points=self.opt_params.subdivision, verbosity=self.opt_params.verbosity) 853 854 # Unpack the results. 855 param_vector, func, iter, warning = results 856 fc = iter 857 gc = 0.0 858 hc = 0.0 859 860 # Return everything. 861 return param_vector, func, iter, fc, gc, hc, warning

862 863 864

865 -class MF_result_command(Result_command):

866 """Class for processing the model-free results.""" 867

868 - def __init__(self, processor, memo_id, param_vector, func, iter, fc, gc, hc, warning, completed):

869 """Set up the class, placing the minimisation results here.""" 870 871 # Execute the base class __init__() method. 872 super(MF_result_command, self).__init__(processor=processor, completed=completed) 873 874 # Store the arguments. 875 self.memo_id = memo_id 876 self.param_vector = param_vector 877 self.func = func 878 self.iter = iter 879 self.fc = fc 880 self.gc = gc 881 self.hc = hc 882 self.warning = warning

883 884

885 - def run(self, processor, memo):

886 """Disassemble the model-free optimisation results. 887 888 @param processor: Unused! 889 @type processor: None 890 @param memo: The model-free memo. 891 @type memo: memo 892 """ 893 894 # Disassemble the results. 895 disassemble_result(param_vector=self.param_vector, func=self.func, iter=self.iter, fc=self.fc, gc=self.gc, hc=self.hc, warning=self.warning, spin=memo.spin, sim_index=memo.sim_index, model_type=memo.model_type, scaling=memo.scaling, scaling_matrix=memo.scaling_matrix)

896

Source Code for Module specific_analyses.model_free.optimisation