specific_fns.frame

54 """Class containing the specific methods of the Frame Order theories.""" 55

56 - def __init__(self):

57 """Initialise the class by placing API_common methods into the API.""" 58 59 # Place methods into the API. 60 self.eliminate = self._eliminate_false 61 self.overfit_deselect = self._overfit_deselect_dummy 62 self.return_conversion_factor = self._return_no_conversion_factor 63 self.return_data_name = self._return_data_name 64 self.return_units = self._return_units_global 65 self.set_param_values = self._set_param_values_global 66 67 # Set up the global parameters. 68 self.GLOBAL_PARAMS.add('ave_pos_alpha', units='rad') 69 self.GLOBAL_PARAMS.add('ave_pos_beta', units='rad') 70 self.GLOBAL_PARAMS.add('ave_pos_gamma', units='rad') 71 self.GLOBAL_PARAMS.add('eigen_alpha', units='rad') 72 self.GLOBAL_PARAMS.add('eigen_beta', units='rad') 73 self.GLOBAL_PARAMS.add('eigen_gamma', units='rad') 74 self.GLOBAL_PARAMS.add('axis_theta', units='rad') 75 self.GLOBAL_PARAMS.add('axis_phi', units='rad') 76 self.GLOBAL_PARAMS.add('cone_theta_x', units='rad') 77 self.GLOBAL_PARAMS.add('cone_theta_y', units='rad') 78 self.GLOBAL_PARAMS.add('cone_theta', units='rad') 79 self.GLOBAL_PARAMS.add('cone_s1') 80 self.GLOBAL_PARAMS.add('cone_sigma_max', units='rad') 81 82 # Set up the spin parameters. 83 self.SPIN_PARAMS.add('heteronuc_type', default='15N') 84 self.SPIN_PARAMS.add('proton_type', default='1H')

85 86

87 - def _assemble_limit_arrays(self):

88 """Assemble and return the limit vectors. 89 90 @return: The lower and upper limit vectors. 91 @rtype: numpy rank-1 array, numpy rank-1 array 92 """ 93 94 # Init. 95 lower = zeros(len(cdp.params), float64) 96 upper = 2.0*pi * ones(len(cdp.params), float64) 97 98 # Return the arrays. 99 return lower, upper

100 101

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

103 """Assemble and return the parameter vector. 104 105 @return: The parameter vector. 106 @rtype: numpy rank-1 array 107 @keyword sim_index: The Monte Carlo simulation index. 108 @type sim_index: int 109 """ 110 111 # Normal values. 112 if sim_index == None: 113 # Initialise the parameter array using the tensor rotation Euler angles (average domain position). 114 if cdp.model in ['free rotor', 'iso cone, torsionless', 'iso cone, free rotor']: 115 param_vect = [cdp.ave_pos_beta, cdp.ave_pos_gamma] 116 else: 117 param_vect = [cdp.ave_pos_alpha, cdp.ave_pos_beta, cdp.ave_pos_gamma] 118 119 # Frame order eigenframe - the full frame. 120 if cdp.model in ['iso cone', 'pseudo-ellipse', 'pseudo-ellipse, torsionless', 'pseudo-ellipse, free rotor']: 121 param_vect.append(cdp.eigen_alpha) 122 param_vect.append(cdp.eigen_beta) 123 param_vect.append(cdp.eigen_gamma) 124 125 # Frame order eigenframe - the isotropic cone axis. 126 elif cdp.model in ['free rotor', 'iso cone, torsionless', 'iso cone, free rotor', 'rotor']: 127 param_vect.append(cdp.axis_theta) 128 param_vect.append(cdp.axis_phi) 129 130 # Cone parameters - pseudo-elliptic cone parameters. 131 if cdp.model in ['pseudo-ellipse', 'pseudo-ellipse, torsionless', 'pseudo-ellipse, free rotor']: 132 param_vect.append(cdp.cone_theta_x) 133 param_vect.append(cdp.cone_theta_y) 134 135 # Cone parameters - single isotropic angle or order parameter. 136 elif cdp.model in ['iso cone', 'iso cone, torsionless']: 137 param_vect.append(cdp.cone_theta) 138 elif cdp.model in ['iso cone, free rotor']: 139 param_vect.append(cdp.cone_s1) 140 141 # Cone parameters - torsion angle. 142 if cdp.model in ['rotor', 'line', 'iso cone', 'pseudo-ellipse']: 143 param_vect.append(cdp.cone_sigma_max) 144 145 # Simulation values. 146 else: 147 # Initialise the parameter array using the tensor rotation Euler angles (average domain position). 148 if cdp.model in ['free rotor', 'iso cone, torsionless', 'iso cone, free rotor']: 149 param_vect = [cdp.ave_pos_beta_sim[sim_index], cdp.ave_pos_gamma_sim[sim_index]] 150 else: 151 param_vect = [cdp.ave_pos_alpha_sim[sim_index], cdp.ave_pos_beta_sim[sim_index], cdp.ave_pos_gamma_sim[sim_index]] 152 153 # Frame order eigenframe - the full frame. 154 if cdp.model in ['iso cone', 'pseudo-ellipse', 'pseudo-ellipse, torsionless', 'pseudo-ellipse, free rotor']: 155 param_vect.append(cdp.eigen_alpha_sim[sim_index]) 156 param_vect.append(cdp.eigen_beta_sim[sim_index]) 157 param_vect.append(cdp.eigen_gamma_sim[sim_index]) 158 159 # Frame order eigenframe - the isotropic cone axis. 160 elif cdp.model in ['free rotor', 'iso cone, torsionless', 'iso cone, free rotor', 'rotor']: 161 param_vect.append(cdp.axis_theta_sim[sim_index]) 162 param_vect.append(cdp.axis_phi_sim[sim_index]) 163 164 # Cone parameters - pseudo-elliptic cone parameters. 165 if cdp.model in ['pseudo-ellipse', 'pseudo-ellipse, torsionless', 'pseudo-ellipse, free rotor']: 166 param_vect.append(cdp.cone_theta_x_sim[sim_index]) 167 param_vect.append(cdp.cone_theta_y_sim[sim_index]) 168 169 # Cone parameters - single isotropic angle or order parameter. 170 elif cdp.model in ['iso cone', 'iso cone, torsionless']: 171 param_vect.append(cdp.cone_theta_sim[sim_index]) 172 elif cdp.model in ['iso cone, free rotor']: 173 param_vect.append(cdp.cone_s1_sim[sim_index]) 174 175 # Cone parameters - torsion angle. 176 if cdp.model in ['rotor', 'line', 'iso cone', 'pseudo-ellipse']: 177 param_vect.append(cdp.cone_sigma_max_sim[sim_index]) 178 179 # Return as a numpy array. 180 return array(param_vect, float64)

181 182

183 - def _back_calc(self):

184 """Back-calculation of the reduced alignment tensor. 185 186 @return: The reduced alignment tensors. 187 @rtype: numpy array 188 """ 189 190 # Get the parameter vector. 191 param_vector = self._assemble_param_vector() 192 193 # Get the data structures for optimisation using the tensors as base data sets. 194 full_tensors, red_tensors, red_tensor_err, full_in_ref_frame = self._minimise_setup_tensors() 195 196 # Set up the optimisation function. 197 target = frame_order.Frame_order(model=cdp.model, full_tensors=full_tensors, red_tensors=red_tensors, red_errors=red_tensor_err, full_in_ref_frame=full_in_ref_frame) 198 199 # Make a single function call. This will cause back calculation and the data will be stored in the class instance. 200 target.func(param_vector) 201 202 # Store the back-calculated tensors. 203 self._store_bc_tensors(target) 204 205 # Return the reduced tensors. 206 return target.red_tensors_bc

207 208

209 - def _cone_pdb(self, size=30.0, file=None, dir=None, inc=40, force=False):

210 """Create a PDB file containing a geometric object representing the Frame Order cone models. 211 212 @param size: The size of the geometric object in Angstroms. 213 @type size: float 214 @param inc: The number of increments for the filling of the cone objects. 215 @type inc: int 216 @param file: The name of the PDB file to create. 217 @type file: str 218 @param dir: The name of the directory to place the PDB file into. 219 @type dir: str 220 @param force: Flag which if set to True will cause any pre-existing file to be 221 overwritten. 222 @type force: bool 223 """ 224 225 # Test if the current data pipe exists. 226 pipes.test() 227 228 # The rigid model cannot be used here. 229 if cdp.model == 'rigid': 230 raise RelaxError("The 'rigid' frame order model has no cone representation.") 231 232 # Test for the necessary data structures. 233 if not hasattr(cdp, 'pivot'): 234 raise RelaxError("The pivot point for the domain motion has not been set.") 235 236 # Negative cone flag. 237 neg_cone = True 238 239 # Monte Carlo simulation flag. 240 sim = False 241 num_sim = 0 242 if hasattr(cdp, 'sim_number'): 243 sim = True 244 num_sim = cdp.sim_number 245 246 # The inversion matrix. 247 inv_mat = -eye(3) 248 249 # Create the structural object. 250 structure = Internal() 251 252 # Create model for the positive and negative images. 253 model = structure.add_model(model=1) 254 if neg_cone: 255 model_neg = structure.add_model(model=2) 256 257 # Create the molecule. 258 structure.add_molecule(name=cdp.model) 259 260 # Alias the molecules. 261 mol = model.mol[0] 262 if neg_cone: 263 mol_neg = model_neg.mol[0] 264 265 266 # The pivot point. 267 ################## 268 269 # Add the pivot point. 270 structure.add_atom(mol_name=cdp.model, pdb_record='HETATM', atom_num=1, atom_name='R', res_name='PIV', res_num=1, pos=cdp.pivot, element='C') 271 272 273 # The central axis. 274 ################### 275 276 # Print out. 277 print("\nGenerating the central axis.") 278 279 # The spherical angles. 280 if cdp.model in ['free rotor', 'iso cone, torsionless', 'iso cone, free rotor', 'rotor']: 281 theta_name = 'axis_theta' 282 phi_name = 'axis_phi' 283 else: 284 theta_name = 'eigen_beta' 285 phi_name = 'eigen_gamma' 286 287 # The axis. 288 axis = zeros(3, float64) 289 spherical_to_cartesian([1.0, getattr(cdp, theta_name), getattr(cdp, phi_name)], axis) 290 print(("Central axis: %s." % axis)) 291 292 # Rotations and inversions. 293 axis_pos = axis 294 axis_neg = dot(inv_mat, axis) 295 296 # Simulation central axis. 297 axis_sim_pos = None 298 axis_sim_neg = None 299 if sim: 300 # Init. 301 axis_sim = zeros((cdp.sim_number, 3), float64) 302 303 # Fill the structure. 304 for i in range(cdp.sim_number): 305 spherical_to_cartesian([1.0, getattr(cdp, theta_name+'_sim')[i], getattr(cdp, phi_name+'_sim')[i]], axis_sim[i]) 306 307 # Inversion. 308 axis_sim_pos = axis_sim 309 axis_sim_neg = transpose(dot(inv_mat, transpose(axis_sim_pos))) 310 311 # Generate the axis vectors. 312 print("\nGenerating the axis vectors.") 313 res_num = generate_vector_residues(mol=mol, vector=axis_pos, atom_name='z-ax', res_name_vect='ZAX', sim_vectors=axis_sim_pos, res_num=2, origin=cdp.pivot, scale=size) 314 315 # The negative. 316 if neg_cone: 317 res_num = generate_vector_residues(mol=mol_neg, vector=axis_neg, atom_name='z-ax', res_name_vect='ZAX', sim_vectors=axis_sim_neg, res_num=2, origin=cdp.pivot, scale=size) 318 319 320 # The x and y axes. 321 ################### 322 323 # Skip models missing the full eigenframe. 324 if cdp.model not in ['free rotor', 'iso cone, torsionless', 'iso cone, free rotor', 'rotor']: 325 # Print out. 326 print("\nGenerating the x and y axes.") 327 328 # The axis system. 329 axes = zeros((3, 3), float64) 330 euler_to_R_zyz(cdp.eigen_alpha, cdp.eigen_beta, cdp.eigen_gamma, axes) 331 print(("Axis system:\n%s" % axes)) 332 333 # Rotations and inversions. 334 axes_pos = axes 335 axes_neg = dot(inv_mat, axes) 336 337 # Simulation central axis. 338 axes_sim_pos = None 339 axes_sim_neg = None 340 if sim: 341 # Init. 342 axes_sim = zeros((3, cdp.sim_number, 3), float64) 343 344 # Fill the structure. 345 for i in range(cdp.sim_number): 346 euler_to_R_zyz(cdp.eigen_alpha_sim[i], cdp.eigen_beta_sim[i], cdp.eigen_gamma_sim[i], axes_sim[:, i]) 347 348 # Rotation and inversion. 349 axes_sim_pos = axes_sim 350 axes_sim_neg = dot(inv_mat, axes_sim_pos) 351 352 # Generate the axis vectors. 353 print("\nGenerating the axis vectors.") 354 label = ['x', 'y'] 355 for i in range(2): 356 # Simulation structures. 357 if sim: 358 axis_sim_pos = axes_sim_pos[:, i] 359 axis_sim_neg = axes_sim_neg[:, i] 360 361 # The vectors. 362 res_num = generate_vector_residues(mol=mol, vector=axes_pos[:, i], atom_name='%s-ax'%label[i], res_name_vect='%sAX'%upper(label[i]), sim_vectors=axis_sim_pos, res_num=res_num+1, origin=cdp.pivot, scale=size) 363 if neg_cone: 364 res_num = generate_vector_residues(mol=mol_neg, vector=axes_neg[:, i], atom_name='%s-ax'%label[i], res_name_vect='%sAX'%upper(label[i]), sim_vectors=axis_sim_neg, res_num=res_num, origin=cdp.pivot, scale=size) 365 366 367 # The cone object. 368 ################## 369 370 # Skip models missing a cone. 371 if cdp.model not in ['rotor', 'free rotor']: 372 # The rotation matrix (rotation from the z-axis to the cone axis). 373 if cdp.model not in ['iso cone, torsionless', 'iso cone, free rotor']: 374 R = axes 375 else: 376 R = zeros((3, 3), float64) 377 two_vect_to_R(array([0, 0, 1], float64), axis, R) 378 379 # Average position rotation. 380 R_pos = R 381 R_neg = dot(inv_mat, R) 382 383 # The pseudo-ellipse cone object. 384 if cdp.model in ['pseudo-ellipse', 'pseudo-ellipse, torsionless', 'pseudo-ellipse, free rotor']: 385 cone = Pseudo_elliptic(cdp.cone_theta_x, cdp.cone_theta_y) 386 387 # The isotropic cone object. 388 else: 389 cone = Iso_cone(cdp.cone_theta) 390 391 # Create the positive and negative cones. 392 create_cone_pdb(mol=mol, cone=cone, start_res=mol.res_num[-1]+1, apex=cdp.pivot, R=R_pos, inc=inc, distribution='regular') 393 394 # The negative. 395 if neg_cone: 396 create_cone_pdb(mol=mol_neg, cone=cone, start_res=mol_neg.res_num[-1]+1, apex=cdp.pivot, R=R_neg, inc=inc, distribution='regular') 397 398 399 # Create the PDB file. 400 ###################### 401 402 # Print out. 403 print("\nGenerating the PDB file.") 404 405 # Write the file. 406 pdb_file = open_write_file(file, dir, force=force) 407 structure.write_pdb(pdb_file) 408 pdb_file.close()

409 410

411 - def _domain_to_pdb(self, domain=None, pdb=None):

412 """Match domains to PDB files. 413 414 @keyword domain: The domain to associate the PDB file to. 415 @type domain: str 416 @keyword pdb: The PDB file to associate the domain to. 417 @type pdb: str 418 """ 419 420 # Check that the domain exists. 421 exists = False 422 for i in range(len(cdp.align_tensors)): 423 if hasattr(cdp.align_tensors[i], 'domain') and domain == cdp.align_tensors[i].domain: 424 exists = True 425 if not exists: 426 raise RelaxError("The domain '%s' cannot be found" % domain) 427 428 # Init if needed. 429 if not hasattr(cdp, 'domain_to_pdb'): 430 cdp.domain_to_pdb = [] 431 432 # Strip the file ending if given. 433 if search('.pdb$', pdb): 434 pdb = pdb[:-4] 435 436 # Add the data. 437 cdp.domain_to_pdb.append([domain, pdb])

438 439

440 - def _grid_row(self, incs, lower, upper, dist_type=None, end_point=True):

441 """Set up a row of the grid search for a given parameter. 442 443 @param incs: The number of increments. 444 @type incs: int 445 @param lower: The lower bounds. 446 @type lower: float 447 @param upper: The upper bounds. 448 @type upper: float 449 @keyword dist_type: The spacing or distribution type between grid nodes. If None, then a linear grid row is returned. If 'acos', then an inverse cos distribution of points is returned (e.g. for uniform sampling in angular space). 450 @type dist_type: None or str 451 @keyword end_point: A flag which if False will cause the end point to be removed. 452 @type end_point: bool 453 @return: The row of the grid. 454 @rtype: list of float 455 """ 456 457 # Init. 458 row = [] 459 460 # Linear grid. 461 if dist_type == None: 462 # Loop over the increments. 463 for i in range(incs): 464 # The row. 465 row.append(lower + i * (upper - lower) / (incs - 1.0)) 466 467 # Inverse cos distribution. 468 elif dist_type == 'acos': 469 # Generate the increment values of v from cos(upper) to cos(lower). 470 v = zeros(incs, float64) 471 val = (cos(lower) - cos(upper)) / (incs - 1.0) 472 for i in range(incs): 473 v[-i-1] = cos(upper) + float(i) * val 474 475 # Generate the distribution. 476 row = arccos(v) 477 478 # Remove the last point. 479 if not end_point: 480 row = row[:-1] 481 482 # Return the row (as a list). 483 return list(row)

484 485

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

487 """Set up the data structures for optimisation using alignment tensors as base data sets. 488 489 @keyword sim_index: The simulation index. This should be None if normal optimisation is 490 desired. 491 @type sim_index: None or int 492 @return: The assembled data structures for using alignment tensors as the base 493 data for optimisation. These include: 494 - full_tensors, the full tensors as concatenated arrays. 495 - red_tensors, the reduced tensors as concatenated arrays. 496 - red_err, the reduced tensor errors as concatenated arrays. 497 @rtype: tuple of 3 numpy nx5D, rank-1 arrays 498 """ 499 500 # Checks. 501 if not hasattr(cdp, 'ref_domain'): 502 raise RelaxError("The reference domain has not been set up.") 503 if not hasattr(cdp.align_tensors, 'reduction'): 504 raise RelaxError("The tensor reductions have not been specified.") 505 for i, tensor in self._tensor_loop(): 506 if not hasattr(tensor, 'domain'): 507 raise RelaxError("The domain that the '%s' tensor is attached to has not been set" % tensor.name) 508 509 # Initialise. 510 n = len(cdp.align_tensors.reduction) 511 full_tensors = zeros(n*5, float64) 512 red_tensors = zeros(n*5, float64) 513 red_err = ones(n*5, float64) * 1e-5 514 full_in_ref_frame = zeros(n, float64) 515 516 # Loop over the full tensors. 517 for i, tensor in self._tensor_loop(red=False): 518 # The full tensor. 519 full_tensors[5*i + 0] = tensor.Axx 520 full_tensors[5*i + 1] = tensor.Ayy 521 full_tensors[5*i + 2] = tensor.Axy 522 full_tensors[5*i + 3] = tensor.Axz 523 full_tensors[5*i + 4] = tensor.Ayz 524 525 # The full tensor corresponds to the frame of reference. 526 if cdp.ref_domain == tensor.domain: 527 full_in_ref_frame[i] = 1 528 529 # Loop over the reduced tensors. 530 for i, tensor in self._tensor_loop(red=True): 531 # The reduced tensor (simulation data). 532 if sim_index != None: 533 red_tensors[5*i + 0] = tensor.Axx_sim[sim_index] 534 red_tensors[5*i + 1] = tensor.Ayy_sim[sim_index] 535 red_tensors[5*i + 2] = tensor.Axy_sim[sim_index] 536 red_tensors[5*i + 3] = tensor.Axz_sim[sim_index] 537 red_tensors[5*i + 4] = tensor.Ayz_sim[sim_index] 538 539 # The reduced tensor. 540 else: 541 red_tensors[5*i + 0] = tensor.Axx 542 red_tensors[5*i + 1] = tensor.Ayy 543 red_tensors[5*i + 2] = tensor.Axy 544 red_tensors[5*i + 3] = tensor.Axz 545 red_tensors[5*i + 4] = tensor.Ayz 546 547 # The reduced tensor errors. 548 if hasattr(tensor, 'Axx_err'): 549 red_err[5*i + 0] = tensor.Axx_err 550 red_err[5*i + 1] = tensor.Ayy_err 551 red_err[5*i + 2] = tensor.Axy_err 552 red_err[5*i + 3] = tensor.Axz_err 553 red_err[5*i + 4] = tensor.Ayz_err 554 555 # Return the data structures. 556 return full_tensors, red_tensors, red_err, full_in_ref_frame

557 558

559 - def _pivot(self, pivot=None):

560 """Set the pivot point for the 2 body motion. 561 562 @param pivot: The pivot point of the two bodies (domains, etc.) in the structural 563 coordinate system. 564 @type pivot: list of num 565 """ 566 567 # Test if the current data pipe exists. 568 pipes.test() 569 570 # Set the pivot point. 571 cdp.pivot = pivot 572 573 # Convert to floats. 574 for i in range(3): 575 cdp.pivot[i] = float(cdp.pivot[i])

576 577

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

579 """Set the reference domain for the frame order, multi-domain models. 580 581 @param ref: The reference domain. 582 @type ref: str 583 """ 584 585 # Test if the current data pipe exists. 586 pipes.test() 587 588 # Test if the reference domain exists. 589 exists = False 590 for tensor_cont in cdp.align_tensors: 591 if hasattr(tensor_cont, 'domain') and tensor_cont.domain == ref: 592 exists = True 593 if not exists: 594 raise RelaxError("The reference domain cannot be found within any of the loaded tensors.") 595 596 # Set the reference domain. 597 cdp.ref_domain = ref 598 599 # Update the model. 600 if hasattr(cdp, 'model'): 601 self._update_model()

602 603

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

605 """Select the Frame Order model. 606 607 @param model: The Frame Order model. This can be one of 'pseudo-ellipse', 'pseudo-ellipse, torsionless', 'pseudo-ellipse, free rotor', 'iso cone', 'iso cone, torsionless', 'iso cone, free rotor', 'line', 'line, torsionless', 'line, free rotor', 'rotor', 'rigid', 'free rotor'. 608 @type model: str 609 """ 610 611 # Test if the current data pipe exists. 612 pipes.test() 613 614 # Test if the model name exists. 615 if not model in ['pseudo-ellipse', 'pseudo-ellipse, torsionless', 'pseudo-ellipse, free rotor', 'iso cone', 'iso cone, torsionless', 'iso cone, free rotor', 'line', 'line, torsionless', 'line, free rotor', 'rotor', 'rigid', 'free rotor']: 616 raise RelaxError("The model name " + repr(model) + " is invalid.") 617 618 # Set the model 619 cdp.model = model 620 621 # Initialise the list of model parameters. 622 cdp.params = [] 623 624 # Update the model. 625 self._update_model()

626 627

628 - def _store_bc_tensors(self, target_fn):

629 """Store the back-calculated reduced alignment tensors. 630 631 @param target_fn: The frame-order target function class. 632 @type target_fn: class instance 633 """ 634 635 # Loop over the reduced tensors. 636 for i, tensor in self._tensor_loop(red=True): 637 # New name. 638 name = tensor.name + ' bc' 639 640 # Initialise the new tensor. 641 align_tensor.init(tensor=name, params=(target_fn.red_tensors_bc[5*i + 0], target_fn.red_tensors_bc[5*i + 1], target_fn.red_tensors_bc[5*i + 2], target_fn.red_tensors_bc[5*i + 3], target_fn.red_tensors_bc[5*i + 4]), param_types=2)

642 643

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

645 """Generator method for looping over the full or reduced tensors. 646 647 @keyword red: A flag which if True causes the reduced tensors to be returned, and if False 648 the full tensors are returned. 649 @type red: bool 650 @return: The tensor index and the tensor. 651 @rtype: (int, AlignTensorData instance) 652 """ 653 654 # Number of tensor pairs. 655 n = len(cdp.align_tensors.reduction) 656 657 # Alias. 658 data = cdp.align_tensors 659 list = data.reduction 660 661 # Full or reduced index. 662 if red: 663 index = 1 664 else: 665 index = 0 666 667 # Loop over the reduction list. 668 for i in range(n): 669 yield i, data[list[i][index]]

670 671

672 - def _update_model(self):

673 """Update the model parameters as necessary.""" 674 675 # Initialise the list of model parameters. 676 if not hasattr(cdp, 'params'): 677 cdp.params = [] 678 679 # Build the parameter name list. 680 if not len(cdp.params): 681 # The tensor rotation, or average domain position. 682 if cdp.model not in ['free rotor', 'iso cone, torsionless', 'iso cone, free rotor']: 683 cdp.params.append('ave_pos_alpha') 684 cdp.params.append('ave_pos_beta') 685 cdp.params.append('ave_pos_gamma') 686 687 # Frame order eigenframe - the full frame. 688 if cdp.model in ['iso cone', 'pseudo-ellipse', 'pseudo-ellipse, torsionless', 'pseudo-ellipse, free rotor']: 689 cdp.params.append('eigen_alpha') 690 cdp.params.append('eigen_beta') 691 cdp.params.append('eigen_gamma') 692 693 # Frame order eigenframe - the isotropic cone axis. 694 elif cdp.model in ['free rotor', 'iso cone, torsionless', 'iso cone, free rotor', 'rotor']: 695 cdp.params.append('axis_theta') 696 cdp.params.append('axis_phi') 697 698 # Cone parameters - pseudo-elliptic cone parameters. 699 if cdp.model in ['pseudo-ellipse', 'pseudo-ellipse, torsionless', 'pseudo-ellipse, free rotor']: 700 cdp.params.append('cone_theta_x') 701 cdp.params.append('cone_theta_y') 702 703 # Cone parameters - single isotropic angle or order parameter. 704 elif cdp.model in ['iso cone', 'iso cone, torsionless']: 705 cdp.params.append('cone_theta') 706 elif cdp.model in ['iso cone, free rotor']: 707 cdp.params.append('cone_s1') 708 709 # Cone parameters - torsion angle. 710 if cdp.model in ['rotor', 'line', 'iso cone', 'pseudo-ellipse']: 711 cdp.params.append('cone_sigma_max') 712 713 # Initialise the parameters in the current data pipe. 714 for param in cdp.params: 715 if not hasattr(cdp, param): 716 setattr(cdp, param, 0.0)

717 718

719 - def _unpack_opt_results(self, results, sim_index=None):

720 """Unpack and store the Frame Order optimisation results. 721 722 @param results: The results tuple returned by the minfx generic_minimise() function. 723 @type results: tuple 724 @param sim_index: The index of the simulation to optimise. This should be None for normal 725 optimisation. 726 @type sim_index: None or int 727 """ 728 729 # Disassemble the results. 730 if len(results) == 4: # Grid search. 731 param_vector, func, iter_count, warning = results 732 f_count = iter_count 733 g_count = 0.0 734 h_count = 0.0 735 else: 736 param_vector, func, iter_count, f_count, g_count, h_count, warning = results 737 738 # Catch infinite chi-squared values. 739 if isInf(func): 740 raise RelaxInfError('chi-squared') 741 742 # Catch chi-squared values of NaN. 743 if isNaN(func): 744 raise RelaxNaNError('chi-squared') 745 746 # Unpack the parameters. 747 ave_pos_alpha, ave_pos_beta, ave_pos_gamma = None, None, None 748 eigen_alpha, eigen_beta, eigen_gamma = None, None, None 749 axis_theta, axis_phi = None, None 750 cone_theta_x, cone_theta_y = None, None 751 cone_theta = None 752 cone_s1 = None 753 cone_sigma_max = None 754 if cdp.model == 'pseudo-ellipse': 755 ave_pos_alpha, ave_pos_beta, ave_pos_gamma, eigen_alpha, eigen_beta, eigen_gamma, cone_theta_x, cone_theta_y, cone_sigma_max = param_vector 756 elif cdp.model in ['pseudo-ellipse, torsionless', 'pseudo-ellipse, free rotor']: 757 ave_pos_alpha, ave_pos_beta, ave_pos_gamma, eigen_alpha, eigen_beta, eigen_gamma, cone_theta_x, cone_theta_y = param_vector 758 elif cdp.model == 'iso cone': 759 ave_pos_alpha, ave_pos_beta, ave_pos_gamma, eigen_alpha, eigen_beta, eigen_gamma, cone_theta, cone_sigma_max = param_vector 760 elif cdp.model in ['iso cone, torsionless']: 761 ave_pos_beta, ave_pos_gamma, axis_theta, axis_phi, cone_theta = param_vector 762 ave_pos_alpha = None 763 elif cdp.model in ['iso cone, free rotor']: 764 ave_pos_beta, ave_pos_gamma, axis_theta, axis_phi, cone_s1 = param_vector 765 ave_pos_alpha = None 766 elif cdp.model == 'line': 767 ave_pos_alpha, ave_pos_beta, ave_pos_gamma, eigen_alpha, eigen_beta, eigen_gamma, cone_theta_x, cone_sigma_max = param_vector 768 elif cdp.model in ['line, torsionless', 'line, free rotor']: 769 ave_pos_alpha, ave_pos_beta, ave_pos_gamma, eigen_alpha, eigen_beta, eigen_gamma, cone_theta_x, cone_sigma_max = param_vector 770 elif cdp.model in ['rotor']: 771 ave_pos_alpha, ave_pos_beta, ave_pos_gamma, axis_theta, axis_phi, cone_sigma_max = param_vector 772 elif cdp.model in ['free rotor']: 773 ave_pos_beta, ave_pos_gamma, axis_theta, axis_phi = param_vector 774 ave_pos_alpha = None 775 elif cdp.model == 'rigid': 776 ave_pos_alpha, ave_pos_beta, ave_pos_gamma = param_vector 777 778 # Monte Carlo simulation data structures. 779 if sim_index != None: 780 # Average position parameters. 781 if ave_pos_alpha != None: 782 cdp.ave_pos_alpha_sim[sim_index] = wrap_angles(ave_pos_alpha, 0.0, 2.0*pi) 783 if ave_pos_beta != None: 784 cdp.ave_pos_beta_sim[sim_index] = wrap_angles(ave_pos_beta, 0.0, 2.0*pi) 785 if ave_pos_gamma != None: 786 cdp.ave_pos_gamma_sim[sim_index] = wrap_angles(ave_pos_gamma, 0.0, 2.0*pi) 787 788 # Eigenframe parameters. 789 if eigen_alpha != None: 790 cdp.eigen_alpha_sim[sim_index] = wrap_angles(eigen_alpha, 0.0, 2.0*pi) 791 if eigen_beta != None: 792 cdp.eigen_beta_sim[sim_index] = wrap_angles(eigen_beta, 0.0, 2.0*pi) 793 if eigen_gamma != None: 794 cdp.eigen_gamma_sim[sim_index] = wrap_angles(eigen_gamma, 0.0, 2.0*pi) 795 if axis_theta != None: 796 cdp.axis_theta_sim[sim_index] = axis_theta 797 if axis_phi != None: 798 cdp.axis_phi_sim[sim_index] = axis_phi 799 800 # Cone parameters. 801 if cone_theta != None: 802 cdp.cone_theta_sim[sim_index] = cone_theta 803 if cone_s1 != None: 804 cdp.cone_s1_sim[sim_index] = cone_s1 805 cdp.cone_theta_sim[sim_index] = order_parameters.iso_cone_S_to_theta(cone_s1) 806 if cone_theta_x != None: 807 cdp.cone_theta_x_sim[sim_index] = cone_theta_x 808 if cone_theta_y != None: 809 cdp.cone_theta_y_sim[sim_index] = cone_theta_y 810 if cone_sigma_max != None: 811 cdp.cone_sigma_max_sim[sim_index] = abs(cone_sigma_max) 812 813 # Optimisation stats. 814 cdp.chi2_sim[sim_index] = func 815 cdp.iter_sim[sim_index] = iter_count 816 cdp.f_count_sim[sim_index] = f_count 817 cdp.g_count_sim[sim_index] = g_count 818 cdp.h_count_sim[sim_index] = h_count 819 cdp.warning_sim[sim_index] = warning 820 821 # Normal data structures. 822 else: 823 # Average position parameters. 824 if ave_pos_alpha != None: 825 cdp.ave_pos_alpha = wrap_angles(ave_pos_alpha, 0.0, 2.0*pi) 826 if ave_pos_beta != None: 827 cdp.ave_pos_beta = wrap_angles(ave_pos_beta, 0.0, 2.0*pi) 828 if ave_pos_gamma != None: 829 cdp.ave_pos_gamma = wrap_angles(ave_pos_gamma, 0.0, 2.0*pi) 830 831 # Eigenframe parameters. 832 if eigen_alpha != None: 833 cdp.eigen_alpha = wrap_angles(eigen_alpha, 0.0, 2.0*pi) 834 if eigen_beta != None: 835 cdp.eigen_beta = wrap_angles(eigen_beta, 0.0, 2.0*pi) 836 if eigen_gamma != None: 837 cdp.eigen_gamma = wrap_angles(eigen_gamma, 0.0, 2.0*pi) 838 if axis_theta != None: 839 cdp.axis_theta = axis_theta 840 if axis_phi != None: 841 cdp.axis_phi = axis_phi 842 843 # Cone parameters. 844 if cone_theta != None: 845 cdp.cone_theta = cone_theta 846 if cone_s1 != None: 847 cdp.cone_s1 = cone_s1 848 cdp.cone_theta = order_parameters.iso_cone_S_to_theta(cone_s1) 849 if cone_theta_x != None: 850 cdp.cone_theta_x = cone_theta_x 851 if cone_theta_y != None: 852 cdp.cone_theta_y = cone_theta_y 853 if cone_sigma_max != None: 854 cdp.cone_sigma_max = abs(cone_sigma_max) 855 856 # Optimisation stats. 857 cdp.chi2 = func 858 cdp.iter = iter_count 859 cdp.f_count = f_count 860 cdp.g_count = g_count 861 cdp.h_count = h_count 862 cdp.warning = warning

863 864

865 - def base_data_loop(self):

866 """Generator method for looping nothing. 867 868 The loop essentially consists of a single element. 869 870 @return: Nothing. 871 @rtype: None 872 """ 873 874 yield None

875 876

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

878 """Calculate the chi-squared value for the current parameter values. 879 880 @keyword spin_id: The spin identification string (unused). 881 @type spin_id: None 882 @keyword verbosity: The amount of information to print. The higher the value, the greater the verbosity. 883 @type verbosity: int 884 @keyword sim_index: The optional MC simulation index (unused). 885 @type sim_index: None or int 886 """ 887 888 # Assemble the parameter vector. 889 param_vector = self._assemble_param_vector() 890 891 # Get the data structures for optimisation using the tensors as base data sets. 892 full_tensors, red_tensors, red_tensor_err, full_in_ref_frame = self._minimise_setup_tensors() 893 894 # Set up the optimisation function. 895 target = frame_order.Frame_order(model=cdp.model, full_tensors=full_tensors, red_tensors=red_tensors, red_errors=red_tensor_err, full_in_ref_frame=full_in_ref_frame) 896 897 # Make a single function call. This will cause back calculation and the data will be stored in the class instance. 898 chi2 = target.func(param_vector) 899 900 # Set the chi2. 901 cdp.chi2 = chi2 902 903 # Store the back-calculated tensors. 904 self._store_bc_tensors(target)

905 906

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

908 """Create the Monte Carlo data by back calculating the reduced tensor data. 909 910 @keyword data_id: Unused. 911 @type data_id: None 912 @return: The Monte Carlo simulation data. 913 @rtype: list of floats 914 """ 915 916 # Back calculate the tensors. 917 red_tensors_bc = self._back_calc() 918 919 # Return the data. 920 return red_tensors_bc

921 922

923 - def data_names(self, set='all', error_names=False, sim_names=False):

924 """Function for returning a list of names of data structures. 925 926 Description 927 =========== 928 929 The names are as follows: 930 931 - 'params', an array of the parameter names associated with the model. 932 - 'chi2', chi-squared value. 933 - 'iter', iterations. 934 - 'f_count', function count. 935 - 'g_count', gradient count. 936 - 'h_count', hessian count. 937 - 'warning', minimisation warning. 938 939 The isotropic cone specific model parameters are: 940 941 - 'axis_theta', the cone axis polar angle (for the isotropic cone model). 942 - 'axis_phi', the cone axis azimuthal angle (for the isotropic cone model). 943 - 'cone_s1', the isotropic cone order parameter. 944 945 946 @keyword set: The set of object names to return. This can be set to 'all' for all 947 names, to 'generic' for generic object names, 'params' for the 948 Frame Order parameter names, or to 'min' for minimisation specific 949 object names. 950 @type set: str 951 @keyword error_names: A flag which if True will add the error object names as well. 952 @type error_names: bool 953 @keyword sim_names: A flag which if True will add the Monte Carlo simulation object 954 names as well. 955 @type sim_names: bool 956 @return: The list of object names. 957 @rtype: list of str 958 """ 959 960 # Initialise. 961 names = [] 962 963 # Generic. 964 if set == 'all' or set == 'generic': 965 names.append('params') 966 967 # The parameter suffix. 968 if error_names: 969 suffix = '_err' 970 elif sim_names: 971 suffix = '_sim' 972 else: 973 suffix = '' 974 975 # Parameters. 976 if (set == 'all' or set == 'params') and hasattr(cdp, 'model'): 977 # Initialise the parameter array using the tensor rotation Euler angles (average domain position). 978 if cdp.model not in ['free rotor', 'iso cone, torsionless', 'iso cone, free rotor']: 979 names.append('ave_pos_alpha%s' % suffix) 980 names.append('ave_pos_beta%s' % suffix) 981 names.append('ave_pos_gamma%s' % suffix) 982 983 # Frame order eigenframe - the full frame. 984 if cdp.model in ['iso cone', 'pseudo-ellipse', 'pseudo-ellipse, torsionless', 'pseudo-ellipse, free rotor']: 985 names.append('eigen_alpha%s' % suffix) 986 names.append('eigen_beta%s' % suffix) 987 names.append('eigen_gamma%s' % suffix) 988 989 # Frame order eigenframe - the isotropic cone axis. 990 elif cdp.model in ['free rotor', 'iso cone, torsionless', 'iso cone, free rotor', 'rotor']: 991 names.append('axis_theta%s' % suffix) 992 names.append('axis_phi%s' % suffix) 993 994 # Cone parameters - pseudo-elliptic cone parameters. 995 if cdp.model in ['pseudo-ellipse', 'pseudo-ellipse, torsionless', 'pseudo-ellipse, free rotor']: 996 names.append('cone_theta_x%s' % suffix) 997 names.append('cone_theta_y%s' % suffix) 998 999 # Cone parameters - single isotropic angle or order parameter. 1000 elif cdp.model in ['iso cone', 'iso cone, torsionless']: 1001 names.append('cone_theta%s' % suffix) 1002 elif cdp.model in ['iso cone, free rotor']: 1003 names.append('cone_s1%s' % suffix) 1004 1005 # Cone parameters - torsion angle. 1006 if cdp.model in ['rotor', 'line', 'iso cone', 'pseudo-ellipse']: 1007 names.append('cone_sigma_max%s' % suffix) 1008 1009 # Minimisation statistics. 1010 if set == 'all' or set == 'min': 1011 names.append('chi2') 1012 names.append('iter') 1013 names.append('f_count') 1014 names.append('g_count') 1015 names.append('h_count') 1016 names.append('warning') 1017 1018 # Return the names. 1019 return names

1020 1021

1022 - def get_param_names(self, model_info=None):

1023 """Return a vector of parameter names. 1024 1025 @keyword model_info: The model index from model_info(). 1026 @type model_info: int 1027 @return: The vector of parameter names. 1028 @rtype: list of str 1029 """ 1030 1031 # Return the parameter object names. 1032 param_names = self.data_names(set='params') 1033 return param_names

1034 1035

1036 - def get_param_values(self, model_info=None, sim_index=None):

1037 """Return a vector of parameter values. 1038 1039 @keyword model_info: The model index from model_info(). This is zero for the global models or equal to the global spin index (which covers the molecule, residue, and spin indices). 1040 @type model_info: int 1041 @keyword sim_index: The Monte Carlo simulation index. 1042 @type sim_index: int 1043 @return: The vector of parameter values. 1044 @rtype: list of str 1045 """ 1046 1047 # Assemble the values and return it. 1048 return self._assemble_param_vector(sim_index=sim_index)

1049 1050

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

1052 """Perform a grid search. 1053 1054 @keyword lower: The lower bounds of the grid search which must be equal to the 1055 number of parameters in the model. 1056 @type lower: list of float 1057 @keyword upper: The upper bounds of the grid search which must be equal to the 1058 number of parameters in the model. 1059 @type upper: list of float 1060 @keyword inc: The increments for each dimension of the space for the grid search. 1061 The number of elements in the array must equal to the number of 1062 parameters in the model. 1063 @type inc: int or list of int 1064 @keyword constraints: If True, constraints are applied during the grid search (eliminating 1065 parts of the grid). If False, no constraints are used. 1066 @type constraints: bool 1067 @keyword verbosity: A flag specifying the amount of information to print. The higher 1068 the value, the greater the verbosity. 1069 @type verbosity: int 1070 @keyword sim_index: The Monte Carlo simulation index. 1071 @type sim_index: None or int 1072 """ 1073 1074 # Test if the Frame Order model has been set up. 1075 if not hasattr(cdp, 'model'): 1076 raise RelaxNoModelError('Frame Order') 1077 1078 # The number of parameters. 1079 n = len(cdp.params) 1080 1081 # If inc is an int, convert it into an array of that value. 1082 if isinstance(inc, int): 1083 incs = [inc]*n 1084 else: 1085 incs = inc 1086 1087 # Sanity check. 1088 if len(incs) != n: 1089 raise RelaxError("The size of the increment list %s does not match the number of parameters in %s." % (incs, cdp.params)) 1090 1091 # Initialise the grid increments structures. 1092 lower_list = lower 1093 upper_list = upper 1094 grid = [] 1095 """This structure is a list of lists. The first dimension corresponds to the model 1096 parameter. The second dimension are the grid node positions.""" 1097 1098 # Generate the grid. 1099 for i in range(n): 1100 # Fixed parameter. 1101 if incs[i] == None: 1102 grid.append(None) 1103 continue 1104 1105 # Reset. 1106 dist_type = None 1107 end_point = True 1108 1109 # Linear angle grid from 0 to one inc before 2pi. 1110 if cdp.params[i] in ['ave_pos_alpha', 'ave_pos_gamma', 'eigen_alpha', 'eigen_gamma', 'axis_phi']: 1111 lower = 0.0 1112 upper = 2*pi * (1.0 - 1.0/incs[i]) 1113 1114 # Arccos grid from 0 to pi. 1115 if cdp.params[i] in ['ave_pos_beta', 'eigen_beta', 'axis_theta']: 1116 # Change the default increment numbers. 1117 if not isinstance(inc, list): 1118 incs[i] = incs[i] / 2 + 1 1119 1120 # The distribution type and end point. 1121 dist_type = 'acos' 1122 end_point = False 1123 1124 # Set the default bounds. 1125 lower = 0.0 1126 upper = pi 1127 1128 # The isotropic cone order parameter. 1129 if cdp.params[i] == 'cone_s1': 1130 lower = -0.125 1131 upper = 1.0 1132 1133 # Linear angle grid from 0 to pi excluding the outer points. 1134 if cdp.params[i] in ['cone_theta', 'cone_theta_x', 'cone_theta_y', 'cone_sigma_max']: 1135 lower = pi * (1.0/incs[i]) 1136 upper = pi * (1.0 - 1.0/incs[i]) 1137 1138 # Over-ride the bounds. 1139 if lower_list: 1140 lower = lower_list[i] 1141 if upper_list: 1142 upper = upper_list[i] 1143 1144 # Append the grid row. 1145 row = self._grid_row(incs[i], lower, upper, dist_type=dist_type, end_point=end_point) 1146 grid.append(row) 1147 1148 # Remove an inc if the end point has been removed. 1149 if not end_point: 1150 incs[i] -= 1 1151 1152 # Total number of points. 1153 total_pts = 1 1154 for i in range(n): 1155 # Fixed parameter. 1156 if grid[i] == None: 1157 continue 1158 1159 total_pts = total_pts * len(grid[i]) 1160 1161 # Check the number. 1162 max_pts = 50e6 1163 if total_pts > max_pts: 1164 raise RelaxError("The total number of grid points '%s' exceeds the maximum of '%s'." % (total_pts, int(max_pts))) 1165 1166 # Build the points array. 1167 pts = zeros((total_pts, n), float64) 1168 indices = zeros(n, int) 1169 for i in range(total_pts): 1170 # Loop over the dimensions. 1171 for j in range(n): 1172 # Fixed parameter. 1173 if grid[j] == None: 1174 # Get the current parameter value. 1175 pts[i, j] = getattr(cdp, cdp.params[j]) 1176 1177 # Add the point coordinate. 1178 else: 1179 pts[i, j] = grid[j][indices[j]] 1180 1181 # Increment the step positions. 1182 for j in range(n): 1183 if incs[j] != None and indices[j] < incs[j]-1: 1184 indices[j] += 1 1185 break # Exit so that the other step numbers are not incremented. 1186 else: 1187 indices[j] = 0 1188 1189 # Minimisation. 1190 self.minimise(min_algor='grid', min_options=pts, constraints=constraints, verbosity=verbosity, sim_index=sim_index)

1191 1192

1193 - def is_spin_param(self, name):

1194 """State that the parameter is not spin specific. 1195 1196 @param name: The name of the parameter. 1197 @type name: str 1198 @return: False. 1199 @rtype: bool 1200 """ 1201 1202 # Not spin specific! 1203 return False

1204 1205

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

1207 """Create bounds for the OpenDX mapping function. 1208 1209 @param param: The name of the parameter to return the lower and upper bounds of. 1210 @type param: str 1211 @param spin_id: The spin identification string (unused). 1212 @type spin_id: None 1213 @return: The upper and lower bounds of the parameter. 1214 @rtype: list of float 1215 """ 1216 1217 # Average position Euler angles. 1218 if param in ['ave_pos_alpha', 'ave_pos_beta', 'ave_pos_gamma']: 1219 return [0.0, 2*pi] 1220 1221 # Axis spherical coordinate theta. 1222 if param == 'axis_theta': 1223 return [0.0, pi] 1224 1225 # Axis spherical coordinate phi. 1226 if param == 'axis_phi': 1227 return [0.0, 2*pi] 1228 1229 # Cone angle. 1230 if param == 'cone_theta': 1231 return [0.0, pi]

1232 1233

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

1235 """Minimisation function. 1236 1237 @param min_algor: The minimisation algorithm to use. 1238 @type min_algor: str 1239 @param min_options: An array of options to be used by the minimisation algorithm. 1240 @type min_options: array of str 1241 @param func_tol: The function tolerance which, when reached, terminates optimisation. 1242 Setting this to None turns of the check. 1243 @type func_tol: None or float 1244 @param grad_tol: The gradient tolerance which, when reached, terminates optimisation. 1245 Setting this to None turns of the check. 1246 @type grad_tol: None or float 1247 @param max_iterations: The maximum number of iterations for the algorithm. 1248 @type max_iterations: int 1249 @param constraints: If True, constraints are used during optimisation. 1250 @type constraints: bool 1251 @param scaling: If True, diagonal scaling is enabled during optimisation to allow 1252 the problem to be better conditioned. 1253 @type scaling: bool 1254 @param verbosity: A flag specifying the amount of information to print. The higher 1255 the value, the greater the verbosity. 1256 @type verbosity: int 1257 @param sim_index: The index of the simulation to optimise. This should be None if 1258 normal optimisation is desired. 1259 @type sim_index: None or int 1260 @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. 1261 @type lower: array of numbers 1262 @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. 1263 @type upper: array of numbers 1264 @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. 1265 @type inc: array of int 1266 """ 1267 1268 # Constraints not implemented yet. 1269 if constraints: 1270 # Turn the constraints off. 1271 constraints = False 1272 1273 # Remove the method of multipliers arg. 1274 if not search('^[Gg]rid', min_algor): 1275 min_algor = min_options[0] 1276 min_options = min_options[1:] 1277 1278 # Throw a warning. 1279 warn(RelaxWarning("Constraints are as of yet not implemented - turning this option off.")) 1280 1281 # Simulated annealing constraints. 1282 lower, upper = self._assemble_limit_arrays() 1283 1284 # Assemble the parameter vector. 1285 param_vector = self._assemble_param_vector() 1286 1287 # Get the data structures for optimisation using the tensors as base data sets. 1288 full_tensors, red_tensors, red_tensor_err, full_in_ref_frame = self._minimise_setup_tensors(sim_index) 1289 1290 # Set up the optimisation function. 1291 target = frame_order.Frame_order(model=cdp.model, full_tensors=full_tensors, red_tensors=red_tensors, red_errors=red_tensor_err, full_in_ref_frame=full_in_ref_frame) 1292 1293 # Grid search. 1294 if search('^[Gg]rid', min_algor): 1295 results = grid_point_array(func=target.func, args=(), points=min_options, verbosity=verbosity) 1296 1297 # Minimisation. 1298 else: 1299 results = generic_minimise(func=target.func, args=(), x0=param_vector, min_algor=min_algor, min_options=min_options, func_tol=func_tol, grad_tol=grad_tol, maxiter=max_iterations, l=lower, u=upper, full_output=True, print_flag=verbosity) 1300 1301 # Unpack the results. 1302 self._unpack_opt_results(results, sim_index) 1303 1304 # Store the back-calculated tensors. 1305 self._store_bc_tensors(target)

1306 1307 1308

1309 - def model_loop(self):

1310 """Dummy generator method. 1311 1312 In this case only a single model per spin system is assumed. Hence the yielded data is the 1313 spin container object. 1314 1315 1316 @return: Information about the model which for this analysis is the spin container. 1317 @rtype: SpinContainer instance 1318 """ 1319 1320 # Don't return anything, just loop once. 1321 yield None

1322 1323

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

1325 """Return the k, n, and chi2 model statistics. 1326 1327 k - number of parameters. 1328 n - number of data points. 1329 chi2 - the chi-squared value. 1330 1331 1332 @keyword model_info: Unused. 1333 @type model_info: None 1334 @keyword spin_id: The spin identification string (unused). 1335 @type spin_id: None 1336 @keyword global_stats: Unused. 1337 @type global_stats: None 1338 @return: The optimisation statistics, in tuple format, of the number of 1339 parameters (k), the number of data points (n), and the chi-squared 1340 value (chi2). 1341 @rtype: tuple of (int, int, float) 1342 """ 1343 1344 # Count the number of parameters. 1345 param_vector = self._assemble_param_vector() 1346 k = len(param_vector) 1347 1348 # The number of data points. 1349 n = len(cdp.align_tensors.reduction) 1350 1351 # The chi2 value. 1352 if not hasattr(cdp, 'chi2'): 1353 raise RelaxError("Statistics are not available, most likely because the model has not been optimised.") 1354 chi2 = cdp.chi2 1355 1356 # Return the data. 1357 return k, n, chi2

1358 1359

1360 - def return_error(self, data_id):

1361 """Return the alignment tensor error structure. 1362 1363 @param data_id: The information yielded by the base_data_loop() generator method. 1364 @type data_id: None 1365 @return: The array of tensor error values. 1366 @rtype: list of float 1367 """ 1368 1369 # Get the tensor data structures. 1370 full_tensors, red_tensors, red_tensor_err, full_in_ref_frame = self._minimise_setup_tensors() 1371 1372 # Return the errors. 1373 return red_tensor_err

1374 1375

1376 - def return_units(self, param):

1377 """Return a string representing the parameters units. 1378 1379 @param param: The name of the parameter to return the units string for. 1380 @type param: str 1381 @return: The parameter units string. 1382 @rtype: str 1383 """ 1384 1385 # Average position Euler angles. 1386 if param in ['ave_pos_alpha', 'ave_pos_beta', 'ave_pos_gamma']: 1387 return 'rad' 1388 1389 # Eigenframe angles. 1390 if param in ['eigen_alpha', 'eigen_beta', 'eigen_gamma', 'axis_theta', 'axis_phi']: 1391 return 'rad' 1392 1393 # Cone angles. 1394 if param in ['cone_theta_x', 'cone_theta_y', 'cone_sigma_max', 'cone_theta']: 1395 return 'rad'

1396 1397

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

1399 """Set the parameter errors. 1400 1401 @param model_info: The model information originating from model_loop() (unused). 1402 @type model_info: None 1403 @param index: The index of the parameter to set the errors for. 1404 @type index: int 1405 @param error: The error value. 1406 @type error: float 1407 """ 1408 1409 # Parameter increment counter. 1410 inc = 0 1411 1412 # Loop over the residue specific parameters. 1413 for param in self.data_names(set='params'): 1414 # Return the parameter array. 1415 if index == inc: 1416 setattr(cdp, param + "_err", error) 1417 1418 # Increment. 1419 inc = inc + 1 1420 1421 # Add some additional parameters. 1422 if cdp.model == 'iso cone, free rotor' and inc == index: 1423 setattr(cdp, 'cone_theta_err', error)

1424 1425 1426

1427 - def set_selected_sim(self, model_info, select_sim):

1428 """Set the simulation selection flag for the spin. 1429 1430 @param model_info: The model information originating from model_loop() (unused). 1431 @type model_info: None 1432 @param select_sim: The selection flag for the simulations. 1433 @type select_sim: bool 1434 """ 1435 1436 # Set the array. 1437 cdp.select_sim = deepcopy(select_sim)

1438 1439

1440 - def sim_init_values(self):

1441 """Initialise the Monte Carlo parameter values.""" 1442 1443 # Get the parameter object names. 1444 param_names = self.data_names(set='params') 1445 1446 # Add some additional parameters. 1447 if cdp.model == 'iso cone, free rotor': 1448 param_names.append('cone_theta') 1449 1450 # Get the minimisation statistic object names. 1451 min_names = self.data_names(set='min') 1452 1453 1454 # Test if Monte Carlo parameter values have already been set. 1455 ############################################################# 1456 1457 # Loop over all the parameter names. 1458 for object_name in param_names: 1459 # Name for the simulation object. 1460 sim_object_name = object_name + '_sim' 1461 1462 # Test if the simulation object already exists. 1463 if hasattr(cdp, sim_object_name): 1464 raise RelaxError("Monte Carlo parameter values have already been set.") 1465 1466 1467 # Set the Monte Carlo parameter values. 1468 ####################################### 1469 1470 # Loop over all the data names. 1471 for object_name in param_names: 1472 # Name for the simulation object. 1473 sim_object_name = object_name + '_sim' 1474 1475 # Create the simulation object. 1476 setattr(cdp, sim_object_name, []) 1477 1478 # Get the simulation object. 1479 sim_object = getattr(cdp, sim_object_name) 1480 1481 # Loop over the simulations. 1482 for j in xrange(cdp.sim_number): 1483 # Copy and append the data. 1484 sim_object.append(deepcopy(getattr(cdp, object_name))) 1485 1486 # Loop over all the minimisation object names. 1487 for object_name in min_names: 1488 # Name for the simulation object. 1489 sim_object_name = object_name + '_sim' 1490 1491 # Create the simulation object. 1492 setattr(cdp, sim_object_name, []) 1493 1494 # Get the simulation object. 1495 sim_object = getattr(cdp, sim_object_name) 1496 1497 # Loop over the simulations. 1498 for j in xrange(cdp.sim_number): 1499 # Copy and append the data. 1500 sim_object.append(deepcopy(getattr(cdp, object_name)))

1501 1502

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

1504 """Pack the Monte Carlo simulation data. 1505 1506 @param data_id: The spin identification string, as yielded by the base_data_loop() generator method. 1507 @type data_id: None 1508 @param sim_data: The Monte Carlo simulation data. 1509 @type sim_data: list of float 1510 """ 1511 1512 # Transpose the data. 1513 sim_data = transpose(sim_data) 1514 1515 # Loop over the reduced tensors. 1516 for i, tensor in self._tensor_loop(red=True): 1517 # Set the reduced tensor simulation data. 1518 tensor.Axx_sim = sim_data[5*i + 0] 1519 tensor.Ayy_sim = sim_data[5*i + 1] 1520 tensor.Axy_sim = sim_data[5*i + 2] 1521 tensor.Axz_sim = sim_data[5*i + 3] 1522 tensor.Ayz_sim = sim_data[5*i + 4]

1523 1524

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

1526 """Return the array of simulation parameter values. 1527 1528 @param model_info: The model information originating from model_loop() (unused). 1529 @type model_info: None 1530 @param index: The index of the parameter to return the array of values for. 1531 @type index: int 1532 """ 1533 1534 # Parameter increment counter. 1535 inc = 0 1536 1537 # Get the parameter object names. 1538 param_names = self.data_names(set='params') 1539 1540 # Loop over the parameters. 1541 for param in param_names: 1542 # Return the parameter array. 1543 if index == inc: 1544 return getattr(cdp, param + "_sim") 1545 1546 # Increment. 1547 inc = inc + 1 1548 1549 # Add some additional parameters. 1550 if cdp.model == 'iso cone, free rotor' and inc == index: 1551 return getattr(cdp, 'cone_theta_sim')

1552 1553

1554 - def sim_return_selected(self, model_info):

1555 """Return the array of selected simulation flags for the spin. 1556 1557 @param model_info: The model information originating from model_loop() (unused). 1558 @type model_info: None 1559 @return: The array of selected simulation flags. 1560 @rtype: list of int 1561 """ 1562 1563 # Return the array. 1564 return cdp.select_sim

Source Code for Module specific_fns.frame_order