Package specific_analyses :: Package relax_disp :: Module parameters
[hide private]
[frames] | no frames]

Source Code for Module specific_analyses.relax_disp.parameters

   1  ############################################################################### 
   2  #                                                                             # 
   3  # Copyright (C) 2003-2008,2013-2014,2019 Edward d'Auvergne                    # 
   4  # Copyright (C) 2006 Chris MacRaild                                           # 
   5  # Copyright (C) 2008-2009 Sebastien Morin                                     # 
   6  # Copyright (C) 2013-2014 Troels E. Linnet                                    # 
   7  #                                                                             # 
   8  # This file is part of the program relax (http://www.nmr-relax.com).          # 
   9  #                                                                             # 
  10  # This program is free software: you can redistribute it and/or modify        # 
  11  # it under the terms of the GNU General Public License as published by        # 
  12  # the Free Software Foundation, either version 3 of the License, or           # 
  13  # (at your option) any later version.                                         # 
  14  #                                                                             # 
  15  # This program is distributed in the hope that it will be useful,             # 
  16  # but WITHOUT ANY WARRANTY; without even the implied warranty of              # 
  17  # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the               # 
  18  # GNU General Public License for more details.                                # 
  19  #                                                                             # 
  20  # You should have received a copy of the GNU General Public License           # 
  21  # along with this program.  If not, see <http://www.gnu.org/licenses/>.       # 
  22  #                                                                             # 
  23  ############################################################################### 
  24   
  25  # Module docstring. 
  26  """Functions relating to the parameters of the relaxation dispersion models.""" 
  27   
  28  # Python module imports. 
  29  from copy import deepcopy 
  30  from numpy import array, float64, median, zeros 
  31  import sys 
  32   
  33  # relax module imports. 
  34  from lib.dispersion.variables import MODEL_LIST_MMQ, \ 
  35      MODEL_M61B, \ 
  36      MODEL_NS_MMQ_3SITE, \ 
  37      MODEL_NS_MMQ_3SITE_LINEAR, \ 
  38      MODEL_NS_R1RHO_3SITE, \ 
  39      MODEL_NS_R1RHO_3SITE_LINEAR,  \ 
  40      PARAMS_R1, \ 
  41      PARAMS_R20, \ 
  42      PARAMS_CHEM_SHIFT_DIFF, \ 
  43      PARAMS_CHEM_SHIFT_DIFF_MMQ, \ 
  44      PARAMS_SPIN 
  45  from lib.errors import RelaxError, RelaxNoSequenceError 
  46  from lib.text.sectioning import subsection 
  47  from pipe_control import pipes 
  48  from pipe_control.mol_res_spin import exists_mol_res_spin_data, return_spin, spin_loop 
  49  from pipe_control.pipes import check_pipe 
  50  from specific_analyses.relax_disp.data import count_spins, generate_r20_key, has_exponential_exp_type, is_r1_optimised, loop_cluster, loop_exp_frq 
  51   
  52   
53 -def assemble_param_vector(spins=None, key=None, sim_index=None):
54 """Assemble the dispersion relaxation dispersion curve fitting parameter vector. 55 56 @keyword spins: The list of spin data containers for the block. 57 @type spins: list of SpinContainer instances 58 @keyword key: The key for the R2eff and I0 parameters. 59 @type key: str or None 60 @keyword sim_index: The optional MC simulation index. 61 @type sim_index: int 62 @return: An array of the parameter values of the dispersion relaxation model. 63 @rtype: numpy float array 64 """ 65 66 # Initialise. 67 param_vector = [] 68 69 # Loop over the parameters of the cluster. 70 for param_name, param_index, spin_index, r20_key in loop_parameters(spins=spins): 71 # Get the value. 72 value = get_value(key=key, spins=spins, sim_index=sim_index, param_name=param_name, spin_index=spin_index, r20_key=r20_key) 73 74 # Add to the vector. 75 param_vector.append(value) 76 77 # Convert all None values to 0.0. 78 for i in range(len(param_vector)): 79 if param_vector[i] == None: 80 param_vector[i] = 0.0 81 82 # Return a numpy array. 83 return array(param_vector, float64)
84 85
86 -def copy(pipe_from=None, pipe_to=None):
87 """Copy dispersion parameters from one data pipe to another. 88 89 For spin clusters, the median of previous values is taken to avoid extreme outliers. 90 91 92 @param pipe_from: The data pipe to copy the value from. This defaults to the current data pipe. 93 @type pipe_from: str 94 @param pipe_to: The data pipe to copy the value to. This defaults to the current data pipe. 95 @type pipe_to: str 96 """ 97 98 # The current data pipe. 99 pipe_orig = pipes.cdp_name() 100 if pipe_from == None: 101 pipe_from = pipe_orig 102 if pipe_to == None: 103 pipe_to = pipe_orig 104 105 # Test that the pipes exist. 106 check_pipe(pipe_from) 107 check_pipe(pipe_to) 108 109 # Test that the pipes are not the same. 110 if pipe_from == pipe_to: 111 raise RelaxError("The source and destination pipes cannot be the same.") 112 113 # Test if the sequence data for pipe_from is loaded. 114 if not exists_mol_res_spin_data(pipe_from): 115 raise RelaxNoSequenceError(pipe_from) 116 117 # Test if the sequence data for pipe_to is loaded. 118 if not exists_mol_res_spin_data(pipe_to): 119 raise RelaxNoSequenceError(pipe_to) 120 121 # Switch to the destination data pipe. 122 pipes.switch(pipe_to) 123 124 # Loop over the clusters. 125 for spin_ids in loop_cluster(): 126 # Initialise some variables. 127 model = None 128 pA = [] 129 pB = [] 130 pC = [] 131 kex = [] 132 kex_AB = [] 133 kex_AC = [] 134 kex_BC = [] 135 k_AB = [] 136 kB = [] 137 kC = [] 138 tex = [] 139 count = 0 140 spins_from = [] 141 spins_to = [] 142 selected_cluster = False 143 144 # Loop over the spins, adding parameters to a list, which in the end will be used to find the median. 145 for id in spin_ids: 146 # Get the spins, then store them. 147 spin_from = return_spin(spin_id=id, pipe=pipe_from) 148 spin_to = return_spin(spin_id=id, pipe=pipe_to) 149 spins_from.append(spin_from) 150 spins_to.append(spin_to) 151 152 # Skip deselected spins. 153 if not spin_from.select or not spin_to.select: 154 continue 155 156 # The first printout. 157 if not selected_cluster: 158 subsection(file=sys.stdout, text="Copying parameters for the spin block %s"%spin_ids, prespace=2) 159 160 # Change the cluster selection flag. 161 selected_cluster = True 162 163 # The model. 164 if not model: 165 model = spin_from.model 166 167 # Check that the models match for all spins of the cluster. 168 if spin_from.model != model: 169 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)) 170 if spin_to.model != model: 171 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)) 172 173 # Sum the source parameters. 174 if 'pA' in spin_from.params: 175 pA.append(spin_from.pA) 176 if 'pB' in spin_from.params: 177 pB.append(spin_from.pB) 178 if 'pC' in spin_from.params: 179 pC.append(spin_from.pC) 180 if 'kex' in spin_from.params: 181 kex.append(spin_from.kex) 182 if 'kex_AB' in spin_from.params: 183 kex_AB.append(spin_from.kex_AB) 184 if 'kex_AC' in spin_from.params: 185 kex_AC.append(spin_from.kex_AC) 186 if 'kex_BC' in spin_from.params: 187 kex_BC.append(spin_from.kex_BC) 188 if 'k_AB' in spin_from.params: 189 k_AB.append(spin_from.k_AB) 190 if 'kB' in spin_from.params: 191 kB.append(spin_from.kB) 192 if 'kC' in spin_from.params: 193 kC.append(spin_from.kC) 194 if 'tex' in spin_from.params: 195 tex.append(spin_from.tex) 196 197 # Increment the spin count. 198 count += 1 199 200 # The cluster is not selected, so move to the next. 201 if not selected_cluster: 202 continue 203 204 # Take median of parameters. 205 if len(pA) > 0: 206 pA = median(pA) 207 print("Median pA value: %.15f" % pA) 208 if len(pB) > 0: 209 pB = median(pB) 210 print("Median pB value: %.15f" % pB) 211 if len(pC) > 0: 212 pC = median(pC) 213 print("Median pC value: %.15f" % pC) 214 if len(kex) > 0: 215 kex = median(kex) 216 print("Median kex value: %.15f" % kex) 217 if len(kex_AB) > 0: 218 kex_AB = median(kex_AB) 219 print("Median k_AB value: %.15f" % kex_AB) 220 if len(kex_AC) > 0: 221 kex_AC = median(kex_AC) 222 print("Median k_AC value: %.15f" % kex_AC) 223 if len(kex_BC) > 0: 224 kex_BC = median(kex_BC) 225 print("Median k_BC value: %.15f" % kex_BC) 226 if len(k_AB) > 0: 227 k_AB = median(k_AB) 228 print("Median k_AB value: %.15f" % k_AB) 229 if len(kB) > 0: 230 kB = median(kB) 231 print("Median kB value: %.15f" % kB) 232 if len(kC) > 0: 233 kC = median(kC) 234 print("Median kC value: %.15f" % kC) 235 if len(tex) > 0: 236 tex = median(tex) 237 print("Median tex value: %.15f" % tex) 238 239 # Loop over the spins, this time copying the parameters. 240 for i in range(len(spin_ids)): 241 # Alias the containers. 242 spin_from = spins_from[i] 243 spin_to = spins_to[i] 244 245 # Skip deselected spins. 246 if not spin_from.select or not spin_to.select: 247 continue 248 249 # The R20 parameters. 250 if 'r2' in spin_from.params: 251 spin_to.r2 = deepcopy(spin_from.r2) 252 if 'r2a' in spin_from.params: 253 spin_to.r2a = deepcopy(spin_from.r2a) 254 if 'r2b' in spin_from.params: 255 spin_to.r2b = deepcopy(spin_from.r2b) 256 257 # The median parameters. 258 if 'pB' in spin_from.params and 'pC' not in spin_from.params: 259 spin_to.pA = pA 260 spin_to.pB = pB 261 spin_to.pC = 1.0 - pA - pB 262 elif 'pA' in spin_from.params: 263 spin_to.pA = pA 264 spin_to.pB = 1.0 - pA 265 if 'kex' in spin_from.params: 266 spin_to.kex = kex 267 if 'kex_AB' in spin_from.params: 268 spin_to.kex_AB = kex_AB 269 if 'kex_AC' in spin_from.params: 270 spin_to.kex_AC = kex_AC 271 if 'kex_BC' in spin_from.params: 272 spin_to.kex_BC = kex_BC 273 if 'k_AB' in spin_from.params: 274 spin_to.k_AB = k_AB 275 if 'kB' in spin_from.params: 276 spin_to.kB = kB 277 if 'kC' in spin_from.params: 278 spin_to.kC = kC 279 if 'tex' in spin_from.params: 280 spin_to.tex = tex 281 282 # All other spin specific parameters. 283 for param in spin_from.params: 284 if param in ['r2', 'r2a', 'r2b', 'pA', 'pB', 'pC', 'kex', 'kex_AB', 'kex_AC', 'kex_BC', 'k_AB', 'kB', 'kC', 'tex']: 285 continue 286 287 # Copy the value. 288 setattr(spin_to, param, deepcopy(getattr(spin_from, param))) 289 290 # Switch back to the original data pipe. 291 pipes.switch(pipe_orig)
292 293
294 -def disassemble_param_vector(param_vector=None, key=None, spins=None, sim_index=None):
295 """Disassemble the relaxation dispersion parameter vector. 296 297 @keyword param_vector: The parameter vector. 298 @type param_vector: numpy array 299 @keyword key: The key for the R2eff and I0 parameters. 300 @type key: str or None 301 @keyword spins: The list of spin data containers for the block. 302 @type spins: list of SpinContainer instances 303 @keyword sim_index: The optional MC simulation index. 304 @type sim_index: int 305 """ 306 307 # Initialise parameters if needed. 308 for spin in spins: 309 # Skip deselected spins. 310 if not spin.select: 311 continue 312 313 # The R2 parameter. 314 if 'r2' in spin.params: 315 if sim_index != None: 316 spin.r2_sim[sim_index] = {} 317 else: 318 spin.r2 = {} 319 320 # The R2A parameter. 321 if 'r2a' in spin.params: 322 if sim_index != None: 323 spin.r2a_sim[sim_index] = {} 324 else: 325 spin.r2a = {} 326 327 # The R2B parameter. 328 if 'r2b' in spin.params: 329 if sim_index != None: 330 spin.r2b_sim[sim_index] = {} 331 else: 332 spin.r2b = {} 333 334 # Loop over the parameters of the cluster, setting the values. 335 for param_name, param_index, spin_index, r20_key in loop_parameters(spins=spins): 336 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)
337 338
339 -def get_param_names(spins=None, full=False):
340 """Generate a list of dispersion parameter names for the given spins. 341 342 @keyword spins: The list of spin data containers for the block. 343 @type spins: list of SpinContainer instances 344 @keyword full: A flag which if True will add the R2 key to the parameter names. 345 @type full: bool 346 """ 347 348 # Initialise the structure. 349 names = [] 350 351 # Loop over the parameters. 352 for param_name, param_index, spin_index, r20_key in loop_parameters(spins=spins): 353 # Set the initial text. 354 param_text = param_name 355 356 # The parameters with additional details. 357 if full and param_name in PARAMS_R20: 358 param_text += " (%s)" % r20_key 359 360 # Append the text. 361 names.append(param_text) 362 363 # Return the structure. 364 return names
365 366
367 -def get_value(key=None, spins=None, sim_index=None, param_name=None, spin_index=None, r20_key=None):
368 """Return the value for the given parameter. 369 370 @keyword key: The key for the R2eff and I0 parameters. 371 @type key: str or None 372 @keyword spins: The list of spin data containers for the block. 373 @type spins: list of SpinContainer instances 374 @keyword sim_index: The optional MC simulation index. 375 @type sim_index: int 376 @keyword param_name: The parameter name. 377 @type param_name: str 378 @keyword spin_index: The spin index (for the cluster). 379 @type spin_index: int 380 @keyword r20_key: The unique R20 parameter key. 381 @type r20_key: str 382 @return: The parameter value. 383 @rtype: float 384 """ 385 386 # Default value of None. 387 value = None 388 389 # Spin specific parameters. 390 if spin_index != None: 391 # Get the simulation value. 392 if sim_index != None: 393 # Get the simulation object. 394 sim_obj = getattr(spins[spin_index], param_name+'_sim') 395 396 # R20 parameter. 397 if r20_key != None: 398 if r20_key in sim_obj[sim_index]: 399 value = sim_obj[sim_index][r20_key] 400 401 # All other parameters. 402 else: 403 value = sim_obj[sim_index] 404 405 # Get the normal value. 406 else: 407 # Get the object. 408 obj = getattr(spins[spin_index], param_name) 409 410 # R20 parameter. 411 if r20_key != None: 412 if r20_key in obj: 413 value = obj[r20_key] 414 415 # All other parameters. 416 else: 417 value = getattr(spins[spin_index], param_name) 418 419 # Cluster specific parameters - use the parameter value from the first spin. 420 else: 421 # Set the simulation value. 422 if sim_index != None: 423 value = getattr(spins[0], param_name+'_sim')[sim_index] 424 425 # Set the normal value. 426 else: 427 value = getattr(spins[0], param_name) 428 429 # The R2eff model parameters. 430 if key != None: 431 if not key in value: 432 value = 0.0 433 else: 434 value = value[key] 435 436 # Return the value. 437 return value
438 439
440 -def linear_constraints(spins=None, scaling_matrix=None):
441 """Set up the relaxation dispersion curve fitting linear constraint matrices A and b. 442 443 Standard notation 444 ================= 445 446 The different constraints used within different models are:: 447 448 0 <= R1 <= 200 449 0 <= R2 <= 200 450 0 <= R2A <= 200 451 0 <= R2B <= 200 452 pC <= pB <= pA <= 1 453 pA >= 0.85 (the skewed condition, pA >> pB) 454 pB >= 0 455 pC >= 0 456 phi_ex >= 0 457 phi_ex_B >= 0 458 phi_ex_C >= 0 459 padw2 >= 0 460 dw >= 0 461 0 <= kex <= 2e6 462 0 <= kB <= 2e6 463 0 <= kC <= 2e6 464 0 <= k_AB <= 100 465 tex >= 0 466 467 468 Matrix notation 469 =============== 470 471 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:: 472 473 | 1 0 0 | | R1 | | 0 | 474 | | | | | | 475 |-1 0 0 | | R1 | | -200 | 476 | | | | | | 477 | 1 0 0 | | R2 | | 0 | 478 | | | | | | 479 |-1 0 0 | | R2 | | -200 | 480 | | | | | | 481 | 1 0 0 | | R2A | | 0 | 482 | | | | | | 483 |-1 0 0 | | R2A | | -200 | 484 | | | | | | 485 | 1 0 0 | | R2B | | 0 | 486 | | | | | | 487 |-1 0 0 | | R2B | | -200 | 488 | | | | | | 489 | 1 0 0 | | pA | | 0.5 | 490 | | | | | | 491 |-1 0 0 | | pA | | -1 | 492 | | | | | | 493 | 1 0 0 | | pA | | 0.85 | 494 | | | | | | 495 | 1 0 0 | | pB | | 0 | 496 | | | | | | 497 | 1 0 0 | | phi_ex | | 0 | 498 | | | | | | 499 | 1 0 0 | . | phi_ex_B | >= | 0 | 500 | | | | | | 501 | 1 0 0 | | phi_ex_C | | 0 | 502 | | | | | | 503 | 1 0 0 | | padw2 | | 0 | 504 | | | | | | 505 | 1 0 0 | | dw | | 0 | 506 | | | | | | 507 | 1 0 0 | | kex | | 0 | 508 | | | | | | 509 |-1 0 0 | | kex | | -2e6 | 510 | | | | | | 511 | 1 0 0 | | kB | | 0 | 512 | | | | | | 513 |-1 0 0 | | kB | | -2e6 | 514 | | | | | | 515 | 1 0 0 | | kC | | 0 | 516 | | | | | | 517 |-1 0 0 | | kC | | -2e6 | 518 | | | | | | 519 | 1 0 0 | | k_AB | | 0 | 520 | | | | | | 521 |-1 0 0 | | k_AB | | 100 | 522 | | | | | | 523 | 1 0 0 | | tex | | 0 | 524 525 526 @keyword spins: The list of spin data containers for the block. 527 @type spins: list of SpinContainer instances 528 @keyword scaling_matrix: The diagonal, square scaling matrix. 529 @type scaling_matrix: numpy diagonal matrix 530 @return: The matrices A and b. 531 @rtype: numpy rank-2 NxM array, numpy rank-1 N array 532 """ 533 534 # Initialisation (0..j..m). 535 A = [] 536 b = [] 537 n = param_num(spins=spins) 538 zero_array = zeros(n, float64) 539 j = 0 540 541 # Loop over the parameters of the cluster. 542 for param_name, param_index, spin_index, r20_key in loop_parameters(spins=spins): 543 # Effective transversal relaxation rate. 544 if param_name == 'r2eff': 545 A.append(zero_array * 0.0) 546 A.append(zero_array * 0.0) 547 A[j][param_index] = 1.0 548 A[j+1][param_index] = -1.0 549 b.append(0.0) 550 b.append(-200.0 / scaling_matrix[param_index, param_index]) 551 j += 2 552 553 # Initial intensity. 554 elif param_name == 'i0': 555 A.append(zero_array * 0.0) 556 A[j][param_index] = 1.0 557 b.append(0.0) 558 j += 1 559 560 # The longitudinal relaxation rates (0 <= r1 <= 200). 561 elif param_name in ['r1']: 562 A.append(zero_array * 0.0) 563 A.append(zero_array * 0.0) 564 A[j][param_index] = 1.0 565 A[j+1][param_index] = -1.0 566 b.append(0.0) 567 b.append(-200.0 / scaling_matrix[param_index, param_index]) 568 j += 2 569 570 # The transversal relaxation rates (0 <= r2 <= 200). 571 elif param_name in PARAMS_R20: 572 A.append(zero_array * 0.0) 573 A.append(zero_array * 0.0) 574 A[j][param_index] = 1.0 575 A[j+1][param_index] = -1.0 576 b.append(0.0) 577 b.append(-200.0 / scaling_matrix[param_index, param_index]) 578 j += 2 579 580 # The pA.pB.dw**2, phi_ex_B, phi_ex_C and pA.dw**2 parameters (phi_ex* >= 0 and padw2 >= 0). 581 elif param_name in ['phi_ex', 'phi_ex_B', 'phi_ex_C', 'padw2']: 582 A.append(zero_array * 0.0) 583 A[j][param_index] = 1.0 584 b.append(0.0) 585 j += 1 586 587 # Chemical exchange difference (dw >= 0). 588 elif param_name in ['dw', 'dw_AB', 'dw_AC', 'dw_BC', 'dwH', 'dwH_AB', 'dwH_AC', 'dwH_BC']: 589 if not spins[0].model in MODEL_LIST_MMQ + [MODEL_NS_R1RHO_3SITE, MODEL_NS_R1RHO_3SITE_LINEAR]: 590 A.append(zero_array * 0.0) 591 A[j][param_index] = 1.0 592 b.append(0.0) 593 j += 1 594 595 # The population of state A and B. 596 elif param_name == 'pA': 597 # First the pA <= 1 constraint (which also rearranged is pB >= 0). 598 A.append(zero_array * 0.0) 599 A[j][param_index] = -1.0 600 b.append(-1.0 / scaling_matrix[param_index, param_index]) 601 j += 1 602 603 # The skewed condition (pA >> pB). 604 if spins[0].model == MODEL_M61B: 605 A.append(zero_array * 0.0) 606 A[j][param_index] = 1.0 607 b.append(0.85 / scaling_matrix[param_index, param_index]) 608 j += 1 609 610 # Otherwise use the pA >= pB constraint. 611 else: 612 A.append(zero_array * 0.0) 613 A[j][param_index] = 1.0 614 b.append(0.5 / scaling_matrix[param_index, param_index]) 615 j += 1 616 617 # The population of state C (the pB parameter is only present when pC exists). 618 elif param_name == 'pB': 619 # Find pA. 620 for param_name2, param_index2, spin_index2, r20_key2 in loop_parameters(spins=spins): 621 if param_name2 == 'pA': 622 # First pC >= 0 (which rearranged is -pA - pB >= -1). 623 A.append(zero_array * 0.0) 624 A[j][param_index2] = -1.0 625 A[j][param_index] = -1.0 626 b.append(-1.0 / scaling_matrix[param_index, param_index]) 627 j += 1 628 629 # Then the pA >= pB constraint. 630 A.append(zero_array * 0.0) 631 A[j][param_index2] = 1.0 632 A[j][param_index] = -1.0 633 b.append(0.0) 634 j += 1 635 break 636 637 # Exchange rates and times (0 <= k <= 2e6). 638 elif param_name in ['kex', 'kex_AB', 'kex_AC', 'kex_BC', 'kB', 'kC']: 639 A.append(zero_array * 0.0) 640 A.append(zero_array * 0.0) 641 A[j][param_index] = 1.0 642 A[j+1][param_index] = -1.0 643 b.append(0.0) 644 b.append(-2e6 / scaling_matrix[param_index, param_index]) 645 j += 2 646 647 # Exchange rates and times (0 <= k_AB <= 100). 648 elif param_name in ['k_AB']: 649 A.append(zero_array * 0.0) 650 A.append(zero_array * 0.0) 651 A[j][param_index] = 1.0 652 A[j+1][param_index] = -1.0 653 b.append(0.0) 654 b.append(-100. / scaling_matrix[param_index, param_index]) 655 j += 2 656 657 # Exchange times (tex >= 0). 658 elif param_name in ['tex']: 659 A.append(zero_array * 0.0) 660 A[j][param_index] = 1.0 661 b.append(0.0) 662 j += 1 663 664 # Convert to numpy data structures. 665 A = array(A, float64) 666 b = array(b, float64) 667 668 # Return the matrices. 669 return A, b
670 671
672 -def loop_parameters(spins=None):
673 """Generator function for looping of the model parameters of the cluster. 674 675 @keyword spins: The list of spin data containers for the block. 676 @type spins: list of SpinContainer instances 677 @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). 678 @rtype: str, int, int, str 679 """ 680 681 # Make sure that the R1 parameter is correctly set up. 682 r1_setup() 683 684 # The parameter index. 685 param_index = -1 686 687 # The R2eff model. 688 if cdp.model_type == 'R2eff': 689 # Loop over the spins. 690 for spin_index in range(len(spins)): 691 # Skip deselected spins. 692 if not spins[spin_index].select: 693 continue 694 695 # One or two parameters. 696 params = ['r2eff'] 697 if has_exponential_exp_type(): 698 params = ['r2eff', 'i0'] 699 700 # Yield the parameters. 701 for param in params: 702 # First increment the indices. 703 param_index += 1 704 705 # Yield the data. 706 yield param, param_index, spin_index, None 707 708 # All other models. 709 else: 710 # First the R1 fit parameter (one per spin per field strength). 711 for spin_index in range(len(spins)): 712 # Skip deselected spins. 713 if not spins[spin_index].select: 714 continue 715 716 # The parameters. 717 for param in PARAMS_R1: 718 if param in spins[spin_index].params: 719 for exp_type, frq in loop_exp_frq(): 720 param_index += 1 721 yield param, param_index, spin_index, generate_r20_key(exp_type=exp_type, frq=frq) 722 723 # Then the R2 parameters (one per spin per field strength). 724 for spin_index in range(len(spins)): 725 # Skip deselected spins. 726 if not spins[spin_index].select: 727 continue 728 729 # The parameters. 730 for param in PARAMS_R20: 731 if param in spins[spin_index].params: 732 for exp_type, frq in loop_exp_frq(): 733 param_index += 1 734 yield param, param_index, spin_index, generate_r20_key(exp_type=exp_type, frq=frq) 735 736 # Then the chemical shift difference parameters (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 for param in PARAMS_CHEM_SHIFT_DIFF: 744 if param in spins[spin_index].params: 745 param_index += 1 746 yield param, param_index, spin_index, None 747 748 # Then a separate block for the proton chemical shift difference parameters for the MQ models (one per spin). 749 for spin_index in range(len(spins)): 750 # Skip deselected spins. 751 if not spins[spin_index].select: 752 continue 753 754 # Yield the data. 755 for param in PARAMS_CHEM_SHIFT_DIFF_MMQ: 756 if param in spins[spin_index].params: 757 param_index += 1 758 yield param, param_index, spin_index, None 759 760 # All other parameters (one per spin cluster). 761 for spin_index in range(len(spins)): 762 # Skip deselected spins. 763 if not spins[spin_index].select: 764 continue 765 766 # The parameters. 767 for param in spins[0].params: 768 if not param in PARAMS_SPIN: 769 param_index += 1 770 yield param, param_index, None, None 771 772 # No more spins. 773 break
774 775
776 -def param_conversion(key=None, spins=None, sim_index=None):
777 """Calculate the non-model parameter values from the current parameter values. 778 779 @keyword key: The key for the R2eff and I0 parameters. 780 @type key: str or None 781 @keyword spins: The list of spin data containers for the block. 782 @type spins: list of SpinContainer instances 783 @keyword sim_index: The optional MC simulation index. 784 @type sim_index: int 785 """ 786 787 # Loop over the parameters of the cluster. 788 for param_name, param_index, spin_index, r20_key in loop_parameters(spins=spins): 789 # Get the value. 790 value = get_value(key=key, spins=spins, sim_index=sim_index, param_name=param_name, spin_index=spin_index, r20_key=r20_key) 791 792 # The pA to pB to pC conversion. 793 if param_name == 'pA': 794 # 3-site exchange with pA and pB optimised. 795 if spins[0].model in [MODEL_NS_MMQ_3SITE, MODEL_NS_MMQ_3SITE_LINEAR, MODEL_NS_R1RHO_3SITE, MODEL_NS_R1RHO_3SITE_LINEAR]: 796 # Get the pB value. 797 pB = get_value(key=key, spins=spins, sim_index=sim_index, param_name='pB', spin_index=spin_index, r20_key=r20_key) 798 799 # Set the pC value. 800 pC = 1.0 - value - pB 801 set_value(value=pC, key=key, spins=spins, sim_index=sim_index, param_name='pC', spin_index=spin_index) 802 803 # 2-site exchange. 804 else: 805 pB = 1.0 - value 806 set_value(value=pB, key=key, spins=spins, sim_index=sim_index, param_name='pB', spin_index=spin_index) 807 808 # The kex to tex conversion. 809 if param_name == 'kex': 810 tex = 1.0 / value 811 set_value(value=tex, key=key, spins=spins, sim_index=sim_index, param_name='tex', spin_index=spin_index) 812 813 # The kex to k_AB and k_BA conversion. 814 if param_name == 'kex' and 'pA' in spins[0].params: 815 # Get pA value. 816 pA = get_value(key=key, spins=spins, sim_index=sim_index, param_name='pA', spin_index=spin_index) 817 818 # Calculate k_AB value and set it. 819 k_AB = value * (1.0 - pA) 820 set_value(value=k_AB, key=key, spins=spins, sim_index=sim_index, param_name='k_AB', spin_index=spin_index) 821 822 # Calculate k_BA value and set it. 823 k_BA = value * pA 824 set_value(value=k_BA, key=key, spins=spins, sim_index=sim_index, param_name='k_BA', spin_index=spin_index) 825 826 # The tex to kex conversion. 827 if param_name == 'tex': 828 kex = 1.0 / value 829 set_value(value=kex, key=key, spins=spins, sim_index=sim_index, param_name='kex', spin_index=spin_index)
830 831
832 -def param_index_to_param_info(index=None, spins=None):
833 """Convert the given parameter array index to parameter identifying information. 834 835 The parameter index will be converted to the parameter name, the relevant spin index in the cluster, and relevant exponential curve key. 836 837 838 @keyword index: The index of the parameter array. 839 @type index: int 840 @keyword spins: The list of spin data containers for the block. 841 @type spins: list of SpinContainer instances 842 @return: The parameter name, the spin index (for the cluster), and the frequency index (for parameters with different values per spectrometer field strength). 843 @rtype: str, int, int 844 """ 845 846 # Loop over the parameters, yielding when a match is found. 847 for param_name, param_index, spin_index, r20_key in loop_parameters(spins=spins): 848 if param_index == index: 849 return param_name, spin_index, r20_key
850 851
852 -def param_num(spins=None):
853 """Determine the number of parameters in the model. 854 855 @keyword spins: The list of spin data containers for the block. 856 @type spins: list of SpinContainer instances 857 @return: The number of model parameters. 858 @rtype: int 859 """ 860 861 # Initialise the number. 862 num = 0 863 864 # The R2eff model. 865 if cdp.model_type == 'R2eff': 866 # Count the selected spins. 867 spin_num = count_spins(spins) 868 869 # Exponential curves (with clustering). 870 if has_exponential_exp_type(): 871 return 2 * spin_num 872 873 # Fixed time period experiments (with clustering). 874 return 1 * spin_num 875 876 # Check the spin cluster. 877 for spin in spins: 878 # Skip deselected spins. 879 if not spin.select: 880 continue 881 882 if len(spin.params) != len(spins[0].params): 883 raise RelaxError("The number of parameters for each spin in the cluster are not the same.") 884 885 # Count the number of R10 parameters. 886 for spin in spins: 887 # Skip deselected spins. 888 if not spin.select: 889 continue 890 891 for i in range(len(spin.params)): 892 if spin.params[i] in PARAMS_R1: 893 for exp_type, frq in loop_exp_frq(): 894 num += 1 895 896 # Count the number of R20 parameters. 897 for spin in spins: 898 # Skip deselected spins. 899 if not spin.select: 900 continue 901 902 for i in range(len(spin.params)): 903 if spin.params[i] in PARAMS_R20: 904 for exp_type, frq in loop_exp_frq(): 905 num += 1 906 907 # Count the number of spin specific parameters for all spins. 908 alread_counted = PARAMS_R1 + PARAMS_R20 909 for spin in spins: 910 # Skip deselected spins. 911 if not spin.select: 912 continue 913 914 for param in spin.params: 915 if param in PARAMS_SPIN and not param in alread_counted: 916 num += 1 917 918 # Count all other parameters, but only for a single spin. 919 alread_counted += PARAMS_SPIN 920 for spin in spins: 921 # Skip deselected spins. 922 if not spin.select: 923 continue 924 925 for param in spin.params: 926 if not param in alread_counted: 927 num += 1 928 break 929 930 # Return the number. 931 return num
932 933
934 -def r1_setup():
935 """Modify the current spin parameter list to either include or exclude the R1 parameter.""" 936 937 # Loop over the spins. 938 for spin, spin_id in spin_loop(return_id=True, skip_desel=True): 939 # No model set up. 940 if not hasattr(spin, 'params') or not hasattr(spin, 'model'): 941 continue 942 943 # Should R1 data be optimised? 944 r1_fit = is_r1_optimised(spin.model) 945 946 # Prepend R1 and add it to the spin container. 947 if r1_fit and 'r1' not in spin.params: 948 spin.params.insert(0, 'r1') 949 if r1_fit and not hasattr(spin, 'r1'): 950 spin.r1 = {} 951 952 # Remove the R1 parameter. 953 if not r1_fit and 'r1' in spin.params: 954 spin.params.remove('r1')
955 956
957 -def set_value(value=None, key=None, spins=None, sim_index=None, param_name=None, spin_index=None, r20_key=None):
958 """Set a specific parameter value. 959 960 @keyword value: The parameter value to set. 961 @type value: float 962 @keyword key: The key for the R2eff and I0 parameters. 963 @type key: str or None 964 @keyword spins: The list of spin data containers for the block. 965 @type spins: list of SpinContainer instances 966 @keyword sim_index: The optional MC simulation index. 967 @type sim_index: int 968 @keyword param_name: The parameter name. 969 @type param_name: str 970 @keyword spin_index: The spin index (for the cluster). 971 @type spin_index: int 972 @keyword r20_key: The unique R20 parameter key. 973 @type r20_key: str 974 """ 975 976 # Spin specific parameters. 977 if spin_index != None: 978 # Set the simulation value. 979 if sim_index != None: 980 # Get the simulation object. 981 obj = getattr(spins[spin_index], param_name+'_sim') 982 983 # R20 parameter. 984 if r20_key != None: 985 obj[sim_index][r20_key] = value 986 987 # All other parameters. 988 else: 989 if key != None: 990 obj[sim_index][key] = value 991 else: 992 obj[sim_index] = value 993 994 # Set the normal value. 995 else: 996 # Get the object. 997 obj = getattr(spins[spin_index], param_name) 998 999 # R20 parameter. 1000 if r20_key != None: 1001 obj[r20_key] = value 1002 1003 # Set the normal value. 1004 else: 1005 if key != None: 1006 obj[key] = value 1007 else: 1008 setattr(spins[spin_index], param_name, value) 1009 1010 # Cluster specific parameters. 1011 else: 1012 # Set the same parameter value for all spins in the cluster. 1013 for spin in spins: 1014 # Skip deselected spins. 1015 if not spin.select: 1016 continue 1017 1018 # Set the simulation value. 1019 if sim_index != None: 1020 sim_obj = getattr(spin, param_name+'_sim') 1021 sim_obj[sim_index] = value 1022 1023 # Set the normal value. 1024 else: 1025 setattr(spin, param_name, value)
1026