specific_analyses.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.set(param='tm', value=param_vector[0], category='sim', sim_index=sim_index) 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.set(param='tm', value=param_vector[0], category='sim', sim_index=sim_index) 89 cdp.diff_tensor.set(param='Da', value=param_vector[1], category='sim', sim_index=sim_index) 90 cdp.diff_tensor.set(param='theta', value=param_vector[2], category='sim', sim_index=sim_index) 91 cdp.diff_tensor.set(param='phi', value=param_vector[3], category='sim', sim_index=sim_index) 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.set(param='tm', value=param_vector[0], category='sim', sim_index=sim_index) 101 cdp.diff_tensor.set(param='Da', value=param_vector[1], category='sim', sim_index=sim_index) 102 cdp.diff_tensor.set(param='Dr', value=param_vector[2], category='sim', sim_index=sim_index) 103 cdp.diff_tensor.set(param='alpha', value=param_vector[3], category='sim', sim_index=sim_index) 104 cdp.diff_tensor.set(param='beta', value=param_vector[4], category='sim', sim_index=sim_index) 105 cdp.diff_tensor.set(param='gamma', value=param_vector[5], category='sim', sim_index=sim_index) 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.set(param='tm', value=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.set(param='tm', value=param_vector[0]) 125 cdp.diff_tensor.set(param='Da', value=param_vector[1]) 126 cdp.diff_tensor.set(param='theta', value=param_vector[2]) 127 cdp.diff_tensor.set(param='phi', value=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.set(param='tm', value=param_vector[0]) 137 cdp.diff_tensor.set(param='Da', value=param_vector[1]) 138 cdp.diff_tensor.set(param='Dr', value=param_vector[2]) 139 cdp.diff_tensor.set(param='alpha', value=param_vector[3]) 140 cdp.diff_tensor.set(param='beta', value=param_vector[4]) 141 cdp.diff_tensor.set(param='gamma', value=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 range(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 # No result. 327 if param_vector == None: 328 return 329 330 # Alias the current data pipe. 331 cdp = pipes.get_pipe() 332 333 # Catch infinite chi-squared values. 334 if isInf(func): 335 raise RelaxInfError('chi-squared') 336 337 # Catch chi-squared values of NaN. 338 if isNaN(func): 339 raise RelaxNaNError('chi-squared') 340 341 # Scaling. 342 if scaling: 343 param_vector = dot(scaling_matrix, param_vector) 344 345 # Check if the chi-squared value is lower. This allows for a parallelised grid search! 346 if sim_index == None: 347 # Get the correct value. 348 chi2 = None 349 if (model_type == 'mf' or model_type == 'local_tm') and hasattr(cdp, 'chi2'): 350 chi2 = spin.chi2 351 if (model_type == 'diff' or model_type == 'all') and hasattr(cdp, 'chi2'): 352 chi2 = cdp.chi2 353 354 # Spin text. 355 spin_text = '' 356 if spin != None and hasattr(spin, '_spin_ids') and len(spin._spin_ids): 357 spin_text = " for the spin '%s'" % spin._spin_ids[0] 358 359 # No improvement. 360 if chi2 != None and func >= chi2: 361 print("Discarding the optimisation results%s, the optimised chi-squared value is higher than the current value (%s >= %s)." % (spin_text, func, chi2)) 362 363 # Exit! 364 return 365 366 # New minimum. 367 else: 368 print("Storing the optimisation results%s, the optimised chi-squared value is lower than the current value (%s < %s)." % (spin_text, func, chi2)) 369 370 # Disassemble the parameter vector. 371 self._disassemble_param_vector(model_type, param_vector=param_vector, spin=spin, sim_index=sim_index) 372 373 # Monte Carlo minimisation statistics. 374 if sim_index != None: 375 # Sequence specific minimisation statistics. 376 if model_type == 'mf' or model_type == 'local_tm': 377 378 # Chi-squared statistic. 379 spin.chi2_sim[sim_index] = func 380 381 # Iterations. 382 spin.iter_sim[sim_index] = iter 383 384 # Function evaluations. 385 spin.f_count_sim[sim_index] = fc 386 387 # Gradient evaluations. 388 spin.g_count_sim[sim_index] = gc 389 390 # Hessian evaluations. 391 spin.h_count_sim[sim_index] = hc 392 393 # Warning. 394 spin.warning_sim[sim_index] = warning 395 396 # Global minimisation statistics. 397 elif model_type == 'diff' or model_type == 'all': 398 # Chi-squared statistic. 399 cdp.chi2_sim[sim_index] = func 400 401 # Iterations. 402 cdp.iter_sim[sim_index] = iter 403 404 # Function evaluations. 405 cdp.f_count_sim[sim_index] = fc 406 407 # Gradient evaluations. 408 cdp.g_count_sim[sim_index] = gc 409 410 # Hessian evaluations. 411 cdp.h_count_sim[sim_index] = hc 412 413 # Warning. 414 cdp.warning_sim[sim_index] = warning 415 416 # Normal statistics. 417 else: 418 # Sequence specific minimisation statistics. 419 if model_type == 'mf' or model_type == 'local_tm': 420 # Chi-squared statistic. 421 spin.chi2 = func 422 423 # Iterations. 424 spin.iter = iter 425 426 # Function evaluations. 427 spin.f_count = fc 428 429 # Gradient evaluations. 430 spin.g_count = gc 431 432 # Hessian evaluations. 433 spin.h_count = hc 434 435 # Warning. 436 spin.warning = warning 437 438 # Global minimisation statistics. 439 elif model_type == 'diff' or model_type == 'all': 440 # Chi-squared statistic. 441 cdp.chi2 = func 442 443 # Iterations. 444 cdp.iter = iter 445 446 # Function evaluations. 447 cdp.f_count = fc 448 449 # Gradient evaluations. 450 cdp.g_count = gc 451 452 # Hessian evaluations. 453 cdp.h_count = hc 454 455 # Warning. 456 cdp.warning = warning

457 458

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

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

532 533

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

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

600 601

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

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

651 652

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

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

988 989

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

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

1201 1202

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

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

1271 1272

1273 - def _reset_min_stats(self):

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

1297 1298

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

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

1483 1484

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

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

1511 1512

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

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

Source Code for Module specific_analyses.model_free.mf_minimise