specific_fns.n_state

60 """Class containing functions for the N-state model.""" 61

62 - def __init__(self):

63 """Initialise the class by placing API_common methods into the API.""" 64 65 # Execute the base class __init__ method. 66 super(N_state_model, self).__init__() 67 68 # Place methods into the API. 69 self.model_loop = self._model_loop_single_global 70 self.overfit_deselect = self._overfit_deselect_dummy 71 self.return_conversion_factor = self._return_no_conversion_factor 72 self.set_selected_sim = self._set_selected_sim_global 73 self.sim_return_selected = self._sim_return_selected_global 74 self.test_grid_ops = self._test_grid_ops_general 75 76 # Set up the spin parameters. 77 self.PARAMS.add('r', scope='spin', units='Angstrom', desc='Bond length', py_type=float, grace_string='Bond length') 78 self.PARAMS.add('csa', scope='spin', units='ppm', desc='CSA value', py_type=float, grace_string='\\qCSA\\Q') 79 self.PARAMS.add('heteronuc_type', scope='spin', default='15N', desc='The heteronucleus type', py_type=str) 80 self.PARAMS.add('proton_type', scope='spin', default='1H', desc='The proton type', py_type=str) 81 82 # Add the minimisation data. 83 self.PARAMS.add_min_data(min_stats_global=False, min_stats_spin=True)

84 85

86 - def _assemble_param_vector(self, sim_index=None):

87 """Assemble all the parameters of the model into a single array. 88 89 @param sim_index: The index of the simulation to optimise. This should be None if 90 normal optimisation is desired. 91 @type sim_index: None or int 92 @return: The parameter vector used for optimisation. 93 @rtype: numpy array 94 """ 95 96 # Test if the model is selected. 97 if not hasattr(cdp, 'model') or not isinstance(cdp.model, str): 98 raise RelaxNoModelError 99 100 # Determine the data type. 101 data_types = self._base_data_types() 102 103 # Initialise the parameter vector. 104 param_vector = [] 105 106 # A RDC or PCS data type requires the alignment tensors to be at the start of the parameter vector (unless the tensors are fixed). 107 if ('rdc' in data_types or 'pcs' in data_types) and not all_tensors_fixed(): 108 # Loop over the alignments, adding the alignment tensor parameters to the parameter vector. 109 for i in xrange(len(cdp.align_tensors)): 110 # No alignment ID, so skip the tensor as it will not be optimised. 111 if cdp.align_tensors[i].name not in cdp.align_ids: 112 continue 113 114 # Fixed tensor. 115 if cdp.align_tensors[i].fixed: 116 continue 117 118 # Add the parameters. 119 param_vector = param_vector + list(cdp.align_tensors[i].A_5D) 120 121 # Monte Carlo simulation data structures. 122 if sim_index != None: 123 # Populations. 124 if cdp.model in ['2-domain', 'population']: 125 probs = cdp.probs_sim[sim_index] 126 127 # Euler angles. 128 if cdp.model == '2-domain': 129 alpha = cdp.alpha_sim[sim_index] 130 beta = cdp.beta_sim[sim_index] 131 gamma = cdp.gamma_sim[sim_index] 132 133 # Normal data structures. 134 else: 135 # Populations. 136 if cdp.model in ['2-domain', 'population']: 137 probs = cdp.probs 138 139 # Euler angles. 140 if cdp.model == '2-domain': 141 alpha = cdp.alpha 142 beta = cdp.beta 143 gamma = cdp.gamma 144 145 # The probabilities (exclude that of state N). 146 if cdp.model in ['2-domain', 'population']: 147 param_vector = param_vector + probs[0:-1] 148 149 # The Euler angles. 150 if cdp.model == '2-domain': 151 for i in xrange(cdp.N): 152 param_vector.append(alpha[i]) 153 param_vector.append(beta[i]) 154 param_vector.append(gamma[i]) 155 156 # The paramagnetic centre. 157 if hasattr(cdp, 'paramag_centre_fixed') and not cdp.paramag_centre_fixed: 158 if not hasattr(cdp, 'paramagnetic_centre'): 159 param_vector.append(0.0) 160 param_vector.append(0.0) 161 param_vector.append(0.0) 162 163 else: 164 param_vector.append(cdp.paramagnetic_centre[0]) 165 param_vector.append(cdp.paramagnetic_centre[1]) 166 param_vector.append(cdp.paramagnetic_centre[2]) 167 168 # Convert all None values to zero (to avoid conversion to NaN). 169 for i in xrange(len(param_vector)): 170 if param_vector[i] == None: 171 param_vector[i] = 0.0 172 173 # Return a numpy arrary. 174 return array(param_vector, float64)

175 176

177 - def _assemble_scaling_matrix(self, data_types=None, scaling=True):

178 """Create and return the scaling matrix. 179 180 @keyword data_types: The base data types used in the optimisation. This list can contain 181 the elements 'rdc', 'pcs' or 'tensor'. 182 @type data_types: list of str 183 @keyword scaling: If False, then the identity matrix will be returned. 184 @type scaling: bool 185 @return: The square and diagonal scaling matrix. 186 @rtype: numpy rank-2 array 187 """ 188 189 # Initialise. 190 scaling_matrix = identity(self._param_num(), float64) 191 192 # Return the identity matrix. 193 if not scaling: 194 return scaling_matrix 195 196 # Starting point of the populations. 197 pop_start = 0 198 if ('rdc' in data_types or 'pcs' in data_types) and not all_tensors_fixed(): 199 # Loop over the alignments. 200 tensor_num = 0 201 for i in range(len(cdp.align_tensors)): 202 # Fixed tensor. 203 if cdp.align_tensors[i].fixed: 204 continue 205 206 # Add the 5 alignment parameters. 207 pop_start = pop_start + 5 208 209 # The alignment parameters. 210 for j in range(5): 211 scaling_matrix[5*tensor_num + j, 5*tensor_num + j] = 1.0 212 213 # Increase the tensor number. 214 tensor_num += 1 215 216 # Loop over the populations, and set the scaling factor. 217 if cdp.model in ['2-domain', 'population']: 218 factor = 0.1 219 for i in xrange(pop_start, pop_start + (cdp.N-1)): 220 scaling_matrix[i, i] = factor 221 222 # The paramagnetic centre. 223 if hasattr(cdp, 'paramag_centre_fixed') and not cdp.paramag_centre_fixed: 224 for i in range(-3, 0): 225 scaling_matrix[i, i] = 1e2 226 227 # Return the matrix. 228 return scaling_matrix

229 230

231 - def _base_data_types(self):

232 """Determine all the base data types. 233 234 The base data types can include:: 235 - 'rdc', residual dipolar couplings. 236 - 'pcs', pseudo-contact shifts. 237 - 'noesy', NOE restraints. 238 - 'tensor', alignment tensors. 239 240 @return: A list of all the base data types. 241 @rtype: list of str 242 """ 243 244 # Array of data types. 245 list = [] 246 247 # RDC search. 248 for spin in spin_loop(): 249 if hasattr(spin, 'rdc'): 250 list.append('rdc') 251 break 252 253 # PCS search. 254 for spin in spin_loop(): 255 if hasattr(spin, 'pcs'): 256 list.append('pcs') 257 break 258 259 # Alignment tensor search. 260 if not ('rdc' in list or 'pcs' in list) and hasattr(cdp, 'align_tensors'): 261 list.append('tensor') 262 263 # NOESY data search. 264 if hasattr(cdp, 'noe_restraints'): 265 list.append('noesy') 266 267 # No data is present. 268 if not list: 269 raise RelaxError("Neither RDC, PCS, NOESY nor alignment tensor data is present.") 270 271 # Return the list. 272 return list

273 274

275 - def _calc_ave_dist(self, atom1, atom2, exp=1):

276 """Calculate the average distances. 277 278 The formula used is:: 279 280 _N_ 281 / 1 \ \ 1/exp 282 <r> = | - > |p1i - p2i|^exp | 283 \ N /__ / 284 i 285 286 where i are the members of the ensemble, N is the total number of structural models, and p1 287 and p2 at the two atom positions. 288 289 290 @param atom1: The atom identification string of the first atom. 291 @type atom1: str 292 @param atom2: The atom identification string of the second atom. 293 @type atom2: str 294 @keyword exp: The exponent used for the averaging, e.g. 1 for linear averaging and -6 for 295 r^-6 NOE averaging. 296 @type exp: int 297 @return: The average distance between the two atoms. 298 @rtype: float 299 """ 300 301 # Get the spin containers. 302 spin1 = return_spin(atom1) 303 spin2 = return_spin(atom2) 304 305 # Loop over each model. 306 num_models = len(spin1.pos) 307 ave_dist = 0.0 308 for i in range(num_models): 309 # Distance to the minus sixth power. 310 dist = norm(spin1.pos[i] - spin2.pos[i]) 311 ave_dist = ave_dist + dist**(exp) 312 313 # Average. 314 ave_dist = ave_dist / num_models 315 316 # The exponent. 317 ave_dist = ave_dist**(1.0/exp) 318 319 # Return the average distance. 320 return ave_dist

321 322

323 - def _CoM(self, pivot_point=None, centre=None):

324 """Centre of mass analysis. 325 326 This function does an analysis of the centre of mass (CoM) of the N states. This includes 327 calculating the order parameter associated with the pivot-CoM vector, and the associated 328 cone of motions. The pivot_point argument must be supplied. If centre is None, then the 329 CoM will be calculated from the selected parts of the loaded structure. Otherwise it will 330 be set to the centre arg. 331 332 @param pivot_point: The pivot point in the structural file(s). 333 @type pivot_point: list of float of length 3 334 @param centre: The optional centre of mass vector. 335 @type centre: list of float of length 3 336 """ 337 338 # Test if the current data pipe exists. 339 pipes.test() 340 341 # Set the pivot point. 342 cdp.pivot_point = pivot_point 343 344 # The centre has been supplied. 345 if centre: 346 cdp.CoM = centre 347 348 # Calculate from the structure file. 349 else: 350 cdp.CoM = generic_fns.structure.mass.centre_of_mass() 351 352 # Calculate the vector between the pivot and CoM points. 353 cdp.pivot_CoM = array(cdp.CoM, float64) - array(cdp.pivot_point, float64) 354 355 # Calculate the unit vector between the pivot and CoM points. 356 unit_vect = cdp.pivot_CoM / norm(cdp.pivot_CoM) 357 358 # Initilise some data structures. 359 R = zeros((3, 3), float64) 360 vectors = zeros((cdp.N, 3), float64) 361 362 # Loop over the N states. 363 for c in xrange(cdp.N): 364 # Generate the rotation matrix. 365 euler_to_R_zyz(cdp.alpha[c], cdp.beta[c], cdp.gamma[c], R) 366 367 # Rotate the unit vector. 368 vectors[c] = dot(R, unit_vect) 369 370 # Multiply by the probability. 371 vectors[c] = vectors[c] * cdp.probs[c] 372 373 # Average of the unit vectors. 374 cdp.ave_unit_pivot_CoM = sum(vectors) 375 376 # The length reduction. 377 cdp.ave_pivot_CoM_red = norm(cdp.ave_unit_pivot_CoM) 378 379 # The aveage pivot-CoM vector. 380 cdp.ave_pivot_CoM = norm(cdp.pivot_CoM) * cdp.ave_unit_pivot_CoM 381 382 # The full length rotated pivot-CoM vector. 383 cdp.full_ave_pivot_CoM = cdp.ave_pivot_CoM / cdp.ave_pivot_CoM_red 384 385 # The cone angle for diffusion on an axially symmetric cone. 386 cdp.theta_diff_on_cone = acos(cdp.ave_pivot_CoM_red) 387 cdp.S_diff_on_cone = (3.0*cos(cdp.theta_diff_on_cone)**2 - 1.0) / 2.0 388 389 # The cone angle and order parameter for diffusion in an axially symmetric cone. 390 cdp.theta_diff_in_cone = acos(2.*cdp.ave_pivot_CoM_red - 1.) 391 cdp.S_diff_in_cone = cos(cdp.theta_diff_in_cone) * (1 + cos(cdp.theta_diff_in_cone)) / 2.0 392 393 # Print out. 394 print(("\n%-40s %-20s" % ("Pivot point:", repr(cdp.pivot_point)))) 395 print(("%-40s %-20s" % ("Moving domain CoM (prior to rotation):", repr(cdp.CoM)))) 396 print(("%-40s %-20s" % ("Pivot-CoM vector", repr(cdp.pivot_CoM)))) 397 print(("%-40s %-20s" % ("Pivot-CoM unit vector:", repr(unit_vect)))) 398 print(("%-40s %-20s" % ("Average of the unit pivot-CoM vectors:", repr(cdp.ave_unit_pivot_CoM)))) 399 print(("%-40s %-20s" % ("Average of the pivot-CoM vector:", repr(cdp.ave_pivot_CoM)))) 400 print(("%-40s %-20s" % ("Full length rotated pivot-CoM vector:", repr(cdp.full_ave_pivot_CoM)))) 401 print(("%-40s %-20s" % ("Length reduction from unity:", repr(cdp.ave_pivot_CoM_red)))) 402 print(("%-40s %.5f rad (%.5f deg)" % ("Cone angle (diffusion on a cone)", cdp.theta_diff_on_cone, cdp.theta_diff_on_cone / (2*pi) *360.))) 403 print(("%-40s S_cone = %.5f (S^2 = %.5f)" % ("S_cone (diffusion on a cone)", cdp.S_diff_on_cone, cdp.S_diff_on_cone**2))) 404 print(("%-40s %.5f rad (%.5f deg)" % ("Cone angle (diffusion in a cone)", cdp.theta_diff_in_cone, cdp.theta_diff_in_cone / (2*pi) *360.))) 405 print(("%-40s S_cone = %.5f (S^2 = %.5f)" % ("S_cone (diffusion in a cone)", cdp.S_diff_in_cone, cdp.S_diff_in_cone**2))) 406 print("\n\n")

407 408

409 - def _cone_pdb(self, cone_type=None, scale=1.0, file=None, dir=None, force=False):

410 """Create a PDB file containing a geometric object representing the various cone models. 411 412 Currently the only cone types supported are 'diff in cone' and 'diff on cone'. 413 414 415 @param cone_type: The type of cone model to represent. 416 @type cone_type: str 417 @param scale: The size of the geometric object is eqaul to the average pivot-CoM 418 vector length multiplied by this scaling factor. 419 @type scale: float 420 @param file: The name of the PDB file to create. 421 @type file: str 422 @param dir: The name of the directory to place the PDB file into. 423 @type dir: str 424 @param force: Flag which if set to True will cause any pre-existing file to be 425 overwritten. 426 @type force: int 427 """ 428 429 # Test if the cone models have been determined. 430 if cone_type == 'diff in cone': 431 if not hasattr(cdp, 'S_diff_in_cone'): 432 raise RelaxError("The diffusion in a cone model has not yet been determined.") 433 elif cone_type == 'diff on cone': 434 if not hasattr(cdp, 'S_diff_on_cone'): 435 raise RelaxError("The diffusion on a cone model has not yet been determined.") 436 else: 437 raise RelaxError("The cone type " + repr(cone_type) + " is unknown.") 438 439 # The number of increments for the filling of the cone objects. 440 inc = 20 441 442 # The rotation matrix. 443 R = zeros((3, 3), float64) 444 two_vect_to_R(array([0, 0, 1], float64), cdp.ave_pivot_CoM/norm(cdp.ave_pivot_CoM), R) 445 446 # The isotropic cone object. 447 if cone_type == 'diff in cone': 448 angle = cdp.theta_diff_in_cone 449 elif cone_type == 'diff on cone': 450 angle = cdp.theta_diff_on_cone 451 cone = Iso_cone(angle) 452 453 # Create the structural object. 454 structure = Internal() 455 456 # Add a structure. 457 structure.add_molecule(name='cone') 458 459 # Alias the single molecule from the single model. 460 mol = structure.structural_data[0].mol[0] 461 462 # Add the pivot point. 463 mol.atom_add(pdb_record='HETATM', atom_num=1, atom_name='R', res_name='PIV', res_num=1, pos=cdp.pivot_point, element='C') 464 465 # Generate the average pivot-CoM vectors. 466 print("\nGenerating the average pivot-CoM vectors.") 467 sim_vectors = None 468 if hasattr(cdp, 'ave_pivot_CoM_sim'): 469 sim_vectors = cdp.ave_pivot_CoM_sim 470 res_num = generic_fns.structure.geometric.generate_vector_residues(mol=mol, vector=cdp.ave_pivot_CoM, atom_name='Ave', res_name_vect='AVE', sim_vectors=sim_vectors, res_num=2, origin=cdp.pivot_point, scale=scale) 471 472 # Generate the cone outer edge. 473 print("\nGenerating the cone outer edge.") 474 cap_start_atom = mol.atom_num[-1]+1 475 generic_fns.structure.geometric.cone_edge(mol=mol, cone=cone, res_name='CON', res_num=3, apex=cdp.pivot_point, R=R, scale=norm(cdp.pivot_CoM), inc=inc) 476 477 # Generate the cone cap, and stitch it to the cone edge. 478 if cone_type == 'diff in cone': 479 print("\nGenerating the cone cap.") 480 cone_start_atom = mol.atom_num[-1]+1 481 generic_fns.structure.geometric.generate_vector_dist(mol=mol, res_name='CON', res_num=3, centre=cdp.pivot_point, R=R, limit_check=cone.limit_check, scale=norm(cdp.pivot_CoM), inc=inc) 482 generic_fns.structure.geometric.stitch_cone_to_edge(mol=mol, cone=cone, dome_start=cone_start_atom, edge_start=cap_start_atom+1, inc=inc) 483 484 # Create the PDB file. 485 print("\nGenerating the PDB file.") 486 pdb_file = open_write_file(file, dir, force=force) 487 structure.write_pdb(pdb_file) 488 pdb_file.close()

489 490

491 - def _disassemble_param_vector(self, param_vector=None, data_types=None, sim_index=None):

492 """Disassemble the parameter vector and place the values into the relevant variables. 493 494 For the 2-domain N-state model, the parameters are stored in the probability and Euler angle data structures. For the population N-state model, only the probabilities are stored. If RDCs are present and alignment tensors are optimised, then these are stored as well. 495 496 @keyword data_types: The base data types used in the optimisation. This list can contain the elements 'rdc', 'pcs' or 'tensor'. 497 @type data_types: list of str 498 @keyword param_vector: The parameter vector returned from optimisation. 499 @type param_vector: numpy array 500 @keyword sim_index: The index of the simulation to optimise. This should be None if normal optimisation is desired. 501 @type sim_index: None or int 502 """ 503 504 # Test if the model is selected. 505 if not hasattr(cdp, 'model') or not isinstance(cdp.model, str): 506 raise RelaxNoModelError 507 508 # Unpack and strip off the alignment tensor parameters. 509 if ('rdc' in data_types or 'pcs' in data_types) and not all_tensors_fixed(): 510 # Loop over the alignments, adding the alignment tensor parameters to the tensor data container. 511 tensor_num = 0 512 for i in xrange(len(cdp.align_tensors)): 513 # No alignment ID, so skip the tensor as it will not be optimised. 514 if cdp.align_tensors[i].name not in cdp.align_ids: 515 continue 516 517 # Fixed tensor. 518 if cdp.align_tensors[i].fixed: 519 continue 520 521 # Normal tensors. 522 if sim_index == None: 523 cdp.align_tensors[i].Axx = param_vector[5*tensor_num] 524 cdp.align_tensors[i].Ayy = param_vector[5*tensor_num+1] 525 cdp.align_tensors[i].Axy = param_vector[5*tensor_num+2] 526 cdp.align_tensors[i].Axz = param_vector[5*tensor_num+3] 527 cdp.align_tensors[i].Ayz = param_vector[5*tensor_num+4] 528 529 # Monte Carlo simulated tensors. 530 else: 531 cdp.align_tensors[i].Axx_sim[sim_index] = param_vector[5*tensor_num] 532 cdp.align_tensors[i].Ayy_sim[sim_index] = param_vector[5*tensor_num+1] 533 cdp.align_tensors[i].Axy_sim[sim_index] = param_vector[5*tensor_num+2] 534 cdp.align_tensors[i].Axz_sim[sim_index] = param_vector[5*tensor_num+3] 535 cdp.align_tensors[i].Ayz_sim[sim_index] = param_vector[5*tensor_num+4] 536 537 # Increase the tensor number. 538 tensor_num += 1 539 540 # Create a new parameter vector without the tensors. 541 param_vector = param_vector[5*tensor_num:] 542 543 # Monte Carlo simulation data structures. 544 if sim_index != None: 545 # Populations. 546 if cdp.model in ['2-domain', 'population']: 547 probs = cdp.probs_sim[sim_index] 548 549 # Euler angles. 550 if cdp.model == '2-domain': 551 alpha = cdp.alpha_sim[sim_index] 552 beta = cdp.beta_sim[sim_index] 553 gamma = cdp.gamma_sim[sim_index] 554 555 # Normal data structures. 556 else: 557 # Populations. 558 if cdp.model in ['2-domain', 'population']: 559 probs = cdp.probs 560 561 # Euler angles. 562 if cdp.model == '2-domain': 563 alpha = cdp.alpha 564 beta = cdp.beta 565 gamma = cdp.gamma 566 567 # The probabilities for states 0 to N-1. 568 if cdp.model in ['2-domain', 'population']: 569 for i in xrange(cdp.N-1): 570 probs[i] = param_vector[i] 571 572 # The probability for state N. 573 probs[-1] = 1 - sum(probs[0:-1]) 574 575 # The Euler angles. 576 if cdp.model == '2-domain': 577 for i in xrange(cdp.N): 578 alpha[i] = param_vector[cdp.N-1 + 3*i] 579 beta[i] = param_vector[cdp.N-1 + 3*i + 1] 580 gamma[i] = param_vector[cdp.N-1 + 3*i + 2] 581 582 # The paramagnetic centre. 583 if hasattr(cdp, 'paramag_centre_fixed') and not cdp.paramag_centre_fixed: 584 # Create the structure if needed. 585 if not hasattr(cdp, 'paramagnetic_centre'): 586 cdp.paramagnetic_centre = zeros(3, float64) 587 588 # The position. 589 cdp.paramagnetic_centre[0] = param_vector[-3] 590 cdp.paramagnetic_centre[1] = param_vector[-2] 591 cdp.paramagnetic_centre[2] = param_vector[-1]

592 593

594 - def _elim_no_prob(self):

595 """Remove all structures or states which have no probability.""" 596 597 # Test if the current data pipe exists. 598 pipes.test() 599 600 # Test if the model is setup. 601 if not hasattr(cdp, 'model'): 602 raise RelaxNoModelError('N-state') 603 604 # Test if there are populations. 605 if not hasattr(cdp, 'probs'): 606 raise RelaxError("The N-state model populations do not exist.") 607 608 # Loop over the structures. 609 i = 0 610 while True: 611 # End condition. 612 if i == cdp.N - 1: 613 break 614 615 # No probability. 616 if cdp.probs[i] < 1e-5: 617 # Remove the probability. 618 cdp.probs.pop(i) 619 620 # Remove the structure. 621 cdp.structure.structural_data.pop(i) 622 623 # Eliminate bond vectors. 624 for spin in spin_loop(): 625 # Position info. 626 if hasattr(spin, 'pos'): 627 spin.pos.pop(i) 628 629 # Vector info. 630 if hasattr(spin, 'xh_vect'): 631 spin.xh_vect.pop(i) 632 if hasattr(spin, 'bond_vect'): 633 spin.bond_vect.pop(i) 634 635 # Update N. 636 cdp.N -= 1 637 638 # Start the loop again without incrementing i. 639 continue 640 641 # Increment i. 642 i += 1

643 644

645 - def _linear_constraints(self, data_types=None, scaling_matrix=None):

646 """Function for setting up the linear constraint matrices A and b. 647 648 Standard notation 649 ================= 650 651 The N-state model constraints are:: 652 653 0 <= pc <= 1, 654 655 where p is the probability and c corresponds to state c. 656 657 658 Matrix notation 659 =============== 660 661 In the notation A.x >= b, where A is an matrix of coefficients, x is an array of parameter 662 values, and b is a vector of scalars, these inequality constraints are:: 663 664 | 1 0 0 | | 0 | 665 | | | | 666 |-1 0 0 | | -1 | 667 | | | | 668 | 0 1 0 | | 0 | 669 | | | p0 | | | 670 | 0 -1 0 | | | | -1 | 671 | | . | p1 | >= | | 672 | 0 0 1 | | | | 0 | 673 | | | p2 | | | 674 | 0 0 -1 | | -1 | 675 | | | | 676 |-1 -1 -1 | | -1 | 677 | | | | 678 | 1 1 1 | | 0 | 679 680 This example is for a 4-state model, the last probability pn is not included as this 681 parameter does not exist (because the sum of pc is equal to 1). The Euler angle parameters 682 have been excluded here but will be included in the returned A and b objects. These 683 parameters simply add columns of zero to the A matrix and have no effect on b. The last two 684 rows correspond to the inequality:: 685 686 0 <= pN <= 1. 687 688 As:: 689 N-1 690 \ 691 pN = 1 - > pc, 692 /__ 693 c=1 694 695 then:: 696 697 -p1 - p2 - ... - p(N-1) >= -1, 698 699 p1 + p2 + ... + p(N-1) >= 0. 700 701 702 @keyword data_types: The base data types used in the optimisation. This list can 703 contain the elements 'rdc', 'pcs' or 'tensor'. 704 @type data_types: list of str 705 @keyword scaling_matrix: The diagonal scaling matrix. 706 @type scaling_matrix: numpy rank-2 square matrix 707 @return: The matrices A and b. 708 @rtype: tuple of len 2 of a numpy rank-2, size NxM matrix and numpy 709 rank-1, size N array 710 """ 711 712 # Starting point of the populations. 713 pop_start = 0 714 if ('rdc' in data_types or 'pcs' in data_types) and not all_tensors_fixed(): 715 # Loop over the alignments. 716 for i in range(len(cdp.align_tensors)): 717 # Fixed tensor. 718 if cdp.align_tensors[i].fixed: 719 continue 720 721 # Add 5 parameters. 722 pop_start += 5 723 724 # Initialisation (0..j..m). 725 A = [] 726 b = [] 727 zero_array = zeros(self._param_num(), float64) 728 i = pop_start 729 j = 0 730 731 # Probability parameters. 732 if cdp.model in ['2-domain', 'population']: 733 # Loop over the prob parameters (N - 1, because the sum of pc is 1). 734 for k in xrange(cdp.N - 1): 735 # 0 <= pc <= 1. 736 A.append(zero_array * 0.0) 737 A.append(zero_array * 0.0) 738 A[j][i] = 1.0 739 A[j+1][i] = -1.0 740 b.append(0.0) 741 b.append(-1.0 / scaling_matrix[i, i]) 742 j = j + 2 743 744 # Increment i. 745 i = i + 1 746 747 # Add the inequalities for pN. 748 A.append(zero_array * 0.0) 749 A.append(zero_array * 0.0) 750 for i in xrange(pop_start, self._param_num()): 751 A[-2][i] = -1.0 752 A[-1][i] = 1.0 753 b.append(-1.0 / scaling_matrix[i, i]) 754 b.append(0.0) 755 756 # Convert to numpy data structures. 757 A = array(A, float64) 758 b = array(b, float64) 759 760 # Return the constraint objects. 761 return A, b

762 763

764 - def _minimise_bc_data(self, model):

765 """Extract and unpack the back calculated data. 766 767 @param model: The instantiated class containing the target function. 768 @type model: class instance 769 """ 770 771 # No alignment tensors, so nothing to do. 772 if not hasattr(cdp, 'align_tensors'): 773 return 774 775 # Loop over each alignment. 776 for i in xrange(len(cdp.align_ids)): 777 # Fixed tensor. 778 if cdp.align_tensors[i].fixed: 779 continue 780 781 # The alignment ID. 782 align_id = cdp.align_ids[i] 783 784 # Data flags 785 rdc_flag = False 786 if hasattr(cdp, 'rdc_ids') and align_id in cdp.rdc_ids: 787 rdc_flag = True 788 pcs_flag = False 789 if hasattr(cdp, 'pcs_ids') and align_id in cdp.pcs_ids: 790 pcs_flag = True 791 792 # Spin loop. 793 data_index = 0 794 for spin in spin_loop(): 795 # Skip deselected spins. 796 if not spin.select: 797 continue 798 799 # Spins with PCS data. 800 if pcs_flag and hasattr(spin, 'pcs'): 801 # Initialise the data structure if necessary. 802 if not hasattr(spin, 'pcs_bc'): 803 spin.pcs_bc = {} 804 805 # Add the back calculated PCS (in ppm). 806 spin.pcs_bc[align_id] = model.deltaij_theta[i, data_index] * 1e6 807 808 # Spins with RDC data. 809 if rdc_flag and hasattr(spin, 'rdc') and (hasattr(spin, 'xh_vect') or hasattr(spin, 'bond_vect')): 810 # Initialise the data structure if necessary. 811 if not hasattr(spin, 'rdc_bc'): 812 spin.rdc_bc = {} 813 814 # Append the back calculated PCS. 815 spin.rdc_bc[align_id] = model.Dij_theta[i, data_index] 816 817 # Increment the spin index if it contains data. 818 if hasattr(spin, 'pcs') or (hasattr(spin, 'rdc') and (hasattr(spin, 'xh_vect') or hasattr(spin, 'bond_vect'))): 819 data_index = data_index + 1

820 821

822 - def _minimise_setup_atomic_pos(self):

823 """Set up the atomic position data structures for optimisation using PCSs and PREs as base data sets. 824 825 @return: The atomic positions (the first index is the spins, the second is the structures, and the third is the atomic coordinates) and the paramagnetic centre. 826 @rtype: numpy rank-3 array, numpy rank-1 array. 827 """ 828 829 # Initialise. 830 atomic_pos = [] 831 832 # Store the atomic positions. 833 for spin in spin_loop(): 834 # Skip deselected spins. 835 if not spin.select: 836 continue 837 838 # Only use spins with alignment/paramagnetic data. 839 if not hasattr(spin, 'pcs') and not hasattr(spin, 'rdc') and not hasattr(spin, 'pre'): 840 continue 841 842 # The position list. 843 if type(spin.pos[0]) in [float, float64]: 844 atomic_pos.append([spin.pos]) 845 else: 846 atomic_pos.append(spin.pos) 847 848 # Convert to numpy objects. 849 atomic_pos = array(atomic_pos, float64) 850 851 # The paramagnetic centre. 852 if hasattr(cdp, 'paramagnetic_centre'): 853 paramag_centre = array(cdp.paramagnetic_centre) 854 else: 855 paramag_centre = zeros(3, float64) 856 857 # Return the data structures. 858 return atomic_pos, paramag_centre

859 860

861 - def _minimise_setup_pcs(self, sim_index=None):

862 """Set up the data structures for optimisation using PCSs as base data sets. 863 864 @keyword sim_index: The index of the simulation to optimise. This should be None if normal optimisation is desired. 865 @type sim_index: None or int 866 @return: The assembled data structures for using PCSs as the base data for optimisation. These include: 867 - the PCS values. 868 - the unit vectors connecting the paramagnetic centre (the electron spin) to 869 - the PCS weight. 870 - the nuclear spin. 871 - the pseudocontact shift constants. 872 @rtype: tuple of (numpy rank-2 array, numpy rank-2 array, numpy rank-2 array, numpy rank-1 array, numpy rank-1 array) 873 """ 874 875 # Data setup tests. 876 if not hasattr(cdp, 'paramagnetic_centre') and (hasattr(cdp, 'paramag_centre_fixed') and cdp.paramag_centre_fixed): 877 raise RelaxError("The paramagnetic centre has not yet been specified.") 878 if not hasattr(cdp, 'temperature'): 879 raise RelaxError("The experimental temperatures have not been set.") 880 if not hasattr(cdp, 'frq'): 881 raise RelaxError("The spectrometer frequencies of the experiments have not been set.") 882 883 # Initialise. 884 pcs = [] 885 pcs_err = [] 886 pcs_weight = [] 887 temp = [] 888 frq = [] 889 890 # The PCS data. 891 for align_id in cdp.align_ids: 892 # No RDC or PCS data, so jump to the next alignment. 893 if (hasattr(cdp, 'rdc_ids') and not align_id in cdp.rdc_ids) and (hasattr(cdp, 'pcs_ids') and not align_id in cdp.pcs_ids): 894 continue 895 896 # Append empty arrays to the PCS structures. 897 pcs.append([]) 898 pcs_err.append([]) 899 pcs_weight.append([]) 900 901 # Get the temperature for the PCS constant. 902 if cdp.temperature.has_key(align_id): 903 temp.append(cdp.temperature[align_id]) 904 else: 905 temp.append(0.0) 906 907 # Get the spectrometer frequency in Tesla units for the PCS constant. 908 if cdp.frq.has_key(align_id): 909 frq.append(cdp.frq[align_id] * 2.0 * pi / g1H) 910 else: 911 frq.append(1e-10) 912 913 # Spin loop. 914 j = 0 915 for spin in spin_loop(): 916 # Skip deselected spins. 917 if not spin.select: 918 continue 919 920 # Skip spins without PCS data. 921 if not hasattr(spin, 'pcs'): 922 # Add rows of None if other alignment data exists. 923 if hasattr(spin, 'rdc'): 924 pcs[-1].append(None) 925 pcs_err[-1].append(None) 926 pcs_weight[-1].append(None) 927 j = j + 1 928 929 # Jump to the next spin. 930 continue 931 932 # Append the PCSs to the list. 933 if align_id in spin.pcs.keys(): 934 if sim_index != None: 935 pcs[-1].append(spin.pcs_sim[align_id][sim_index]) 936 else: 937 pcs[-1].append(spin.pcs[align_id]) 938 else: 939 pcs[-1].append(None) 940 941 # Append the PCS errors. 942 if hasattr(spin, 'pcs_err') and align_id in spin.pcs_err.keys(): 943 pcs_err[-1].append(spin.pcs_err[align_id]) 944 else: 945 pcs_err[-1].append(None) 946 947 # Append the weight. 948 if hasattr(spin, 'pcs_weight') and align_id in spin.pcs_weight.keys(): 949 pcs_weight[-1].append(spin.pcs_weight[align_id]) 950 else: 951 pcs_weight[-1].append(1.0) 952 953 # Spin index. 954 j = j + 1 955 956 # Convert to numpy objects. 957 pcs = array(pcs, float64) 958 pcs_err = array(pcs_err, float64) 959 pcs_weight = array(pcs_weight, float64) 960 961 # Convert the PCS from ppm to no units. 962 pcs = pcs * 1e-6 963 pcs_err = pcs_err * 1e-6 964 965 # Return the data structures. 966 return pcs, pcs_err, pcs_weight, temp, frq

967 968

969 - def _minimise_setup_rdcs(self, sim_index=None):

970 """Set up the data structures for optimisation using RDCs as base data sets. 971 972 @keyword sim_index: The index of the simulation to optimise. This should be None if normal optimisation is desired. 973 @type sim_index: None or int 974 @return: The assembled data structures for using RDCs as the base data for optimisation. These include: 975 - rdc, the RDC values. 976 - rdc_err, the RDC errors. 977 - rdc_weight, the RDC weights. 978 - vectors, the heteronucleus to proton vectors. 979 - rdc_const, the dipolar constants. 980 @rtype: tuple of (numpy rank-2 array, numpy rank-2 array, numpy rank-2 array) 981 """ 982 983 # Initialise. 984 rdc = [] 985 rdc_err = [] 986 rdc_weight = [] 987 unit_vect = [] 988 rdc_const = [] 989 990 # The unit vectors and RDC constants. 991 for spin, spin_id in spin_loop(return_id=True): 992 # Skip deselected spins. 993 if not spin.select: 994 continue 995 996 # Only use spins with RDC data. 997 if not hasattr(spin, 'rdc'): 998 # Add rows of None if other alignment data exists. 999 if hasattr(spin, 'pcs'): 1000 unit_vect.append(None) 1001 rdc_const.append(None) 1002 1003 # Jump to the next spin. 1004 continue 1005 1006 # RDC data exists but the XH bond vectors are missing? 1007 if not hasattr(spin, 'members') and not hasattr(spin, 'xh_vect') and not hasattr(spin, 'bond_vect'): 1008 # Throw a warning. 1009 warn(RelaxWarning("RDC data exists but the XH bond vectors are missing, skipping spin %s." % spin_id)) 1010 1011 # Add rows of None if other data exists. 1012 if hasattr(spin, 'pcs'): 1013 unit_vect.append(None) 1014 rdc_const.append(None) 1015 1016 # Jump to the next spin. 1017 continue 1018 1019 # Pseudo-atom set up. 1020 if hasattr(spin, 'members'): 1021 # Skip non-Me groups. 1022 if len(spin.members) != 3: 1023 warn(RelaxWarning("Only methyl group pseudo atoms are supported due to their fast rotation, skipping spin %s." % spin_id)) 1024 continue 1025 1026 # The summed vector. 1027 vect = zeros(3, float64) 1028 for i in range(3): 1029 # Get the spin. 1030 spin_i = return_spin(spin.members[i]) 1031 1032 # Add the bond vector. 1033 if hasattr(spin_i, 'xh_vect'): 1034 obj = getattr(spin_i, 'xh_vect') 1035 else: 1036 obj = getattr(spin_i, 'bond_vect') 1037 vect = vect + obj 1038 1039 # Normalise. 1040 vect = vect / norm(vect) 1041 1042 # Normal spin set up. 1043 else: 1044 # Append the RDC and XH vectors to the lists. 1045 if hasattr(spin, 'xh_vect'): 1046 vect = getattr(spin, 'xh_vect') 1047 else: 1048 vect = getattr(spin, 'bond_vect') 1049 1050 # Add the bond vectors. 1051 if isinstance(vect[0], float): 1052 unit_vect.append([vect]) 1053 else: 1054 unit_vect.append(vect) 1055 1056 # Checks. 1057 if not hasattr(spin, 'heteronuc_type'): 1058 raise RelaxSpinTypeError 1059 if not hasattr(spin, 'proton_type'): 1060 raise RelaxProtonTypeError 1061 if not hasattr(spin, 'r'): 1062 raise RelaxNoValueError("bond length") 1063 1064 # Gyromagnetic ratios. 1065 gx = return_gyromagnetic_ratio(spin.heteronuc_type) 1066 gh = return_gyromagnetic_ratio(spin.proton_type) 1067 1068 # Calculate the RDC dipolar constant (in Hertz, and the 3 comes from the alignment tensor), and append it to the list. 1069 rdc_const.append(3.0/(2.0*pi) * dipolar_constant(gx, gh, spin.r)) 1070 1071 # Fix the unit vector data structure. 1072 num = None 1073 for i in range(len(unit_vect)): 1074 # Number of vectors. 1075 if num == None: 1076 if unit_vect[i] != None: 1077 num = len(unit_vect[i]) 1078 continue 1079 1080 # Check. 1081 if unit_vect[i] != None and len(unit_vect[i]) != num: 1082 raise RelaxError, "The number of bond vectors for all spins do no match:\n%s" % unit_vect 1083 1084 # Missing unit vectors. 1085 if num == None: 1086 raise RelaxError, "No bond vectors could be found." 1087 1088 # Update None entries. 1089 for i in range(len(unit_vect)): 1090 if unit_vect[i] == None: 1091 unit_vect[i] = [[None, None, None]]*num 1092 1093 # The RDC data. 1094 for align_id in cdp.align_ids: 1095 # No RDC or PCS data, so jump to the next alignment. 1096 if (hasattr(cdp, 'rdc_ids') and not align_id in cdp.rdc_ids) and (hasattr(cdp, 'pcs_ids') and not align_id in cdp.pcs_ids): 1097 continue 1098 1099 # Append empty arrays to the RDC structures. 1100 rdc.append([]) 1101 rdc_err.append([]) 1102 rdc_weight.append([]) 1103 1104 # Spin loop. 1105 for spin in spin_loop(): 1106 # Skip deselected spins. 1107 if not spin.select: 1108 continue 1109 1110 # Skip spins without RDC data or XH bond vectors. 1111 if not hasattr(spin, 'rdc') or (not hasattr(spin, 'members') and not hasattr(spin, 'xh_vect') and not hasattr(spin, 'bond_vect')): 1112 # Add rows of None if other alignment data exists. 1113 if hasattr(spin, 'pcs'): 1114 rdc[-1].append(None) 1115 rdc_err[-1].append(None) 1116 rdc_weight[-1].append(None) 1117 1118 # Jump to the next spin. 1119 continue 1120 1121 # Defaults of None. 1122 value = None 1123 error = None 1124 1125 # Pseudo-atom set up. 1126 if hasattr(spin, 'members') and align_id in spin.rdc.keys(): 1127 # Skip non-Me groups. 1128 if len(spin.members) != 3: 1129 continue 1130 1131 # The RDC for the Me-pseudo spin where: 1132 # <D> = -1/3 Dpar. 1133 # See Verdier, et al., JMR, 2003, 163, 353-359. 1134 if sim_index != None: 1135 value = -3.0 * spin.rdc_sim[align_id][sim_index] 1136 else: 1137 value = -3.0 * spin.rdc[align_id] 1138 1139 # The error. 1140 if hasattr(spin, 'rdc_err') and align_id in spin.rdc_err.keys(): 1141 error = -3.0 * spin.rdc_err[align_id] 1142 1143 # Normal spin set up. 1144 elif align_id in spin.rdc.keys(): 1145 # The RDC. 1146 if sim_index != None: 1147 value = spin.rdc_sim[align_id][sim_index] 1148 else: 1149 value = spin.rdc[align_id] 1150 1151 # The error. 1152 if hasattr(spin, 'rdc_err') and align_id in spin.rdc_err.keys(): 1153 error = spin.rdc_err[align_id] 1154 1155 # Append the RDCs to the list. 1156 rdc[-1].append(value) 1157 1158 # Append the RDC errors. 1159 rdc_err[-1].append(error) 1160 1161 # Append the weight. 1162 if hasattr(spin, 'rdc_weight') and align_id in spin.rdc_weight.keys(): 1163 rdc_weight[-1].append(spin.rdc_weight[align_id]) 1164 else: 1165 rdc_weight[-1].append(1.0) 1166 1167 # Convert to numpy objects. 1168 rdc = array(rdc, float64) 1169 rdc_err = array(rdc_err, float64) 1170 rdc_weight = array(rdc_weight, float64) 1171 unit_vect = array(unit_vect, float64) 1172 rdc_const = array(rdc_const, float64) 1173 1174 # Return the data structures. 1175 return rdc, rdc_err, rdc_weight, unit_vect, rdc_const

1176 1177

1178 - def _minimise_setup_tensors(self, sim_index=None):

1179 """Set up the data structures for optimisation using alignment tensors as base data sets. 1180 1181 @keyword sim_index: The index of the simulation to optimise. This should be None if 1182 normal optimisation is desired. 1183 @type sim_index: None or int 1184 @return: The assembled data structures for using alignment tensors as the base 1185 data for optimisation. These include: 1186 - full_tensors, the data of the full alignment tensors. 1187 - red_tensor_elem, the tensors as concatenated rank-1 5D arrays. 1188 - red_tensor_err, the tensor errors as concatenated rank-1 5D 1189 arrays. 1190 - full_in_ref_frame, flags specifying if the tensor in the reference 1191 frame is the full or reduced tensor. 1192 @rtype: tuple of (list, numpy rank-1 array, numpy rank-1 array, numpy rank-1 1193 array) 1194 """ 1195 1196 # Initialise. 1197 n = len(cdp.align_tensors.reduction) 1198 full_tensors = zeros(n*5, float64) 1199 red_tensors = zeros(n*5, float64) 1200 red_err = ones(n*5, float64) * 1e-5 1201 full_in_ref_frame = zeros(n, float64) 1202 1203 # Loop over the full tensors. 1204 for i, tensor in self._tensor_loop(red=False): 1205 # The full tensor. 1206 full_tensors[5*i + 0] = tensor.Axx 1207 full_tensors[5*i + 1] = tensor.Ayy 1208 full_tensors[5*i + 2] = tensor.Axy 1209 full_tensors[5*i + 3] = tensor.Axz 1210 full_tensors[5*i + 4] = tensor.Ayz 1211 1212 # The full tensor corresponds to the frame of reference. 1213 if cdp.ref_domain == tensor.domain: 1214 full_in_ref_frame[i] = 1 1215 1216 # Loop over the reduced tensors. 1217 for i, tensor in self._tensor_loop(red=True): 1218 # The reduced tensor (simulation data). 1219 if sim_index != None: 1220 red_tensors[5*i + 0] = tensor.Axx_sim[sim_index] 1221 red_tensors[5*i + 1] = tensor.Ayy_sim[sim_index] 1222 red_tensors[5*i + 2] = tensor.Axy_sim[sim_index] 1223 red_tensors[5*i + 3] = tensor.Axz_sim[sim_index] 1224 red_tensors[5*i + 4] = tensor.Ayz_sim[sim_index] 1225 1226 # The reduced tensor. 1227 else: 1228 red_tensors[5*i + 0] = tensor.Axx 1229 red_tensors[5*i + 1] = tensor.Ayy 1230 red_tensors[5*i + 2] = tensor.Axy 1231 red_tensors[5*i + 3] = tensor.Axz 1232 red_tensors[5*i + 4] = tensor.Ayz 1233 1234 # The reduced tensor errors. 1235 if hasattr(tensor, 'Axx_err'): 1236 red_err[5*i + 0] = tensor.Axx_err 1237 red_err[5*i + 1] = tensor.Ayy_err 1238 red_err[5*i + 2] = tensor.Axy_err 1239 red_err[5*i + 3] = tensor.Axz_err 1240 red_err[5*i + 4] = tensor.Ayz_err 1241 1242 # Return the data structures. 1243 return full_tensors, red_tensors, red_err, full_in_ref_frame

1244 1245

1246 - def _minimise_setup_fixed_tensors(self):

1247 """Set up the data structures for the fixed alignment tensors. 1248 1249 @return: The assembled data structures for the fixed alignment tensors. 1250 @rtype: numpy rank-1 array. 1251 """ 1252 1253 # Initialise. 1254 n = num_tensors(skip_fixed=False) - num_tensors(skip_fixed=True) 1255 tensors = zeros(n*5, float64) 1256 1257 # Nothing to do. 1258 if n == 0: 1259 return None 1260 1261 # Loop over the tensors. 1262 index = 0 1263 for i in range(len(cdp.align_tensors)): 1264 # Skip unfixed tensors. 1265 if not cdp.align_tensors[i].fixed: 1266 continue 1267 1268 # The real tensors. 1269 tensors[5*index + 0] = cdp.align_tensors[i].Axx 1270 tensors[5*index + 1] = cdp.align_tensors[i].Ayy 1271 tensors[5*index + 2] = cdp.align_tensors[i].Axy 1272 tensors[5*index + 3] = cdp.align_tensors[i].Axz 1273 tensors[5*index + 4] = cdp.align_tensors[i].Ayz 1274 1275 # Increment the index. 1276 index += 1 1277 1278 # Return the data structure. 1279 return tensors

1280 1281

1282 - def _num_data_points(self):

1283 """Determine the number of data points used in the model. 1284 1285 @return: The number, n, of data points in the model. 1286 @rtype: int 1287 """ 1288 1289 # Determine the data type. 1290 data_types = self._base_data_types() 1291 1292 # Init. 1293 n = 0 1294 1295 # Spin loop. 1296 for spin in spin_loop(): 1297 # Skip deselected spins. 1298 if not spin.select: 1299 continue 1300 1301 # RDC data (skipping array elements set to None). 1302 if 'rdc' in data_types: 1303 if hasattr(spin, 'rdc'): 1304 for rdc in spin.rdc: 1305 if isinstance(rdc, float): 1306 n = n + 1 1307 1308 # PCS data (skipping array elements set to None). 1309 if 'pcs' in data_types: 1310 if hasattr(spin, 'pcs'): 1311 for pcs in spin.pcs: 1312 if isinstance(pcs, float): 1313 n = n + 1 1314 1315 # Alignment tensors. 1316 if 'tensor' in data_types: 1317 n = n + 5*len(cdp.align_tensors) 1318 1319 # Return the value. 1320 return n

1321 1322

1323 - def _number_of_states(self, N=None):

1324 """Set the number of states in the N-state model. 1325 1326 @param N: The number of states. 1327 @type N: int 1328 """ 1329 1330 # Test if the current data pipe exists. 1331 pipes.test() 1332 1333 # Test if the model is setup. 1334 if not hasattr(cdp, 'model'): 1335 raise RelaxNoModelError('N-state') 1336 1337 # Set the value of N. 1338 cdp.N = N 1339 1340 # Update the model. 1341 self._update_model()

1342 1343

1344 - def _param_model_index(self, param):

1345 """Return the N-state model index for the given parameter. 1346 1347 @param param: The N-state model parameter. 1348 @type param: str 1349 @return: The N-state model index. 1350 @rtype: str 1351 """ 1352 1353 # Probability. 1354 if search('^p[0-9]*$', param): 1355 return int(param[1:]) 1356 1357 # Alpha Euler angle. 1358 if search('^alpha', param): 1359 return int(param[5:]) 1360 1361 # Beta Euler angle. 1362 if search('^beta', param): 1363 return int(param[4:]) 1364 1365 # Gamma Euler angle. 1366 if search('^gamma', param): 1367 return int(param[5:]) 1368 1369 # Model independent parameter. 1370 return None

1371 1372

1373 - def _param_num(self):

1374 """Determine the number of parameters in the model. 1375 1376 @return: The number of model parameters. 1377 @rtype: int 1378 """ 1379 1380 # Determine the data type. 1381 data_types = self._base_data_types() 1382 1383 # Init. 1384 num = 0 1385 1386 # Alignment tensor params. 1387 if ('rdc' in data_types or 'pcs' in data_types) and not all_tensors_fixed(): 1388 # Loop over the alignments. 1389 for i in xrange(len(cdp.align_tensors)): 1390 # No alignment ID, so skip the tensor as it is not part of the parameter set. 1391 if cdp.align_tensors[i].name not in cdp.align_ids: 1392 continue 1393 1394 # Fixed tensor. 1395 if cdp.align_tensors[i].fixed: 1396 continue 1397 1398 # Add 5 tensor parameters. 1399 num += 5 1400 1401 # Populations. 1402 if cdp.model in ['2-domain', 'population']: 1403 num = num + (cdp.N - 1) 1404 1405 # Euler angles. 1406 if cdp.model == '2-domain': 1407 num = num + 3*cdp.N 1408 1409 # The paramagnetic centre. 1410 if hasattr(cdp, 'paramag_centre_fixed') and not cdp.paramag_centre_fixed: 1411 num = num + 3 1412 1413 # Return the param number. 1414 return num

1415 1416

1417 - def _ref_domain(self, ref=None):

1418 """Set the reference domain for the '2-domain' N-state model. 1419 1420 @param ref: The reference domain. 1421 @type ref: str 1422 """ 1423 1424 # Test if the current data pipe exists. 1425 pipes.test() 1426 1427 # Test if the model is setup. 1428 if not hasattr(cdp, 'model'): 1429 raise RelaxNoModelError('N-state') 1430 1431 # Test that the correct model is set. 1432 if cdp.model != '2-domain': 1433 raise RelaxError("Setting the reference domain is only possible for the '2-domain' N-state model.") 1434 1435 # Test if the reference domain exists. 1436 exists = False 1437 for tensor_cont in cdp.align_tensors: 1438 if tensor_cont.domain == ref: 1439 exists = True 1440 if not exists: 1441 raise RelaxError("The reference domain cannot be found within any of the loaded tensors.") 1442 1443 # Set the reference domain. 1444 cdp.ref_domain = ref 1445 1446 # Update the model. 1447 self._update_model()

1448 1449

1450 - def _select_model(self, model=None):

1451 """Select the N-state model type. 1452 1453 @param model: The N-state model type. Can be one of '2-domain', 'population', or 'fixed'. 1454 @type model: str 1455 """ 1456 1457 # Test if the current data pipe exists. 1458 pipes.test() 1459 1460 # Test if the model name exists. 1461 if not model in ['2-domain', 'population', 'fixed']: 1462 raise RelaxError("The model name " + repr(model) + " is invalid.") 1463 1464 # Test if the model is setup. 1465 if hasattr(cdp, 'model'): 1466 warn(RelaxWarning("The N-state model has already been set up. Switching from model '%s' to '%s'." % (cdp.model, model))) 1467 1468 # Set the model 1469 cdp.model = model 1470 1471 # Initialise the list of model parameters. 1472 cdp.params = [] 1473 1474 # Update the model. 1475 self._update_model()

1476 1477

1478 - def _target_fn_setup(self, sim_index=None, scaling=True):

1479 """Initialise the target function for optimisation or direct calculation. 1480 1481 @param sim_index: The index of the simulation to optimise. This should be None if normal optimisation is desired. 1482 @type sim_index: None or int 1483 @param scaling: If True, diagonal scaling is enabled during optimisation to allow the problem to be better conditioned. 1484 @type scaling: bool 1485 """ 1486 1487 # Test if the N-state model has been set up. 1488 if not hasattr(cdp, 'model'): 1489 raise RelaxNoModelError('N-state') 1490 1491 # '2-domain' model setup tests. 1492 if cdp.model == '2-domain': 1493 # The number of states. 1494 if not hasattr(cdp, 'N'): 1495 raise RelaxError("The number of states has not been set.") 1496 1497 # The reference domain. 1498 if not hasattr(cdp, 'ref_domain'): 1499 raise RelaxError("The reference domain has not been set.") 1500 1501 # Update the model parameters if necessary. 1502 self._update_model() 1503 1504 # Create the initial parameter vector. 1505 param_vector = self._assemble_param_vector(sim_index=sim_index) 1506 1507 # Determine if alignment tensors or RDCs are to be used. 1508 data_types = self._base_data_types() 1509 1510 # The probabilities. 1511 probs = None 1512 if hasattr(cdp, 'probs'): 1513 probs = cdp.probs 1514 1515 # Diagonal scaling. 1516 scaling_matrix = None 1517 if len(param_vector): 1518 scaling_matrix = self._assemble_scaling_matrix(data_types=data_types, scaling=scaling) 1519 param_vector = dot(inv(scaling_matrix), param_vector) 1520 1521 # Get the data structures for optimisation using the tensors as base data sets. 1522 full_tensors, red_tensor_elem, red_tensor_err, full_in_ref_frame = None, None, None, None 1523 if 'tensor' in data_types: 1524 full_tensors, red_tensor_elem, red_tensor_err, full_in_ref_frame = self._minimise_setup_tensors(sim_index=sim_index) 1525 1526 # Get the data structures for optimisation using PCSs as base data sets. 1527 pcs, pcs_err, pcs_weight, temp, frq = None, None, None, None, None 1528 if 'pcs' in data_types: 1529 pcs, pcs_err, pcs_weight, temp, frq = self._minimise_setup_pcs(sim_index=sim_index) 1530 1531 # Get the data structures for optimisation using RDCs as base data sets. 1532 rdcs, rdc_err, rdc_weight, xh_vect, rdc_dj = None, None, None, None, None 1533 if 'rdc' in data_types: 1534 rdcs, rdc_err, rdc_weight, xh_vect, rdc_dj = self._minimise_setup_rdcs(sim_index=sim_index) 1535 1536 # Get the fixed tensors. 1537 fixed_tensors = None 1538 if 'rdc' in data_types or 'pcs' in data_types: 1539 full_tensors = self._minimise_setup_fixed_tensors() 1540 1541 # The flag list. 1542 fixed_tensors = [] 1543 for i in range(len(cdp.align_tensors)): 1544 if cdp.align_tensors[i].fixed: 1545 fixed_tensors.append(True) 1546 else: 1547 fixed_tensors.append(False) 1548 1549 # Get the atomic_positions. 1550 atomic_pos, paramag_centre, centre_fixed = None, None, True 1551 if 'pcs' in data_types or 'pre' in data_types: 1552 atomic_pos, paramag_centre = self._minimise_setup_atomic_pos() 1553 1554 # Optimisation of the centre. 1555 if hasattr(cdp, 'paramag_centre_fixed'): 1556 centre_fixed = cdp.paramag_centre_fixed 1557 1558 # Set up the class instance containing the target function. 1559 model = N_state_opt(model=cdp.model, N=cdp.N, init_params=param_vector, probs=probs, full_tensors=full_tensors, red_data=red_tensor_elem, red_errors=red_tensor_err, full_in_ref_frame=full_in_ref_frame, fixed_tensors=fixed_tensors, pcs=pcs, rdcs=rdcs, pcs_errors=pcs_err, rdc_errors=rdc_err, pcs_weights=pcs_weight, rdc_weights=rdc_weight, xh_vect=xh_vect, temp=temp, frq=frq, dip_const=rdc_dj, atomic_pos=atomic_pos, paramag_centre=paramag_centre, scaling_matrix=scaling_matrix, centre_fixed=centre_fixed) 1560 1561 # Return the data. 1562 return model, param_vector, data_types, scaling_matrix

1563 1564

1565 - def _tensor_loop(self, red=False):

1566 """Generator method for looping over the full or reduced tensors. 1567 1568 @keyword red: A flag which if True causes the reduced tensors to be returned, and if False 1569 the full tensors are returned. 1570 @type red: bool 1571 @return: The tensor index and the tensor. 1572 @rtype: (int, AlignTensorData instance) 1573 """ 1574 1575 # Number of tensor pairs. 1576 n = len(cdp.align_tensors.reduction) 1577 1578 # Alias. 1579 data = cdp.align_tensors 1580 list = data.reduction 1581 1582 # Full or reduced index. 1583 if red: 1584 index = 1 1585 else: 1586 index = 0 1587 1588 # Loop over the reduction list. 1589 for i in range(n): 1590 yield i, data[list[i][index]]

1591 1592

1593 - def _update_model(self):

1594 """Update the model parameters as necessary.""" 1595 1596 # Initialise the list of model parameters. 1597 if not hasattr(cdp, 'params'): 1598 cdp.params = [] 1599 1600 # Determine the number of states (loaded as structural models), if not already set. 1601 if not hasattr(cdp, 'N'): 1602 # Set the number. 1603 if hasattr(cdp, 'structure'): 1604 cdp.N = cdp.structure.num_models() 1605 1606 # Otherwise return as the rest cannot be updated without N. 1607 else: 1608 return 1609 1610 # Set up the parameter arrays. 1611 if not cdp.params: 1612 # Add the probability or population weight parameters. 1613 if cdp.model in ['2-domain', 'population']: 1614 for i in xrange(cdp.N-1): 1615 cdp.params.append('p' + repr(i)) 1616 1617 # Add the Euler angle parameters. 1618 if cdp.model == '2-domain': 1619 for i in xrange(cdp.N): 1620 cdp.params.append('alpha' + repr(i)) 1621 cdp.params.append('beta' + repr(i)) 1622 cdp.params.append('gamma' + repr(i)) 1623 1624 # Initialise the probability and Euler angle arrays. 1625 if cdp.model in ['2-domain', 'population']: 1626 if not hasattr(cdp, 'probs'): 1627 cdp.probs = [None] * cdp.N 1628 if cdp.model == '2-domain': 1629 if not hasattr(cdp, 'alpha'): 1630 cdp.alpha = [None] * cdp.N 1631 if not hasattr(cdp, 'beta'): 1632 cdp.beta = [None] * cdp.N 1633 if not hasattr(cdp, 'gamma'): 1634 cdp.gamma = [None] * cdp.N 1635 1636 # Determine the data type. 1637 data_types = self._base_data_types() 1638 1639 # Set up tensors for each alignment. 1640 if hasattr(cdp, 'align_ids'): 1641 for id in cdp.align_ids: 1642 # No tensors initialised. 1643 if not hasattr(cdp, 'align_tensors'): 1644 generic_fns.align_tensor.init(tensor=id, params=[0.0, 0.0, 0.0, 0.0, 0.0]) 1645 1646 # Find if the tensor corresponding to the id exists. 1647 exists = False 1648 for tensor in cdp.align_tensors: 1649 if id == tensor.name: 1650 exists = True 1651 1652 # Initialise the tensor. 1653 if not exists: 1654 generic_fns.align_tensor.init(tensor=id, params=[0.0, 0.0, 0.0, 0.0, 0.0])

1655 1656

1657 - def base_data_loop(self):

1658 """Loop over the base data of the spins - RDCs, PCSs, and NOESY data. 1659 1660 This loop iterates for each data point (RDC, PCS, NOESY) for each spin, returning the identification information. 1661 1662 @return: A list of the spin ID string, the data type ('rdc', 'pcs', 'noesy'), and the alignment ID if required. 1663 @rtype: list of str 1664 """ 1665 1666 # Loop over the spins. 1667 for spin, spin_id in spin_loop(return_id=True): 1668 # Re-initialise the data structure. 1669 base_ids = [spin_id, None, None] 1670 1671 # Skip deselected spins. 1672 if not spin.select: 1673 continue 1674 1675 # RDC data. 1676 if hasattr(spin, 'rdc'): 1677 base_ids[1] = 'rdc' 1678 1679 # Loop over the alignment IDs. 1680 for id in cdp.rdc_ids: 1681 # Add the ID. 1682 base_ids[2] = id 1683 1684 # Yield the set. 1685 yield base_ids 1686 1687 # PCS data. 1688 if hasattr(spin, 'pcs'): 1689 base_ids[1] = 'pcs' 1690 1691 # Loop over the alignment IDs. 1692 for id in cdp.pcs_ids: 1693 # Add the ID. 1694 base_ids[2] = id 1695 1696 # Yield the set. 1697 yield base_ids 1698 1699 # NOESY data. 1700 if hasattr(spin, 'noesy'): 1701 base_ids[1] = 'noesy' 1702 1703 # Loop over the alignment IDs. 1704 for id in cdp.noesy_ids: 1705 # Add the ID. 1706 base_ids[2] = id 1707 1708 # Yield the set. 1709 yield base_ids

1710 1711

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

1713 """Calculation function. 1714 1715 Currently this function simply calculates the NOESY flat-bottom quadratic energy potential, 1716 if NOE restraints are available. 1717 1718 @keyword spin_id: The spin identification string (unused). 1719 @type spin_id: None or str 1720 @keyword verbosity: The amount of information to print. The higher the value, the greater the verbosity. 1721 @type verbosity: int 1722 @keyword sim_index: The MC simulation index (unused). 1723 @type sim_index: None 1724 """ 1725 1726 # Set up the target function for direct calculation. 1727 model, param_vector, data_types, scaling_matrix = self._target_fn_setup() 1728 1729 # Calculate the chi-squared value. 1730 if model: 1731 # Make a function call. 1732 chi2 = model.func(param_vector) 1733 1734 # Store the global chi-squared value. 1735 cdp.chi2 = chi2 1736 1737 # Store the back-calculated data. 1738 self._minimise_bc_data(model) 1739 1740 # Calculate the RDC Q-factors. 1741 if 'rdc' in data_types: 1742 rdc.q_factors() 1743 1744 # Calculate the PCS Q-factors. 1745 if 'pcs' in data_types: 1746 pcs.q_factors() 1747 1748 # NOE potential. 1749 if hasattr(cdp, 'noe_restraints'): 1750 # Init some numpy arrays. 1751 num_restraints = len(cdp.noe_restraints) 1752 dist = zeros(num_restraints, float64) 1753 pot = zeros(num_restraints, float64) 1754 lower = zeros(num_restraints, float64) 1755 upper = zeros(num_restraints, float64) 1756 1757 # Loop over the NOEs. 1758 for i in range(num_restraints): 1759 # Create arrays of the NOEs. 1760 lower[i] = cdp.noe_restraints[i][2] 1761 upper[i] = cdp.noe_restraints[i][3] 1762 1763 # Calculate the average distances, using -6 power averaging. 1764 dist[i] = self._calc_ave_dist(cdp.noe_restraints[i][0], cdp.noe_restraints[i][1], exp=-6) 1765 1766 # Calculate the quadratic potential. 1767 quad_pot(dist, pot, lower, upper) 1768 1769 # Store the distance and potential information. 1770 cdp.ave_dist = [] 1771 cdp.quad_pot = [] 1772 for i in range(num_restraints): 1773 cdp.ave_dist.append([cdp.noe_restraints[i][0], cdp.noe_restraints[i][1], dist[i]]) 1774 cdp.quad_pot.append([cdp.noe_restraints[i][0], cdp.noe_restraints[i][1], pot[i]])

1775 1776

1777 - def create_mc_data(self, data_id=None):

1778 """Create the Monte Carlo data by back-calculation. 1779 1780 @keyword data_id: The list of spin ID, data type, and alignment ID, as yielded by the base_data_loop() generator method. 1781 @type data_id: str 1782 @return: The Monte Carlo Ri data. 1783 @rtype: list of floats 1784 """ 1785 1786 # Initialise the MC data structure. 1787 mc_data = [] 1788 1789 # Get the spin container and global spin index. 1790 spin = return_spin(data_id[0]) 1791 1792 # RDC data. 1793 if data_id[1] == 'rdc' and hasattr(spin, 'rdc'): 1794 # Does back-calculated data exist? 1795 if not hasattr(spin, 'rdc_bc'): 1796 self.calculate() 1797 1798 # The data. 1799 if not hasattr(spin, 'rdc_bc') or not spin.rdc_bc.has_key(data_id[2]): 1800 data = None 1801 else: 1802 data = spin.rdc_bc[data_id[2]] 1803 1804 # Append the data. 1805 mc_data.append(data) 1806 1807 # PCS data. 1808 elif data_id[1] == 'pcs' and hasattr(spin, 'pcs'): 1809 # Does back-calculated data exist? 1810 if not hasattr(spin, 'pcs_bc'): 1811 self.calculate() 1812 1813 # The data. 1814 if not hasattr(spin, 'pcs_bc') or not spin.pcs_bc.has_key(data_id[2]): 1815 data = None 1816 else: 1817 data = spin.pcs_bc[data_id[2]] 1818 1819 # Append the data. 1820 mc_data.append(data) 1821 1822 # NOESY data. 1823 elif data_id[1] == 'noesy' and hasattr(spin, 'noesy'): 1824 # Does back-calculated data exist? 1825 if not hasattr(spin, 'noesy_bc'): 1826 self.calculate() 1827 1828 # Append the data. 1829 mc_data.append(spin.noesy_bc) 1830 1831 # Return the data. 1832 return mc_data

1833 1834 1835 default_value_doc = Desc_container("N-state model default values") 1836 _table = uf_tables.add_table(label="table: N-state default values", caption="N-state model default values.") 1837 _table.add_headings(["Data type", "Object name", "Value"]) 1838 _table.add_row(["Probabilities", "'p0', 'p1', 'p2', ..., 'pN'", "1/N"]) 1839 _table.add_row(["Euler angle alpha", "'alpha0', 'alpha1', ...", "(c+1) * pi / (N+1)"]) 1840 _table.add_row(["Euler angle beta", "'beta0', 'beta1', ...", "(c+1) * pi / (N+1)"]) 1841 _table.add_row(["Euler angle gamma", "'gamma0', 'gamma1', ...", "(c+1) * pi / (N+1)"]) 1842 default_value_doc.add_table(_table.label) 1843 default_value_doc.add_paragraph("In this table, N is the total number of states and c is the index of a given state ranging from 0 to N-1. The default probabilities are all set to be equal whereas the angles are given a range of values so that no 2 states are equal at the start of optimisation.") 1844 default_value_doc.add_paragraph("Note that setting the probability for state N will do nothing as it is equal to one minus all the other probabilities.") 1845

1846 - def default_value(self, param):

1847 """The default N-state model parameter values. 1848 1849 @param param: The N-state model parameter. 1850 @type param: str 1851 @return: The default value. 1852 @rtype: float 1853 """ 1854 1855 # Split the parameter into its base name and index. 1856 name = self.return_data_name(param) 1857 index = self._param_model_index(param) 1858 1859 # The number of states as a float. 1860 N = float(cdp.N) 1861 1862 # Probability. 1863 if name == 'probs': 1864 return 1.0 / N 1865 1866 # Euler angles. 1867 elif name == 'alpha' or name == 'beta' or name == 'gamma': 1868 return (float(index)+1) * pi / (N+1.0)

1869 1870

1871 - def grid_search(self, lower=None, upper=None, inc=None, constraints=False, verbosity=0, sim_index=None):

1872 """The grid search function. 1873 1874 @param lower: The lower bounds of the grid search which must be equal to the number of 1875 parameters in the model. 1876 @type lower: array of numbers 1877 @param upper: The upper bounds of the grid search which must be equal to the number of 1878 parameters in the model. 1879 @type upper: array of numbers 1880 @param inc: The increments for each dimension of the space for the grid search. The 1881 number of elements in the array must equal to the number of parameters 1882 in the model. 1883 @type inc: array of int 1884 @param constraints: If True, constraints are applied during the grid search (elinating parts 1885 of the grid). If False, no constraints are used. 1886 @type constraints: bool 1887 @param verbosity: A flag specifying the amount of information to print. The higher the 1888 value, the greater the verbosity. 1889 @type verbosity: int 1890 """ 1891 1892 # Test if the N-state model has been set up. 1893 if not hasattr(cdp, 'model'): 1894 raise RelaxNoModelError('N-state') 1895 1896 # The number of parameters. 1897 n = self._param_num() 1898 1899 # Make sure that the length of the parameter array is > 0. 1900 if n == 0: 1901 print("Cannot run a grid search on a model with zero parameters, skipping the grid search.") 1902 return 1903 1904 # Test the grid search options. 1905 self.test_grid_ops(lower=lower, upper=upper, inc=inc, n=n) 1906 1907 # If inc is a single int, convert it into an array of that value. 1908 if isinstance(inc, int): 1909 inc = [inc]*n 1910 1911 # Setup the default bounds. 1912 if not lower: 1913 # Init. 1914 lower = [] 1915 upper = [] 1916 1917 # Loop over the parameters. 1918 for i in range(n): 1919 # i is in the parameter array. 1920 if i < len(cdp.params): 1921 # Probabilities (default values). 1922 if search('^p', cdp.params[i]): 1923 lower.append(0.0) 1924 upper.append(1.0) 1925 1926 # Angles (default values). 1927 if search('^alpha', cdp.params[i]) or search('^gamma', cdp.params[i]): 1928 lower.append(0.0) 1929 upper.append(2*pi) 1930 elif search('^beta', cdp.params[i]): 1931 lower.append(0.0) 1932 upper.append(pi) 1933 1934 # The paramagnetic centre. 1935 elif hasattr(cdp, 'paramag_centre_fixed') and not cdp.paramag_centre_fixed and (n - i) <= 3: 1936 lower.append(-100) 1937 upper.append(100) 1938 1939 # Otherwise this must be an alignment tensor component. 1940 else: 1941 lower.append(-1e-3) 1942 upper.append(1e-3) 1943 1944 # Minimisation. 1945 self.minimise(min_algor='grid', lower=lower, upper=upper, inc=inc, constraints=constraints, verbosity=verbosity, sim_index=sim_index)

1946 1947

1948 - def is_spin_param(self, name):

1949 """Determine whether the given parameter is spin specific. 1950 1951 @param name: The name of the parameter. 1952 @type name: str 1953 @return: False 1954 @rtype: bool 1955 """ 1956 1957 # Spin specific parameters. 1958 if name in ['r', 'heteronuc_type', 'proton_type']: 1959 return True 1960 1961 # All other parameters are global. 1962 return False

1963 1964

1965 - def map_bounds(self, param, spin_id=None):

1966 """Create bounds for the OpenDX mapping function. 1967 1968 @param param: The name of the parameter to return the lower and upper bounds of. 1969 @type param: str 1970 @param spin_id: The spin identification string (unused). 1971 @type spin_id: None 1972 @return: The upper and lower bounds of the parameter. 1973 @rtype: list of float 1974 """ 1975 1976 # Paramagnetic centre. 1977 if search('^paramag_[xyz]$', param): 1978 return [-100.0, 100.0]

1979 1980

1981 - def minimise(self, min_algor=None, min_options=None, func_tol=None, grad_tol=None, max_iterations=None, constraints=False, scaling=True, verbosity=0, sim_index=None, lower=None, upper=None, inc=None):

1982 """Minimisation function. 1983 1984 @param min_algor: The minimisation algorithm to use. 1985 @type min_algor: str 1986 @param min_options: An array of options to be used by the minimisation algorithm. 1987 @type min_options: array of str 1988 @param func_tol: The function tolerance which, when reached, terminates optimisation. Setting this to None turns of the check. 1989 @type func_tol: None or float 1990 @param grad_tol: The gradient tolerance which, when reached, terminates optimisation. Setting this to None turns of the check. 1991 @type grad_tol: None or float 1992 @param max_iterations: The maximum number of iterations for the algorithm. 1993 @type max_iterations: int 1994 @param constraints: If True, constraints are used during optimisation. 1995 @type constraints: bool 1996 @param scaling: If True, diagonal scaling is enabled during optimisation to allow the problem to be better conditioned. 1997 @type scaling: bool 1998 @param verbosity: A flag specifying the amount of information to print. The higher the value, the greater the verbosity. 1999 @type verbosity: int 2000 @param sim_index: The index of the simulation to optimise. This should be None if normal optimisation is desired. 2001 @type sim_index: None or int 2002 @keyword lower: The lower bounds of the grid search which must be equal to the number of parameters in the model. This optional argument is only used when doing a grid search. 2003 @type lower: array of numbers 2004 @keyword upper: The upper bounds of the grid search which must be equal to the number of parameters in the model. This optional argument is only used when doing a grid search. 2005 @type upper: array of numbers 2006 @keyword inc: The increments for each dimension of the space for the grid search. The number of elements in the array must equal to the number of parameters in the model. This argument is only used when doing a grid search. 2007 @type inc: array of int 2008 """ 2009 2010 # Set up the target function for direct calculation. 2011 model, param_vector, data_types, scaling_matrix = self._target_fn_setup(sim_index=sim_index, scaling=scaling) 2012 2013 # Nothing to do! 2014 if not len(param_vector): 2015 warn(RelaxWarning("The model has no parameters, minimisation cannot be performed.")) 2016 return 2017 2018 # Right, constraints cannot be used for the 'fixed' model. 2019 if constraints and cdp.model == 'fixed': 2020 warn(RelaxWarning("Turning constraints off. These cannot be used for the 'fixed' model.")) 2021 constraints = False 2022 2023 # Pop out the Method of Multipliers algorithm. 2024 if min_algor == 'Method of Multipliers': 2025 min_algor = min_options[0] 2026 min_options = min_options[1:] 2027 2028 # And constraints absolutely must be used for the 'population' model. 2029 if not constraints and cdp.model == 'population': 2030 warn(RelaxWarning("Turning constraints on. These absolutely must be used for the 'population' model.")) 2031 constraints = True 2032 2033 # Add the Method of Multipliers algorithm. 2034 min_options = (min_algor,) + min_options 2035 min_algor = 'Method of Multipliers' 2036 2037 # Linear constraints. 2038 if constraints: 2039 A, b = self._linear_constraints(data_types=data_types, scaling_matrix=scaling_matrix) 2040 else: 2041 A, b = None, None 2042 2043 # Grid search. 2044 if search('^[Gg]rid', min_algor): 2045 # Scaling. 2046 if scaling: 2047 for i in xrange(len(param_vector)): 2048 lower[i] = lower[i] / scaling_matrix[i, i] 2049 upper[i] = upper[i] / scaling_matrix[i, i] 2050 2051 # The search. 2052 results = grid(func=model.func, args=(), num_incs=inc, lower=lower, upper=upper, A=A, b=b, verbosity=verbosity) 2053 2054 # Unpack the results. 2055 param_vector, func, iter_count, warning = results 2056 f_count = iter_count 2057 g_count = 0.0 2058 h_count = 0.0 2059 2060 # Minimisation. 2061 else: 2062 results = generic_minimise(func=model.func, dfunc=model.dfunc, d2func=model.d2func, args=(), x0=param_vector, min_algor=min_algor, min_options=min_options, func_tol=func_tol, grad_tol=grad_tol, maxiter=max_iterations, A=A, b=b, full_output=1, print_flag=verbosity) 2063 2064 # Unpack the results. 2065 if results == None: 2066 return 2067 param_vector, func, iter_count, f_count, g_count, h_count, warning = results 2068 2069 # Catch infinite chi-squared values. 2070 if isInf(func): 2071 raise RelaxInfError('chi-squared') 2072 2073 # Catch chi-squared values of NaN. 2074 if isNaN(func): 2075 raise RelaxNaNError('chi-squared') 2076 2077 # Make a last function call to update the back-calculated RDC and PCS structures to the optimal values. 2078 chi2 = model.func(param_vector) 2079 2080 # Scaling. 2081 if scaling: 2082 param_vector = dot(scaling_matrix, param_vector) 2083 2084 # Disassemble the parameter vector. 2085 self._disassemble_param_vector(param_vector=param_vector, data_types=data_types, sim_index=sim_index) 2086 2087 # Monte Carlo minimisation statistics. 2088 if sim_index != None: 2089 # Chi-squared statistic. 2090 cdp.chi2_sim[sim_index] = func 2091 2092 # Iterations. 2093 cdp.iter_sim[sim_index] = iter_count 2094 2095 # Function evaluations. 2096 cdp.f_count_sim[sim_index] = f_count 2097 2098 # Gradient evaluations. 2099 cdp.g_count_sim[sim_index] = g_count 2100 2101 # Hessian evaluations. 2102 cdp.h_count_sim[sim_index] = h_count 2103 2104 # Warning. 2105 cdp.warning_sim[sim_index] = warning 2106 2107 # Normal statistics. 2108 else: 2109 # Chi-squared statistic. 2110 cdp.chi2 = func 2111 2112 # Iterations. 2113 cdp.iter = iter_count 2114 2115 # Function evaluations. 2116 cdp.f_count = f_count 2117 2118 # Gradient evaluations. 2119 cdp.g_count = g_count 2120 2121 # Hessian evaluations. 2122 cdp.h_count = h_count 2123 2124 # Warning. 2125 cdp.warning = warning 2126 2127 # Statistical analysis. 2128 if sim_index == None and ('rdc' in data_types or 'pcs' in data_types): 2129 # Get the final back calculated data (for the Q-factor and 2130 self._minimise_bc_data(model) 2131 2132 # Calculate the RDC Q-factors. 2133 if 'rdc' in data_types: 2134 rdc.q_factors() 2135 2136 # Calculate the PCS Q-factors. 2137 if 'pcs' in data_types: 2138 pcs.q_factors()

2139 2140

2141 - def model_statistics(self, model_info=None, spin_id=None, global_stats=None):

2142 """Return the k, n, and chi2 model statistics. 2143 2144 k - number of parameters. 2145 n - number of data points. 2146 chi2 - the chi-squared value. 2147 2148 2149 @keyword model_info: The data returned from model_loop() (unused). 2150 @type model_info: None 2151 @keyword spin_id: The spin identification string. This is ignored in the N-state model. 2152 @type spin_id: None or str 2153 @keyword global_stats: A parameter which determines if global or local statistics are returned. For the N-state model, this argument is ignored. 2154 @type global_stats: None or bool 2155 @return: The optimisation statistics, in tuple format, of the number of parameters (k), the number of data points (n), and the chi-squared value (chi2). 2156 @rtype: tuple of (int, int, float) 2157 """ 2158 2159 # Return the values. 2160 return self._param_num(), self._num_data_points(), cdp.chi2

2161 2162 2163 return_data_name_doc = Desc_container("N-state model data type string matching patterns") 2164 _table = uf_tables.add_table(label="table: N-state data type patterns", caption="N-state model data type string matching patterns.") 2165 _table.add_headings(["Data type", "Object name", "Patterns"]) 2166 _table.add_row(["Probabilities", "'probs'", "'p0', 'p1', 'p2', ..., 'pN'"]) 2167 _table.add_row(["Euler angle alpha", "'alpha'", "'alpha0', 'alpha1', ..."]) 2168 _table.add_row(["Euler angle beta", "'beta'", "'beta0', 'beta1', ..."]) 2169 _table.add_row(["Euler angle gamma", "'gamma'", "'gamma0', 'gamma1', ..."]) 2170 _table.add_row(["Bond length", "'r'", "'^r$' or '[Bb]ond[ -_][Ll]ength'"]) 2171 _table.add_row(["Heteronucleus type", "'heteronuc_type'", "'^[Hh]eteronucleus$'"]) 2172 _table.add_row(["Proton type", "'proton_type'", "'^[Pp]roton$'"]) 2173 return_data_name_doc.add_table(_table.label) 2174 return_data_name_doc.add_paragraph("The objects corresponding to the object names are lists (or arrays) with each element corrsponding to each state.") 2175

2176 - def return_data_name(self, param):

2177 """Return a unique identifying string for the N-state model parameter. 2178 2179 @param param: The N-state model parameter. 2180 @type param: str 2181 @return: The unique parameter identifying string. 2182 @rtype: str 2183 """ 2184 2185 # Probability. 2186 if search('^p[0-9]*$', param): 2187 return 'probs' 2188 2189 # Alpha Euler angle. 2190 if search('^alpha', param): 2191 return 'alpha' 2192 2193 # Beta Euler angle. 2194 if search('^beta', param): 2195 return 'beta' 2196 2197 # Gamma Euler angle. 2198 if search('^gamma', param): 2199 return 'gamma' 2200 2201 # Bond length. 2202 if search('^r$', param) or search('[Bb]ond[ -_][Ll]ength', param): 2203 return 'r' 2204 2205 # Heteronucleus type. 2206 if param == 'heteronuc_type': 2207 return 'heteronuc_type' 2208 2209 # Proton type. 2210 if param == 'proton_type': 2211 return 'proton_type' 2212 2213 # Paramagnetic centre. 2214 if search('^paramag_[xyz]$', param): 2215 return param

2216 2217

2218 - def return_error(self, data_id=None):

2219 """Create and return the spin specific Monte Carlo Ri error structure. 2220 2221 @keyword data_id: The list of spin ID, data type, and alignment ID, as yielded by the base_data_loop() generator method. 2222 @type data_id: str 2223 @return: The Monte Carlo simulation data errors. 2224 @rtype: list of floats 2225 """ 2226 2227 # Initialise the MC data structure. 2228 mc_errors = [] 2229 2230 # Get the spin container and global spin index. 2231 spin = return_spin(data_id[0]) 2232 2233 # Skip deselected spins. 2234 if not spin.select: 2235 return 2236 2237 # RDC data. 2238 if data_id[1] == 'rdc' and hasattr(spin, 'rdc'): 2239 # Do errors exist? 2240 if not hasattr(spin, 'rdc_err'): 2241 raise RelaxError("The RDC errors are missing for spin '%s'." % data_id[0]) 2242 2243 # Append the data. 2244 mc_errors.append(spin.rdc_err[data_id[2]]) 2245 2246 # PCS data. 2247 elif data_id[1] == 'pcs' and hasattr(spin, 'pcs'): 2248 # Do errors exist? 2249 if not hasattr(spin, 'pcs_err'): 2250 raise RelaxError("The PCS errors are missing for spin '%s'." % data_id[0]) 2251 2252 # Append the data. 2253 mc_errors.append(spin.pcs_err[data_id[2]]) 2254 2255 # NOESY data. 2256 elif hasattr(spin, 'noesy'): 2257 # Do errors exist? 2258 if not hasattr(spin, 'noesy_err'): 2259 raise RelaxError("The NOESY errors are missing for spin '%s'." % data_id[0]) 2260 2261 # Append the data. 2262 mc_errors.append(spin.noesy_err) 2263 2264 # Return the errors. 2265 return mc_errors

2266 2267

2268 - def return_grace_string(self, param):

2269 """Return the Grace string representation of the parameter. 2270 2271 This is used for axis labelling. 2272 2273 @param param: The specific analysis parameter. 2274 @type param: str 2275 @return: The Grace string representation of the parameter. 2276 @rtype: str 2277 """ 2278 2279 # The measured PCS. 2280 if param == 'pcs': 2281 return "Measured PCS" 2282 2283 # The back-calculated PCS. 2284 if param == 'pcs_bc': 2285 return "Back-calculated PCS" 2286 2287 # The measured RDC. 2288 if param == 'rdc': 2289 return "Measured RDC" 2290 2291 # The back-calculated RDC. 2292 if param == 'rdc_bc': 2293 return "Back-calculated RDC"

2294 2295

2296 - def return_units(self, param):

2297 """Return a string representing the parameters units. 2298 2299 @param param: The name of the parameter to return the units string for. 2300 @type param: str 2301 @return: The parameter units string. 2302 @rtype: str 2303 """ 2304 2305 # PCSs. 2306 if param == 'pcs' or param == 'pcs_bc': 2307 return 'ppm' 2308 2309 # RDCs. 2310 if param == 'rdc' or param == 'rdc_bc': 2311 return 'Hz'

2312 2313 2314 set_doc = Desc_container("N-state model set details") 2315 set_doc.add_paragraph("Setting parameters for the N-state model is a little different from the other type of analyses as each state has a set of parameters with the same names as the other states. To set the parameters for a specific state c (ranging from 0 for the first to N-1 for the last, the number c should be added to the end of the parameter name. So the Euler angle gamma of the third state is specified using the string 'gamma2'.") 2316 2317

2318 - def set_error(self, model_info, index, error):

2319 """Set the parameter errors. 2320 2321 @param model_info: The global model index originating from model_loop(). 2322 @type model_info: int 2323 @param index: The index of the parameter to set the errors for. 2324 @type index: int 2325 @param error: The error value. 2326 @type error: float 2327 """ 2328 2329 # Align parameters. 2330 names = ['Axx', 'Ayy', 'Axy', 'Axz', 'Ayz'] 2331 2332 # Alignment tensor parameters. 2333 if index < len(cdp.align_ids)*5: 2334 # The tensor and parameter index. 2335 param_index = index % 5 2336 tensor_index = (index - index % 5) / 5 2337 2338 # Set the error. 2339 tensor = return_tensor(index=tensor_index, skip_fixed=True) 2340 return setattr(tensor, names[param_index]+'_err', error)

2341 2342

2343 - def set_param_values(self, param=None, value=None, spin_id=None, force=True):

2344 """Set the N-state model parameter values. 2345 2346 @keyword param: The parameter name list. 2347 @type param: list of str 2348 @keyword value: The parameter value list. 2349 @type value: list 2350 @keyword spin_id: The spin identification string (unused). 2351 @type spin_id: None 2352 @keyword force: A flag which if True will cause current values to be overwritten. If False, a RelaxError will raised if the parameter value is already set. 2353 @type force: bool 2354 """ 2355 2356 # Checks. 2357 arg_check.is_str_list(param, 'parameter name') 2358 arg_check.is_list(value, 'parameter value') 2359 2360 # Loop over the parameters. 2361 for i in range(len(param)): 2362 # Get the object's name. 2363 obj_name = self.return_data_name(param[i]) 2364 2365 # Is the parameter is valid? 2366 if not obj_name: 2367 raise RelaxError("The parameter '%s' is not valid for this data pipe type." % param[i]) 2368 2369 # Set the indexed parameter. 2370 if obj_name in ['probs', 'alpha', 'beta', 'gamma']: 2371 # The index. 2372 index = self._param_model_index(param[i]) 2373 2374 # Set. 2375 obj = getattr(cdp, obj_name) 2376 obj[index] = value[i] 2377 2378 # The paramagnetic centre. 2379 if search('^paramag_[xyz]$', obj_name): 2380 # Init. 2381 if not hasattr(cdp, 'paramagnetic_centre'): 2382 cdp.paramagnetic_centre = zeros(3, float64) 2383 2384 # Set the coordinate. 2385 if obj_name == 'paramag_x': 2386 index = 0 2387 elif obj_name == 'paramag_y': 2388 index = 1 2389 else: 2390 index = 2 2391 2392 # Set the value in Angstrom. 2393 cdp.paramagnetic_centre[index] = value[i] 2394 2395 # Set the spin parameters. 2396 else: 2397 for spin in spin_loop(spin_id): 2398 setattr(spin, obj_name, value[i])

2399 2400

2401 - def sim_init_values(self):

2402 """Initialise the Monte Carlo parameter values.""" 2403 2404 # Get the minimisation statistic object names. 2405 min_names = self.data_names(set='min') 2406 2407 # Alignments. 2408 if hasattr(cdp, 'align_tensors'): 2409 # The parameter names. 2410 names = ['Axx', 'Ayy', 'Axy', 'Axz', 'Ayz'] 2411 2412 # Loop over the alignments, adding the alignment tensor parameters to the tensor data container. 2413 for i in xrange(len(cdp.align_tensors)): 2414 # Fixed tensor. 2415 if cdp.align_tensors[i].fixed: 2416 continue 2417 2418 # Loop over all the parameter names. 2419 for object_name in names: 2420 # Name for the simulation object. 2421 sim_object_name = object_name + '_sim' 2422 2423 # Create the simulation object. 2424 setattr(cdp.align_tensors[i], sim_object_name, []) 2425 2426 # Get the simulation object. 2427 sim_object = getattr(cdp.align_tensors[i], sim_object_name) 2428 2429 # Set the initial simulation values to the optimised tensor parameter values. 2430 for j in xrange(cdp.sim_number): 2431 sim_object.append(getattr(cdp.align_tensors[i], object_name)) 2432 2433 # Loop over all the minimisation object names. 2434 for object_name in min_names: 2435 # Name for the simulation object. 2436 sim_object_name = object_name + '_sim' 2437 2438 # Create the simulation object. 2439 setattr(cdp, sim_object_name, []) 2440 2441 # Get the simulation object. 2442 sim_object = getattr(cdp, sim_object_name) 2443 2444 # Loop over the simulations. 2445 for j in xrange(cdp.sim_number): 2446 # Append None to fill the structure. 2447 sim_object.append(None)

2448 2449

2450 - def sim_pack_data(self, data_id, sim_data):

2451 """Pack the Monte Carlo simulation data. 2452 2453 @keyword data_id: The list of spin ID, data type, and alignment ID, as yielded by the base_data_loop() generator method. 2454 @type data_id: list of str 2455 @param sim_data: The Monte Carlo simulation data. 2456 @type sim_data: list of float 2457 """ 2458 2459 # Get the spin container. 2460 spin = return_spin(data_id[0]) 2461 2462 # RDC data. 2463 if data_id[1] == 'rdc' and hasattr(spin, 'rdc'): 2464 # Initialise. 2465 if not hasattr(spin, 'rdc_sim'): 2466 spin.rdc_sim = {} 2467 2468 # Store the data structure. 2469 spin.rdc_sim[data_id[2]] = [] 2470 for i in range(cdp.sim_number): 2471 spin.rdc_sim[data_id[2]].append(sim_data[i][0]) 2472 2473 # PCS data. 2474 if data_id[1] == 'pcs' and hasattr(spin, 'pcs'): 2475 # Initialise. 2476 if not hasattr(spin, 'pcs_sim'): 2477 spin.pcs_sim = {} 2478 2479 # Store the data structure. 2480 spin.pcs_sim[data_id[2]] = [] 2481 for i in range(cdp.sim_number): 2482 spin.pcs_sim[data_id[2]].append(sim_data[i][0]) 2483 2484 # NOESY data. 2485 if data_id[1] == 'noesy' and hasattr(spin, 'noesy'): 2486 # Store the data structure. 2487 spin.noesy_sim = [] 2488 for i in range(cdp.sim_number): 2489 spin.noesy_sim[data_id[2]].append(sim_data[i][0])

2490 2491

2492 - def sim_return_param(self, model_info, index):

2493 """Return the array of simulation parameter values. 2494 2495 @param model_info: The global model index originating from model_loop(). 2496 @type model_info: int 2497 @param index: The index of the parameter to return the array of values for. 2498 @type index: int 2499 @return: The array of simulation parameter values. 2500 @rtype: list of float 2501 """ 2502 2503 # Align parameters. 2504 names = ['Axx', 'Ayy', 'Axy', 'Axz', 'Ayz'] 2505 2506 # Alignment tensor parameters. 2507 if index < num_tensors(skip_fixed=True)*5: 2508 # The tensor and parameter index. 2509 param_index = index % 5 2510 tensor_index = (index - index % 5) / 5 2511 2512 # Return the simulation parameter array. 2513 tensor = return_tensor(index=tensor_index, skip_fixed=True) 2514 return getattr(tensor, names[param_index]+'_sim')

Source Code for Module specific_fns.n_state_model