specific_fns.model_free.mf

56 """Class containing functions specific to model-free optimisation.""" 57

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

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

296 297

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

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

456 457

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

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

531 532

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

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

599 600

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

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

650 651

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

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

987 988

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

990 """Set up all the data required for minimisation. 991 992 @param data_store: A data storage container. 993 @type data_store: class instance 994 @param min_algor: The minimisation algorithm to use. 995 @type min_algor: str 996 @param num_data_sets: The number of data sets. 997 @type num_data_sets: int 998 @param min_options: The minimisation options array. 999 @type min_options: list 1000 @keyword spin: The spin data container. 1001 @type spin: SpinContainer instance 1002 @keyword sim_index: The optional MC simulation index. 1003 @type sim_index: int 1004 @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) 1005 @rtype: tuple 1006 """ 1007 1008 # Initialise the data structures for the model-free function. 1009 data_store.ri_data = [] 1010 data_store.ri_data_err = [] 1011 data_store.equations = [] 1012 data_store.param_types = [] 1013 data_store.param_values = None 1014 data_store.r = [] 1015 data_store.csa = [] 1016 data_store.num_frq = [] 1017 data_store.frq = [] 1018 data_store.num_ri = [] 1019 data_store.remap_table = [] 1020 data_store.noe_r1_table = [] 1021 data_store.ri_types = [] 1022 data_store.gx = [] 1023 data_store.gh = [] 1024 data_store.num_params = [] 1025 data_store.xh_unit_vectors = [] 1026 if data_store.model_type == 'local_tm': 1027 data_store.mf_params = [] 1028 elif data_store.model_type == 'diff': 1029 data_store.param_values = [] 1030 1031 # Set up the data for the back_calc function. 1032 if min_algor == 'back_calc': 1033 # The spin data. 1034 data_store.ri_data = [0.0] 1035 data_store.ri_data_err = [0.000001] 1036 data_store.equations = [spin.equation] 1037 data_store.param_types = [spin.params] 1038 data_store.csa = [spin.csa] 1039 data_store.num_frq = [1] 1040 data_store.frq = [[min_options[3]]] 1041 data_store.num_ri = [1] 1042 data_store.remap_table = [[0]] 1043 data_store.noe_r1_table = [[None]] 1044 data_store.ri_types = [[min_options[2]]] 1045 data_store.gx = [return_gyromagnetic_ratio(spin.isotope)] 1046 1047 # The interatomic data. 1048 interatoms = return_interatom_list(data_store.spin_id) 1049 for i in range(len(interatoms)): 1050 # No relaxation mechanism. 1051 if not interatoms[i].dipole_pair: 1052 continue 1053 1054 # The surrounding spins. 1055 if data_store.spin_id != interatoms[i].spin_id1: 1056 spin_id2 = interatoms[i].spin_id1 1057 else: 1058 spin_id2 = interatoms[i].spin_id2 1059 spin2 = return_spin(spin_id2) 1060 1061 # The data. 1062 data_store.r = [interatoms[i].r] 1063 data_store.gh = [return_gyromagnetic_ratio(spin2.isotope)] 1064 if data_store.model_type != 'local_tm' and cdp.diff_tensor.type != 'sphere': 1065 data_store.xh_unit_vectors = [interatoms[i].vector] 1066 else: 1067 data_store.xh_unit_vectors = [None] 1068 1069 # Count the number of model-free parameters for the spin index. 1070 data_store.num_params = [len(spin.params)] 1071 1072 # Loop over the number of data sets. 1073 for j in range(num_data_sets): 1074 # Set the spin index and get the spin, if not already set. 1075 if data_store.model_type == 'diff' or data_store.model_type == 'all': 1076 spin_index = j 1077 spin, data_store.spin_id = return_spin_from_index(global_index=spin_index, return_spin_id=True) 1078 1079 # Skip deselected spins. 1080 if not spin.select: 1081 continue 1082 1083 # Skip spins where there is no data or errors. 1084 if not hasattr(spin, 'ri_data') or not hasattr(spin, 'ri_data_err'): 1085 continue 1086 1087 # Make sure that the errors are strictly positive numbers. 1088 for ri_id in cdp.ri_ids: 1089 # Skip missing data. 1090 if not ri_id in spin.ri_data_err: 1091 continue 1092 1093 # Alias. 1094 err = spin.ri_data_err[ri_id] 1095 1096 # Checks. 1097 if err != None and err == 0.0: 1098 raise RelaxError("Zero error for spin '%s' for the relaxation data ID '%s', minimisation not possible." % (data_store.spin_id, ri_id)) 1099 elif err != None and err < 0.0: 1100 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)) 1101 1102 # The relaxation data optimisation structures. 1103 data = self._relax_data_opt_structs(spin, sim_index=sim_index) 1104 1105 # Append the data. 1106 data_store.ri_data.append(data[0]) 1107 data_store.ri_data_err.append(data[1]) 1108 data_store.num_frq.append(data[2]) 1109 data_store.num_ri.append(data[3]) 1110 data_store.ri_types.append(data[4]) 1111 data_store.frq.append(data[5]) 1112 data_store.remap_table.append(data[6]) 1113 data_store.noe_r1_table.append(data[7]) 1114 if sim_index == None or data_store.model_type == 'diff': 1115 data_store.csa.append(spin.csa) 1116 else: 1117 data_store.csa.append(spin.csa_sim[sim_index]) 1118 1119 # Repackage the spin data. 1120 data_store.equations.append(spin.equation) 1121 data_store.param_types.append(spin.params) 1122 data_store.gx.append(return_gyromagnetic_ratio(spin.isotope)) 1123 1124 # Repackage the interatomic data. 1125 interatoms = return_interatom_list(data_store.spin_id) 1126 for i in range(len(interatoms)): 1127 # No relaxation mechanism. 1128 if not interatoms[i].dipole_pair: 1129 continue 1130 1131 # The surrounding spins. 1132 if data_store.spin_id != interatoms[i].spin_id1: 1133 spin_id2 = interatoms[i].spin_id1 1134 else: 1135 spin_id2 = interatoms[i].spin_id2 1136 spin2 = return_spin(spin_id2) 1137 1138 # The data. 1139 data_store.gh.append(return_gyromagnetic_ratio(spin2.isotope)) 1140 if sim_index == None or data_store.model_type == 'diff' or not hasattr(interatoms[i], 'r_sim'): 1141 data_store.r.append(interatoms[i].r) 1142 else: 1143 data_store.r.append(interatoms[i].r_sim[sim_index]) 1144 1145 # Vectors. 1146 if data_store.model_type != 'local_tm' and cdp.diff_tensor.type != 'sphere': 1147 # Check that this is a single vector! 1148 if arg_check.is_num_list(interatoms[i].vector[0], raise_error=False): 1149 raise RelaxMultiVectorError(data_store.spin_id) 1150 1151 # Store the vector. 1152 data_store.xh_unit_vectors.append(interatoms[i].vector) 1153 1154 # No vector. 1155 else: 1156 data_store.xh_unit_vectors.append(None) 1157 1158 # Stop - only one mechanism is current supported. 1159 break 1160 1161 # Model-free parameter values. 1162 if data_store.model_type == 'local_tm': 1163 pass 1164 1165 # Count the number of model-free parameters for the spin index. 1166 data_store.num_params.append(len(spin.params)) 1167 1168 # Repackage the parameter values for minimising just the diffusion tensor parameters. 1169 if data_store.model_type == 'diff': 1170 data_store.param_values.append(self._assemble_param_vector(model_type='mf')) 1171 1172 # Convert to numpy arrays. 1173 for k in range(len(data_store.ri_data)): 1174 data_store.ri_data[k] = array(data_store.ri_data[k], float64) 1175 data_store.ri_data_err[k] = array(data_store.ri_data_err[k], float64) 1176 1177 # Diffusion tensor type. 1178 if data_store.model_type == 'local_tm': 1179 data_store.diff_type = 'sphere' 1180 else: 1181 data_store.diff_type = cdp.diff_tensor.type 1182 1183 # Package the diffusion tensor parameters. 1184 data_store.diff_params = None 1185 if data_store.model_type == 'mf': 1186 # Spherical diffusion. 1187 if data_store.diff_type == 'sphere': 1188 data_store.diff_params = [cdp.diff_tensor.tm] 1189 1190 # Spheroidal diffusion. 1191 elif data_store.diff_type == 'spheroid': 1192 data_store.diff_params = [cdp.diff_tensor.tm, cdp.diff_tensor.Da, cdp.diff_tensor.theta, cdp.diff_tensor.phi] 1193 1194 # Ellipsoidal diffusion. 1195 elif data_store.diff_type == 'ellipsoid': 1196 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] 1197 elif min_algor == 'back_calc' and data_store.model_type == 'local_tm': 1198 # Spherical diffusion. 1199 data_store.diff_params = [spin.local_tm]

1200 1201

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

1203 """Package the relaxation data into the data structures used for optimisation. 1204 1205 @param spin: The spin container to extract the data from. 1206 @type spin: SpinContainer instance 1207 @keyword sim_index: The optional MC simulation index. 1208 @type sim_index: int 1209 @return: The structures ri_data, ri_data_err, num_frq, num_ri, ri_ids, frq, remap_table, noe_r1_table. 1210 @rtype: tuple 1211 """ 1212 1213 # Initialise the data. 1214 ri_data = [] 1215 ri_data_err = [] 1216 ri_labels = [] 1217 frq = [] 1218 remap_table = [] 1219 noe_r1_table = [] 1220 1221 # Loop over the relaxation data. 1222 for ri_id in cdp.ri_ids: 1223 # Skip missing data. 1224 if ri_id not in spin.ri_data: 1225 continue 1226 1227 # The Rx data. 1228 if sim_index == None: 1229 data = spin.ri_data[ri_id] 1230 else: 1231 data = spin.ri_data_sim[ri_id][sim_index] 1232 1233 # The errors. 1234 err = spin.ri_data_err[ri_id] 1235 1236 # Missing data, so don't add it. 1237 if data == None or err == None: 1238 continue 1239 1240 # Append the data and error. 1241 ri_data.append(data) 1242 ri_data_err.append(err) 1243 1244 # The labels. 1245 ri_labels.append(cdp.ri_type[ri_id]) 1246 1247 # The frequencies. 1248 if cdp.frq[ri_id] not in frq: 1249 frq.append(cdp.frq[ri_id]) 1250 1251 # The remap table. 1252 remap_table.append(frq.index(cdp.frq[ri_id])) 1253 1254 # The NOE to R1 mapping table. 1255 noe_r1_table.append(None) 1256 1257 # The number of data sets. 1258 num_ri = len(ri_data) 1259 1260 # Fill the NOE to R1 mapping table. 1261 for i in range(num_ri): 1262 # If the data corresponds to 'NOE', try to find if the corresponding R1 data. 1263 if cdp.ri_type[cdp.ri_ids[i]] == 'NOE': 1264 for j in range(num_ri): 1265 if cdp.ri_type[cdp.ri_ids[j]] == 'R1' and cdp.frq[cdp.ri_ids[i]] == cdp.frq[cdp.ri_ids[j]]: 1266 noe_r1_table[i] = j 1267 1268 # Return the structures. 1269 return ri_data, ri_data_err, len(frq), num_ri, ri_labels, frq, remap_table, noe_r1_table

1270 1271

1272 - def _reset_min_stats(self):

1273 """Reset all the minimisation statistics. 1274 1275 All global and spin specific values will be set to None. 1276 """ 1277 1278 # Global stats. 1279 if hasattr(cdp, 'chi2'): 1280 cdp.chi2 = None 1281 cdp.iter = None 1282 cdp.f_count = None 1283 cdp.g_count = None 1284 cdp.h_count = None 1285 cdp.warning = None 1286 1287 # Spin specific stats. 1288 for spin in spin_loop(): 1289 if hasattr(spin, 'chi2'): 1290 spin.chi2 = None 1291 spin.iter = None 1292 spin.f_count = None 1293 spin.g_count = None 1294 spin.h_count = None 1295 spin.warning = None

1296 1297

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

1299 """Calculation of the model-free chi-squared value. 1300 1301 @keyword spin_id: The spin identification string. 1302 @type spin_id: str 1303 @keyword verbosity: The amount of information to print. The higher the value, the greater the verbosity. 1304 @type verbosity: int 1305 @keyword sim_index: The optional MC simulation index. 1306 @type sim_index: int 1307 """ 1308 1309 # Test if sequence data is loaded. 1310 if not exists_mol_res_spin_data(): 1311 raise RelaxNoSequenceError 1312 1313 # Determine the model type. 1314 model_type = self._determine_model_type() 1315 1316 # Test if diffusion tensor data exists. 1317 if model_type != 'local_tm' and not diff_data_exists(): 1318 raise RelaxNoTensorError('diffusion') 1319 1320 # Test if the PDB file has been loaded. 1321 if model_type != 'local_tm' and cdp.diff_tensor.type != 'sphere' and not hasattr(cdp, 'structure'): 1322 raise RelaxNoPdbError 1323 1324 # Loop over the spins. 1325 for spin, id in spin_loop(spin_id, return_id=True): 1326 # Skip deselected spins. 1327 if not spin.select: 1328 continue 1329 1330 # Test if the model-free model has been setup. 1331 if not spin.model: 1332 raise RelaxNoModelError 1333 1334 # Test if the nuclear isotope type has been set. 1335 if not hasattr(spin, 'isotope'): 1336 raise RelaxSpinTypeError 1337 1338 # Test if the model-free parameter values exist. 1339 unset_param = self._are_mf_params_set(spin) 1340 if unset_param != None: 1341 raise RelaxNoValueError(unset_param) 1342 1343 # Test if the CSA value has been set. 1344 if not hasattr(spin, 'csa') or spin.csa == None: 1345 raise RelaxNoValueError("CSA") 1346 1347 # Test the interatomic data. 1348 interatoms = return_interatom_list(spin_id) 1349 for interatom in interatoms: 1350 # No relaxation mechanism. 1351 if not interatom.dipole_pair: 1352 continue 1353 1354 # Test if unit vectors exist. 1355 if model_type != 'local_tm' and cdp.diff_tensor.type != 'sphere' and not hasattr(interatom, 'vector'): 1356 raise RelaxNoVectorsError 1357 1358 # Test if multiple unit vectors exist. 1359 if model_type != 'local_tm' and cdp.diff_tensor.type != 'sphere' and hasattr(interatom, 'vector') and arg_check.is_num_list(interatom.vector[0], raise_error=False): 1360 raise RelaxMultiVectorError 1361 1362 # The interacting spin. 1363 if spin_id != interatom.spin_id1: 1364 spin_id2 = interatom.spin_id1 1365 else: 1366 spin_id2 = interatom.spin_id2 1367 spin2 = return_spin(spin_id2) 1368 1369 # Test if the nuclear isotope type has been set. 1370 if not hasattr(spin2, 'isotope'): 1371 raise RelaxSpinTypeError 1372 1373 # Test if the interatomic distance has been set. 1374 if not hasattr(interatom, 'r') or interatom.r == None: 1375 raise RelaxNoValueError("interatomic distance", spin_id=spin_id, spin_id2=spin_id2) 1376 1377 # Skip spins where there is no data or errors. 1378 if not hasattr(spin, 'ri_data') or not hasattr(spin, 'ri_data_err'): 1379 continue 1380 1381 # Make sure that the errors are strictly positive numbers. 1382 for ri_id in cdp.ri_ids: 1383 # Alias. 1384 err = spin.ri_data_err[ri_id] 1385 1386 # Checks. 1387 if err != None and err == 0.0: 1388 raise RelaxError("Zero error for spin '%s' for the relaxation data ID '%s', minimisation not possible." % (id, ri_id)) 1389 elif err != None and err < 0.0: 1390 raise RelaxError("Negative error of %s for spin '%s' for the relaxation data ID '%s', minimisation not possible." % (err, id, ri_id)) 1391 1392 # Create the initial parameter vector. 1393 param_vector = self._assemble_param_vector(spin=spin, sim_index=sim_index) 1394 1395 # The relaxation data optimisation structures. 1396 data = self._relax_data_opt_structs(spin, sim_index=sim_index) 1397 1398 # The spin data. 1399 ri_data = [array(data[0])] 1400 ri_data_err = [array(data[1])] 1401 num_frq = [data[2]] 1402 num_ri = [data[3]] 1403 ri_labels = [data[4]] 1404 frq = [data[5]] 1405 remap_table = [data[6]] 1406 noe_r1_table = [data[7]] 1407 gx = [return_gyromagnetic_ratio(spin.isotope)] 1408 if sim_index == None: 1409 csa = [spin.csa] 1410 else: 1411 csa = [spin.csa_sim[sim_index]] 1412 1413 # The interatomic data. 1414 interatoms = return_interatom_list(spin_id) 1415 for i in range(len(interatoms)): 1416 # No relaxation mechanism. 1417 if not interatoms[i].dipole_pair: 1418 continue 1419 1420 # The surrounding spins. 1421 if spin_id != interatoms[i].spin_id1: 1422 spin_id2 = interatoms[i].spin_id1 1423 else: 1424 spin_id2 = interatoms[i].spin_id2 1425 spin2 = return_spin(spin_id2) 1426 1427 # The data. 1428 if sim_index == None: 1429 r = [interatoms[i].r] 1430 else: 1431 r = [interatoms[i].r_sim[sim_index]] 1432 1433 # Vectors. 1434 if model_type != 'local_tm' and cdp.diff_tensor.type != 'sphere': 1435 xh_unit_vectors = [interatoms[i].vector] 1436 else: 1437 xh_unit_vectors = [None] 1438 1439 # Gyromagnetic ratios. 1440 gh = [return_gyromagnetic_ratio(spin2.isotope)] 1441 1442 # Count the number of model-free parameters for the residue index. 1443 num_params = [len(spin.params)] 1444 1445 # Repackage the parameter values as a local model (ignore if the diffusion tensor is not fixed). 1446 param_values = [self._assemble_param_vector(model_type='mf')] 1447 1448 # Package the diffusion tensor parameters. 1449 if model_type == 'local_tm': 1450 diff_params = [spin.local_tm] 1451 diff_type = 'sphere' 1452 else: 1453 # Diff type. 1454 diff_type = cdp.diff_tensor.type 1455 1456 # Spherical diffusion. 1457 if diff_type == 'sphere': 1458 diff_params = [cdp.diff_tensor.tm] 1459 1460 # Spheroidal diffusion. 1461 elif diff_type == 'spheroid': 1462 diff_params = [cdp.diff_tensor.tm, cdp.diff_tensor.Da, cdp.diff_tensor.theta, cdp.diff_tensor.phi] 1463 1464 # Ellipsoidal diffusion. 1465 elif diff_type == 'ellipsoid': 1466 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] 1467 1468 # Initialise the model-free function. 1469 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) 1470 1471 # Chi-squared calculation. 1472 try: 1473 chi2 = mf.func(param_vector) 1474 except OverflowError: 1475 chi2 = 1e200 1476 1477 # Global chi-squared value. 1478 if model_type == 'all' or model_type == 'diff': 1479 cdp.chi2 = cdp.chi2 + chi2 1480 else: 1481 spin.chi2 = chi2

1482 1483

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

1485 """The model-free grid search function. 1486 1487 @keyword lower: The lower bounds of the grid search which must be equal to the 1488 number of parameters in the model. 1489 @type lower: array of numbers 1490 @keyword upper: The upper bounds of the grid search which must be equal to the 1491 number of parameters in the model. 1492 @type upper: array of numbers 1493 @keyword inc: The increments for each dimension of the space for the grid search. 1494 The number of elements in the array must equal to the number of 1495 parameters in the model. 1496 @type inc: array of int 1497 @keyword constraints: If True, constraints are applied during the grid search (eliminating 1498 parts of the grid). If False, no constraints are used. 1499 @type constraints: bool 1500 @keyword verbosity: A flag specifying the amount of information to print. The higher 1501 the value, the greater the verbosity. 1502 @type verbosity: int 1503 @keyword sim_index: The index of the simulation to apply the grid search to. If None, 1504 the normal model is optimised. 1505 @type sim_index: int 1506 """ 1507 1508 # Minimisation. 1509 self.minimise(min_algor='grid', lower=lower, upper=upper, inc=inc, constraints=constraints, verbosity=verbosity, sim_index=sim_index)

1510 1511

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

1513 """Model-free minimisation function. 1514 1515 Three categories of models exist for which the approach to minimisation is different. These 1516 are: 1517 1518 Single spin optimisations: The 'mf' and 'local_tm' model types which are the 1519 model-free parameters for single spins, optionally with a local tm parameter. These 1520 models have no global parameters. 1521 1522 Diffusion tensor optimisations: The 'diff' diffusion tensor model type. No spin 1523 specific parameters exist. 1524 1525 Optimisation of everything: The 'all' model type consisting of all model-free and all 1526 diffusion tensor parameters. 1527 1528 1529 @keyword min_algor: The minimisation algorithm to use. 1530 @type min_algor: str 1531 @keyword min_options: An array of options to be used by the minimisation algorithm. 1532 @type min_options: array of str 1533 @keyword func_tol: The function tolerance which, when reached, terminates optimisation. 1534 Setting this to None turns of the check. 1535 @type func_tol: None or float 1536 @keyword grad_tol: The gradient tolerance which, when reached, terminates optimisation. 1537 Setting this to None turns of the check. 1538 @type grad_tol: None or float 1539 @keyword max_iterations: The maximum number of iterations for the algorithm. 1540 @type max_iterations: int 1541 @keyword constraints: If True, constraints are used during optimisation. 1542 @type constraints: bool 1543 @keyword scaling: If True, diagonal scaling is enabled during optimisation to allow 1544 the problem to be better conditioned. 1545 @type scaling: bool 1546 @keyword verbosity: The amount of information to print. The higher the value, the 1547 greater the verbosity. 1548 @type verbosity: int 1549 @keyword sim_index: The index of the simulation to optimise. This should be None if 1550 normal optimisation is desired. 1551 @type sim_index: None or int 1552 @keyword lower: The lower bounds of the grid search which must be equal to the 1553 number of parameters in the model. This optional argument is only 1554 used when doing a grid search. 1555 @type lower: array of numbers 1556 @keyword upper: The upper bounds of the grid search which must be equal to the 1557 number of parameters in the model. This optional argument is only 1558 used when doing a grid search. 1559 @type upper: array of numbers 1560 @keyword inc: The increments for each dimension of the space for the grid search. 1561 The number of elements in the array must equal to the number of 1562 parameters in the model. This argument is only used when doing a 1563 grid search. 1564 @type inc: array of int 1565 """ 1566 1567 # Test if sequence data is loaded. 1568 if not exists_mol_res_spin_data(): 1569 raise RelaxNoSequenceError 1570 1571 # Test if the model-free model has been setup, and that the nuclear isotope types have been set. 1572 for spin in spin_loop(): 1573 # Skip deselected spins. 1574 if not spin.select: 1575 continue 1576 1577 # Not setup. 1578 if not spin.model: 1579 raise RelaxNoModelError 1580 1581 # Test if the nuclear isotope type has been set. 1582 if not hasattr(spin, 'isotope'): 1583 raise RelaxSpinTypeError 1584 1585 # Reset the minimisation statistics. 1586 if sim_index == None and min_algor != 'back_calc': 1587 self._reset_min_stats() 1588 1589 # Containers for the model-free data and optimisation parameters. 1590 data_store = Data_container() 1591 opt_params = Data_container() 1592 1593 # Add the imported parameters. 1594 data_store.h_bar = h_bar 1595 data_store.mu0 = mu0 1596 opt_params.min_algor = min_algor 1597 opt_params.min_options = min_options 1598 opt_params.func_tol = func_tol 1599 opt_params.grad_tol = grad_tol 1600 opt_params.max_iterations = max_iterations 1601 1602 # Add the keyword args. 1603 opt_params.verbosity = verbosity 1604 1605 # Determine the model type. 1606 data_store.model_type = self._determine_model_type() 1607 if not data_store.model_type: 1608 return 1609 1610 # Model type for the back-calculate function. 1611 if min_algor == 'back_calc' and data_store.model_type != 'local_tm': 1612 data_store.model_type = 'mf' 1613 1614 # Test if diffusion tensor data exists. 1615 if data_store.model_type != 'local_tm' and not diffusion_tensor.diff_data_exists(): 1616 raise RelaxNoTensorError('diffusion') 1617 1618 # Tests for the PDB file and unit vectors. 1619 if data_store.model_type != 'local_tm' and cdp.diff_tensor.type != 'sphere': 1620 # Test if the structure file has been loaded. 1621 if not hasattr(cdp, 'structure'): 1622 raise RelaxNoPdbError 1623 1624 # Test if unit vectors exist. 1625 for spin, spin_id in spin_loop(return_id=True): 1626 # Skip deselected spins. 1627 if not spin.select: 1628 continue 1629 1630 # Get the interatomic data container. 1631 interatoms = return_interatom_list(spin_id) 1632 1633 # Unit vectors. 1634 for i in range(len(interatoms)): 1635 # No relaxation mechanism. 1636 if not interatoms[i].dipole_pair: 1637 continue 1638 1639 # Check for the vectors. 1640 if not hasattr(interatoms[i], 'vector'): 1641 raise RelaxNoVectorsError 1642 1643 # Test if the model-free parameter values are set for minimising diffusion tensor parameters by themselves. 1644 if data_store.model_type == 'diff': 1645 # Loop over the sequence. 1646 for spin in spin_loop(): 1647 unset_param = self._are_mf_params_set(spin) 1648 if unset_param != None: 1649 raise RelaxNoValueError(unset_param) 1650 1651 # Print out. 1652 if verbosity >= 1: 1653 if data_store.model_type == 'mf': 1654 print("Only the model-free parameters for single spins will be used.") 1655 elif data_store.model_type == 'local_mf': 1656 print("Only a local tm value together with the model-free parameters for single spins will be used.") 1657 elif data_store.model_type == 'diff': 1658 print("Only diffusion tensor parameters will be used.") 1659 elif data_store.model_type == 'all': 1660 print("The diffusion tensor parameters together with the model-free parameters for all spins will be used.") 1661 1662 # Test if the CSA and interatomic distances have been set. 1663 for spin, spin_id in spin_loop(return_id=True): 1664 # Skip deselected spins. 1665 if not spin.select: 1666 continue 1667 1668 # CSA value. 1669 if not hasattr(spin, 'csa') or spin.csa == None: 1670 raise RelaxNoValueError("CSA") 1671 1672 # Get the interatomic data container. 1673 interatoms = return_interatom_list(spin_id) 1674 1675 # Interatomic distances. 1676 count = 0 1677 for i in range(len(interatoms)): 1678 # No relaxation mechanism. 1679 if not interatoms[i].dipole_pair: 1680 continue 1681 1682 # Check for the distances. 1683 if not hasattr(interatoms[i], 'r') or interatoms[i].r == None: 1684 raise RelaxNoValueError("interatomic distance", spin_id=interatoms[i].spin_id1, spin_id2=interatoms[i].spin_id2) 1685 1686 # Count the number of interactions. 1687 count += 1 1688 1689 # Too many interactions. 1690 if count > 1: 1691 raise RelaxError("The spin '%s' has %s dipolar relaxation interactions defined, but only a maximum of one is currently supported." % (spin_id, count)) 1692 1693 # Number of spins, minimisation instances, and data sets for each model type. 1694 if data_store.model_type == 'mf' or data_store.model_type == 'local_tm': 1695 num_instances = count_spins(skip_desel=False) 1696 num_data_sets = 1 1697 data_store.num_spins = 1 1698 elif data_store.model_type == 'diff' or data_store.model_type == 'all': 1699 num_instances = 1 1700 num_data_sets = count_spins(skip_desel=False) 1701 data_store.num_spins = count_spins() 1702 1703 # Number of spins, minimisation instances, and data sets for the back-calculate function. 1704 if min_algor == 'back_calc': 1705 num_instances = 1 1706 num_data_sets = 0 1707 data_store.num_spins = 1 1708 1709 # Get the Processor box singleton (it contains the Processor instance) and alias the Processor. 1710 processor_box = Processor_box() 1711 processor = processor_box.processor 1712 1713 # Loop over the minimisation instances. 1714 ####################################### 1715 1716 for i in range(num_instances): 1717 # Get the spin container if required. 1718 if data_store.model_type == 'diff' or data_store.model_type == 'all': 1719 spin_index = None 1720 spin, data_store.spin_id = None, None 1721 elif min_algor == 'back_calc': 1722 spin_index = opt_params.min_options[0] 1723 spin, data_store.spin_id = return_spin_from_index(global_index=spin_index, return_spin_id=True) 1724 else: 1725 spin_index = i 1726 spin, data_store.spin_id = return_spin_from_index(global_index=spin_index, return_spin_id=True) 1727 1728 # Individual spin stuff. 1729 if spin and (data_store.model_type == 'mf' or data_store.model_type == 'local_tm') and not min_algor == 'back_calc': 1730 # Skip deselected spins. 1731 if not spin.select: 1732 continue 1733 1734 # Skip spins missing relaxation data or errors. 1735 if not hasattr(spin, 'ri_data') or not hasattr(spin, 'ri_data_err'): 1736 continue 1737 1738 # Skip spins missing the dipolar interaction. 1739 if spin and (data_store.model_type == 'mf' or data_store.model_type == 'local_tm'): 1740 interatoms = return_interatom_list(data_store.spin_id) 1741 if not len(interatoms): 1742 continue 1743 1744 # Parameter vector and diagonal scaling. 1745 if min_algor == 'back_calc': 1746 # Create the initial parameter vector. 1747 opt_params.param_vector = self._assemble_param_vector(spin=spin, model_type=data_store.model_type) 1748 1749 # Diagonal scaling. 1750 data_store.scaling_matrix = None 1751 1752 else: 1753 # Create the initial parameter vector. 1754 opt_params.param_vector = self._assemble_param_vector(spin=spin, sim_index=sim_index) 1755 1756 # The number of parameters. 1757 num_params = len(opt_params.param_vector) 1758 1759 # Diagonal scaling. 1760 data_store.scaling_matrix = self._assemble_scaling_matrix(num_params, model_type=data_store.model_type, spin=spin, scaling=scaling) 1761 if len(data_store.scaling_matrix): 1762 opt_params.param_vector = dot(inv(data_store.scaling_matrix), opt_params.param_vector) 1763 1764 # Configure the grid search. 1765 opt_params.inc, opt_params.lower, opt_params.upper = None, None, None 1766 if match('^[Gg]rid', min_algor): 1767 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) 1768 1769 # Scaling of values for the set function. 1770 if match('^[Ss]et', min_algor): 1771 opt_params.min_options = dot(inv(data_store.scaling_matrix), opt_params.min_options) 1772 1773 # Linear constraints. 1774 if constraints: 1775 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) 1776 else: 1777 opt_params.A, opt_params.b = None, None 1778 1779 # Get the data for minimisation. 1780 self._minimise_data_setup(data_store, min_algor, num_data_sets, opt_params.min_options, spin=spin, sim_index=sim_index) 1781 1782 # Setup the minimisation algorithm when constraints are present. 1783 if constraints and not match('^[Gg]rid', min_algor): 1784 algor = opt_params.min_options[0] 1785 else: 1786 algor = min_algor 1787 1788 # Initialise the function to minimise (for back-calculation and LM minimisation). 1789 if min_algor == 'back_calc' or match('[Ll][Mm]$', algor) or match('[Ll]evenburg-[Mm]arquardt$', algor): 1790 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) 1791 1792 # Levenberg-Marquardt minimisation. 1793 if match('[Ll][Mm]$', algor) or match('[Ll]evenburg-[Mm]arquardt$', algor): 1794 # Total number of ri. 1795 number_ri = 0 1796 for k in range(len(ri_data_err)): 1797 number_ri = number_ri + len(ri_data_err[k]) 1798 1799 # Reconstruct the error data structure. 1800 lm_error = zeros(number_ri, float64) 1801 index = 0 1802 for k in range(len(ri_data_err)): 1803 lm_error[index:index+len(ri_data_err[k])] = ri_data_err[k] 1804 index = index + len(ri_data_err[k]) 1805 1806 opt_params.min_options = opt_params.min_options + (self.mf.lm_dri, lm_error) 1807 1808 # Back-calculation. 1809 if min_algor == 'back_calc': 1810 return self.mf.calc_ri() 1811 1812 # Parallelised grid search for the diffusion parameter space. 1813 if match('^[Gg]rid', min_algor) and data_store.model_type == 'diff': 1814 # Print out. 1815 print("Parallelised diffusion tensor grid search.") 1816 1817 # Loop over each grid subdivision. 1818 for subdivision in grid_split(divisions=processor.processor_size(), lower=opt_params.lower, upper=opt_params.upper, inc=opt_params.inc): 1819 # Set the points. 1820 opt_params.subdivision = subdivision 1821 1822 # Grid search initialisation. 1823 command = MF_grid_command() 1824 1825 # Pass in the data and optimisation parameters. 1826 command.store_data(deepcopy(data_store), deepcopy(opt_params)) 1827 1828 # Set up the model-free memo and add it to the processor queue. 1829 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) 1830 processor.add_to_queue(command, memo) 1831 1832 # Execute the queued elements. 1833 processor.run_queue() 1834 1835 # Exit this method. 1836 return 1837 1838 # Normal grid search (command initialisation). 1839 if search('^[Gg]rid', min_algor): 1840 command = MF_grid_command() 1841 1842 # Minimisation of all other model types (command initialisation). 1843 else: 1844 command = MF_minimise_command() 1845 1846 # Pass in the data and optimisation parameters. 1847 command.store_data(deepcopy(data_store), deepcopy(opt_params)) 1848 1849 # Set up the model-free memo and add it to the processor queue. 1850 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) 1851 processor.add_to_queue(command, memo) 1852 1853 # Execute the queued elements. 1854 processor.run_queue()

Source Code for Module specific_fns.model_free.mf_minimise