specific_fns.model_free.mf

53 """Class containing functions specific to model-free optimisation.""" 54

55 - def _disassemble_param_vector(self, model_type, param_vector=None, spin=None, spin_id=None, sim_index=None):

56 """Disassemble the model-free parameter vector. 57 58 @param model_type: The model-free model type. This must be one of 'mf', 'local_tm', 59 'diff', or 'all'. 60 @type model_type: str 61 @keyword param_vector: The model-free parameter vector. 62 @type param_vector: numpy array 63 @keyword spin: The spin data container. If this argument is supplied, then the spin_id 64 argument will be ignored. 65 @type spin: SpinContainer instance 66 @keyword spin_id: The spin identification string. 67 @type spin_id: str 68 @keyword sim_index: The optional MC simulation index. 69 @type sim_index: int 70 """ 71 72 # Initialise. 73 param_index = 0 74 75 # Diffusion tensor parameters of the Monte Carlo simulations. 76 if sim_index != None and (model_type == 'diff' or model_type == 'all'): 77 # Spherical diffusion. 78 if cdp.diff_tensor.type == 'sphere': 79 # Sim values. 80 cdp.diff_tensor.tm_sim[sim_index] = param_vector[0] 81 82 # Parameter index. 83 param_index = param_index + 1 84 85 # Spheroidal diffusion. 86 elif cdp.diff_tensor.type == 'spheroid': 87 # Sim values. 88 cdp.diff_tensor.tm_sim[sim_index] = param_vector[0] 89 cdp.diff_tensor.Da_sim[sim_index] = param_vector[1] 90 cdp.diff_tensor.theta_sim[sim_index] = param_vector[2] 91 cdp.diff_tensor.phi_sim[sim_index] = param_vector[3] 92 diffusion_tensor.fold_angles(sim_index=sim_index) 93 94 # Parameter index. 95 param_index = param_index + 4 96 97 # Ellipsoidal diffusion. 98 elif cdp.diff_tensor.type == 'ellipsoid': 99 # Sim values. 100 cdp.diff_tensor.tm_sim[sim_index] = param_vector[0] 101 cdp.diff_tensor.Da_sim[sim_index] = param_vector[1] 102 cdp.diff_tensor.Dr_sim[sim_index] = param_vector[2] 103 cdp.diff_tensor.alpha_sim[sim_index] = param_vector[3] 104 cdp.diff_tensor.beta_sim[sim_index] = param_vector[4] 105 cdp.diff_tensor.gamma_sim[sim_index] = param_vector[5] 106 diffusion_tensor.fold_angles(sim_index=sim_index) 107 108 # Parameter index. 109 param_index = param_index + 6 110 111 # Diffusion tensor parameters. 112 elif model_type == 'diff' or model_type == 'all': 113 # Spherical diffusion. 114 if cdp.diff_tensor.type == 'sphere': 115 # Values. 116 cdp.diff_tensor.tm = param_vector[0] 117 118 # Parameter index. 119 param_index = param_index + 1 120 121 # Spheroidal diffusion. 122 elif cdp.diff_tensor.type == 'spheroid': 123 # Values. 124 cdp.diff_tensor.tm = param_vector[0] 125 cdp.diff_tensor.Da = param_vector[1] 126 cdp.diff_tensor.theta = param_vector[2] 127 cdp.diff_tensor.phi = param_vector[3] 128 diffusion_tensor.fold_angles() 129 130 # Parameter index. 131 param_index = param_index + 4 132 133 # Ellipsoidal diffusion. 134 elif cdp.diff_tensor.type == 'ellipsoid': 135 # Values. 136 cdp.diff_tensor.tm = param_vector[0] 137 cdp.diff_tensor.Da = param_vector[1] 138 cdp.diff_tensor.Dr = param_vector[2] 139 cdp.diff_tensor.alpha = param_vector[3] 140 cdp.diff_tensor.beta = param_vector[4] 141 cdp.diff_tensor.gamma = param_vector[5] 142 diffusion_tensor.fold_angles() 143 144 # Parameter index. 145 param_index = param_index + 6 146 147 # Model-free parameters. 148 if model_type != 'diff': 149 # The loop. 150 if spin: 151 loop = [spin] 152 else: 153 loop = spin_loop(spin_id) 154 155 # Loop over the spins. 156 for spin in loop: 157 # Skip deselected spins. 158 if not spin.select: 159 continue 160 161 # Loop over the model-free parameters. 162 for j in xrange(len(spin.params)): 163 # Local tm. 164 if spin.params[j] == 'local_tm': 165 if sim_index == None: 166 spin.local_tm = param_vector[param_index] 167 else: 168 spin.local_tm_sim[sim_index] = param_vector[param_index] 169 170 # S2. 171 elif spin.params[j] == 's2': 172 if sim_index == None: 173 spin.s2 = param_vector[param_index] 174 else: 175 spin.s2_sim[sim_index] = param_vector[param_index] 176 177 # S2f. 178 elif spin.params[j] == 's2f': 179 if sim_index == None: 180 spin.s2f = param_vector[param_index] 181 else: 182 spin.s2f_sim[sim_index] = param_vector[param_index] 183 184 # S2s. 185 elif spin.params[j] == 's2s': 186 if sim_index == None: 187 spin.s2s = param_vector[param_index] 188 else: 189 spin.s2s_sim[sim_index] = param_vector[param_index] 190 191 # te. 192 elif spin.params[j] == 'te': 193 if sim_index == None: 194 spin.te = param_vector[param_index] 195 else: 196 spin.te_sim[sim_index] = param_vector[param_index] 197 198 # tf. 199 elif spin.params[j] == 'tf': 200 if sim_index == None: 201 spin.tf = param_vector[param_index] 202 else: 203 spin.tf_sim[sim_index] = param_vector[param_index] 204 205 # ts. 206 elif spin.params[j] == 'ts': 207 if sim_index == None: 208 spin.ts = param_vector[param_index] 209 else: 210 spin.ts_sim[sim_index] = param_vector[param_index] 211 212 # Rex. 213 elif spin.params[j] == 'rex': 214 if sim_index == None: 215 spin.rex = param_vector[param_index] 216 else: 217 spin.rex_sim[sim_index] = param_vector[param_index] 218 219 # r. 220 elif spin.params[j] == 'r': 221 if sim_index == None: 222 spin.r = param_vector[param_index] 223 else: 224 spin.r_sim[sim_index] = param_vector[param_index] 225 226 # CSA. 227 elif spin.params[j] == 'csa': 228 if sim_index == None: 229 spin.csa = param_vector[param_index] 230 else: 231 spin.csa_sim[sim_index] = param_vector[param_index] 232 233 # Unknown parameter. 234 else: 235 raise RelaxError("Unknown parameter.") 236 237 # Increment the parameter index. 238 param_index = param_index + 1 239 240 # Calculate all order parameters after unpacking the vector. 241 if model_type != 'diff': 242 # The loop. 243 if spin: 244 loop = [spin] 245 else: 246 loop = spin_loop(spin_id) 247 248 # Loop over the spins. 249 for spin in loop: 250 # Skip deselected residues. 251 if not spin.select: 252 continue 253 254 # Normal values. 255 if sim_index == None: 256 # S2. 257 if 's2' not in spin.params and 's2f' in spin.params and 's2s' in spin.params: 258 spin.s2 = spin.s2f * spin.s2s 259 260 # S2f. 261 if 's2f' not in spin.params and 's2' in spin.params and 's2s' in spin.params: 262 if spin.s2s == 0.0: 263 spin.s2f = 1e99 264 else: 265 spin.s2f = spin.s2 / spin.s2s 266 267 # S2s. 268 if 's2s' not in spin.params and 's2' in spin.params and 's2f' in spin.params: 269 if spin.s2f == 0.0: 270 spin.s2s = 1e99 271 else: 272 spin.s2s = spin.s2 / spin.s2f 273 274 # Simulation values. 275 else: 276 # S2. 277 if 's2' not in spin.params and 's2f' in spin.params and 's2s' in spin.params: 278 spin.s2_sim[sim_index] = spin.s2f_sim[sim_index] * spin.s2s_sim[sim_index] 279 280 # S2f. 281 if 's2f' not in spin.params and 's2' in spin.params and 's2s' in spin.params: 282 if spin.s2s_sim[sim_index] == 0.0: 283 spin.s2f_sim[sim_index] = 1e99 284 else: 285 spin.s2f_sim[sim_index] = spin.s2_sim[sim_index] / spin.s2s_sim[sim_index] 286 287 # S2s. 288 if 's2s' not in spin.params and 's2' in spin.params and 's2f' in spin.params: 289 if spin.s2f_sim[sim_index] == 0.0: 290 spin.s2s_sim[sim_index] = 1e99 291 else: 292 spin.s2s_sim[sim_index] = spin.s2_sim[sim_index] / spin.s2f_sim[sim_index]

293 294

295 - def _disassemble_result(self, 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):

296 """Disassemble the optimisation results. 297 298 @keyword param_vector: The model-free parameter vector. 299 @type param_vector: numpy array 300 @keyword func: The optimised chi-squared value. 301 @type func: float 302 @keyword iter: The number of optimisation steps required to find the minimum. 303 @type iter: int 304 @keyword fc: The function count. 305 @type fc: int 306 @keyword gc: The gradient count. 307 @type gc: int 308 @keyword hc: The Hessian count. 309 @type hc: int 310 @keyword warning: Any optimisation warnings. 311 @type warning: str or None 312 @keyword spin: The spin container. 313 @type spin: SpinContainer instance or None 314 @keyword sim_index: The Monte Carlo simulation index. 315 @type sim_index: int or None 316 @keyword model_type: The model-free model type, one of 'mf', 'local_tm', 'diff', or 317 'all'. 318 @type model_type: str 319 @keyword scaling: If True, diagonal scaling is enabled during optimisation to 320 allow the problem to be better conditioned. 321 @type scaling: bool 322 @keyword scaling_matrix: The diagonal, square scaling matrix. 323 @type scaling_matrix: numpy diagonal matrix 324 """ 325 326 # Alias the current data pipe. 327 cdp = pipes.get_pipe() 328 329 # Catch infinite chi-squared values. 330 if isInf(func): 331 raise RelaxInfError('chi-squared') 332 333 # Catch chi-squared values of NaN. 334 if isNaN(func): 335 raise RelaxNaNError('chi-squared') 336 337 # Scaling. 338 if scaling: 339 param_vector = dot(scaling_matrix, param_vector) 340 341 # Check if the chi-squared value is lower. This allows for a parallelised grid search! 342 if sim_index == None: 343 # Get the correct value. 344 chi2 = None 345 if (model_type == 'mf' or model_type == 'local_tm') and hasattr(cdp, 'chi2'): 346 chi2 = spin.chi2 347 if (model_type == 'diff' or model_type == 'all') and hasattr(cdp, 'chi2'): 348 chi2 = cdp.chi2 349 350 # No improvement. 351 if chi2 != None and func >= chi2: 352 print("Discarding the optimisation results, the optimised chi-squared value is higher than the current value (%s >= %s)." % (func, chi2)) 353 354 # Exit! 355 return 356 357 # New minimum. 358 else: 359 print("Storing the optimisation results, the optimised chi-squared value is lower than the current value (%s < %s)." % (func, chi2)) 360 361 # Disassemble the parameter vector. 362 self._disassemble_param_vector(model_type, param_vector=param_vector, spin=spin, sim_index=sim_index) 363 364 # Monte Carlo minimisation statistics. 365 if sim_index != None: 366 # Sequence specific minimisation statistics. 367 if model_type == 'mf' or model_type == 'local_tm': 368 369 # Chi-squared statistic. 370 spin.chi2_sim[sim_index] = func 371 372 # Iterations. 373 spin.iter_sim[sim_index] = iter 374 375 # Function evaluations. 376 spin.f_count_sim[sim_index] = fc 377 378 # Gradient evaluations. 379 spin.g_count_sim[sim_index] = gc 380 381 # Hessian evaluations. 382 spin.h_count_sim[sim_index] = hc 383 384 # Warning. 385 spin.warning_sim[sim_index] = warning 386 387 # Global minimisation statistics. 388 elif model_type == 'diff' or model_type == 'all': 389 # Chi-squared statistic. 390 cdp.chi2_sim[sim_index] = func 391 392 # Iterations. 393 cdp.iter_sim[sim_index] = iter 394 395 # Function evaluations. 396 cdp.f_count_sim[sim_index] = fc 397 398 # Gradient evaluations. 399 cdp.g_count_sim[sim_index] = gc 400 401 # Hessian evaluations. 402 cdp.h_count_sim[sim_index] = hc 403 404 # Warning. 405 cdp.warning_sim[sim_index] = warning 406 407 # Normal statistics. 408 else: 409 # Sequence specific minimisation statistics. 410 if model_type == 'mf' or model_type == 'local_tm': 411 # Chi-squared statistic. 412 spin.chi2 = func 413 414 # Iterations. 415 spin.iter = iter 416 417 # Function evaluations. 418 spin.f_count = fc 419 420 # Gradient evaluations. 421 spin.g_count = gc 422 423 # Hessian evaluations. 424 spin.h_count = hc 425 426 # Warning. 427 spin.warning = warning 428 429 # Global minimisation statistics. 430 elif model_type == 'diff' or model_type == 'all': 431 # Chi-squared statistic. 432 cdp.chi2 = func 433 434 # Iterations. 435 cdp.iter = iter 436 437 # Function evaluations. 438 cdp.f_count = fc 439 440 # Gradient evaluations. 441 cdp.g_count = gc 442 443 # Hessian evaluations. 444 cdp.h_count = hc 445 446 # Warning. 447 cdp.warning = warning

448 449

450 - def _grid_search_config(self, num_params, spin=None, spin_id=None, lower=None, upper=None, inc=None, scaling_matrix=None, verbosity=1):

451 """Configure the grid search. 452 453 @param num_params: The number of parameters in the model. 454 @type num_params: int 455 @keyword spin: The spin data container. 456 @type spin: SpinContainer instance 457 @keyword spin_id: The spin identification string. 458 @type spin_id: str 459 @keyword lower: The lower bounds of the grid search which must be equal to the 460 number of parameters in the model. 461 @type lower: array of numbers 462 @keyword upper: The upper bounds of the grid search which must be equal to the 463 number of parameters in the model. 464 @type upper: array of numbers 465 @keyword inc: The increments for each dimension of the space for the grid 466 search. The number of elements in the array must equal to the 467 number of parameters in the model. 468 @type inc: array of int 469 @keyword scaling_matrix: The diagonal, square scaling matrix. 470 @type scaling_matrix: numpy diagonal matrix 471 @keyword verbosity: A flag specifying the amount of information to print. The 472 higher the value, the greater the verbosity. 473 @type verbosity: int 474 """ 475 476 # Test the grid search options. 477 self.test_grid_ops(lower=lower, upper=upper, inc=inc, n=num_params) 478 479 # If inc is a single int, convert it into an array of that value. 480 if isinstance(inc, int): 481 inc = [inc]*num_params 482 483 # Set up the default bounds. 484 if not lower: 485 # Init. 486 lower = [] 487 upper = [] 488 489 # Determine the model type. 490 model_type = self._determine_model_type() 491 492 # Minimisation options for diffusion tensor parameters. 493 if model_type == 'diff' or model_type == 'all': 494 # Get the diffusion tensor specific configuration. 495 self._grid_search_diff_bounds(lower, upper) 496 497 # Model-free parameters (residue specific parameters). 498 if model_type != 'diff': 499 # The loop. 500 if spin: 501 loop = [spin] 502 else: 503 loop = spin_loop(spin_id) 504 505 # Loop over the spins. 506 for spin in loop: 507 # Skip deselected residues. 508 if not spin.select: 509 continue 510 511 # Get the spin specific configuration. 512 self._grid_search_spin_bounds(spin, lower, upper) 513 514 # Diagonal scaling of minimisation options. 515 lower_new = [] 516 upper_new = [] 517 for i in xrange(num_params): 518 lower_new.append(lower[i] / scaling_matrix[i, i]) 519 upper_new.append(upper[i] / scaling_matrix[i, i]) 520 521 # Return the minimisation options. 522 return inc, lower_new, upper_new

523 524

525 - def _grid_search_diff_bounds(self, lower, upper):

526 """Set up the default grid search bounds the diffusion tensor. 527 528 This method appends the default bounds to the lower and upper lists. 529 530 @param lower: The lower bound list to append to. 531 @type lower: list 532 @param upper: The upper bound list to append to. 533 @type upper: list 534 """ 535 536 # Spherical diffusion {tm}. 537 if cdp.diff_tensor.type == 'sphere': 538 lower.append(1.0 * 1e-9) 539 upper.append(12.0 * 1e-9) 540 541 # Spheroidal diffusion {tm, Da, theta, phi}. 542 if cdp.diff_tensor.type == 'spheroid': 543 # tm. 544 lower.append(1.0 * 1e-9) 545 upper.append(12.0 * 1e-9) 546 547 # Da. 548 if cdp.diff_tensor.spheroid_type == 'prolate': 549 lower.append(0.0) 550 upper.append(1e7) 551 elif cdp.diff_tensor.spheroid_type == 'oblate': 552 lower.append(-1e7) 553 upper.append(0.0) 554 else: 555 lower.append(-1e7) 556 upper.append(1e7) 557 558 # theta. 559 lower.append(0.0) 560 upper.append(pi) 561 562 # phi. 563 lower.append(0.0) 564 upper.append(pi) 565 566 # Ellipsoidal diffusion {tm, Da, Dr, alpha, beta, gamma}. 567 elif cdp.diff_tensor.type == 'ellipsoid': 568 # tm. 569 lower.append(1.0 * 1e-9) 570 upper.append(12.0 * 1e-9) 571 572 # Da. 573 lower.append(0.0) 574 upper.append(1e7) 575 576 # Dr. 577 lower.append(0.0) 578 upper.append(1.0) 579 580 # alpha. 581 lower.append(0.0) 582 upper.append(pi) 583 584 # beta. 585 lower.append(0.0) 586 upper.append(pi) 587 588 # gamma. 589 lower.append(0.0) 590 upper.append(pi)

591 592

593 - def _grid_search_spin_bounds(self, spin, lower, upper):

594 """Set up the default grid search bounds for a single spin. 595 596 This method appends the default bounds to the lower and upper lists. The ordering of the 597 lists in min_options matches that of the params list in the spin container. 598 599 @param spin: A SpinContainer object. 600 @type spin: class instance 601 @param lower: The lower bound list to append to. 602 @type lower: list 603 @param upper: The upper bound list to append to. 604 @type upper: list 605 """ 606 607 # Loop over the model-free parameters. 608 for i in xrange(len(spin.params)): 609 # Local tm. 610 if spin.params[i] == 'local_tm': 611 lower.append(1.0 * 1e-9) 612 upper.append(12.0 * 1e-9) 613 614 # {S2, S2f, S2s}. 615 elif match('s2', spin.params[i]): 616 lower.append(0.0) 617 upper.append(1.0) 618 619 # {te, tf, ts}. 620 elif match('t', spin.params[i]): 621 lower.append(0.0) 622 upper.append(500.0 * 1e-12) 623 624 # Rex. 625 elif spin.params[i] == 'rex': 626 lower.append(0.0) 627 upper.append(5.0 / (2.0 * pi * cdp.frq[cdp.ri_ids[0]])**2) 628 629 # Bond length. 630 elif spin.params[i] == 'r': 631 lower.append(1.0 * 1e-10) 632 upper.append(1.05 * 1e-10) 633 634 # CSA. 635 elif spin.params[i] == 'csa': 636 lower.append(-120 * 1e-6) 637 upper.append(-200 * 1e-6) 638 639 # Unknown option. 640 else: 641 raise RelaxError("Unknown model-free parameter.")

642 643

644 - def _linear_constraints(self, num_params, model_type=None, spin=None, spin_id=None, scaling_matrix=None):

645 """Set up the model-free linear constraint matrices A and b. 646 647 Standard notation 648 ================= 649 650 The order parameter constraints are:: 651 652 0 <= S2 <= 1 653 0 <= S2f <= 1 654 0 <= S2s <= 1 655 656 By substituting the formula S2 = S2f.S2s into the above inequalities, the additional two 657 inequalities can be derived:: 658 659 S2 <= S2f 660 S2 <= S2s 661 662 Correlation time constraints are:: 663 664 te >= 0 665 tf >= 0 666 ts >= 0 667 668 tf <= ts 669 670 te, tf, ts <= 2 * tm 671 672 Additional constraints used include:: 673 674 Rex >= 0 675 0.9e-10 <= r <= 2e-10 676 -300e-6 <= CSA <= 0 677 678 679 Rearranged notation 680 =================== 681 682 The above inequality constraints can be rearranged into:: 683 684 S2 >= 0 685 -S2 >= -1 686 S2f >= 0 687 -S2f >= -1 688 S2s >= 0 689 -S2s >= -1 690 S2f - S2 >= 0 691 S2s - S2 >= 0 692 te >= 0 693 tf >= 0 694 ts >= 0 695 ts - tf >= 0 696 Rex >= 0 697 r >= 0.9e-10 698 -r >= -2e-10 699 CSA >= -300e-6 700 -CSA >= 0 701 702 703 Matrix notation 704 =============== 705 706 In the notation A.x >= b, where A is an matrix of coefficients, x is an array of parameter 707 values, and b is a vector of scalars, these inequality constraints are:: 708 709 | 1 0 0 0 0 0 0 0 0 | | 0 | 710 | | | | 711 |-1 0 0 0 0 0 0 0 0 | | -1 | 712 | | | | 713 | 0 1 0 0 0 0 0 0 0 | | 0 | 714 | | | | 715 | 0 -1 0 0 0 0 0 0 0 | | -1 | 716 | | | | 717 | 0 0 1 0 0 0 0 0 0 | | S2 | | 0 | 718 | | | | | | 719 | 0 0 -1 0 0 0 0 0 0 | | S2f | | -1 | 720 | | | | | | 721 |-1 1 0 0 0 0 0 0 0 | | S2s | | 0 | 722 | | | | | | 723 |-1 0 1 0 0 0 0 0 0 | | te | | 0 | 724 | | | | | | 725 | 0 0 0 1 0 0 0 0 0 | . | tf | >= | 0 | 726 | | | | | | 727 | 0 0 0 0 1 0 0 0 0 | | ts | | 0 | 728 | | | | | | 729 | 0 0 0 0 0 1 0 0 0 | | Rex | | 0 | 730 | | | | | | 731 | 0 0 0 0 -1 1 0 0 0 | | r | | 0 | 732 | | | | | | 733 | 0 0 0 0 0 0 1 0 0 | | CSA | | 0 | 734 | | | | 735 | 0 0 0 0 0 0 0 1 0 | | 0.9e-10 | 736 | | | | 737 | 0 0 0 0 0 0 0 -1 0 | | -2e-10 | 738 | | | | 739 | 0 0 0 0 0 0 0 0 1 | | -300e-6 | 740 | | | | 741 | 0 0 0 0 0 0 0 0 -1 | | 0 | 742 743 744 @param num_params: The number of parameters in the model. 745 @type num_params: int 746 @keyword model_type: The model type, one of 'all', 'diff', 'mf', or 'local_tm'. 747 @type model_type: str 748 @keyword spin: The spin data container. If this argument is supplied, then the 749 spin_id argument will be ignored. 750 @type spin: SpinContainer instance 751 @keyword spin_id: The spin identification string. 752 @type spin_id: str 753 @keyword scaling_matrix: The diagonal, square scaling matrix. 754 @type scaling_matrix: numpy diagonal matrix 755 """ 756 757 # Upper limit flag for correlation times. 758 upper_time_limit = 1 759 760 # Initialisation (0..j..m). 761 A = [] 762 b = [] 763 zero_array = zeros(num_params, float64) 764 i = 0 765 j = 0 766 767 # Diffusion tensor parameters. 768 if model_type != 'mf' and diffusion_tensor.diff_data_exists(): 769 # Spherical diffusion. 770 if cdp.diff_tensor.type == 'sphere': 771 # 0 <= tm <= 200 ns. 772 A.append(zero_array * 0.0) 773 A.append(zero_array * 0.0) 774 A[j][i] = 1.0 775 A[j+1][i] = -1.0 776 b.append(0.0 / scaling_matrix[i, i]) 777 b.append(-200.0 * 1e-9 / scaling_matrix[i, i]) 778 i = i + 1 779 j = j + 2 780 781 # Spheroidal diffusion. 782 elif cdp.diff_tensor.type == 'spheroid': 783 # 0 <= tm <= 200 ns. 784 A.append(zero_array * 0.0) 785 A.append(zero_array * 0.0) 786 A[j][i] = 1.0 787 A[j+1][i] = -1.0 788 b.append(0.0 / scaling_matrix[i, i]) 789 b.append(-200.0 * 1e-9 / scaling_matrix[i, i]) 790 i = i + 1 791 j = j + 2 792 793 # Prolate diffusion, Da >= 0. 794 if cdp.diff_tensor.spheroid_type == 'prolate': 795 A.append(zero_array * 0.0) 796 A[j][i] = 1.0 797 b.append(0.0 / scaling_matrix[i, i]) 798 i = i + 1 799 j = j + 1 800 801 # Add two to i for the theta and phi parameters. 802 i = i + 2 803 804 # Oblate diffusion, Da <= 0. 805 elif cdp.diff_tensor.spheroid_type == 'oblate': 806 A.append(zero_array * 0.0) 807 A[j][i] = -1.0 808 b.append(0.0 / scaling_matrix[i, i]) 809 i = i + 1 810 j = j + 1 811 812 # Add two to i for the theta and phi parameters. 813 i = i + 2 814 815 else: 816 # Add three to i for the Da, theta and phi parameters. 817 i = i + 3 818 819 # Ellipsoidal diffusion. 820 elif cdp.diff_tensor.type == 'ellipsoid': 821 # 0 <= tm <= 200 ns. 822 A.append(zero_array * 0.0) 823 A.append(zero_array * 0.0) 824 A[j][i] = 1.0 825 A[j+1][i] = -1.0 826 b.append(0.0 / scaling_matrix[i, i]) 827 b.append(-200.0 * 1e-9 / scaling_matrix[i, i]) 828 i = i + 1 829 j = j + 2 830 831 # Da >= 0. 832 A.append(zero_array * 0.0) 833 A[j][i] = 1.0 834 b.append(0.0 / scaling_matrix[i, i]) 835 i = i + 1 836 j = j + 1 837 838 # 0 <= Dr <= 1. 839 A.append(zero_array * 0.0) 840 A.append(zero_array * 0.0) 841 A[j][i] = 1.0 842 A[j+1][i] = -1.0 843 b.append(0.0 / scaling_matrix[i, i]) 844 b.append(-1.0 / scaling_matrix[i, i]) 845 i = i + 1 846 j = j + 2 847 848 # Add three to i for the alpha, beta, and gamma parameters. 849 i = i + 3 850 851 # Model-free parameters. 852 if model_type != 'diff': 853 # The loop. 854 if spin: 855 loop = [spin] 856 else: 857 loop = spin_loop(spin_id) 858 859 # Loop over the spins. 860 for spin in loop: 861 # Skip deselected spins. 862 if not spin.select: 863 continue 864 865 # Save current value of i. 866 old_i = i 867 868 # Loop over the model-free parameters. 869 for l in xrange(len(spin.params)): 870 # Local tm. 871 if spin.params[l] == 'local_tm': 872 if upper_time_limit: 873 # 0 <= tm <= 200 ns. 874 A.append(zero_array * 0.0) 875 A.append(zero_array * 0.0) 876 A[j][i] = 1.0 877 A[j+1][i] = -1.0 878 b.append(0.0 / scaling_matrix[i, i]) 879 b.append(-200.0 * 1e-9 / scaling_matrix[i, i]) 880 j = j + 2 881 else: 882 # 0 <= tm. 883 A.append(zero_array * 0.0) 884 A[j][i] = 1.0 885 b.append(0.0 / scaling_matrix[i, i]) 886 j = j + 1 887 888 # Order parameters {S2, S2f, S2s}. 889 elif match('s2', spin.params[l]): 890 # 0 <= S2 <= 1. 891 A.append(zero_array * 0.0) 892 A.append(zero_array * 0.0) 893 A[j][i] = 1.0 894 A[j+1][i] = -1.0 895 b.append(0.0 / scaling_matrix[i, i]) 896 b.append(-1.0 / scaling_matrix[i, i]) 897 j = j + 2 898 899 # S2 <= S2f and S2 <= S2s. 900 if spin.params[l] == 's2': 901 for m in xrange(len(spin.params)): 902 if spin.params[m] == 's2f' or spin.params[m] == 's2s': 903 A.append(zero_array * 0.0) 904 A[j][i] = -1.0 905 A[j][old_i+m] = 1.0 906 b.append(0.0) 907 j = j + 1 908 909 # Correlation times {te, tf, ts}. 910 elif match('t[efs]', spin.params[l]): 911 # te, tf, ts >= 0. 912 A.append(zero_array * 0.0) 913 A[j][i] = 1.0 914 b.append(0.0 / scaling_matrix[i, i]) 915 j = j + 1 916 917 # tf <= ts. 918 if spin.params[l] == 'ts': 919 for m in xrange(len(spin.params)): 920 if spin.params[m] == 'tf': 921 A.append(zero_array * 0.0) 922 A[j][i] = 1.0 923 A[j][old_i+m] = -1.0 924 b.append(0.0) 925 j = j + 1 926 927 # te, tf, ts <= 2 * tm. (tf not needed because tf <= ts). 928 if upper_time_limit: 929 if not spin.params[l] == 'tf': 930 if model_type == 'mf': 931 A.append(zero_array * 0.0) 932 A[j][i] = -1.0 933 b.append(-2.0 * cdp.diff_tensor.tm / scaling_matrix[i, i]) 934 else: 935 A.append(zero_array * 0.0) 936 A[j][0] = 2.0 937 A[j][i] = -1.0 938 b.append(0.0) 939 940 j = j + 1 941 942 # Rex. 943 elif spin.params[l] == 'rex': 944 A.append(zero_array * 0.0) 945 A[j][i] = 1.0 946 b.append(0.0 / scaling_matrix[i, i]) 947 j = j + 1 948 949 # Bond length. 950 elif spin.params[l] == 'r': 951 # 0.9e-10 <= r <= 2e-10. 952 A.append(zero_array * 0.0) 953 A.append(zero_array * 0.0) 954 A[j][i] = 1.0 955 A[j+1][i] = -1.0 956 b.append(0.9e-10 / scaling_matrix[i, i]) 957 b.append(-2e-10 / scaling_matrix[i, i]) 958 j = j + 2 959 960 # CSA. 961 elif spin.params[l] == 'csa': 962 # -300e-6 <= CSA <= 0. 963 A.append(zero_array * 0.0) 964 A.append(zero_array * 0.0) 965 A[j][i] = 1.0 966 A[j+1][i] = -1.0 967 b.append(-300e-6 / scaling_matrix[i, i]) 968 b.append(0.0 / scaling_matrix[i, i]) 969 j = j + 2 970 971 # Increment i. 972 i = i + 1 973 974 # Convert to numpy data structures. 975 A = array(A, int8) 976 b = array(b, float64) 977 978 return A, b

979 980

981 - def _minimise_data_setup(self, data_store, min_algor, num_data_sets, min_options, spin=None, sim_index=None):

982 """Set up all the data required for minimisation. 983 984 @param data_store: A data storage container. 985 @type data_store: class instance 986 @param min_algor: The minimisation algorithm to use. 987 @type min_algor: str 988 @param num_data_sets: The number of data sets. 989 @type num_data_sets: int 990 @param min_options: The minimisation options array. 991 @type min_options: list 992 @keyword spin: The spin data container. 993 @type spin: SpinContainer instance 994 @keyword sim_index: The optional MC simulation index. 995 @type sim_index: int 996 @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) 997 @rtype: tuple 998 """ 999 1000 # Initialise the data structures for the model-free function. 1001 data_store.ri_data = [] 1002 data_store.ri_data_err = [] 1003 data_store.equations = [] 1004 data_store.param_types = [] 1005 data_store.param_values = None 1006 data_store.r = [] 1007 data_store.csa = [] 1008 data_store.num_frq = [] 1009 data_store.frq = [] 1010 data_store.num_ri = [] 1011 data_store.remap_table = [] 1012 data_store.noe_r1_table = [] 1013 data_store.ri_types = [] 1014 data_store.gx = [] 1015 data_store.gh = [] 1016 data_store.num_params = [] 1017 data_store.xh_unit_vectors = [] 1018 if data_store.model_type == 'local_tm': 1019 data_store.mf_params = [] 1020 elif data_store.model_type == 'diff': 1021 data_store.param_values = [] 1022 1023 # Set up the data for the back_calc function. 1024 if min_algor == 'back_calc': 1025 # The data. 1026 data_store.ri_data = [0.0] 1027 data_store.ri_data_err = [0.000001] 1028 data_store.equations = [spin.equation] 1029 data_store.param_types = [spin.params] 1030 data_store.r = [spin.r] 1031 data_store.csa = [spin.csa] 1032 data_store.num_frq = [1] 1033 data_store.frq = [[min_options[3]]] 1034 data_store.num_ri = [1] 1035 data_store.remap_table = [[0]] 1036 data_store.noe_r1_table = [[None]] 1037 data_store.ri_types = [[min_options[2]]] 1038 data_store.gx = [return_gyromagnetic_ratio(spin.heteronuc_type)] 1039 data_store.gh = [return_gyromagnetic_ratio(spin.proton_type)] 1040 if data_store.model_type != 'local_tm' and cdp.diff_tensor.type != 'sphere': 1041 data_store.xh_unit_vectors = [spin.xh_vect] 1042 else: 1043 data_store.xh_unit_vectors = [None] 1044 1045 # Count the number of model-free parameters for the spin index. 1046 data_store.num_params = [len(spin.params)] 1047 1048 # Loop over the number of data sets. 1049 for j in xrange(num_data_sets): 1050 # Set the spin index and get the spin, if not already set. 1051 if data_store.model_type == 'diff' or data_store.model_type == 'all': 1052 spin_index = j 1053 spin, data_store.spin_id = return_spin_from_index(global_index=spin_index, return_spin_id=True) 1054 1055 # Skip deselected spins. 1056 if not spin.select: 1057 continue 1058 1059 # Skip spins where there is no data or errors. 1060 if not hasattr(spin, 'ri_data') or not hasattr(spin, 'ri_data_err'): 1061 continue 1062 1063 # Make sure that the errors are strictly positive numbers. 1064 for ri_id in cdp.ri_ids: 1065 # Alias. 1066 err = spin.ri_data_err[ri_id] 1067 1068 # Checks. 1069 if err != None and err == 0.0: 1070 raise RelaxError("Zero error for spin '%s' for the relaxation data ID '%s', minimisation not possible." % (errid)) 1071 elif err != None and err < 0.0: 1072 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)) 1073 1074 # The relaxation data optimisation structures. 1075 data = self._relax_data_opt_structs(spin, sim_index=sim_index) 1076 1077 # Append the data. 1078 data_store.ri_data.append(data[0]) 1079 data_store.ri_data_err.append(data[1]) 1080 data_store.num_frq.append(data[2]) 1081 data_store.num_ri.append(data[3]) 1082 data_store.ri_types.append(data[4]) 1083 data_store.frq.append(data[5]) 1084 data_store.remap_table.append(data[6]) 1085 data_store.noe_r1_table.append(data[7]) 1086 1087 # Repackage the data. 1088 data_store.equations.append(spin.equation) 1089 data_store.param_types.append(spin.params) 1090 data_store.gx.append(return_gyromagnetic_ratio(spin.heteronuc_type)) 1091 data_store.gh.append(return_gyromagnetic_ratio(spin.proton_type)) 1092 if sim_index == None or data_store.model_type == 'diff': 1093 data_store.r.append(spin.r) 1094 data_store.csa.append(spin.csa) 1095 else: 1096 data_store.r.append(spin.r_sim[sim_index]) 1097 data_store.csa.append(spin.csa_sim[sim_index]) 1098 1099 # Model-free parameter values. 1100 if data_store.model_type == 'local_tm': 1101 pass 1102 1103 # Vectors. 1104 if data_store.model_type != 'local_tm' and cdp.diff_tensor.type != 'sphere': 1105 data_store.xh_unit_vectors.append(spin.xh_vect) 1106 else: 1107 data_store.xh_unit_vectors.append(None) 1108 1109 # Count the number of model-free parameters for the spin index. 1110 data_store.num_params.append(len(spin.params)) 1111 1112 # Repackage the parameter values for minimising just the diffusion tensor parameters. 1113 if data_store.model_type == 'diff': 1114 data_store.param_values.append(self._assemble_param_vector(model_type='mf')) 1115 1116 # Convert to numpy arrays. 1117 for k in xrange(len(data_store.ri_data)): 1118 data_store.ri_data[k] = array(data_store.ri_data[k], float64) 1119 data_store.ri_data_err[k] = array(data_store.ri_data_err[k], float64) 1120 1121 # Diffusion tensor type. 1122 if data_store.model_type == 'local_tm': 1123 data_store.diff_type = 'sphere' 1124 else: 1125 data_store.diff_type = cdp.diff_tensor.type 1126 1127 # Package the diffusion tensor parameters. 1128 data_store.diff_params = None 1129 if data_store.model_type == 'mf': 1130 # Spherical diffusion. 1131 if data_store.diff_type == 'sphere': 1132 data_store.diff_params = [cdp.diff_tensor.tm] 1133 1134 # Spheroidal diffusion. 1135 elif data_store.diff_type == 'spheroid': 1136 data_store.diff_params = [cdp.diff_tensor.tm, cdp.diff_tensor.Da, cdp.diff_tensor.theta, cdp.diff_tensor.phi] 1137 1138 # Ellipsoidal diffusion. 1139 elif data_store.diff_type == 'ellipsoid': 1140 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] 1141 elif min_algor == 'back_calc' and data_store.model_type == 'local_tm': 1142 # Spherical diffusion. 1143 data_store.diff_params = [spin.local_tm]

1144 1145

1146 - def _relax_data_opt_structs(self, spin, sim_index=None):

1147 """Package the relaxation data into the data structures used for optimisation. 1148 1149 @param spin: The spin container to extract the data from. 1150 @type spin: SpinContainer instance 1151 @keyword sim_index: The optional MC simulation index. 1152 @type sim_index: int 1153 @return: The structures ri_data, ri_data_err, num_frq, num_ri, ri_ids, frq, remap_table, noe_r1_table. 1154 @rtype: tuple 1155 """ 1156 1157 # Initialise the data. 1158 ri_data = [] 1159 ri_data_err = [] 1160 ri_labels = [] 1161 frq = [] 1162 remap_table = [] 1163 noe_r1_table = [] 1164 1165 # Loop over the relaxation data. 1166 for ri_id in cdp.ri_ids: 1167 # The Rx data. 1168 if sim_index == None: 1169 data = spin.ri_data[ri_id] 1170 else: 1171 data = spin.ri_data_sim[ri_id][sim_index] 1172 1173 # The errors. 1174 err = spin.ri_data_err[ri_id] 1175 1176 # Missing data, so don't add it. 1177 if data == None or err == None: 1178 continue 1179 1180 # Append the data and error. 1181 ri_data.append(data) 1182 ri_data_err.append(err) 1183 1184 # The labels. 1185 ri_labels.append(cdp.ri_type[ri_id]) 1186 1187 # The frequencies. 1188 if cdp.frq[ri_id] not in frq: 1189 frq.append(cdp.frq[ri_id]) 1190 1191 # The remap table. 1192 remap_table.append(frq.index(cdp.frq[ri_id])) 1193 1194 # The NOE to R1 mapping table. 1195 noe_r1_table.append(None) 1196 1197 # The number of data sets. 1198 num_ri = len(ri_data) 1199 1200 # Fill the NOE to R1 mapping table. 1201 for i in range(num_ri): 1202 # If the data corresponds to 'NOE', try to find if the corresponding R1 data. 1203 if cdp.ri_type[cdp.ri_ids[i]] == 'NOE': 1204 for j in range(num_ri): 1205 if cdp.ri_type[cdp.ri_ids[j]] == 'R1' and cdp.frq[cdp.ri_ids[i]] == cdp.frq[cdp.ri_ids[j]]: 1206 noe_r1_table[i] = j 1207 1208 # Return the structures. 1209 return ri_data, ri_data_err, len(frq), num_ri, ri_labels, frq, remap_table, noe_r1_table

1210 1211

1212 - def _reset_min_stats(self):

1213 """Reset all the minimisation statistics. 1214 1215 All global and spin specific values will be set to None. 1216 """ 1217 1218 # Global stats. 1219 if hasattr(cdp, 'chi2'): 1220 cdp.chi2 = None 1221 cdp.iter = None 1222 cdp.f_count = None 1223 cdp.g_count = None 1224 cdp.h_count = None 1225 cdp.warning = None 1226 1227 # Spin specific stats. 1228 for spin in spin_loop(): 1229 if hasattr(spin, 'chi2'): 1230 spin.chi2 = None 1231 spin.iter = None 1232 spin.f_count = None 1233 spin.g_count = None 1234 spin.h_count = None 1235 spin.warning = None

1236 1237

1238 - def calculate(self, spin_id=None, verbosity=1, sim_index=None):

1239 """Calculation of the model-free chi-squared value. 1240 1241 @keyword spin_id: The spin identification string. 1242 @type spin_id: str 1243 @keyword verbosity: The amount of information to print. The higher the value, the greater the verbosity. 1244 @type verbosity: int 1245 @keyword sim_index: The optional MC simulation index. 1246 @type sim_index: int 1247 """ 1248 1249 # Test if sequence data is loaded. 1250 if not exists_mol_res_spin_data(): 1251 raise RelaxNoSequenceError 1252 1253 # Determine the model type. 1254 model_type = self._determine_model_type() 1255 1256 # Test if diffusion tensor data exists. 1257 if model_type != 'local_tm' and not diff_data_exists(): 1258 raise RelaxNoTensorError('diffusion') 1259 1260 # Test if the PDB file has been loaded. 1261 if model_type != 'local_tm' and cdp.diff_tensor.type != 'sphere' and not hasattr(cdp, 'structure'): 1262 raise RelaxNoPdbError 1263 1264 # Loop over the spins. 1265 for spin, id in spin_loop(spin_id, return_id=True): 1266 # Skip deselected spins. 1267 if not spin.select: 1268 continue 1269 1270 # Test if the model-free model has been setup. 1271 if not spin.model: 1272 raise RelaxNoModelError 1273 1274 # Test if unit vectors exist. 1275 if model_type != 'local_tm' and cdp.diff_tensor.type != 'sphere' and not hasattr(spin, 'xh_vect'): 1276 raise RelaxNoVectorsError 1277 1278 # Test if multiple unit vectors exist. 1279 if model_type != 'local_tm' and cdp.diff_tensor.type != 'sphere' and hasattr(spin, 'xh_vect') and arg_check.is_num_list(spin.xh_vect[0], raise_error=False): 1280 raise RelaxMultiVectorError 1281 1282 # Test if the spin type has been set. 1283 if not hasattr(spin, 'heteronuc_type'): 1284 raise RelaxSpinTypeError 1285 1286 # Test if the type attached proton has been set. 1287 if not hasattr(spin, 'proton_type'): 1288 raise RelaxProtonTypeError 1289 1290 # Test if the model-free parameter values exist. 1291 unset_param = self._are_mf_params_set(spin) 1292 if unset_param != None: 1293 raise RelaxNoValueError(unset_param) 1294 1295 # Test if the CSA value has been set. 1296 if not hasattr(spin, 'csa') or spin.csa == None: 1297 raise RelaxNoValueError("CSA") 1298 1299 # Test if the bond length value has been set. 1300 if not hasattr(spin, 'r') or spin.r == None: 1301 raise RelaxNoValueError("bond length") 1302 1303 # Skip spins where there is no data or errors. 1304 if not hasattr(spin, 'ri_data') or not hasattr(spin, 'ri_data_err'): 1305 continue 1306 1307 # Make sure that the errors are strictly positive numbers. 1308 for ri_id in cdp.ri_ids: 1309 # Alias. 1310 err = spin.ri_data_err[ri_id] 1311 1312 # Checks. 1313 if err != None and err == 0.0: 1314 raise RelaxError("Zero error for spin '%s' for the relaxation data ID '%s', minimisation not possible." % (id, ri_id)) 1315 elif err != None and err < 0.0: 1316 raise RelaxError("Negative error of %s for spin '%s' for the relaxation data ID '%s', minimisation not possible." % (err, id, ri_id)) 1317 1318 # Create the initial parameter vector. 1319 param_vector = self._assemble_param_vector(spin=spin, sim_index=sim_index) 1320 1321 # The relaxation data optimisation structures. 1322 data = self._relax_data_opt_structs(spin, sim_index=sim_index) 1323 1324 # Append the data. 1325 ri_data = [array(data[0])] 1326 ri_data_err = [array(data[1])] 1327 num_frq = [data[2]] 1328 num_ri = [data[3]] 1329 ri_labels = [data[4]] 1330 frq = [data[5]] 1331 remap_table = [data[6]] 1332 noe_r1_table = [data[7]] 1333 1334 # Repackage the spin. 1335 if sim_index == None: 1336 r = [spin.r] 1337 csa = [spin.csa] 1338 else: 1339 r = [spin.r_sim[sim_index]] 1340 csa = [spin.csa_sim[sim_index]] 1341 1342 # Vectors. 1343 if model_type != 'local_tm' and cdp.diff_tensor.type != 'sphere': 1344 xh_unit_vectors = [spin.xh_vect] 1345 else: 1346 xh_unit_vectors = [None] 1347 1348 # Gyromagnetic ratios. 1349 gx = [return_gyromagnetic_ratio(spin.heteronuc_type)] 1350 gh = [return_gyromagnetic_ratio(spin.proton_type)] 1351 1352 # Count the number of model-free parameters for the residue index. 1353 num_params = [len(spin.params)] 1354 1355 # Repackage the parameter values as a local model (ignore if the diffusion tensor is not fixed). 1356 param_values = [self._assemble_param_vector(model_type='mf')] 1357 1358 # Package the diffusion tensor parameters. 1359 if model_type == 'local_tm': 1360 diff_params = [spin.local_tm] 1361 diff_type = 'sphere' 1362 else: 1363 # Diff type. 1364 diff_type = cdp.diff_tensor.type 1365 1366 # Spherical diffusion. 1367 if diff_type == 'sphere': 1368 diff_params = [cdp.diff_tensor.tm] 1369 1370 # Spheroidal diffusion. 1371 elif diff_type == 'spheroid': 1372 diff_params = [cdp.diff_tensor.tm, cdp.diff_tensor.Da, cdp.diff_tensor.theta, cdp.diff_tensor.phi] 1373 1374 # Ellipsoidal diffusion. 1375 elif diff_type == 'ellipsoid': 1376 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] 1377 1378 # Initialise the model-free function. 1379 mf = Mf(init_params=param_vector, model_type='mf', diff_type=diff_type, diff_params=diff_params, num_spins=1, equations=[spin.equation], param_types=[spin.params], param_values=param_values, relax_data=ri_data, errors=ri_data_err, bond_length=r, csa=csa, num_frq=num_frq, frq=frq, num_ri=num_ri, remap_table=remap_table, noe_r1_table=noe_r1_table, ri_labels=ri_labels, gx=gx, gh=gh, h_bar=h_bar, mu0=mu0, num_params=num_params, vectors=xh_unit_vectors) 1380 1381 # Chi-squared calculation. 1382 try: 1383 chi2 = mf.func(param_vector) 1384 except OverflowError: 1385 chi2 = 1e200 1386 1387 # Global chi-squared value. 1388 if model_type == 'all' or model_type == 'diff': 1389 cdp.chi2 = cdp.chi2 + chi2 1390 else: 1391 spin.chi2 = chi2

1392 1393

1394 - def grid_search(self, lower=None, upper=None, inc=None, constraints=True, verbosity=1, sim_index=None):

1395 """The model-free grid search function. 1396 1397 @keyword lower: The lower bounds of the grid search which must be equal to the 1398 number of parameters in the model. 1399 @type lower: array of numbers 1400 @keyword upper: The upper bounds of the grid search which must be equal to the 1401 number of parameters in the model. 1402 @type upper: array of numbers 1403 @keyword inc: The increments for each dimension of the space for the grid search. 1404 The number of elements in the array must equal to the number of 1405 parameters in the model. 1406 @type inc: array of int 1407 @keyword constraints: If True, constraints are applied during the grid search (eliminating 1408 parts of the grid). If False, no constraints are used. 1409 @type constraints: bool 1410 @keyword verbosity: A flag specifying the amount of information to print. The higher 1411 the value, the greater the verbosity. 1412 @type verbosity: int 1413 @keyword sim_index: The index of the simulation to apply the grid search to. If None, 1414 the normal model is optimised. 1415 @type sim_index: int 1416 """ 1417 1418 # Minimisation. 1419 self.minimise(min_algor='grid', lower=lower, upper=upper, inc=inc, constraints=constraints, verbosity=verbosity, sim_index=sim_index)

1420 1421

1422 - def minimise(self, min_algor=None, min_options=None, func_tol=None, grad_tol=None, max_iterations=None, constraints=False, scaling=True, verbosity=0, sim_index=None, lower=None, upper=None, inc=None):

1423 """Model-free minimisation function. 1424 1425 Three categories of models exist for which the approach to minimisation is different. These 1426 are: 1427 1428 Single spin optimisations: The 'mf' and 'local_tm' model types which are the 1429 model-free parameters for single spins, optionally with a local tm parameter. These 1430 models have no global parameters. 1431 1432 Diffusion tensor optimisations: The 'diff' diffusion tensor model type. No spin 1433 specific parameters exist. 1434 1435 Optimisation of everything: The 'all' model type consisting of all model-free and all 1436 diffusion tensor parameters. 1437 1438 1439 @keyword min_algor: The minimisation algorithm to use. 1440 @type min_algor: str 1441 @keyword min_options: An array of options to be used by the minimisation algorithm. 1442 @type min_options: array of str 1443 @keyword func_tol: The function tolerance which, when reached, terminates optimisation. 1444 Setting this to None turns of the check. 1445 @type func_tol: None or float 1446 @keyword grad_tol: The gradient tolerance which, when reached, terminates optimisation. 1447 Setting this to None turns of the check. 1448 @type grad_tol: None or float 1449 @keyword max_iterations: The maximum number of iterations for the algorithm. 1450 @type max_iterations: int 1451 @keyword constraints: If True, constraints are used during optimisation. 1452 @type constraints: bool 1453 @keyword scaling: If True, diagonal scaling is enabled during optimisation to allow 1454 the problem to be better conditioned. 1455 @type scaling: bool 1456 @keyword verbosity: The amount of information to print. The higher the value, the 1457 greater the verbosity. 1458 @type verbosity: int 1459 @keyword sim_index: The index of the simulation to optimise. This should be None if 1460 normal optimisation is desired. 1461 @type sim_index: None or int 1462 @keyword lower: The lower bounds of the grid search which must be equal to the 1463 number of parameters in the model. This optional argument is only 1464 used when doing a grid search. 1465 @type lower: array of numbers 1466 @keyword upper: The upper bounds of the grid search which must be equal to the 1467 number of parameters in the model. This optional argument is only 1468 used when doing a grid search. 1469 @type upper: array of numbers 1470 @keyword inc: The increments for each dimension of the space for the grid search. 1471 The number of elements in the array must equal to the number of 1472 parameters in the model. This argument is only used when doing a 1473 grid search. 1474 @type inc: array of int 1475 """ 1476 1477 # Test if sequence data is loaded. 1478 if not exists_mol_res_spin_data(): 1479 raise RelaxNoSequenceError 1480 1481 # Test if the model-free model has been setup, and that the heteronucleus and attached proton type have been set. 1482 for spin in spin_loop(): 1483 # Skip deselected spins. 1484 if not spin.select: 1485 continue 1486 1487 # Not setup. 1488 if not spin.model: 1489 raise RelaxNoModelError 1490 1491 # Test if the spin type has been set. 1492 if not hasattr(spin, 'heteronuc_type'): 1493 raise RelaxSpinTypeError 1494 1495 # Test if the type attached proton has been set. 1496 if not hasattr(spin, 'proton_type'): 1497 raise RelaxProtonTypeError 1498 1499 # Reset the minimisation statistics. 1500 if sim_index == None and min_algor != 'back_calc': 1501 self._reset_min_stats() 1502 1503 # Containers for the model-free data and optimisation parameters. 1504 data_store = Data_container() 1505 opt_params = Data_container() 1506 1507 # Add the imported parameters. 1508 data_store.h_bar = h_bar 1509 data_store.mu0 = mu0 1510 opt_params.min_algor = min_algor 1511 opt_params.min_options = min_options 1512 opt_params.func_tol = func_tol 1513 opt_params.grad_tol = grad_tol 1514 opt_params.max_iterations = max_iterations 1515 1516 # Add the keyword args. 1517 opt_params.verbosity = verbosity 1518 1519 # Determine the model type. 1520 data_store.model_type = self._determine_model_type() 1521 if not data_store.model_type: 1522 return 1523 1524 # Model type for the back-calculate function. 1525 if min_algor == 'back_calc' and data_store.model_type != 'local_tm': 1526 data_store.model_type = 'mf' 1527 1528 # Test if diffusion tensor data exists. 1529 if data_store.model_type != 'local_tm' and not diffusion_tensor.diff_data_exists(): 1530 raise RelaxNoTensorError('diffusion') 1531 1532 # Tests for the PDB file and unit vectors. 1533 if data_store.model_type != 'local_tm' and cdp.diff_tensor.type != 'sphere': 1534 # Test if the structure file has been loaded. 1535 if not hasattr(cdp, 'structure'): 1536 raise RelaxNoPdbError 1537 1538 # Test if unit vectors exist. 1539 for spin in spin_loop(): 1540 # Skip deselected spins. 1541 if not spin.select: 1542 continue 1543 1544 # Unit vector. 1545 if not hasattr(spin, 'xh_vect'): 1546 raise RelaxNoVectorsError 1547 1548 # Test if the model-free parameter values are set for minimising diffusion tensor parameters by themselves. 1549 if data_store.model_type == 'diff': 1550 # Loop over the sequence. 1551 for spin in spin_loop(): 1552 unset_param = self._are_mf_params_set(spin) 1553 if unset_param != None: 1554 raise RelaxNoValueError(unset_param) 1555 1556 # Print out. 1557 if verbosity >= 1: 1558 if data_store.model_type == 'mf': 1559 print("Only the model-free parameters for single spins will be used.") 1560 elif data_store.model_type == 'local_mf': 1561 print("Only a local tm value together with the model-free parameters for single spins will be used.") 1562 elif data_store.model_type == 'diff': 1563 print("Only diffusion tensor parameters will be used.") 1564 elif data_store.model_type == 'all': 1565 print("The diffusion tensor parameters together with the model-free parameters for all spins will be used.") 1566 1567 # Test if the CSA and bond length values have been set. 1568 for spin in spin_loop(): 1569 # Skip deselected spins. 1570 if not spin.select: 1571 continue 1572 1573 # CSA value. 1574 if not hasattr(spin, 'csa') or spin.csa == None: 1575 raise RelaxNoValueError("CSA") 1576 1577 # Bond length value. 1578 if not hasattr(spin, 'r') or spin.r == None: 1579 raise RelaxNoValueError("bond length") 1580 1581 # Number of spins, minimisation instances, and data sets for each model type. 1582 if data_store.model_type == 'mf' or data_store.model_type == 'local_tm': 1583 num_instances = count_spins(skip_desel=False) 1584 num_data_sets = 1 1585 data_store.num_spins = 1 1586 elif data_store.model_type == 'diff' or data_store.model_type == 'all': 1587 num_instances = 1 1588 num_data_sets = count_spins(skip_desel=False) 1589 data_store.num_spins = count_spins() 1590 1591 # Number of spins, minimisation instances, and data sets for the back-calculate function. 1592 if min_algor == 'back_calc': 1593 num_instances = 1 1594 num_data_sets = 0 1595 data_store.num_spins = 1 1596 1597 # Get the Processor box singleton (it contains the Processor instance) and alias the Processor. 1598 processor_box = Processor_box() 1599 processor = processor_box.processor 1600 1601 # Loop over the minimisation instances. 1602 ####################################### 1603 1604 for i in xrange(num_instances): 1605 # Get the spin container if required. 1606 if data_store.model_type == 'diff' or data_store.model_type == 'all': 1607 spin_index = None 1608 spin, data_store.spin_id = None, None 1609 elif min_algor == 'back_calc': 1610 spin_index = opt_params.min_options[0] 1611 spin, data_store.spin_id = return_spin_from_index(global_index=spin_index, return_spin_id=True) 1612 else: 1613 spin_index = i 1614 spin, data_store.spin_id = return_spin_from_index(global_index=spin_index, return_spin_id=True) 1615 1616 # Individual spin stuff. 1617 if spin and (data_store.model_type == 'mf' or data_store.model_type == 'local_tm') and not min_algor == 'back_calc': 1618 # Skip deselected spins. 1619 if not spin.select: 1620 continue 1621 1622 # Skip spins missing relaxation data or errors. 1623 if not hasattr(spin, 'ri_data') or not hasattr(spin, 'ri_data_err'): 1624 continue 1625 1626 # Parameter vector and diagonal scaling. 1627 if min_algor == 'back_calc': 1628 # Create the initial parameter vector. 1629 opt_params.param_vector = self._assemble_param_vector(spin=spin, model_type=data_store.model_type) 1630 1631 # Diagonal scaling. 1632 data_store.scaling_matrix = None 1633 1634 else: 1635 # Create the initial parameter vector. 1636 opt_params.param_vector = self._assemble_param_vector(spin=spin, sim_index=sim_index) 1637 1638 # The number of parameters. 1639 num_params = len(opt_params.param_vector) 1640 1641 # Diagonal scaling. 1642 data_store.scaling_matrix = self._assemble_scaling_matrix(num_params, model_type=data_store.model_type, spin=spin, scaling=scaling) 1643 if len(data_store.scaling_matrix): 1644 opt_params.param_vector = dot(inv(data_store.scaling_matrix), opt_params.param_vector) 1645 1646 # Configure the grid search. 1647 opt_params.inc, opt_params.lower, opt_params.upper = None, None, None 1648 if match('^[Gg]rid', min_algor): 1649 opt_params.inc, opt_params.lower, opt_params.upper = self._grid_search_config(num_params, spin=spin, lower=lower, upper=upper, inc=inc, scaling_matrix=data_store.scaling_matrix) 1650 1651 # Scaling of values for the set function. 1652 if match('^[Ss]et', min_algor): 1653 opt_params.min_options = dot(inv(data_store.scaling_matrix), opt_params.min_options) 1654 1655 # Linear constraints. 1656 if constraints: 1657 opt_params.A, opt_params.b = self._linear_constraints(num_params, model_type=data_store.model_type, spin=spin, scaling_matrix=data_store.scaling_matrix) 1658 else: 1659 opt_params.A, opt_params.b = None, None 1660 1661 # Get the data for minimisation. 1662 self._minimise_data_setup(data_store, min_algor, num_data_sets, opt_params.min_options, spin=spin, sim_index=sim_index) 1663 1664 # Setup the minimisation algorithm when constraints are present. 1665 if constraints and not match('^[Gg]rid', min_algor): 1666 algor = opt_params.min_options[0] 1667 else: 1668 algor = min_algor 1669 1670 # Initialise the function to minimise (for back-calculation and LM minimisation). 1671 if min_algor == 'back_calc' or match('[Ll][Mm]$', algor) or match('[Ll]evenburg-[Mm]arquardt$', algor): 1672 self.mf = Mf(init_params=opt_params.param_vector, model_type=data_store.model_type, diff_type=data_store.diff_type, diff_params=data_store.diff_params, scaling_matrix=data_store.scaling_matrix, num_spins=data_store.num_spins, equations=data_store.equations, param_types=data_store.param_types, param_values=data_store.param_values, relax_data=data_store.ri_data, errors=data_store.ri_data_err, bond_length=data_store.r, csa=data_store.csa, num_frq=data_store.num_frq, frq=data_store.frq, num_ri=data_store.num_ri, remap_table=data_store.remap_table, noe_r1_table=data_store.noe_r1_table, ri_labels=data_store.ri_types, gx=data_store.gx, gh=data_store.gh, h_bar=data_store.h_bar, mu0=data_store.mu0, num_params=data_store.num_params, vectors=data_store.xh_unit_vectors) 1673 1674 # Levenberg-Marquardt minimisation. 1675 if match('[Ll][Mm]$', algor) or match('[Ll]evenburg-[Mm]arquardt$', algor): 1676 # Total number of ri. 1677 number_ri = 0 1678 for k in xrange(len(ri_data_err)): 1679 number_ri = number_ri + len(ri_data_err[k]) 1680 1681 # Reconstruct the error data structure. 1682 lm_error = zeros(number_ri, float64) 1683 index = 0 1684 for k in xrange(len(ri_data_err)): 1685 lm_error[index:index+len(ri_data_err[k])] = ri_data_err[k] 1686 index = index + len(ri_data_err[k]) 1687 1688 opt_params.min_options = opt_params.min_options + (self.mf.lm_dri, lm_error) 1689 1690 # Back-calculation. 1691 if min_algor == 'back_calc': 1692 return self.mf.calc_ri() 1693 1694 # Parallelised grid search for the diffusion parameter space. 1695 if match('^[Gg]rid', min_algor) and data_store.model_type == 'diff': 1696 # Print out. 1697 print("Parallelised diffusion tensor grid search.") 1698 1699 # Loop over each grid subdivision. 1700 for subdivision in grid_split(divisions=processor.processor_size(), lower=opt_params.lower, upper=opt_params.upper, inc=opt_params.inc): 1701 # Set the points. 1702 opt_params.subdivision = subdivision 1703 1704 # Grid search initialisation. 1705 command = MF_grid_command() 1706 1707 # Pass in the data and optimisation parameters. 1708 command.store_data(deepcopy(data_store), deepcopy(opt_params)) 1709 1710 # Set up the model-free memo and add it to the processor queue. 1711 memo = MF_memo(model_free=self, model_type=data_store.model_type, spin=spin, sim_index=sim_index, scaling=scaling, scaling_matrix=data_store.scaling_matrix) 1712 processor.add_to_queue(command, memo) 1713 1714 # Execute the queued elements. 1715 processor.run_queue() 1716 1717 # Exit this method. 1718 return 1719 1720 # Normal grid search (command initialisation). 1721 if search('^[Gg]rid', min_algor): 1722 command = MF_grid_command() 1723 1724 # Minimisation of all other model types (command initialisation). 1725 else: 1726 command = MF_minimise_command() 1727 1728 # Pass in the data and optimisation parameters. 1729 command.store_data(deepcopy(data_store), deepcopy(opt_params)) 1730 1731 # Set up the model-free memo and add it to the processor queue. 1732 memo = MF_memo(model_free=self, model_type=data_store.model_type, spin=spin, sim_index=sim_index, scaling=scaling, scaling_matrix=data_store.scaling_matrix) 1733 processor.add_to_queue(command, memo) 1734 1735 # Execute the queued elements. 1736 processor.run_queue()

Source Code for Module specific_fns.model_free.mf_minimise