specific_analyses.relax_disp.parameters

1 ############################################################################### 2 # # 3 # Copyright (C) 2004-2014 Edward d'Auvergne # 4 # Copyright (C) 2009 Sebastien Morin # 5 # Copyright (C) 2014 Troels E. Linnet # 6 # # 7 # This file is part of the program relax (http://www.nmr-relax.com). # 8 # # 9 # This program is free software: you can redistribute it and/or modify # 10 # it under the terms of the GNU General Public License as published by # 11 # the Free Software Foundation, either version 3 of the License, or # 12 # (at your option) any later version. # 13 # # 14 # This program is distributed in the hope that it will be useful, # 15 # but WITHOUT ANY WARRANTY; without even the implied warranty of # 16 # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the # 17 # GNU General Public License for more details. # 18 # # 19 # You should have received a copy of the GNU General Public License # 20 # along with this program. If not, see <http://www.gnu.org/licenses/>. # 21 # # 22 ############################################################################### 23 24 # Module docstring. 25 """Functions relating to the parameters of the relaxation dispersion models.""" 26 27 # Python module imports. 28 from copy import deepcopy 29 from numpy import array, float64, median, zeros 30 import sys 31 32 # relax module imports. 33 from lib.dispersion.variables import MODEL_LIST_MMQ, MODEL_M61B, MODEL_NS_MMQ_3SITE, MODEL_NS_MMQ_3SITE_LINEAR, MODEL_NS_R1RHO_3SITE, MODEL_NS_R1RHO_3SITE_LINEAR, PARAMS_R20 34 from lib.errors import RelaxError, RelaxNoSequenceError 35 from lib.text.sectioning import subsection 36 from pipe_control import pipes 37 from pipe_control.mol_res_spin import exists_mol_res_spin_data, return_spin, spin_loop 38 from pipe_control.pipes import check_pipe 39 from specific_analyses.relax_disp.data import count_spins, generate_r20_key, has_exponential_exp_type, is_r1_optimised, loop_cluster, loop_exp_frq 40 41

42 -def assemble_param_vector(spins=None, key=None, sim_index=None):

43 """Assemble the dispersion relaxation dispersion curve fitting parameter vector. 44 45 @keyword spins: The list of spin data containers for the block. 46 @type spins: list of SpinContainer instances 47 @keyword key: The key for the R2eff and I0 parameters. 48 @type key: str or None 49 @keyword sim_index: The optional MC simulation index. 50 @type sim_index: int 51 @return: An array of the parameter values of the dispersion relaxation model. 52 @rtype: numpy float array 53 """ 54 55 # Initialise. 56 param_vector = [] 57 58 # Loop over the parameters of the cluster. 59 for param_name, param_index, spin_index, r20_key in loop_parameters(spins=spins): 60 # Get the value. 61 value = get_value(key=key, spins=spins, sim_index=sim_index, param_name=param_name, spin_index=spin_index, r20_key=r20_key) 62 63 # Add to the vector. 64 param_vector.append(value) 65 66 # Convert all None values to 0.0. 67 for i in range(len(param_vector)): 68 if param_vector[i] == None: 69 param_vector[i] = 0.0 70 71 # Return a numpy array. 72 return array(param_vector, float64)

73 74

75 -def copy(pipe_from=None, pipe_to=None):

76 """Copy dispersion parameters from one data pipe to another, taking the median of previous values to a start value for clusters. 77 Taking the median prevent averaging extreme outliers. 78 79 @param pipe_from: The data pipe to copy the value from. This defaults to the current data pipe. 80 @type pipe_from: str 81 @param pipe_to: The data pipe to copy the value to. This defaults to the current data pipe. 82 @type pipe_to: str 83 """ 84 85 # The current data pipe. 86 pipe_orig = pipes.cdp_name() 87 if pipe_from == None: 88 pipe_from = pipe_orig 89 if pipe_to == None: 90 pipe_to = pipe_orig 91 92 # Test that the pipes exist. 93 check_pipe(pipe_from) 94 check_pipe(pipe_to) 95 96 # Test that the pipes are not the same. 97 if pipe_from == pipe_to: 98 raise RelaxError("The source and destination pipes cannot be the same.") 99 100 # Test if the sequence data for pipe_from is loaded. 101 if not exists_mol_res_spin_data(pipe_from): 102 raise RelaxNoSequenceError(pipe_from) 103 104 # Test if the sequence data for pipe_to is loaded. 105 if not exists_mol_res_spin_data(pipe_to): 106 raise RelaxNoSequenceError(pipe_to) 107 108 # Switch to the destination data pipe. 109 pipes.switch(pipe_to) 110 111 # Loop over the clusters. 112 for spin_ids in loop_cluster(): 113 # Initialise some variables. 114 model = None 115 pA = [] 116 pB = [] 117 pC = [] 118 kex = [] 119 kex_AB = [] 120 kex_AC = [] 121 kex_BC = [] 122 k_AB = [] 123 kB = [] 124 kC = [] 125 tex = [] 126 count = 0 127 spins_from = [] 128 spins_to = [] 129 selected_cluster = False 130 131 # Loop over the spins, adding parameters to a list, which in the end will be used to find the median. 132 for id in spin_ids: 133 # Get the spins, then store them. 134 spin_from = return_spin(id, pipe=pipe_from) 135 spin_to = return_spin(id, pipe=pipe_to) 136 spins_from.append(spin_from) 137 spins_to.append(spin_to) 138 139 # Skip deselected spins. 140 if not spin_from.select or not spin_to.select: 141 continue 142 143 # The first printout. 144 if not selected_cluster: 145 subsection(file=sys.stdout, text="Copying parameters for the spin block %s"%spin_ids, prespace=2) 146 147 # Change the cluster selection flag. 148 selected_cluster = True 149 150 # The model. 151 if not model: 152 model = spin_from.model 153 154 # Check that the models match for all spins of the cluster. 155 if spin_from.model != model: 156 raise RelaxError("The model '%s' of spin '%s' from the source data pipe does not match the '%s' model of previous spins of the cluster." % (spin_from.model, id, model)) 157 if spin_to.model != model: 158 raise RelaxError("The model '%s' of spin '%s' from the destination data pipe does not match the '%s' model of previous spins of the cluster." % (spin_from.model, id, model)) 159 160 # Sum the source parameters. 161 if 'pA' in spin_from.params: 162 pA.append(spin_from.pA) 163 if 'pB' in spin_from.params: 164 pB.append(spin_from.pB) 165 if 'pC' in spin_from.params: 166 pC.append(spin_from.pC) 167 if 'kex' in spin_from.params: 168 kex.append(spin_from.kex) 169 if 'kex_AB' in spin_from.params: 170 kex_AB.append(spin_from.kex_AB) 171 if 'kex_AC' in spin_from.params: 172 kex_AC.append(spin_from.kex_AC) 173 if 'kex_BC' in spin_from.params: 174 kex_BC.append(spin_from.kex_BC) 175 if 'k_AB' in spin_from.params: 176 k_AB.append(spin_from.k_AB) 177 if 'kB' in spin_from.params: 178 kB.append(spin_from.kB) 179 if 'kC' in spin_from.params: 180 kC.append(spin_from.kC) 181 if 'tex' in spin_from.params: 182 tex.append(spin_from.tex) 183 184 # Increment the spin count. 185 count += 1 186 187 # The cluster is not selected, so move to the next. 188 if not selected_cluster: 189 continue 190 191 # Take median of parameters. 192 if len(pA) > 0: 193 pA = median(pA) 194 print("Median pA value: %.15f" % pA) 195 if len(pB) > 0: 196 pB = median(pB) 197 print("Median pB value: %.15f" % pB) 198 if len(pC) > 0: 199 pC = median(pC) 200 print("Median pC value: %.15f" % pC) 201 if len(kex) > 0: 202 kex = median(kex) 203 print("Median kex value: %.15f" % kex) 204 if len(kex_AB) > 0: 205 kex_AB = median(kex_AB) 206 print("Median k_AB value: %.15f" % kex_AB) 207 if len(kex_AC) > 0: 208 kex_AC = median(kex_AC) 209 print("Median k_AC value: %.15f" % kex_AC) 210 if len(kex_BC) > 0: 211 kex_BC = median(kex_BC) 212 print("Median k_BC value: %.15f" % kex_BC) 213 if len(k_AB) > 0: 214 k_AB = median(k_AB) 215 print("Median k_AB value: %.15f" % k_AB) 216 if len(kB) > 0: 217 kB = median(kB) 218 print("Median kB value: %.15f" % kB) 219 if len(kC) > 0: 220 kC = median(kC) 221 print("Median kC value: %.15f" % kC) 222 if len(tex) > 0: 223 tex = median(tex) 224 print("Median tex value: %.15f" % tex) 225 226 # Loop over the spins, this time copying the parameters. 227 for i in range(len(spin_ids)): 228 # Alias the containers. 229 spin_from = spins_from[i] 230 spin_to = spins_to[i] 231 232 # Skip deselected spins. 233 if not spin_from.select or not spin_to.select: 234 continue 235 236 # The R20 parameters. 237 if 'r2' in spin_from.params: 238 spin_to.r2 = deepcopy(spin_from.r2) 239 if 'r2a' in spin_from.params: 240 spin_to.r2a = deepcopy(spin_from.r2a) 241 if 'r2b' in spin_from.params: 242 spin_to.r2b = deepcopy(spin_from.r2b) 243 244 # The median parameters. 245 if 'pB' in spin_from.params and 'pC' not in spin_from.params: 246 spin_to.pA = pA 247 spin_to.pB = pB 248 spin_to.pC = 1.0 - pA - pB 249 elif 'pA' in spin_from.params: 250 spin_to.pA = pA 251 spin_to.pB = 1.0 - pA 252 if 'kex' in spin_from.params: 253 spin_to.kex = kex 254 if 'kex_AB' in spin_from.params: 255 spin_to.kex_AB = kex_AB 256 if 'kex_AC' in spin_from.params: 257 spin_to.kex_AC = kex_AC 258 if 'kex_BC' in spin_from.params: 259 spin_to.kex_BC = kex_BC 260 if 'k_AB' in spin_from.params: 261 spin_to.k_AB = k_AB 262 if 'kB' in spin_from.params: 263 spin_to.kB = kB 264 if 'kC' in spin_from.params: 265 spin_to.kC = kC 266 if 'tex' in spin_from.params: 267 spin_to.tex = tex 268 269 # All other spin specific parameters. 270 for param in spin_from.params: 271 if param in ['r2', 'r2a', 'r2b', 'pA', 'pB', 'pC', 'kex', 'kex_AB', 'kex_AC', 'kex_BC', 'k_AB', 'kB', 'kC', 'tex']: 272 continue 273 274 # Copy the value. 275 setattr(spin_to, param, deepcopy(getattr(spin_from, param))) 276 277 # Switch back to the original data pipe. 278 pipes.switch(pipe_orig)

279 280

281 -def disassemble_param_vector(param_vector=None, key=None, spins=None, sim_index=None):

282 """Disassemble the parameter vector. 283 284 @keyword param_vector: The parameter vector. 285 @type param_vector: numpy array 286 @keyword key: The key for the R2eff and I0 parameters. 287 @type key: str or None 288 @keyword spins: The list of spin data containers for the block. 289 @type spins: list of SpinContainer instances 290 @keyword sim_index: The optional MC simulation index. 291 @type sim_index: int 292 """ 293 294 # Initialise parameters if needed. 295 for spin in spins: 296 # Skip deselected spins. 297 if not spin.select: 298 continue 299 300 # The R2 parameter. 301 if 'r2' in spin.params: 302 if sim_index != None: 303 spin.r2_sim[sim_index] = {} 304 else: 305 spin.r2 = {} 306 307 # The R2A parameter. 308 if 'r2a' in spin.params: 309 if sim_index != None: 310 spin.r2a_sim[sim_index] = {} 311 else: 312 spin.r2a = {} 313 314 # The R2B parameter. 315 if 'r2b' in spin.params: 316 if sim_index != None: 317 spin.r2b_sim[sim_index] = {} 318 else: 319 spin.r2b = {} 320 321 # Loop over the parameters of the cluster, setting the values. 322 for param_name, param_index, spin_index, r20_key in loop_parameters(spins=spins): 323 set_value(value=param_vector[param_index], key=key, spins=spins, sim_index=sim_index, param_name=param_name, spin_index=spin_index, r20_key=r20_key)

324 325

326 -def get_param_names(spins=None, full=False):

327 """Generate a list of dispersion parameter names for the given spins. 328 329 @keyword spins: The list of spin data containers for the block. 330 @type spins: list of SpinContainer instances 331 @keyword full: A flag which if True will add the R2 key to the parameter names. 332 @type full: bool 333 """ 334 335 # Initialise the structure. 336 names = [] 337 338 # Loop over the parameters. 339 for param_name, param_index, spin_index, r20_key in loop_parameters(spins=spins): 340 # Set the initial text. 341 param_text = param_name 342 343 # The parameters with additional details. 344 if full and param_name in PARAMS_R20: 345 param_text += " (%s)" % r20_key 346 347 # Append the text. 348 names.append(param_text) 349 350 # Return the structure. 351 return names

352 353

354 -def get_value(key=None, spins=None, sim_index=None, param_name=None, spin_index=None, r20_key=None):

355 """Return the value for the given parameter. 356 357 @keyword key: The key for the R2eff and I0 parameters. 358 @type key: str or None 359 @keyword spins: The list of spin data containers for the block. 360 @type spins: list of SpinContainer instances 361 @keyword sim_index: The optional MC simulation index. 362 @type sim_index: int 363 @keyword param_name: The parameter name. 364 @type param_name: str 365 @keyword spin_index: The spin index (for the cluster). 366 @type spin_index: int 367 @keyword r20_key: The unique R20 parameter key. 368 @type r20_key: str 369 @return: The parameter value. 370 @rtype: float 371 """ 372 373 # Default value of None. 374 value = None 375 376 # Spin specific parameters. 377 if spin_index != None: 378 # Get the simulation value. 379 if sim_index != None: 380 # Get the simulation object. 381 sim_obj = getattr(spins[spin_index], param_name+'_sim') 382 383 # R20 parameter. 384 if r20_key != None: 385 if r20_key in sim_obj[sim_index]: 386 value = sim_obj[sim_index][r20_key] 387 388 # All other parameters. 389 else: 390 value = sim_obj[sim_index] 391 392 # Get the normal value. 393 else: 394 # Get the object. 395 obj = getattr(spins[spin_index], param_name) 396 397 # R20 parameter. 398 if r20_key != None: 399 if r20_key in obj: 400 value = obj[r20_key] 401 402 # All other parameters. 403 else: 404 value = getattr(spins[spin_index], param_name) 405 406 # Cluster specific parameters - use the parameter value from the first spin. 407 else: 408 # Set the simulation value. 409 if sim_index != None: 410 value = getattr(spins[0], param_name+'_sim')[sim_index] 411 412 # Set the normal value. 413 else: 414 value = getattr(spins[0], param_name) 415 416 # The R2eff model parameters. 417 if key != None: 418 if not key in value: 419 value = 0.0 420 else: 421 value = value[key] 422 423 # Return the value. 424 return value

425 426

427 -def linear_constraints(spins=None, scaling_matrix=None):

428 """Set up the relaxation dispersion curve fitting linear constraint matrices A and b. 429 430 Standard notation 431 ================= 432 433 The different constraints used within different models are:: 434 435 0 <= R1 <= 200 436 0 <= R2 <= 200 437 0 <= R2A <= 200 438 0 <= R2B <= 200 439 pC <= pB <= pA <= 1 440 pA >= 0.85 (the skewed condition, pA >> pB) 441 pB >= 0 442 pC >= 0 443 phi_ex >= 0 444 phi_ex_B >= 0 445 phi_ex_C >= 0 446 padw2 >= 0 447 dw >= 0 448 0 <= kex <= 2e6 449 0 <= kB <= 2e6 450 0 <= kC <= 2e6 451 0 <= k_AB <= 100 452 tex >= 0 453 454 455 Matrix notation 456 =============== 457 458 In the notation A.x >= b, where A is a matrix of coefficients, x is an array of parameter values, and b is a vector of scalars, these inequality constraints are:: 459 460 | 1 0 0 | | R1 | | 0 | 461 | | | | | | 462 |-1 0 0 | | R1 | | -200 | 463 | | | | | | 464 | 1 0 0 | | R2 | | 0 | 465 | | | | | | 466 |-1 0 0 | | R2 | | -200 | 467 | | | | | | 468 | 1 0 0 | | R2A | | 0 | 469 | | | | | | 470 |-1 0 0 | | R2A | | -200 | 471 | | | | | | 472 | 1 0 0 | | R2B | | 0 | 473 | | | | | | 474 |-1 0 0 | | R2B | | -200 | 475 | | | | | | 476 | 1 0 0 | | pA | | 0.5 | 477 | | | | | | 478 |-1 0 0 | | pA | | -1 | 479 | | | | | | 480 | 1 0 0 | | pA | | 0.85 | 481 | | | | | | 482 | 1 0 0 | | pB | | 0 | 483 | | | | | | 484 | 1 0 0 | | phi_ex | | 0 | 485 | | | | | | 486 | 1 0 0 | . | phi_ex_B | >= | 0 | 487 | | | | | | 488 | 1 0 0 | | phi_ex_C | | 0 | 489 | | | | | | 490 | 1 0 0 | | padw2 | | 0 | 491 | | | | | | 492 | 1 0 0 | | dw | | 0 | 493 | | | | | | 494 | 1 0 0 | | kex | | 0 | 495 | | | | | | 496 |-1 0 0 | | kex | | -2e6 | 497 | | | | | | 498 | 1 0 0 | | kB | | 0 | 499 | | | | | | 500 |-1 0 0 | | kB | | -2e6 | 501 | | | | | | 502 | 1 0 0 | | kC | | 0 | 503 | | | | | | 504 |-1 0 0 | | kC | | -2e6 | 505 | | | | | | 506 | 1 0 0 | | k_AB | | 0 | 507 | | | | | | 508 |-1 0 0 | | k_AB | | 100 | 509 | | | | | | 510 | 1 0 0 | | tex | | 0 | 511 512 513 @keyword spins: The list of spin data containers for the block. 514 @type spins: list of SpinContainer instances 515 @keyword scaling_matrix: The diagonal, square scaling matrix. 516 @type scaling_matrix: numpy diagonal matrix 517 @return: The matrices A and b. 518 @rtype: numpy rank-2 NxM array, numpy rank-1 N array 519 """ 520 521 # Initialisation (0..j..m). 522 A = [] 523 b = [] 524 n = param_num(spins=spins) 525 zero_array = zeros(n, float64) 526 j = 0 527 528 # Loop over the parameters of the cluster. 529 for param_name, param_index, spin_index, r20_key in loop_parameters(spins=spins): 530 # Effective transversal relaxation rate. 531 if param_name == 'r2eff': 532 A.append(zero_array * 0.0) 533 A.append(zero_array * 0.0) 534 A[j][param_index] = 1.0 535 A[j+1][param_index] = -1.0 536 b.append(0.0) 537 b.append(-200.0 / scaling_matrix[param_index, param_index]) 538 j += 2 539 540 # Initial intensity. 541 elif param_name == 'i0': 542 A.append(zero_array * 0.0) 543 A[j][param_index] = 1.0 544 b.append(0.0) 545 j += 1 546 547 # The longitudinal relaxation rates (0 <= r1 <= 200). 548 elif param_name in ['r1']: 549 A.append(zero_array * 0.0) 550 A.append(zero_array * 0.0) 551 A[j][param_index] = 1.0 552 A[j+1][param_index] = -1.0 553 b.append(0.0) 554 b.append(-200.0 / scaling_matrix[param_index, param_index]) 555 j += 2 556 557 # The transversal relaxation rates (0 <= r2 <= 200). 558 elif param_name in PARAMS_R20: 559 A.append(zero_array * 0.0) 560 A.append(zero_array * 0.0) 561 A[j][param_index] = 1.0 562 A[j+1][param_index] = -1.0 563 b.append(0.0) 564 b.append(-200.0 / scaling_matrix[param_index, param_index]) 565 j += 2 566 567 # The pA.pB.dw**2, phi_ex_B, phi_ex_C and pA.dw**2 parameters (phi_ex* >= 0 and padw2 >= 0). 568 elif param_name in ['phi_ex', 'phi_ex_B', 'phi_ex_C', 'padw2']: 569 A.append(zero_array * 0.0) 570 A[j][param_index] = 1.0 571 b.append(0.0) 572 j += 1 573 574 # Chemical exchange difference (dw >= 0). 575 elif param_name in ['dw', 'dw_AB', 'dw_AC', 'dw_BC', 'dwH', 'dwH_AB', 'dwH_AC', 'dwH_BC']: 576 if not spins[0].model in MODEL_LIST_MMQ + [MODEL_NS_R1RHO_3SITE, MODEL_NS_R1RHO_3SITE_LINEAR]: 577 A.append(zero_array * 0.0) 578 A[j][param_index] = 1.0 579 b.append(0.0) 580 j += 1 581 582 # The population of state A and B. 583 elif param_name == 'pA': 584 # First the pA <= 1 constraint (which also rearranged is pB >= 0). 585 A.append(zero_array * 0.0) 586 A[j][param_index] = -1.0 587 b.append(-1.0 / scaling_matrix[param_index, param_index]) 588 j += 1 589 590 # The skewed condition (pA >> pB). 591 if spins[0].model == MODEL_M61B: 592 A.append(zero_array * 0.0) 593 A[j][param_index] = 1.0 594 b.append(0.85 / scaling_matrix[param_index, param_index]) 595 j += 1 596 597 # Otherwise use the pA >= pB constraint. 598 else: 599 A.append(zero_array * 0.0) 600 A[j][param_index] = 1.0 601 b.append(0.5 / scaling_matrix[param_index, param_index]) 602 j += 1 603 604 # The population of state C (the pB parameter is only present when pC exists). 605 elif param_name == 'pB': 606 # Find pA. 607 for param_name2, param_index2, spin_index2, r20_key2 in loop_parameters(spins=spins): 608 if param_name2 == 'pA': 609 # First pC >= 0 (which rearranged is -pA - pB >= -1). 610 A.append(zero_array * 0.0) 611 A[j][param_index2] = -1.0 612 A[j][param_index] = -1.0 613 b.append(-1.0 / scaling_matrix[param_index, param_index]) 614 j += 1 615 616 # Then the pA >= pC constraint. 617 A.append(zero_array * 0.0) 618 A[j][param_index2] = 1.0 619 A[j][param_index] = -1.0 620 b.append(0.0) 621 j += 1 622 break 623 624 # Exchange rates and times (0 <= k <= 2e6). 625 elif param_name in ['kex', 'kex_AB', 'kex_AC', 'kex_BC', 'kB', 'kC']: 626 A.append(zero_array * 0.0) 627 A.append(zero_array * 0.0) 628 A[j][param_index] = 1.0 629 A[j+1][param_index] = -1.0 630 b.append(0.0) 631 b.append(-2e6 / scaling_matrix[param_index, param_index]) 632 j += 2 633 634 # Exchange rates and times (0 <= k_AB <= 100). 635 elif param_name in ['k_AB']: 636 A.append(zero_array * 0.0) 637 A.append(zero_array * 0.0) 638 A[j][param_index] = 1.0 639 A[j+1][param_index] = -1.0 640 b.append(0.0) 641 b.append(-100. / scaling_matrix[param_index, param_index]) 642 j += 2 643 644 # Exchange times (tex >= 0). 645 elif param_name in ['tex']: 646 A.append(zero_array * 0.0) 647 A[j][param_index] = 1.0 648 b.append(0.0) 649 j += 1 650 651 # Convert to numpy data structures. 652 A = array(A, float64) 653 b = array(b, float64) 654 655 # Return the matrices. 656 return A, b

657 658

659 -def loop_parameters(spins=None):

660 """Generator function for looping of the parameters of the cluster. 661 662 @keyword spins: The list of spin data containers for the block. 663 @type spins: list of SpinContainer instances 664 @return: The parameter name, the parameter index (for the parameter vector), the spin index (for the cluster), and the R20 parameter key (for R20, R20A, and R20B parameters stored as dictionaries). 665 @rtype: str, int, int, str 666 """ 667 668 # Make sure that the R1 parameter is correctly set up. 669 r1_setup() 670 671 # The parameter index. 672 param_index = -1 673 674 # The R2eff model. 675 if cdp.model_type == 'R2eff': 676 # Loop over the spins. 677 for spin_index in range(len(spins)): 678 # Skip deselected spins. 679 if not spins[spin_index].select: 680 continue 681 682 # Yield the two parameters. 683 params = ['r2eff', 'i0'] 684 for i in range(2): 685 # First increment the indices. 686 param_index += 1 687 688 # Yield the data. 689 yield params[i], param_index, spin_index, None 690 691 # All other models. 692 else: 693 # First the R1 fit parameter (one per spin per field strength). 694 for spin_index in range(len(spins)): 695 # Skip deselected spins. 696 if not spins[spin_index].select: 697 continue 698 699 # The R1 parameter. 700 if 'r1' in spins[spin_index].params: 701 for exp_type, frq in loop_exp_frq(): 702 param_index += 1 703 yield 'r1', param_index, spin_index, generate_r20_key(exp_type=exp_type, frq=frq) 704 705 # Then the R2 parameters (one per spin per field strength). 706 for spin_index in range(len(spins)): 707 # Skip deselected spins. 708 if not spins[spin_index].select: 709 continue 710 711 # The R2 parameter. 712 if 'r2' in spins[spin_index].params: 713 for exp_type, frq in loop_exp_frq(): 714 param_index += 1 715 yield 'r2', param_index, spin_index, generate_r20_key(exp_type=exp_type, frq=frq) 716 717 # The R2A parameter. 718 if 'r2a' in spins[spin_index].params: 719 for exp_type, frq in loop_exp_frq(): 720 param_index += 1 721 yield 'r2a', param_index, spin_index, generate_r20_key(exp_type=exp_type, frq=frq) 722 723 # The R2B parameter. 724 if 'r2b' in spins[spin_index].params: 725 for exp_type, frq in loop_exp_frq(): 726 param_index += 1 727 yield 'r2b', param_index, spin_index, generate_r20_key(exp_type=exp_type, frq=frq) 728 729 # Then the chemical shift difference parameters 'phi_ex', 'phi_ex_B', 'phi_ex_C', 'padw2', 'dw', 'dw_AB', 'dw_BC', 'dw_AB' (one per spin). 730 for spin_index in range(len(spins)): 731 # Skip deselected spins. 732 if not spins[spin_index].select: 733 continue 734 735 # Yield the data. 736 if 'phi_ex' in spins[spin_index].params: 737 param_index += 1 738 yield 'phi_ex', param_index, spin_index, None 739 if 'phi_ex_B' in spins[spin_index].params: 740 param_index += 1 741 yield 'phi_ex_B', param_index, spin_index, None 742 if 'phi_ex_C' in spins[spin_index].params: 743 param_index += 1 744 yield 'phi_ex_C', param_index, spin_index, None 745 if 'padw2' in spins[spin_index].params: 746 param_index += 1 747 yield 'padw2', param_index, spin_index, None 748 if 'dw' in spins[spin_index].params: 749 param_index += 1 750 yield 'dw', param_index, spin_index, None 751 if 'dw_AB' in spins[spin_index].params: 752 param_index += 1 753 yield 'dw_AB', param_index, spin_index, None 754 if 'dw_BC' in spins[spin_index].params: 755 param_index += 1 756 yield 'dw_BC', param_index, spin_index, None 757 if 'dw_AC' in spins[spin_index].params: 758 param_index += 1 759 yield 'dw_AC', param_index, spin_index, None 760 761 # Then a separate block for the proton chemical shift difference parameters for the MQ models (one per spin). 762 for spin_index in range(len(spins)): 763 # Skip deselected spins. 764 if not spins[spin_index].select: 765 continue 766 767 if 'dwH' in spins[spin_index].params: 768 param_index += 1 769 yield 'dwH', param_index, spin_index, None 770 if 'dwH_AB' in spins[spin_index].params: 771 param_index += 1 772 yield 'dwH_AB', param_index, spin_index, None 773 if 'dwH_BC' in spins[spin_index].params: 774 param_index += 1 775 yield 'dwH_BC', param_index, spin_index, None 776 if 'dwH_AC' in spins[spin_index].params: 777 param_index += 1 778 yield 'dwH_AC', param_index, spin_index, None 779 780 # All other parameters (one per spin cluster). 781 for param in spins[0].params: 782 if not param in ['r1', 'r2', 'r2a', 'r2b', 'phi_ex', 'phi_ex_B', 'phi_ex_C', 'padw2', 'dw', 'dw_AB', 'dw_BC', 'dw_AB', 'dwH', 'dwH_AB', 'dwH_BC', 'dwH_AB']: 783 param_index += 1 784 yield param, param_index, None, None

785 786

787 -def param_conversion(key=None, spins=None, sim_index=None):

788 """Convert Disassemble the parameter vector. 789 790 @keyword key: The key for the R2eff and I0 parameters. 791 @type key: str or None 792 @keyword spins: The list of spin data containers for the block. 793 @type spins: list of SpinContainer instances 794 @keyword sim_index: The optional MC simulation index. 795 @type sim_index: int 796 """ 797 798 # Loop over the parameters of the cluster. 799 for param_name, param_index, spin_index, r20_key in loop_parameters(spins=spins): 800 # Get the value. 801 value = get_value(key=key, spins=spins, sim_index=sim_index, param_name=param_name, spin_index=spin_index, r20_key=r20_key) 802 803 # The pA to pB to pC conversion. 804 if param_name == 'pA': 805 # 3-site exchange. 806 if spins[0].model in [MODEL_NS_MMQ_3SITE, MODEL_NS_MMQ_3SITE_LINEAR]: 807 # Get the pB value. 808 pB = get_value(key=key, spins=spins, sim_index=sim_index, param_name='pB', spin_index=spin_index, r20_key=r20_key) 809 810 # Set the pC value. 811 pC = 1.0 - value - pB 812 set_value(value=pC, key=key, spins=spins, sim_index=sim_index, param_name='pC', spin_index=spin_index) 813 814 # 2-site exchange. 815 else: 816 pB = 1.0 - value 817 set_value(value=pB, key=key, spins=spins, sim_index=sim_index, param_name='pB', spin_index=spin_index) 818 819 # The kex to tex conversion. 820 if param_name == 'kex': 821 tex = 1.0 / value 822 set_value(value=tex, key=key, spins=spins, sim_index=sim_index, param_name='tex', spin_index=spin_index) 823 824 # The kex to k_AB and k_BA conversion. 825 if param_name == 'kex' and 'pA' in spins[0].params: 826 # Get pA value. 827 pA = get_value(key=key, spins=spins, sim_index=sim_index, param_name='pA', spin_index=spin_index) 828 829 # Calculate k_AB value and set it. 830 k_AB = value * (1.0 - pA) 831 set_value(value=k_AB, key=key, spins=spins, sim_index=sim_index, param_name='k_AB', spin_index=spin_index) 832 833 # Calculate k_BA value and set it. 834 k_BA = value * pA 835 set_value(value=k_BA, key=key, spins=spins, sim_index=sim_index, param_name='k_BA', spin_index=spin_index) 836 837 # The tex to kex conversion. 838 if param_name == 'tex': 839 kex = 1.0 / value 840 set_value(value=kex, key=key, spins=spins, sim_index=sim_index, param_name='kex', spin_index=spin_index)

841 842

843 -def param_index_to_param_info(index=None, spins=None):

844 """Convert the given parameter array index to parameter identifying information. 845 846 The parameter index will be converted to the parameter name, the relevant spin index in the cluster, and relevant exponential curve key. 847 848 @keyword index: The index of the parameter array. 849 @type index: int 850 @keyword spins: The list of spin data containers for the block. 851 @type spins: list of SpinContainer instances 852 @return: The parameter name, the spin index (for the cluster), and the frequency index (for parameters with different values per spectrometer field strength). 853 @rtype: str, int, int 854 """ 855 856 # Loop over the parameters, yielding when a match is found. 857 for param_name, param_index, spin_index, r20_key in loop_parameters(spins=spins): 858 if param_index == index: 859 return param_name, spin_index, r20_key

860 861

862 -def param_num(spins=None):

863 """Determine the number of parameters in the model. 864 865 @keyword spins: The list of spin data containers for the block. 866 @type spins: list of SpinContainer instances 867 @return: The number of model parameters. 868 @rtype: int 869 """ 870 871 # Initialise the number. 872 num = 0 873 874 # The R2eff model. 875 if cdp.model_type == 'R2eff': 876 # Count the selected spins. 877 spin_num = count_spins(spins) 878 879 # Exponential curves (with clustering). 880 if has_exponential_exp_type(): 881 return 2 * spin_num 882 883 # Fixed time period experiments (with clustering). 884 return 1 * spin_num 885 886 # Check the spin cluster. 887 for spin in spins: 888 # Skip deselected spins. 889 if not spin.select: 890 continue 891 892 if len(spin.params) != len(spins[0].params): 893 raise RelaxError("The number of parameters for each spin in the cluster are not the same.") 894 895 # Count the number of R10 parameters. 896 for spin in spins: 897 # Skip deselected spins. 898 if not spin.select: 899 continue 900 901 for i in range(len(spin.params)): 902 if spin.params[i] in ['r1']: 903 for exp_type, frq in loop_exp_frq(): 904 num += 1 905 906 # Count the number of R20 parameters. 907 for spin in spins: 908 # Skip deselected spins. 909 if not spin.select: 910 continue 911 912 for i in range(len(spin.params)): 913 if spin.params[i] in PARAMS_R20: 914 for exp_type, frq in loop_exp_frq(): 915 num += 1 916 917 # Count the number of spin specific parameters for all spins. 918 spin_params = ['phi_ex', 'phi_ex_B', 'phi_ex_C', 'padw2', 'dw', 'dwH'] 919 for spin in spins: 920 # Skip deselected spins. 921 if not spin.select: 922 continue 923 924 for i in range(len(spin.params)): 925 if spin.params[i] in spin_params: 926 num += 1 927 928 # Count all other parameters, but only for a single spin. 929 all_params = ['r1'] + PARAMS_R20 + spin_params 930 for spin in spins: 931 # Skip deselected spins. 932 if not spin.select: 933 continue 934 935 for i in range(len(spin.params)): 936 if not spin.params[i] in all_params: 937 num += 1 938 break 939 940 # Return the number. 941 return num

942 943

944 -def r1_setup():

945 """Modify the current spin parameter list to either include or exclude the R1 parameter.""" 946 947 # Loop over the spins. 948 for spin, spin_id in spin_loop(return_id=True, skip_desel=True): 949 # No model set up. 950 if not hasattr(spin, 'params') or not hasattr(spin, 'model'): 951 continue 952 953 # Should R1 data be optimised? 954 r1_fit = is_r1_optimised(spin.model) 955 956 # Prepend R1 and add it to the spin container. 957 if r1_fit and 'r1' not in spin.params: 958 spin.params.insert(0, 'r1') 959 if r1_fit and not hasattr(spin, 'r1'): 960 spin.r1 = {} 961 962 # Remove the R1 parameter. 963 if not r1_fit and 'r1' in spin.params: 964 spin.params.remove('r1')

965 966

967 -def set_value(value=None, key=None, spins=None, sim_index=None, param_name=None, spin_index=None, r20_key=None):

968 """Return the value for the given parameter. 969 970 @keyword value: The parameter value to set. 971 @type value: float 972 @keyword key: The key for the R2eff and I0 parameters. 973 @type key: str or None 974 @keyword spins: The list of spin data containers for the block. 975 @type spins: list of SpinContainer instances 976 @keyword sim_index: The optional MC simulation index. 977 @type sim_index: int 978 @keyword param_name: The parameter name. 979 @type param_name: str 980 @keyword spin_index: The spin index (for the cluster). 981 @type spin_index: int 982 @keyword r20_key: The unique R20 parameter key. 983 @type r20_key: str 984 """ 985 986 # Spin specific parameters. 987 if spin_index != None: 988 # Set the simulation value. 989 if sim_index != None: 990 # Get the simulation object. 991 obj = getattr(spins[spin_index], param_name+'_sim') 992 993 # R20 parameter. 994 if r20_key != None: 995 obj[sim_index][r20_key] = value 996 997 # All other parameters. 998 else: 999 if key != None: 1000 obj[sim_index][key] = value 1001 else: 1002 obj[sim_index] = value 1003 1004 # Set the normal value. 1005 else: 1006 # Get the object. 1007 obj = getattr(spins[spin_index], param_name) 1008 1009 # R20 parameter. 1010 if r20_key != None: 1011 obj[r20_key] = value 1012 1013 # Set the normal value. 1014 else: 1015 if key != None: 1016 obj[key] = value 1017 else: 1018 setattr(spins[spin_index], param_name, value) 1019 1020 # Cluster specific parameters. 1021 else: 1022 # Set the same parameter value for all spins in the cluster. 1023 for spin in spins: 1024 # Skip deselected spins. 1025 if not spin.select: 1026 continue 1027 1028 # Set the simulation value. 1029 if sim_index != None: 1030 sim_obj = getattr(spin, param_name+'_sim') 1031 sim_obj[sim_index] = value 1032 1033 # Set the normal value. 1034 else: 1035 setattr(spin, param_name, value)

1036

Source Code for Module specific_analyses.relax_disp.parameters