specific_analyses.relax_disp.parameters

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

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

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

71 72

73 -def assemble_scaling_matrix(spins=None, key=None, scaling=True):

74 """Create and return the scaling matrix. 75 76 @keyword spins: The list of spin data containers for the block. 77 @type spins: list of SpinContainer instances 78 @keyword key: The key for the R2eff and I0 parameters. 79 @type key: str or None 80 @keyword scaling: A flag which if False will cause the identity matrix to be returned. 81 @type scaling: bool 82 @return: The diagonal and square scaling matrix. 83 @rtype: numpy diagonal matrix 84 """ 85 86 # Initialise. 87 scaling_matrix = identity(param_num(spins=spins), float64) 88 i = 0 89 param_index = 0 90 91 # No diagonal scaling. 92 if not scaling: 93 return scaling_matrix 94 95 # Loop over the parameters of the cluster. 96 for param_name, param_index, spin_index, r20_key in loop_parameters(spins=spins): 97 # Transversal relaxation rate scaling. 98 if param_name in ['r2', 'r2a', 'r2b']: 99 scaling_matrix[param_index, param_index] = 10 100 101 # The pA.pB.dw**2, phi_ex_B, phi_ex_C and pA.dw**2 parameters. 102 elif param_name in ['phi_ex', 'phi_ex_B', 'phi_ex_C', 'padw2']: 103 scaling_matrix[param_index, param_index] = 1 104 105 # Chemical shift difference between states A and B scaling. 106 elif param_name in ['dw', 'dw_AB', 'dw_AC', 'dw_BC', 'dwH', 'dwH_AB', 'dwH_AC', 'dwH_BC']: 107 scaling_matrix[param_index, param_index] = 1 108 109 # The population of state X. 110 elif param_name in ['pA', 'pB', 'pC']: 111 scaling_matrix[param_index, param_index] = 1 112 113 # Exchange rate scaling. 114 elif param_name in ['kex', 'kex_AB', 'kex_AC', 'kex_BC', 'k_AB', 'kB', 'kC']: 115 scaling_matrix[param_index, param_index] = 10000 116 117 # Time of exchange scaling. 118 elif param_name == 'tex': 119 scaling_matrix[param_index, param_index] = 1e-4 120 121 # Return the scaling matrix. 122 return scaling_matrix

123 124

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

126 """Copy dispersion parameters from one data pipe to another, taking the median of previous values to a start value for clusters. 127 Taking the median prevent averaging extreme outliers. 128 129 @param pipe_from: The data pipe to copy the value from. This defaults to the current data pipe. 130 @type pipe_from: str 131 @param pipe_to: The data pipe to copy the value to. This defaults to the current data pipe. 132 @type pipe_to: str 133 """ 134 135 # The current data pipe. 136 pipe_orig = pipes.cdp_name() 137 if pipe_from == None: 138 pipe_from = pipe_orig 139 if pipe_to == None: 140 pipe_to = pipe_orig 141 142 # Test that the pipes exist. 143 pipes.test(pipe_from) 144 pipes.test(pipe_to) 145 146 # Test that the pipes are not the same. 147 if pipe_from == pipe_to: 148 raise RelaxError("The source and destination pipes cannot be the same.") 149 150 # Test if the sequence data for pipe_from is loaded. 151 if not exists_mol_res_spin_data(pipe_from): 152 raise RelaxNoSequenceError(pipe_from) 153 154 # Test if the sequence data for pipe_to is loaded. 155 if not exists_mol_res_spin_data(pipe_to): 156 raise RelaxNoSequenceError(pipe_to) 157 158 # Switch to the destination data pipe. 159 pipes.switch(pipe_to) 160 161 # Loop over the clusters. 162 for spin_ids in loop_cluster(): 163 # Initialise some variables. 164 model = None 165 pA = [] 166 pB = [] 167 pC = [] 168 kex = [] 169 kex_AB = [] 170 kex_AC = [] 171 kex_BC = [] 172 k_AB = [] 173 kB = [] 174 kC = [] 175 tex = [] 176 count = 0 177 spins_from = [] 178 spins_to = [] 179 selected_cluster = False 180 181 # Loop over the spins, adding parameters to a list, which in the end will be used to find the median. 182 for id in spin_ids: 183 # Get the spins, then store them. 184 spin_from = return_spin(id, pipe=pipe_from) 185 spin_to = return_spin(id, pipe=pipe_to) 186 spins_from.append(spin_from) 187 spins_to.append(spin_to) 188 189 # Skip deselected spins. 190 if not spin_from.select or not spin_to.select: 191 continue 192 193 # The first printout. 194 if not selected_cluster: 195 subsection(file=sys.stdout, text="Copying parameters for the spin block %s"%spin_ids, prespace=2) 196 197 # Change the cluster selection flag. 198 selected_cluster = True 199 200 # The model. 201 if not model: 202 model = spin_from.model 203 204 # Check that the models match for all spins of the cluster. 205 if spin_from.model != model: 206 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)) 207 if spin_to.model != model: 208 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)) 209 210 # Sum the source parameters. 211 if 'pA' in spin_from.params: 212 pA.append(spin_from.pA) 213 if 'pB' in spin_from.params: 214 pB.append(spin_from.pB) 215 if 'pC' in spin_from.params: 216 pC.append(spin_from.pC) 217 if 'kex' in spin_from.params: 218 kex.append(spin_from.kex) 219 if 'kex_AB' in spin_from.params: 220 kex_AB.append(spin_from.kex_AB) 221 if 'kex_AC' in spin_from.params: 222 kex_AC.append(spin_from.kex_AC) 223 if 'kex_BC' in spin_from.params: 224 kex_BC.append(spin_from.kex_BC) 225 if 'k_AB' in spin_from.params: 226 k_AB.append(spin_from.k_AB) 227 if 'kB' in spin_from.params: 228 kB.append(spin_from.kB) 229 if 'kC' in spin_from.params: 230 kC.append(spin_from.kC) 231 if 'tex' in spin_from.params: 232 tex.append(spin_from.tex) 233 234 # Increment the spin count. 235 count += 1 236 237 # The cluster is not selected, so move to the next. 238 if not selected_cluster: 239 continue 240 241 # Take median of parameters. 242 if len(pA) > 0: 243 pA = median(pA) 244 print("Median pA value: %.15f" % pA) 245 if len(pB) > 0: 246 pB = median(pB) 247 print("Median pB value: %.15f" % pB) 248 if len(pC) > 0: 249 pC = median(pC) 250 print("Median pC value: %.15f" % pC) 251 if len(kex) > 0: 252 kex = median(kex) 253 print("Median kex value: %.15f" % kex) 254 if len(kex_AB) > 0: 255 kex_AB = median(kex_AB) 256 print("Median k_AB value: %.15f" % kex_AB) 257 if len(kex_AC) > 0: 258 kex_AC = median(kex_AC) 259 print("Median k_AC value: %.15f" % kex_AC) 260 if len(kex_BC) > 0: 261 kex_BC = median(kex_BC) 262 print("Median k_BC value: %.15f" % kex_BC) 263 if len(k_AB) > 0: 264 k_AB = median(k_AB) 265 print("Median k_AB value: %.15f" % k_AB) 266 if len(kB) > 0: 267 kB = median(kB) 268 print("Median kB value: %.15f" % kB) 269 if len(kC) > 0: 270 kC = median(kC) 271 print("Median kC value: %.15f" % kC) 272 if len(tex) > 0: 273 tex = median(tex) 274 print("Median tex value: %.15f" % tex) 275 276 # Loop over the spins, this time copying the parameters. 277 for i in range(len(spin_ids)): 278 # Alias the containers. 279 spin_from = spins_from[i] 280 spin_to = spins_to[i] 281 282 # Skip deselected spins. 283 if not spin_from.select or not spin_to.select: 284 continue 285 286 # The R20 parameters. 287 if 'r2' in spin_from.params: 288 spin_to.r2 = deepcopy(spin_from.r2) 289 if 'r2a' in spin_from.params: 290 spin_to.r2a = deepcopy(spin_from.r2a) 291 if 'r2b' in spin_from.params: 292 spin_to.r2b = deepcopy(spin_from.r2b) 293 294 # The median parameters. 295 if 'pB' in spin_from.params and 'pC' not in spin_from.params: 296 spin_to.pA = pA 297 spin_to.pB = pB 298 spin_to.pC = 1.0 - pA - pB 299 elif 'pA' in spin_from.params: 300 spin_to.pA = pA 301 spin_to.pB = 1.0 - pA 302 if 'kex' in spin_from.params: 303 spin_to.kex = kex 304 if 'kex_AB' in spin_from.params: 305 spin_to.kex_AB = kex_AB 306 if 'kex_AC' in spin_from.params: 307 spin_to.kex_AC = kex_AC 308 if 'kex_BC' in spin_from.params: 309 spin_to.kex_BC = kex_BC 310 if 'k_AB' in spin_from.params: 311 spin_to.k_AB = k_AB 312 if 'kB' in spin_from.params: 313 spin_to.kB = kB 314 if 'kC' in spin_from.params: 315 spin_to.kC = kC 316 if 'tex' in spin_from.params: 317 spin_to.tex = tex 318 319 # All other spin specific parameters. 320 for param in spin_from.params: 321 if param in ['r2', 'pA', 'pB', 'pC', 'kex', 'kex_AB', 'kex_AC', 'kex_BC', 'k_AB', 'kB', 'kC', 'tex']: 322 continue 323 324 # Copy the value. 325 setattr(spin_to, param, deepcopy(getattr(spin_from, param))) 326 327 # Switch back to the original data pipe. 328 pipes.switch(pipe_orig)

329 330

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

332 """Disassemble the parameter vector. 333 334 @keyword param_vector: The parameter vector. 335 @type param_vector: numpy array 336 @keyword key: The key for the R2eff and I0 parameters. 337 @type key: str or None 338 @keyword spins: The list of spin data containers for the block. 339 @type spins: list of SpinContainer instances 340 @keyword sim_index: The optional MC simulation index. 341 @type sim_index: int 342 """ 343 344 # Initialise parameters if needed. 345 for spin in spins: 346 # Skip deselected spins. 347 if not spin.select: 348 continue 349 350 # The R2 parameter. 351 if 'r2' in spin.params: 352 if sim_index != None: 353 spin.r2_sim[sim_index] = {} 354 else: 355 spin.r2 = {} 356 357 # The R2A parameter. 358 if 'r2a' in spin.params: 359 if sim_index != None: 360 spin.r2a_sim[sim_index] = {} 361 else: 362 spin.r2a = {} 363 364 # The R2B parameter. 365 if 'r2b' in spin.params: 366 if sim_index != None: 367 spin.r2b_sim[sim_index] = {} 368 else: 369 spin.r2b = {} 370 371 # Loop over the parameters of the cluster, setting the values. 372 for param_name, param_index, spin_index, r20_key in loop_parameters(spins=spins): 373 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)

374 375

376 -def get_param_names(spins=None):

377 """Generate a list of dispersion parameter names for the given spins. 378 379 @keyword spins: The list of spin data containers for the block. 380 @type spins: list of SpinContainer instances 381 """ 382 383 # Initialise the structure. 384 names = [] 385 386 # Loop over the parameters. 387 for param_name, param_index, spin_index, r20_key in loop_parameters(spins=spins): 388 # Set the initial text. 389 param_text = param_name 390 391 # The parameters with additional details. 392 if param_name in ['r2', 'r2a', 'r2b']: 393 param_text += " (%s)" % r20_key 394 395 # Append the text. 396 names.append(param_text) 397 398 # Return the structure. 399 return names

400 401

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

403 """Return the value for the given parameter. 404 405 @keyword key: The key for the R2eff and I0 parameters. 406 @type key: str or None 407 @keyword spins: The list of spin data containers for the block. 408 @type spins: list of SpinContainer instances 409 @keyword sim_index: The optional MC simulation index. 410 @type sim_index: int 411 @keyword param_name: The parameter name. 412 @type param_name: str 413 @keyword spin_index: The spin index (for the cluster). 414 @type spin_index: int 415 @keyword r20_key: The unique R20 parameter key. 416 @type r20_key: str 417 @return: The parameter value. 418 @rtype: float 419 """ 420 421 # Default value of 0.0. 422 value = 0.0 423 424 # Spin specific parameters. 425 if spin_index != None: 426 # Get the simulation value. 427 if sim_index != None: 428 # Get the simulation object. 429 sim_obj = getattr(spins[spin_index], param_name+'_sim') 430 431 # R20 parameter. 432 if r20_key != None: 433 if r20_key in sim_obj[sim_index].keys(): 434 value = sim_obj[sim_index][r20_key] 435 436 # All other parameters. 437 else: 438 value = sim_obj[sim_index] 439 440 # Get the normal value. 441 else: 442 # Get the object. 443 obj = getattr(spins[spin_index], param_name) 444 445 # R20 parameter. 446 if r20_key != None: 447 if r20_key in obj.keys(): 448 value = obj[r20_key] 449 450 # All other parameters. 451 else: 452 value = getattr(spins[spin_index], param_name) 453 454 # Cluster specific parameters - use the parameter value from the first spin. 455 else: 456 # Set the simulation value. 457 if sim_index != None: 458 value = getattr(spins[0], param_name+'_sim')[sim_index] 459 460 # Set the normal value. 461 else: 462 value = getattr(spins[0], param_name) 463 464 # The R2eff model parameters. 465 if key != None: 466 if not key in value.keys(): 467 value = 0.0 468 else: 469 value = value[key] 470 471 # Return the value. 472 return value

473 474

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

476 """Set up the relaxation dispersion curve fitting linear constraint matrices A and b. 477 478 Standard notation 479 ================= 480 481 The different constraints used within different models are:: 482 483 0 <= R2 <= 200 484 0 <= R2A <= 200 485 0 <= R2B <= 200 486 pC <= pB <= pA <= 1 487 pA >= 0.85 (the skewed condition, pA >> pB) 488 pB >= 0 489 pC >= 0 490 phi_ex >= 0 491 phi_ex_B >= 0 492 phi_ex_C >= 0 493 padw2 >= 0 494 dw >= 0 495 0 <= kex <= 2e6 496 0 <= k_AB <= 2e6 497 0 <= kB <= 2e6 498 0 <= kC <= 2e6 499 tex >= 0 500 k_AB >= 0 501 502 503 Matrix notation 504 =============== 505 506 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:: 507 508 | 1 0 0 | | R2 | | 0 | 509 | | | | | | 510 |-1 0 0 | | R2 | | -200 | 511 | | | | | | 512 | 1 0 0 | | R2A | | 0 | 513 | | | | | | 514 |-1 0 0 | | R2A | | -200 | 515 | | | | | | 516 | 1 0 0 | | R2B | | 0 | 517 | | | | | | 518 |-1 0 0 | | R2B | | -200 | 519 | | | | | | 520 | 1 0 0 | | pA | | 0.5 | 521 | | | | | | 522 |-1 0 0 | | pA | | -1 | 523 | | | | | | 524 | 1 0 0 | | pA | | 0.85 | 525 | | | | | | 526 | 1 0 0 | | pB | | 0 | 527 | | | | | | 528 | 1 0 0 | | phi_ex | | 0 | 529 | | | | | | 530 | 1 0 0 | . | phi_ex_B | >= | 0 | 531 | | | | | | 532 | 1 0 0 | | phi_ex_C | | 0 | 533 | | | | | | 534 | 1 0 0 | | padw2 | | 0 | 535 | | | | | | 536 | 1 0 0 | | dw | | 0 | 537 | | | | | | 538 | 1 0 0 | | kex | | 0 | 539 | | | | | | 540 |-1 0 0 | | kex | | -2e6 | 541 | | | | | | 542 | 1 0 0 | | kB | | 0 | 543 | | | | | | 544 |-1 0 0 | | kB | | -2e6 | 545 | | | | | | 546 | 1 0 0 | | kC | | 0 | 547 | | | | | | 548 |-1 0 0 | | kC | | -2e6 | 549 | | | | | | 550 | 1 0 0 | | tex | | 0 | 551 | | | | | | 552 | 1 0 0 | | k_AB | | 0 | 553 554 555 @keyword spins: The list of spin data containers for the block. 556 @type spins: list of SpinContainer instances 557 @keyword scaling_matrix: The diagonal, square scaling matrix. 558 @type scaling_matrix: numpy diagonal matrix 559 @return: The matrices A and b. 560 @rtype: numpy rank-2 NxM array, numpy rank-1 N array 561 """ 562 563 # Initialisation (0..j..m). 564 A = [] 565 b = [] 566 n = param_num(spins=spins) 567 zero_array = zeros(n, float64) 568 j = 0 569 570 # Loop over the parameters of the cluster. 571 for param_name, param_index, spin_index, r20_key in loop_parameters(spins=spins): 572 # Effective transversal relaxation rate. 573 if param_name == 'r2eff': 574 A.append(zero_array * 0.0) 575 A.append(zero_array * 0.0) 576 A[j][param_index] = 1.0 577 A[j+1][param_index] = -1.0 578 b.append(0.0) 579 b.append(-200.0 / scaling_matrix[param_index, param_index]) 580 j += 2 581 582 # Initial intensity. 583 elif param_name == 'i0': 584 A.append(zero_array * 0.0) 585 A[j][param_index] = 1.0 586 b.append(0.0) 587 j += 1 588 589 # The transversal relaxation rates (0 <= r2 <= 200). 590 elif param_name in ['r2', 'r2a', 'r2b']: 591 A.append(zero_array * 0.0) 592 A.append(zero_array * 0.0) 593 A[j][param_index] = 1.0 594 A[j+1][param_index] = -1.0 595 b.append(0.0) 596 b.append(-200.0 / scaling_matrix[param_index, param_index]) 597 j += 2 598 599 # The pA.pB.dw**2, phi_ex_B, phi_ex_C and pA.dw**2 parameters (phi_ex* >= 0 and padw2 >= 0). 600 elif param_name in ['phi_ex', 'phi_ex_B', 'phi_ex_C', 'padw2']: 601 A.append(zero_array * 0.0) 602 A[j][param_index] = 1.0 603 b.append(0.0) 604 j += 1 605 606 # Chemical exchange difference (dw >= 0). 607 elif param_name in ['dw', 'dw_AB', 'dw_AC', 'dw_BC', 'dwH', 'dwH_AB', 'dwH_AC', 'dwH_BC']: 608 if not spins[0].model in MODEL_LIST_MMQ + [MODEL_NS_R1RHO_3SITE, MODEL_NS_R1RHO_3SITE_LINEAR]: 609 A.append(zero_array * 0.0) 610 A[j][param_index] = 1.0 611 b.append(0.0) 612 j += 1 613 614 # The population of state A and B. 615 elif param_name == 'pA': 616 # First the pA <= 1 constraint (which also rearranged is pB >= 0). 617 A.append(zero_array * 0.0) 618 A[j][param_index] = -1.0 619 b.append(-1.0 / scaling_matrix[param_index, param_index]) 620 j += 1 621 622 # The skewed condition (pA >> pB). 623 if spins[0].model == MODEL_M61B: 624 A.append(zero_array * 0.0) 625 A[j][param_index] = 1.0 626 b.append(0.85 / scaling_matrix[param_index, param_index]) 627 j += 1 628 629 # Otherwise use the pA >= pB constraint. 630 else: 631 A.append(zero_array * 0.0) 632 A[j][param_index] = 1.0 633 b.append(0.5 / scaling_matrix[param_index, param_index]) 634 j += 1 635 636 # The population of state C (the pB parameter is only present when pC exists). 637 elif param_name == 'pB': 638 # Find pA. 639 for param_name2, param_index2, spin_index2, r20_key2 in loop_parameters(spins=spins): 640 if param_name2 == 'pA': 641 # First pC >= 0 (which rearranged is -pA - pB >= -1). 642 A.append(zero_array * 0.0) 643 A[j][param_index2] = -1.0 644 A[j][param_index] = -1.0 645 b.append(-1.0 / scaling_matrix[param_index, param_index]) 646 j += 1 647 648 # Then the pA >= pC constraint. 649 A.append(zero_array * 0.0) 650 A[j][param_index2] = 1.0 651 A[j][param_index] = -1.0 652 b.append(0.0) 653 j += 1 654 break 655 656 # Exchange rates and times (0 <= k <= 2e6). 657 elif param_name in ['kex', 'kex_AB', 'kex_AC', 'kex_BC', 'k_AB', 'kB', 'kC']: 658 A.append(zero_array * 0.0) 659 A.append(zero_array * 0.0) 660 A[j][param_index] = 1.0 661 A[j+1][param_index] = -1.0 662 b.append(0.0) 663 b.append(-2e6 / scaling_matrix[param_index, param_index]) 664 j += 2 665 666 # Exchange times (tex >= 0). 667 elif param_name in ['tex']: 668 A.append(zero_array * 0.0) 669 A[j][param_index] = 1.0 670 b.append(0.0) 671 j += 1 672 673 # Convert to numpy data structures. 674 A = array(A, float64) 675 b = array(b, float64) 676 677 # Return the matrices. 678 return A, b

679 680

681 -def loop_parameters(spins=None):

682 """Generator function for looping of the parameters of the cluster. 683 684 @keyword spins: The list of spin data containers for the block. 685 @type spins: list of SpinContainer instances 686 @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). 687 @rtype: str, int, int, str 688 """ 689 690 # The parameter index. 691 param_index = -1 692 693 # The R2eff model. 694 if cdp.model_type == 'R2eff': 695 # Loop over the spins. 696 for spin_index in range(len(spins)): 697 # Skip deselected spins. 698 if not spins[spin_index].select: 699 continue 700 701 # Yield the two parameters. 702 params = ['r2eff', 'i0'] 703 for i in range(2): 704 # First increment the indices. 705 param_index += 1 706 707 # Yield the data. 708 yield params[i], param_index, spin_index, None 709 710 # All other models. 711 else: 712 # First the R2 parameters (one per spin per field strength). 713 for spin_index in range(len(spins)): 714 # Skip deselected spins. 715 if not spins[spin_index].select: 716 continue 717 718 # The R2 parameter. 719 if 'r2' in spins[0].params: 720 for exp_type, frq in loop_exp_frq(): 721 param_index += 1 722 yield 'r2', param_index, spin_index, generate_r20_key(exp_type=exp_type, frq=frq) 723 724 # The R2A parameter. 725 if 'r2a' in spins[0].params: 726 for exp_type, frq in loop_exp_frq(): 727 param_index += 1 728 yield 'r2a', param_index, spin_index, generate_r20_key(exp_type=exp_type, frq=frq) 729 730 # The R2B parameter. 731 if 'r2b' in spins[0].params: 732 for exp_type, frq in loop_exp_frq(): 733 param_index += 1 734 yield 'r2b', param_index, spin_index, generate_r20_key(exp_type=exp_type, frq=frq) 735 736 # 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). 737 for spin_index in range(len(spins)): 738 # Skip deselected spins. 739 if not spins[spin_index].select: 740 continue 741 742 # Yield the data. 743 if 'phi_ex' in spins[spin_index].params: 744 param_index += 1 745 yield 'phi_ex', param_index, spin_index, None 746 if 'phi_ex_B' in spins[spin_index].params: 747 param_index += 1 748 yield 'phi_ex_B', param_index, spin_index, None 749 if 'phi_ex_C' in spins[spin_index].params: 750 param_index += 1 751 yield 'phi_ex_C', param_index, spin_index, None 752 if 'padw2' in spins[spin_index].params: 753 param_index += 1 754 yield 'padw2', param_index, spin_index, None 755 if 'dw' in spins[spin_index].params: 756 param_index += 1 757 yield 'dw', param_index, spin_index, None 758 if 'dw_AB' in spins[spin_index].params: 759 param_index += 1 760 yield 'dw_AB', param_index, spin_index, None 761 if 'dw_BC' in spins[spin_index].params: 762 param_index += 1 763 yield 'dw_BC', param_index, spin_index, None 764 if 'dw_AC' in spins[spin_index].params: 765 param_index += 1 766 yield 'dw_AC', param_index, spin_index, None 767 768 # Then a separate block for the proton chemical shift difference parameters for the MQ models (one per spin). 769 for spin_index in range(len(spins)): 770 # Skip deselected spins. 771 if not spins[spin_index].select: 772 continue 773 774 if 'dwH' in spins[spin_index].params: 775 param_index += 1 776 yield 'dwH', param_index, spin_index, None 777 if 'dwH_AB' in spins[spin_index].params: 778 param_index += 1 779 yield 'dwH_AB', param_index, spin_index, None 780 if 'dwH_BC' in spins[spin_index].params: 781 param_index += 1 782 yield 'dwH_BC', param_index, spin_index, None 783 if 'dwH_AC' in spins[spin_index].params: 784 param_index += 1 785 yield 'dwH_AC', param_index, spin_index, None 786 787 # All other parameters (one per spin cluster). 788 for param in spins[0].params: 789 if not param in ['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']: 790 param_index += 1 791 yield param, param_index, None, None

792 793

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

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

848 849

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

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

867 868

869 -def param_num(spins=None):

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

938 939

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

941 """Return the value for the given parameter. 942 943 @keyword value: The parameter value to set. 944 @type value: float 945 @keyword key: The key for the R2eff and I0 parameters. 946 @type key: str or None 947 @keyword spins: The list of spin data containers for the block. 948 @type spins: list of SpinContainer instances 949 @keyword sim_index: The optional MC simulation index. 950 @type sim_index: int 951 @keyword param_name: The parameter name. 952 @type param_name: str 953 @keyword spin_index: The spin index (for the cluster). 954 @type spin_index: int 955 @keyword r20_key: The unique R20 parameter key. 956 @type r20_key: str 957 """ 958 959 # Spin specific parameters. 960 if spin_index != None: 961 # Set the simulation value. 962 if sim_index != None: 963 # Get the simulation object. 964 obj = getattr(spins[spin_index], param_name+'_sim') 965 966 # R20 parameter. 967 if r20_key != None: 968 obj[sim_index][r20_key] = value 969 970 # All other parameters. 971 else: 972 if key != None: 973 obj[sim_index][key] = value 974 else: 975 obj[sim_index] = value 976 977 # Set the normal value. 978 else: 979 # Get the object. 980 obj = getattr(spins[spin_index], param_name) 981 982 # R20 parameter. 983 if r20_key != None: 984 obj[r20_key] = value 985 986 # Set the normal value. 987 else: 988 if key != None: 989 obj[key] = value 990 else: 991 setattr(spins[spin_index], param_name, value) 992 993 # Cluster specific parameters. 994 else: 995 # Set the same parameter value for all spins in the cluster. 996 for spin in spins: 997 # Skip deselected spins. 998 if not spin.select: 999 continue 1000 1001 # Set the simulation value. 1002 if sim_index != None: 1003 sim_obj = getattr(spin, param_name+'_sim') 1004 sim_obj[sim_index] = value 1005 1006 # Set the normal value. 1007 else: 1008 setattr(spin, param_name, value)

1009

Source Code for Module specific_analyses.relax_disp.parameters