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 # Check that this is a single vector! 1106 if arg_check.is_num_list(spin.xh_vect[0], raise_error=False): 1107 raise RelaxMultiVectorError(data_store.spin_id) 1108 1109 # Store the vector. 1110 data_store.xh_unit_vectors.append(spin.xh_vect) 1111 1112 # No vector. 1113 else: 1114 data_store.xh_unit_vectors.append(None) 1115 1116 # Count the number of model-free parameters for the spin index. 1117 data_store.num_params.append(len(spin.params)) 1118 1119 # Repackage the parameter values for minimising just the diffusion tensor parameters. 1120 if data_store.model_type == 'diff': 1121 data_store.param_values.append(self._assemble_param_vector(model_type='mf')) 1122 1123 # Convert to numpy arrays. 1124 for k in xrange(len(data_store.ri_data)): 1125 data_store.ri_data[k] = array(data_store.ri_data[k], float64) 1126 data_store.ri_data_err[k] = array(data_store.ri_data_err[k], float64) 1127 1128 # Diffusion tensor type. 1129 if data_store.model_type == 'local_tm': 1130 data_store.diff_type = 'sphere' 1131 else: 1132 data_store.diff_type = cdp.diff_tensor.type 1133 1134 # Package the diffusion tensor parameters. 1135 data_store.diff_params = None 1136 if data_store.model_type == 'mf': 1137 # Spherical diffusion. 1138 if data_store.diff_type == 'sphere': 1139 data_store.diff_params = [cdp.diff_tensor.tm] 1140 1141 # Spheroidal diffusion. 1142 elif data_store.diff_type == 'spheroid': 1143 data_store.diff_params = [cdp.diff_tensor.tm, cdp.diff_tensor.Da, cdp.diff_tensor.theta, cdp.diff_tensor.phi] 1144 1145 # Ellipsoidal diffusion. 1146 elif data_store.diff_type == 'ellipsoid': 1147 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] 1148 elif min_algor == 'back_calc' and data_store.model_type == 'local_tm': 1149 # Spherical diffusion. 1150 data_store.diff_params = [spin.local_tm]

1151 1152

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

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

1217 1218

1219 - def _reset_min_stats(self):

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

1243 1244

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

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

1399 1400

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

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

1427 1428

1429 - 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):

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

Source Code for Module specific_fns.model_free.mf_minimise