specific_fns.model_free.mf

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

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

56 """Disassemble the model-free parameter vector. 57 58 @param model_type: The model-free model type. This must be one of 'mf', 'local_tm', 59 'diff', or 'all'. 60 @type model_type: str 61 @keyword param_vector: The model-free parameter vector. 62 @type param_vector: numpy array 63 @keyword spin: The spin data container. If this argument is supplied, then the spin_id 64 argument will be ignored. 65 @type spin: SpinContainer instance 66 @keyword spin_id: The spin identification string. 67 @type spin_id: str 68 @keyword sim_index: The optional MC simulation index. 69 @type sim_index: int 70 """ 71 72 # Initialise. 73 param_index = 0 74 75 # Diffusion tensor parameters of the Monte Carlo simulations. 76 if sim_index != None and (model_type == 'diff' or model_type == 'all'): 77 # Spherical diffusion. 78 if cdp.diff_tensor.type == 'sphere': 79 # Sim values. 80 cdp.diff_tensor.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 # Alias the current data pipe. 327 cdp = pipes.get_pipe() 328 329 # Catch infinite chi-squared values. 330 if isInf(func): 331 raise RelaxInfError('chi-squared') 332 333 # Catch chi-squared values of NaN. 334 if isNaN(func): 335 raise RelaxNaNError('chi-squared') 336 337 # Scaling. 338 if scaling: 339 param_vector = dot(scaling_matrix, param_vector) 340 341 # Check if the chi-squared value is lower. This allows for a parallelised grid search! 342 if sim_index == None: 343 # Get the correct value. 344 chi2 = None 345 if (model_type == 'mf' or model_type == 'local_tm') and hasattr(cdp, 'chi2'): 346 chi2 = spin.chi2 347 if (model_type == 'diff' or model_type == 'all') and hasattr(cdp, 'chi2'): 348 chi2 = cdp.chi2 349 350 # Spin text. 351 spin_text = '' 352 if spin != None and hasattr(spin, '_spin_ids') and len(spin._spin_ids): 353 spin_text = " for the spin '%s'" % spin._spin_ids[0] 354 355 # No improvement. 356 if chi2 != None and func >= chi2: 357 print("Discarding the optimisation results%s, the optimised chi-squared value is higher than the current value (%s >= %s)." % (spin_text, func, chi2)) 358 359 # Exit! 360 return 361 362 # New minimum. 363 else: 364 print("Storing the optimisation results%s, the optimised chi-squared value is lower than the current value (%s < %s)." % (spin_text, func, chi2)) 365 366 # Disassemble the parameter vector. 367 self._disassemble_param_vector(model_type, param_vector=param_vector, spin=spin, sim_index=sim_index) 368 369 # Monte Carlo minimisation statistics. 370 if sim_index != None: 371 # Sequence specific minimisation statistics. 372 if model_type == 'mf' or model_type == 'local_tm': 373 374 # Chi-squared statistic. 375 spin.chi2_sim[sim_index] = func 376 377 # Iterations. 378 spin.iter_sim[sim_index] = iter 379 380 # Function evaluations. 381 spin.f_count_sim[sim_index] = fc 382 383 # Gradient evaluations. 384 spin.g_count_sim[sim_index] = gc 385 386 # Hessian evaluations. 387 spin.h_count_sim[sim_index] = hc 388 389 # Warning. 390 spin.warning_sim[sim_index] = warning 391 392 # Global minimisation statistics. 393 elif model_type == 'diff' or model_type == 'all': 394 # Chi-squared statistic. 395 cdp.chi2_sim[sim_index] = func 396 397 # Iterations. 398 cdp.iter_sim[sim_index] = iter 399 400 # Function evaluations. 401 cdp.f_count_sim[sim_index] = fc 402 403 # Gradient evaluations. 404 cdp.g_count_sim[sim_index] = gc 405 406 # Hessian evaluations. 407 cdp.h_count_sim[sim_index] = hc 408 409 # Warning. 410 cdp.warning_sim[sim_index] = warning 411 412 # Normal statistics. 413 else: 414 # Sequence specific minimisation statistics. 415 if model_type == 'mf' or model_type == 'local_tm': 416 # Chi-squared statistic. 417 spin.chi2 = func 418 419 # Iterations. 420 spin.iter = iter 421 422 # Function evaluations. 423 spin.f_count = fc 424 425 # Gradient evaluations. 426 spin.g_count = gc 427 428 # Hessian evaluations. 429 spin.h_count = hc 430 431 # Warning. 432 spin.warning = warning 433 434 # Global minimisation statistics. 435 elif model_type == 'diff' or model_type == 'all': 436 # Chi-squared statistic. 437 cdp.chi2 = func 438 439 # Iterations. 440 cdp.iter = iter 441 442 # Function evaluations. 443 cdp.f_count = fc 444 445 # Gradient evaluations. 446 cdp.g_count = gc 447 448 # Hessian evaluations. 449 cdp.h_count = hc 450 451 # Warning. 452 cdp.warning = warning

453 454

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

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

528 529

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

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

596 597

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

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

647 648

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

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

984 985

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

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

1193 1194

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

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

1259 1260

1261 - def _reset_min_stats(self):

1262 """Reset all the minimisation statistics. 1263 1264 All global and spin specific values will be set to None. 1265 """ 1266 1267 # Global stats. 1268 if hasattr(cdp, 'chi2'): 1269 cdp.chi2 = None 1270 cdp.iter = None 1271 cdp.f_count = None 1272 cdp.g_count = None 1273 cdp.h_count = None 1274 cdp.warning = None 1275 1276 # Spin specific stats. 1277 for spin in spin_loop(): 1278 if hasattr(spin, 'chi2'): 1279 spin.chi2 = None 1280 spin.iter = None 1281 spin.f_count = None 1282 spin.g_count = None 1283 spin.h_count = None 1284 spin.warning = None

1285 1286

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

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

1471 1472

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

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

1499 1500

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

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

Source Code for Module specific_fns.model_free.mf_minimise