Package specific_fns :: Module frame_order
[hide private]
[frames] | no frames]

Source Code for Module specific_fns.frame_order

   1  ############################################################################### 
   2  #                                                                             # 
   3  # Copyright (C) 2009-2013 Edward d'Auvergne                                   # 
   4  #                                                                             # 
   5  # This file is part of the program relax (http://www.nmr-relax.com).          # 
   6  #                                                                             # 
   7  # This program is free software: you can redistribute it and/or modify        # 
   8  # it under the terms of the GNU General Public License as published by        # 
   9  # the Free Software Foundation, either version 3 of the License, or           # 
  10  # (at your option) any later version.                                         # 
  11  #                                                                             # 
  12  # This program is distributed in the hope that it will be useful,             # 
  13  # but WITHOUT ANY WARRANTY; without even the implied warranty of              # 
  14  # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the               # 
  15  # GNU General Public License for more details.                                # 
  16  #                                                                             # 
  17  # You should have received a copy of the GNU General Public License           # 
  18  # along with this program.  If not, see <http://www.gnu.org/licenses/>.       # 
  19  #                                                                             # 
  20  ############################################################################### 
  21   
  22  # Module docstring. 
  23  """Analysis specific code for the Frame Order theory of domain dynamics.""" 
  24   
  25  # Python module imports. 
  26  from copy import deepcopy 
  27  from math import cos, pi 
  28  from minfx.generic import generic_minimise 
  29  from minfx.grid import grid_point_array 
  30  from numpy import arccos, array, dot, eye, float64, ones, transpose, zeros 
  31  from re import search 
  32  from warnings import warn 
  33   
  34  # relax module imports. 
  35  from float import isNaN, isInf 
  36  from generic_fns import align_tensor, pipes 
  37  from generic_fns.angles import wrap_angles 
  38  from generic_fns.structure.cones import Iso_cone, Pseudo_elliptic 
  39  from generic_fns.structure.geometric import create_cone_pdb, generate_vector_dist, generate_vector_residues 
  40  from generic_fns.structure.internal import Internal 
  41  from maths_fns import frame_order, order_parameters 
  42  from maths_fns.coord_transform import spherical_to_cartesian 
  43  from maths_fns.rotation_matrix import euler_to_R_zyz, two_vect_to_R 
  44  from relax_errors import RelaxError, RelaxInfError, RelaxModelError, RelaxNaNError, RelaxNoModelError 
  45  from relax_io import open_write_file 
  46  from relax_warnings import RelaxWarning, RelaxDeselectWarning 
  47  from specific_fns.api_base import API_base 
  48  from specific_fns.api_common import API_common 
  49   
  50   
51 -class Frame_order(API_base, API_common):
52 """Class containing the specific methods of the Frame Order theories.""" 53
54 - def __init__(self):
55 """Initialise the class by placing API_common methods into the API.""" 56 57 # Execute the base class __init__ method. 58 super(Frame_order, self).__init__() 59 60 # Place methods into the API. 61 self.eliminate = self._eliminate_false 62 self.overfit_deselect = self._overfit_deselect_dummy 63 self.return_conversion_factor = self._return_no_conversion_factor 64 self.set_param_values = self._set_param_values_global 65 66 # Set up the global parameters. 67 self.PARAMS.add('ave_pos_alpha', scope='global', units='rad', desc='The average position alpha Euler angle', py_type=float, set='params', err=True, sim=True) 68 self.PARAMS.add('ave_pos_beta', scope='global', units='rad', desc='The average position beta Euler angle', py_type=float, set='params', err=True, sim=True) 69 self.PARAMS.add('ave_pos_gamma', scope='global', units='rad', desc='The average position gamma Euler angle', py_type=float, set='params', err=True, sim=True) 70 self.PARAMS.add('eigen_alpha', scope='global', units='rad', desc='The Eigenframe alpha Euler angle', py_type=float, set='params', err=True, sim=True) 71 self.PARAMS.add('eigen_beta', scope='global', units='rad', desc='The Eigenframe beta Euler angle', py_type=float, set='params', err=True, sim=True) 72 self.PARAMS.add('eigen_gamma', scope='global', units='rad', desc='The Eigenframe gamma Euler angle', py_type=float, set='params', err=True, sim=True) 73 self.PARAMS.add('axis_theta', scope='global', units='rad', desc='The cone axis polar angle (for the isotropic cone model)', py_type=float, set='params', err=True, sim=True) 74 self.PARAMS.add('axis_phi', scope='global', units='rad', desc='The cone axis azimuthal angle (for the isotropic cone model)', py_type=float, set='params', err=True, sim=True) 75 self.PARAMS.add('cone_theta_x', scope='global', units='rad', desc='The pseudo-ellipse cone opening half-angle for the x-axis', py_type=float, set='params', err=True, sim=True) 76 self.PARAMS.add('cone_theta_y', scope='global', units='rad', desc='The pseudo-ellipse cone opening half-angle for the y-axis', py_type=float, set='params', err=True, sim=True) 77 self.PARAMS.add('cone_theta', scope='global', units='rad', desc='The isotropic cone opening half-angle', py_type=float, set='params', err=True, sim=True) 78 self.PARAMS.add('cone_s1', scope='global', units='', desc='The isotropic cone order parameter', py_type=float, set='params', err=True, sim=True) 79 self.PARAMS.add('cone_sigma_max', scope='global', units='rad', desc='The torsion angle', py_type=float, set='params', err=True, sim=True) 80 self.PARAMS.add('params', scope='global', desc='The model parameters', py_type=list) 81 82 # Add minimisation structures. 83 self.PARAMS.add_min_data(min_stats_global=True)
84 85
86 - def _assemble_limit_arrays(self):
87 """Assemble and return the limit vectors. 88 89 @return: The lower and upper limit vectors. 90 @rtype: numpy rank-1 array, numpy rank-1 array 91 """ 92 93 # Init. 94 lower = zeros(len(cdp.params), float64) 95 upper = 2.0*pi * ones(len(cdp.params), float64) 96 97 # Return the arrays. 98 return lower, upper
99 100
101 - def _assemble_param_vector(self, sim_index=None):
102 """Assemble and return the parameter vector. 103 104 @return: The parameter vector. 105 @rtype: numpy rank-1 array 106 @keyword sim_index: The Monte Carlo simulation index. 107 @type sim_index: int 108 """ 109 110 # Normal values. 111 if sim_index == None: 112 # Initialise the parameter array using the tensor rotation Euler angles (average domain position). 113 if cdp.model in ['free rotor', 'iso cone, torsionless', 'iso cone, free rotor']: 114 param_vect = [cdp.ave_pos_beta, cdp.ave_pos_gamma] 115 else: 116 param_vect = [cdp.ave_pos_alpha, cdp.ave_pos_beta, cdp.ave_pos_gamma] 117 118 # Frame order eigenframe - the full frame. 119 if cdp.model in ['iso cone', 'pseudo-ellipse', 'pseudo-ellipse, torsionless', 'pseudo-ellipse, free rotor']: 120 param_vect.append(cdp.eigen_alpha) 121 param_vect.append(cdp.eigen_beta) 122 param_vect.append(cdp.eigen_gamma) 123 124 # Frame order eigenframe - the isotropic cone axis. 125 elif cdp.model in ['free rotor', 'iso cone, torsionless', 'iso cone, free rotor', 'rotor']: 126 param_vect.append(cdp.axis_theta) 127 param_vect.append(cdp.axis_phi) 128 129 # Cone parameters - pseudo-elliptic cone parameters. 130 if cdp.model in ['pseudo-ellipse', 'pseudo-ellipse, torsionless', 'pseudo-ellipse, free rotor']: 131 param_vect.append(cdp.cone_theta_x) 132 param_vect.append(cdp.cone_theta_y) 133 134 # Cone parameters - single isotropic angle or order parameter. 135 elif cdp.model in ['iso cone', 'iso cone, torsionless']: 136 param_vect.append(cdp.cone_theta) 137 elif cdp.model in ['iso cone, free rotor']: 138 param_vect.append(cdp.cone_s1) 139 140 # Cone parameters - torsion angle. 141 if cdp.model in ['rotor', 'line', 'iso cone', 'pseudo-ellipse']: 142 param_vect.append(cdp.cone_sigma_max) 143 144 # Simulation values. 145 else: 146 # Initialise the parameter array using the tensor rotation Euler angles (average domain position). 147 if cdp.model in ['free rotor', 'iso cone, torsionless', 'iso cone, free rotor']: 148 param_vect = [cdp.ave_pos_beta_sim[sim_index], cdp.ave_pos_gamma_sim[sim_index]] 149 else: 150 param_vect = [cdp.ave_pos_alpha_sim[sim_index], cdp.ave_pos_beta_sim[sim_index], cdp.ave_pos_gamma_sim[sim_index]] 151 152 # Frame order eigenframe - the full frame. 153 if cdp.model in ['iso cone', 'pseudo-ellipse', 'pseudo-ellipse, torsionless', 'pseudo-ellipse, free rotor']: 154 param_vect.append(cdp.eigen_alpha_sim[sim_index]) 155 param_vect.append(cdp.eigen_beta_sim[sim_index]) 156 param_vect.append(cdp.eigen_gamma_sim[sim_index]) 157 158 # Frame order eigenframe - the isotropic cone axis. 159 elif cdp.model in ['free rotor', 'iso cone, torsionless', 'iso cone, free rotor', 'rotor']: 160 param_vect.append(cdp.axis_theta_sim[sim_index]) 161 param_vect.append(cdp.axis_phi_sim[sim_index]) 162 163 # Cone parameters - pseudo-elliptic cone parameters. 164 if cdp.model in ['pseudo-ellipse', 'pseudo-ellipse, torsionless', 'pseudo-ellipse, free rotor']: 165 param_vect.append(cdp.cone_theta_x_sim[sim_index]) 166 param_vect.append(cdp.cone_theta_y_sim[sim_index]) 167 168 # Cone parameters - single isotropic angle or order parameter. 169 elif cdp.model in ['iso cone', 'iso cone, torsionless']: 170 param_vect.append(cdp.cone_theta_sim[sim_index]) 171 elif cdp.model in ['iso cone, free rotor']: 172 param_vect.append(cdp.cone_s1_sim[sim_index]) 173 174 # Cone parameters - torsion angle. 175 if cdp.model in ['rotor', 'line', 'iso cone', 'pseudo-ellipse']: 176 param_vect.append(cdp.cone_sigma_max_sim[sim_index]) 177 178 # Return as a numpy array. 179 return array(param_vect, float64)
180 181
182 - def _back_calc(self):
183 """Back-calculation of the reduced alignment tensor. 184 185 @return: The reduced alignment tensors. 186 @rtype: numpy array 187 """ 188 189 # Get the parameter vector. 190 param_vector = self._assemble_param_vector() 191 192 # Get the data structures for optimisation using the tensors as base data sets. 193 full_tensors, red_tensors, red_tensor_err, full_in_ref_frame = self._minimise_setup_tensors() 194 195 # Set up the optimisation function. 196 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) 197 198 # Make a single function call. This will cause back calculation and the data will be stored in the class instance. 199 target.func(param_vector) 200 201 # Store the back-calculated tensors. 202 self._store_bc_tensors(target) 203 204 # Return the reduced tensors. 205 return target.red_tensors_bc
206 207
208 - def _cone_pdb(self, size=30.0, file=None, dir=None, inc=40, force=False):
209 """Create a PDB file containing a geometric object representing the Frame Order cone models. 210 211 @param size: The size of the geometric object in Angstroms. 212 @type size: float 213 @param inc: The number of increments for the filling of the cone objects. 214 @type inc: int 215 @param file: The name of the PDB file to create. 216 @type file: str 217 @param dir: The name of the directory to place the PDB file into. 218 @type dir: str 219 @param force: Flag which if set to True will cause any pre-existing file to be 220 overwritten. 221 @type force: bool 222 """ 223 224 # Test if the current data pipe exists. 225 pipes.test() 226 227 # The rigid model cannot be used here. 228 if cdp.model == 'rigid': 229 raise RelaxError("The 'rigid' frame order model has no cone representation.") 230 231 # Test for the necessary data structures. 232 if not hasattr(cdp, 'pivot'): 233 raise RelaxError("The pivot point for the domain motion has not been set.") 234 235 # Negative cone flag. 236 neg_cone = True 237 238 # Monte Carlo simulation flag. 239 sim = False 240 num_sim = 0 241 if hasattr(cdp, 'sim_number'): 242 sim = True 243 num_sim = cdp.sim_number 244 245 # The inversion matrix. 246 inv_mat = -eye(3) 247 248 # Create the structural object. 249 structure = Internal() 250 251 # Create model for the positive and negative images. 252 model = structure.add_model(model=1) 253 if neg_cone: 254 model_neg = structure.add_model(model=2) 255 256 # Create the molecule. 257 structure.add_molecule(name=cdp.model) 258 259 # Alias the molecules. 260 mol = model.mol[0] 261 if neg_cone: 262 mol_neg = model_neg.mol[0] 263 264 265 # The pivot point. 266 ################## 267 268 # Add the pivot point. 269 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') 270 271 272 # The central axis. 273 ################### 274 275 # Print out. 276 print("\nGenerating the central axis.") 277 278 # The spherical angles. 279 if cdp.model in ['free rotor', 'iso cone, torsionless', 'iso cone, free rotor', 'rotor']: 280 theta_name = 'axis_theta' 281 phi_name = 'axis_phi' 282 else: 283 theta_name = 'eigen_beta' 284 phi_name = 'eigen_gamma' 285 286 # The axis. 287 axis = zeros(3, float64) 288 spherical_to_cartesian([1.0, getattr(cdp, theta_name), getattr(cdp, phi_name)], axis) 289 print("Central axis: %s." % axis) 290 291 # Rotations and inversions. 292 axis_pos = axis 293 axis_neg = dot(inv_mat, axis) 294 295 # Simulation central axis. 296 axis_sim_pos = None 297 axis_sim_neg = None 298 if sim: 299 # Init. 300 axis_sim = zeros((cdp.sim_number, 3), float64) 301 302 # Fill the structure. 303 for i in range(cdp.sim_number): 304 spherical_to_cartesian([1.0, getattr(cdp, theta_name+'_sim')[i], getattr(cdp, phi_name+'_sim')[i]], axis_sim[i]) 305 306 # Inversion. 307 axis_sim_pos = axis_sim 308 axis_sim_neg = transpose(dot(inv_mat, transpose(axis_sim_pos))) 309 310 # Generate the axis vectors. 311 print("\nGenerating the axis vectors.") 312 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) 313 314 # The negative. 315 if neg_cone: 316 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) 317 318 319 # The x and y axes. 320 ################### 321 322 # Skip models missing the full eigenframe. 323 if cdp.model not in ['free rotor', 'iso cone, torsionless', 'iso cone, free rotor', 'rotor']: 324 # Print out. 325 print("\nGenerating the x and y axes.") 326 327 # The axis system. 328 axes = zeros((3, 3), float64) 329 euler_to_R_zyz(cdp.eigen_alpha, cdp.eigen_beta, cdp.eigen_gamma, axes) 330 print("Axis system:\n%s" % axes) 331 332 # Rotations and inversions. 333 axes_pos = axes 334 axes_neg = dot(inv_mat, axes) 335 336 # Simulation central axis. 337 axes_sim_pos = None 338 axes_sim_neg = None 339 if sim: 340 # Init. 341 axes_sim = zeros((3, cdp.sim_number, 3), float64) 342 343 # Fill the structure. 344 for i in range(cdp.sim_number): 345 euler_to_R_zyz(cdp.eigen_alpha_sim[i], cdp.eigen_beta_sim[i], cdp.eigen_gamma_sim[i], axes_sim[:, i]) 346 347 # Rotation and inversion. 348 axes_sim_pos = axes_sim 349 axes_sim_neg = dot(inv_mat, axes_sim_pos) 350 351 # Generate the axis vectors. 352 print("\nGenerating the axis vectors.") 353 label = ['x', 'y'] 354 for i in range(2): 355 # Simulation structures. 356 if sim: 357 axis_sim_pos = axes_sim_pos[:, i] 358 axis_sim_neg = axes_sim_neg[:, i] 359 360 # The vectors. 361 res_num = generate_vector_residues(mol=mol, vector=axes_pos[:, i], atom_name='%s-ax'%label[i], res_name_vect='%sAX'%label[i].upper(), sim_vectors=axis_sim_pos, res_num=res_num+1, origin=cdp.pivot, scale=size) 362 if neg_cone: 363 res_num = generate_vector_residues(mol=mol_neg, vector=axes_neg[:, i], atom_name='%s-ax'%label[i], res_name_vect='%sAX'%label[i].upper(), sim_vectors=axis_sim_neg, res_num=res_num, origin=cdp.pivot, scale=size) 364 365 366 # The cone object. 367 ################## 368 369 # Skip models missing a cone. 370 if cdp.model not in ['rotor', 'free rotor']: 371 # The rotation matrix (rotation from the z-axis to the cone axis). 372 if cdp.model not in ['iso cone, torsionless', 'iso cone, free rotor']: 373 R = axes 374 else: 375 R = zeros((3, 3), float64) 376 two_vect_to_R(array([0, 0, 1], float64), axis, R) 377 378 # Average position rotation. 379 R_pos = R 380 R_neg = dot(inv_mat, R) 381 382 # The pseudo-ellipse cone object. 383 if cdp.model in ['pseudo-ellipse', 'pseudo-ellipse, torsionless', 'pseudo-ellipse, free rotor']: 384 cone = Pseudo_elliptic(cdp.cone_theta_x, cdp.cone_theta_y) 385 386 # The isotropic cone object. 387 else: 388 cone = Iso_cone(cdp.cone_theta) 389 390 # Create the positive and negative cones. 391 create_cone_pdb(mol=mol, cone=cone, start_res=mol.res_num[-1]+1, apex=cdp.pivot, R=R_pos, inc=inc, distribution='regular') 392 393 # The negative. 394 if neg_cone: 395 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') 396 397 398 # Create the PDB file. 399 ###################### 400 401 # Print out. 402 print("\nGenerating the PDB file.") 403 404 # Write the file. 405 pdb_file = open_write_file(file, dir, force=force) 406 structure.write_pdb(pdb_file) 407 pdb_file.close()
408 409
410 - def _domain_to_pdb(self, domain=None, pdb=None):
411 """Match domains to PDB files. 412 413 @keyword domain: The domain to associate the PDB file to. 414 @type domain: str 415 @keyword pdb: The PDB file to associate the domain to. 416 @type pdb: str 417 """ 418 419 # Check that the domain exists. 420 exists = False 421 for i in range(len(cdp.align_tensors)): 422 if hasattr(cdp.align_tensors[i], 'domain') and domain == cdp.align_tensors[i].domain: 423 exists = True 424 if not exists: 425 raise RelaxError("The domain '%s' cannot be found" % domain) 426 427 # Init if needed. 428 if not hasattr(cdp, 'domain_to_pdb'): 429 cdp.domain_to_pdb = [] 430 431 # Strip the file ending if given. 432 if search('.pdb$', pdb): 433 pdb = pdb[:-4] 434 435 # Add the data. 436 cdp.domain_to_pdb.append([domain, pdb])
437 438
439 - def _grid_row(self, incs, lower, upper, dist_type=None, end_point=True):
440 """Set up a row of the grid search for a given parameter. 441 442 @param incs: The number of increments. 443 @type incs: int 444 @param lower: The lower bounds. 445 @type lower: float 446 @param upper: The upper bounds. 447 @type upper: float 448 @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). 449 @type dist_type: None or str 450 @keyword end_point: A flag which if False will cause the end point to be removed. 451 @type end_point: bool 452 @return: The row of the grid. 453 @rtype: list of float 454 """ 455 456 # Init. 457 row = [] 458 459 # Linear grid. 460 if dist_type == None: 461 # Loop over the increments. 462 for i in range(incs): 463 # The row. 464 row.append(lower + i * (upper - lower) / (incs - 1.0)) 465 466 # Inverse cos distribution. 467 elif dist_type == 'acos': 468 # Generate the increment values of v from cos(upper) to cos(lower). 469 v = zeros(incs, float64) 470 val = (cos(lower) - cos(upper)) / (incs - 1.0) 471 for i in range(incs): 472 v[-i-1] = cos(upper) + float(i) * val 473 474 # Generate the distribution. 475 row = arccos(v) 476 477 # Remove the last point. 478 if not end_point: 479 row = row[:-1] 480 481 # Return the row (as a list). 482 return list(row)
483 484
485 - def _minimise_setup_tensors(self, sim_index=None):
486 """Set up the data structures for optimisation using alignment tensors as base data sets. 487 488 @keyword sim_index: The simulation index. This should be None if normal optimisation is 489 desired. 490 @type sim_index: None or int 491 @return: The assembled data structures for using alignment tensors as the base 492 data for optimisation. These include: 493 - full_tensors, the full tensors as concatenated arrays. 494 - red_tensors, the reduced tensors as concatenated arrays. 495 - red_err, the reduced tensor errors as concatenated arrays. 496 @rtype: tuple of 3 numpy nx5D, rank-1 arrays 497 """ 498 499 # Checks. 500 if not hasattr(cdp, 'ref_domain'): 501 raise RelaxError("The reference domain has not been set up.") 502 if not hasattr(cdp.align_tensors, 'reduction'): 503 raise RelaxError("The tensor reductions have not been specified.") 504 for i, tensor in self._tensor_loop(): 505 if not hasattr(tensor, 'domain'): 506 raise RelaxError("The domain that the '%s' tensor is attached to has not been set" % tensor.name) 507 508 # Initialise. 509 n = len(cdp.align_tensors.reduction) 510 full_tensors = zeros(n*5, float64) 511 red_tensors = zeros(n*5, float64) 512 red_err = ones(n*5, float64) * 1e-5 513 full_in_ref_frame = zeros(n, float64) 514 515 # Loop over the full tensors. 516 for i, tensor in self._tensor_loop(red=False): 517 # The full tensor. 518 full_tensors[5*i + 0] = tensor.Axx 519 full_tensors[5*i + 1] = tensor.Ayy 520 full_tensors[5*i + 2] = tensor.Axy 521 full_tensors[5*i + 3] = tensor.Axz 522 full_tensors[5*i + 4] = tensor.Ayz 523 524 # The full tensor corresponds to the frame of reference. 525 if cdp.ref_domain == tensor.domain: 526 full_in_ref_frame[i] = 1 527 528 # Loop over the reduced tensors. 529 for i, tensor in self._tensor_loop(red=True): 530 # The reduced tensor (simulation data). 531 if sim_index != None: 532 red_tensors[5*i + 0] = tensor.Axx_sim[sim_index] 533 red_tensors[5*i + 1] = tensor.Ayy_sim[sim_index] 534 red_tensors[5*i + 2] = tensor.Axy_sim[sim_index] 535 red_tensors[5*i + 3] = tensor.Axz_sim[sim_index] 536 red_tensors[5*i + 4] = tensor.Ayz_sim[sim_index] 537 538 # The reduced tensor. 539 else: 540 red_tensors[5*i + 0] = tensor.Axx 541 red_tensors[5*i + 1] = tensor.Ayy 542 red_tensors[5*i + 2] = tensor.Axy 543 red_tensors[5*i + 3] = tensor.Axz 544 red_tensors[5*i + 4] = tensor.Ayz 545 546 # The reduced tensor errors. 547 if hasattr(tensor, 'Axx_err'): 548 red_err[5*i + 0] = tensor.Axx_err 549 red_err[5*i + 1] = tensor.Ayy_err 550 red_err[5*i + 2] = tensor.Axy_err 551 red_err[5*i + 3] = tensor.Axz_err 552 red_err[5*i + 4] = tensor.Ayz_err 553 554 # Return the data structures. 555 return full_tensors, red_tensors, red_err, full_in_ref_frame
556 557
558 - def _pivot(self, pivot=None):
559 """Set the pivot point for the 2 body motion. 560 561 @param pivot: The pivot point of the two bodies (domains, etc.) in the structural 562 coordinate system. 563 @type pivot: list of num 564 """ 565 566 # Test if the current data pipe exists. 567 pipes.test() 568 569 # Set the pivot point. 570 cdp.pivot = pivot 571 572 # Convert to floats. 573 for i in range(3): 574 cdp.pivot[i] = float(cdp.pivot[i])
575 576
577 - def _ref_domain(self, ref=None):
578 """Set the reference domain for the frame order, multi-domain models. 579 580 @param ref: The reference domain. 581 @type ref: str 582 """ 583 584 # Test if the current data pipe exists. 585 pipes.test() 586 587 # Test if the reference domain exists. 588 exists = False 589 for tensor_cont in cdp.align_tensors: 590 if hasattr(tensor_cont, 'domain') and tensor_cont.domain == ref: 591 exists = True 592 if not exists: 593 raise RelaxError("The reference domain cannot be found within any of the loaded tensors.") 594 595 # Set the reference domain. 596 cdp.ref_domain = ref 597 598 # Update the model. 599 if hasattr(cdp, 'model'): 600 self._update_model()
601 602
603 - def _select_model(self, model=None):
604 """Select the Frame Order model. 605 606 @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'. 607 @type model: str 608 """ 609 610 # Test if the current data pipe exists. 611 pipes.test() 612 613 # Test if the model name exists. 614 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']: 615 raise RelaxError("The model name " + repr(model) + " is invalid.") 616 617 # Set the model 618 cdp.model = model 619 620 # Initialise the list of model parameters. 621 cdp.params = [] 622 623 # Update the model. 624 self._update_model()
625 626
627 - def _store_bc_tensors(self, target_fn):
628 """Store the back-calculated reduced alignment tensors. 629 630 @param target_fn: The frame-order target function class. 631 @type target_fn: class instance 632 """ 633 634 # Loop over the reduced tensors. 635 for i, tensor in self._tensor_loop(red=True): 636 # New name. 637 name = tensor.name + ' bc' 638 639 # Initialise the new tensor. 640 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)
641 642
643 - def _tensor_loop(self, red=False):
644 """Generator method for looping over the full or reduced tensors. 645 646 @keyword red: A flag which if True causes the reduced tensors to be returned, and if False 647 the full tensors are returned. 648 @type red: bool 649 @return: The tensor index and the tensor. 650 @rtype: (int, AlignTensorData instance) 651 """ 652 653 # Number of tensor pairs. 654 n = len(cdp.align_tensors.reduction) 655 656 # Alias. 657 data = cdp.align_tensors 658 list = data.reduction 659 660 # Full or reduced index. 661 if red: 662 index = 1 663 else: 664 index = 0 665 666 # Loop over the reduction list. 667 for i in range(n): 668 yield i, data[list[i][index]]
669 670
671 - def _update_model(self):
672 """Update the model parameters as necessary.""" 673 674 # Initialise the list of model parameters. 675 if not hasattr(cdp, 'params'): 676 cdp.params = [] 677 678 # Build the parameter name list. 679 if not len(cdp.params): 680 # The tensor rotation, or average domain position. 681 if cdp.model not in ['free rotor', 'iso cone, torsionless', 'iso cone, free rotor']: 682 cdp.params.append('ave_pos_alpha') 683 cdp.params.append('ave_pos_beta') 684 cdp.params.append('ave_pos_gamma') 685 686 # Frame order eigenframe - the full frame. 687 if cdp.model in ['iso cone', 'pseudo-ellipse', 'pseudo-ellipse, torsionless', 'pseudo-ellipse, free rotor']: 688 cdp.params.append('eigen_alpha') 689 cdp.params.append('eigen_beta') 690 cdp.params.append('eigen_gamma') 691 692 # Frame order eigenframe - the isotropic cone axis. 693 elif cdp.model in ['free rotor', 'iso cone, torsionless', 'iso cone, free rotor', 'rotor']: 694 cdp.params.append('axis_theta') 695 cdp.params.append('axis_phi') 696 697 # Cone parameters - pseudo-elliptic cone parameters. 698 if cdp.model in ['pseudo-ellipse', 'pseudo-ellipse, torsionless', 'pseudo-ellipse, free rotor']: 699 cdp.params.append('cone_theta_x') 700 cdp.params.append('cone_theta_y') 701 702 # Cone parameters - single isotropic angle or order parameter. 703 elif cdp.model in ['iso cone', 'iso cone, torsionless']: 704 cdp.params.append('cone_theta') 705 elif cdp.model in ['iso cone, free rotor']: 706 cdp.params.append('cone_s1') 707 708 # Cone parameters - torsion angle. 709 if cdp.model in ['rotor', 'line', 'iso cone', 'pseudo-ellipse']: 710 cdp.params.append('cone_sigma_max') 711 712 # Initialise the parameters in the current data pipe. 713 for param in cdp.params: 714 if not hasattr(cdp, param): 715 setattr(cdp, param, 0.0)
716 717
718 - def _unpack_opt_results(self, results, sim_index=None):
719 """Unpack and store the Frame Order optimisation results. 720 721 @param results: The results tuple returned by the minfx generic_minimise() function. 722 @type results: tuple 723 @param sim_index: The index of the simulation to optimise. This should be None for normal 724 optimisation. 725 @type sim_index: None or int 726 """ 727 728 # Disassemble the results. 729 if len(results) == 4: # Grid search. 730 param_vector, func, iter_count, warning = results 731 f_count = iter_count 732 g_count = 0.0 733 h_count = 0.0 734 else: 735 param_vector, func, iter_count, f_count, g_count, h_count, warning = results 736 737 # Catch infinite chi-squared values. 738 if isInf(func): 739 raise RelaxInfError('chi-squared') 740 741 # Catch chi-squared values of NaN. 742 if isNaN(func): 743 raise RelaxNaNError('chi-squared') 744 745 # Unpack the parameters. 746 ave_pos_alpha, ave_pos_beta, ave_pos_gamma = None, None, None 747 eigen_alpha, eigen_beta, eigen_gamma = None, None, None 748 axis_theta, axis_phi = None, None 749 cone_theta_x, cone_theta_y = None, None 750 cone_theta = None 751 cone_s1 = None 752 cone_sigma_max = None 753 if cdp.model == 'pseudo-ellipse': 754 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 755 elif cdp.model in ['pseudo-ellipse, torsionless', 'pseudo-ellipse, free rotor']: 756 ave_pos_alpha, ave_pos_beta, ave_pos_gamma, eigen_alpha, eigen_beta, eigen_gamma, cone_theta_x, cone_theta_y = param_vector 757 elif cdp.model == 'iso cone': 758 ave_pos_alpha, ave_pos_beta, ave_pos_gamma, eigen_alpha, eigen_beta, eigen_gamma, cone_theta, cone_sigma_max = param_vector 759 elif cdp.model in ['iso cone, torsionless']: 760 ave_pos_beta, ave_pos_gamma, axis_theta, axis_phi, cone_theta = param_vector 761 ave_pos_alpha = None 762 elif cdp.model in ['iso cone, free rotor']: 763 ave_pos_beta, ave_pos_gamma, axis_theta, axis_phi, cone_s1 = param_vector 764 ave_pos_alpha = None 765 elif cdp.model == 'line': 766 ave_pos_alpha, ave_pos_beta, ave_pos_gamma, eigen_alpha, eigen_beta, eigen_gamma, cone_theta_x, cone_sigma_max = param_vector 767 elif cdp.model in ['line, torsionless', 'line, free rotor']: 768 ave_pos_alpha, ave_pos_beta, ave_pos_gamma, eigen_alpha, eigen_beta, eigen_gamma, cone_theta_x, cone_sigma_max = param_vector 769 elif cdp.model in ['rotor']: 770 ave_pos_alpha, ave_pos_beta, ave_pos_gamma, axis_theta, axis_phi, cone_sigma_max = param_vector 771 elif cdp.model in ['free rotor']: 772 ave_pos_beta, ave_pos_gamma, axis_theta, axis_phi = param_vector 773 ave_pos_alpha = None 774 elif cdp.model == 'rigid': 775 ave_pos_alpha, ave_pos_beta, ave_pos_gamma = param_vector 776 777 # Monte Carlo simulation data structures. 778 if sim_index != None: 779 # Average position parameters. 780 if ave_pos_alpha != None: 781 cdp.ave_pos_alpha_sim[sim_index] = wrap_angles(ave_pos_alpha, 0.0, 2.0*pi) 782 if ave_pos_beta != None: 783 cdp.ave_pos_beta_sim[sim_index] = wrap_angles(ave_pos_beta, 0.0, 2.0*pi) 784 if ave_pos_gamma != None: 785 cdp.ave_pos_gamma_sim[sim_index] = wrap_angles(ave_pos_gamma, 0.0, 2.0*pi) 786 787 # Eigenframe parameters. 788 if eigen_alpha != None: 789 cdp.eigen_alpha_sim[sim_index] = wrap_angles(eigen_alpha, 0.0, 2.0*pi) 790 if eigen_beta != None: 791 cdp.eigen_beta_sim[sim_index] = wrap_angles(eigen_beta, 0.0, 2.0*pi) 792 if eigen_gamma != None: 793 cdp.eigen_gamma_sim[sim_index] = wrap_angles(eigen_gamma, 0.0, 2.0*pi) 794 if axis_theta != None: 795 cdp.axis_theta_sim[sim_index] = axis_theta 796 if axis_phi != None: 797 cdp.axis_phi_sim[sim_index] = axis_phi 798 799 # Cone parameters. 800 if cone_theta != None: 801 cdp.cone_theta_sim[sim_index] = cone_theta 802 if cone_s1 != None: 803 cdp.cone_s1_sim[sim_index] = cone_s1 804 cdp.cone_theta_sim[sim_index] = order_parameters.iso_cone_S_to_theta(cone_s1) 805 if cone_theta_x != None: 806 cdp.cone_theta_x_sim[sim_index] = cone_theta_x 807 if cone_theta_y != None: 808 cdp.cone_theta_y_sim[sim_index] = cone_theta_y 809 if cone_sigma_max != None: 810 cdp.cone_sigma_max_sim[sim_index] = abs(cone_sigma_max) 811 812 # Optimisation stats. 813 cdp.chi2_sim[sim_index] = func 814 cdp.iter_sim[sim_index] = iter_count 815 cdp.f_count_sim[sim_index] = f_count 816 cdp.g_count_sim[sim_index] = g_count 817 cdp.h_count_sim[sim_index] = h_count 818 cdp.warning_sim[sim_index] = warning 819 820 # Normal data structures. 821 else: 822 # Average position parameters. 823 if ave_pos_alpha != None: 824 cdp.ave_pos_alpha = wrap_angles(ave_pos_alpha, 0.0, 2.0*pi) 825 if ave_pos_beta != None: 826 cdp.ave_pos_beta = wrap_angles(ave_pos_beta, 0.0, 2.0*pi) 827 if ave_pos_gamma != None: 828 cdp.ave_pos_gamma = wrap_angles(ave_pos_gamma, 0.0, 2.0*pi) 829 830 # Eigenframe parameters. 831 if eigen_alpha != None: 832 cdp.eigen_alpha = wrap_angles(eigen_alpha, 0.0, 2.0*pi) 833 if eigen_beta != None: 834 cdp.eigen_beta = wrap_angles(eigen_beta, 0.0, 2.0*pi) 835 if eigen_gamma != None: 836 cdp.eigen_gamma = wrap_angles(eigen_gamma, 0.0, 2.0*pi) 837 if axis_theta != None: 838 cdp.axis_theta = axis_theta 839 if axis_phi != None: 840 cdp.axis_phi = axis_phi 841 842 # Cone parameters. 843 if cone_theta != None: 844 cdp.cone_theta = cone_theta 845 if cone_s1 != None: 846 cdp.cone_s1 = cone_s1 847 cdp.cone_theta = order_parameters.iso_cone_S_to_theta(cone_s1) 848 if cone_theta_x != None: 849 cdp.cone_theta_x = cone_theta_x 850 if cone_theta_y != None: 851 cdp.cone_theta_y = cone_theta_y 852 if cone_sigma_max != None: 853 cdp.cone_sigma_max = abs(cone_sigma_max) 854 855 # Optimisation stats. 856 cdp.chi2 = func 857 cdp.iter = iter_count 858 cdp.f_count = f_count 859 cdp.g_count = g_count 860 cdp.h_count = h_count 861 cdp.warning = warning
862 863
864 - def base_data_loop(self):
865 """Generator method for looping nothing. 866 867 The loop essentially consists of a single element. 868 869 @return: Nothing. 870 @rtype: None 871 """ 872 873 yield None
874 875
876 - def calculate(self, spin_id=None, verbosity=1, sim_index=None):
877 """Calculate the chi-squared value for the current parameter values. 878 879 @keyword spin_id: The spin identification string (unused). 880 @type spin_id: None 881 @keyword verbosity: The amount of information to print. The higher the value, the greater the verbosity. 882 @type verbosity: int 883 @keyword sim_index: The optional MC simulation index (unused). 884 @type sim_index: None or int 885 """ 886 887 # Assemble the parameter vector. 888 param_vector = self._assemble_param_vector() 889 890 # Get the data structures for optimisation using the tensors as base data sets. 891 full_tensors, red_tensors, red_tensor_err, full_in_ref_frame = self._minimise_setup_tensors() 892 893 # Set up the optimisation function. 894 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) 895 896 # Make a single function call. This will cause back calculation and the data will be stored in the class instance. 897 chi2 = target.func(param_vector) 898 899 # Set the chi2. 900 cdp.chi2 = chi2 901 902 # Store the back-calculated tensors. 903 self._store_bc_tensors(target)
904 905
906 - def create_mc_data(self, data_id=None):
907 """Create the Monte Carlo data by back calculating the reduced tensor data. 908 909 @keyword data_id: Unused. 910 @type data_id: None 911 @return: The Monte Carlo simulation data. 912 @rtype: list of floats 913 """ 914 915 # Back calculate the tensors. 916 red_tensors_bc = self._back_calc() 917 918 # Return the data. 919 return red_tensors_bc
920 921
922 - def get_param_names(self, model_info=None):
923 """Return a vector of parameter names. 924 925 @keyword model_info: The model index from model_info(). 926 @type model_info: int 927 @return: The vector of parameter names. 928 @rtype: list of str 929 """ 930 931 # First update the model, if needed. 932 self._update_model() 933 934 # Return the parameter list object. 935 return cdp.params
936 937
938 - def get_param_values(self, model_info=None, sim_index=None):
939 """Return a vector of parameter values. 940 941 @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). 942 @type model_info: int 943 @keyword sim_index: The Monte Carlo simulation index. 944 @type sim_index: int 945 @return: The vector of parameter values. 946 @rtype: list of str 947 """ 948 949 # Assemble the values and return it. 950 return self._assemble_param_vector(sim_index=sim_index)
951 952
953 - def grid_search(self, lower=None, upper=None, inc=None, constraints=False, verbosity=0, sim_index=None):
954 """Perform a grid search. 955 956 @keyword lower: The lower bounds of the grid search which must be equal to the 957 number of parameters in the model. 958 @type lower: list of float 959 @keyword upper: The upper bounds of the grid search which must be equal to the 960 number of parameters in the model. 961 @type upper: list of float 962 @keyword inc: The increments for each dimension of the space for the grid search. 963 The number of elements in the array must equal to the number of 964 parameters in the model. 965 @type inc: int or list of int 966 @keyword constraints: If True, constraints are applied during the grid search (eliminating 967 parts of the grid). If False, no constraints are used. 968 @type constraints: bool 969 @keyword verbosity: A flag specifying the amount of information to print. The higher 970 the value, the greater the verbosity. 971 @type verbosity: int 972 @keyword sim_index: The Monte Carlo simulation index. 973 @type sim_index: None or int 974 """ 975 976 # Test if the Frame Order model has been set up. 977 if not hasattr(cdp, 'model'): 978 raise RelaxNoModelError('Frame Order') 979 980 # The number of parameters. 981 n = len(cdp.params) 982 983 # If inc is an int, convert it into an array of that value. 984 if isinstance(inc, int): 985 incs = [inc]*n 986 else: 987 incs = inc 988 989 # Sanity check. 990 if len(incs) != n: 991 raise RelaxError("The size of the increment list %s does not match the number of parameters in %s." % (incs, cdp.params)) 992 993 # Initialise the grid increments structures. 994 lower_list = lower 995 upper_list = upper 996 grid = [] 997 """This structure is a list of lists. The first dimension corresponds to the model 998 parameter. The second dimension are the grid node positions.""" 999 1000 # Generate the grid. 1001 for i in range(n): 1002 # Fixed parameter. 1003 if incs[i] == None: 1004 grid.append(None) 1005 continue 1006 1007 # Reset. 1008 dist_type = None 1009 end_point = True 1010 1011 # Linear angle grid from 0 to one inc before 2pi. 1012 if cdp.params[i] in ['ave_pos_alpha', 'ave_pos_gamma', 'eigen_alpha', 'eigen_gamma', 'axis_phi']: 1013 lower = 0.0 1014 upper = 2*pi * (1.0 - 1.0/incs[i]) 1015 1016 # Arccos grid from 0 to pi. 1017 if cdp.params[i] in ['ave_pos_beta', 'eigen_beta', 'axis_theta']: 1018 # Change the default increment numbers. 1019 if not isinstance(inc, list): 1020 incs[i] = int(incs[i] / 2) + 1 1021 1022 # The distribution type and end point. 1023 dist_type = 'acos' 1024 end_point = False 1025 1026 # Set the default bounds. 1027 lower = 0.0 1028 upper = pi 1029 1030 # The isotropic cone order parameter. 1031 if cdp.params[i] == 'cone_s1': 1032 lower = -0.125 1033 upper = 1.0 1034 1035 # Linear angle grid from 0 to pi excluding the outer points. 1036 if cdp.params[i] in ['cone_theta', 'cone_theta_x', 'cone_theta_y', 'cone_sigma_max']: 1037 lower = pi * (1.0/incs[i]) 1038 upper = pi * (1.0 - 1.0/incs[i]) 1039 1040 # Over-ride the bounds. 1041 if lower_list: 1042 lower = lower_list[i] 1043 if upper_list: 1044 upper = upper_list[i] 1045 1046 # Append the grid row. 1047 row = self._grid_row(incs[i], lower, upper, dist_type=dist_type, end_point=end_point) 1048 grid.append(row) 1049 1050 # Remove an inc if the end point has been removed. 1051 if not end_point: 1052 incs[i] -= 1 1053 1054 # Total number of points. 1055 total_pts = 1 1056 for i in range(n): 1057 # Fixed parameter. 1058 if grid[i] == None: 1059 continue 1060 1061 total_pts = total_pts * len(grid[i]) 1062 1063 # Check the number. 1064 max_pts = 50e6 1065 if total_pts > max_pts: 1066 raise RelaxError("The total number of grid points '%s' exceeds the maximum of '%s'." % (total_pts, int(max_pts))) 1067 1068 # Build the points array. 1069 pts = zeros((total_pts, n), float64) 1070 indices = zeros(n, int) 1071 for i in range(total_pts): 1072 # Loop over the dimensions. 1073 for j in range(n): 1074 # Fixed parameter. 1075 if grid[j] == None: 1076 # Get the current parameter value. 1077 pts[i, j] = getattr(cdp, cdp.params[j]) 1078 1079 # Add the point coordinate. 1080 else: 1081 pts[i, j] = grid[j][indices[j]] 1082 1083 # Increment the step positions. 1084 for j in range(n): 1085 if incs[j] != None and indices[j] < incs[j]-1: 1086 indices[j] += 1 1087 break # Exit so that the other step numbers are not incremented. 1088 else: 1089 indices[j] = 0 1090 1091 # Minimisation. 1092 self.minimise(min_algor='grid', min_options=pts, constraints=constraints, verbosity=verbosity, sim_index=sim_index)
1093 1094
1095 - def is_spin_param(self, name):
1096 """State that the parameter is not spin specific. 1097 1098 @param name: The name of the parameter. 1099 @type name: str 1100 @return: False. 1101 @rtype: bool 1102 """ 1103 1104 # Not spin specific! 1105 return False
1106 1107
1108 - def map_bounds(self, param, spin_id=None):
1109 """Create bounds for the OpenDX mapping function. 1110 1111 @param param: The name of the parameter to return the lower and upper bounds of. 1112 @type param: str 1113 @param spin_id: The spin identification string (unused). 1114 @type spin_id: None 1115 @return: The upper and lower bounds of the parameter. 1116 @rtype: list of float 1117 """ 1118 1119 # Average position Euler angles. 1120 if param in ['ave_pos_alpha', 'ave_pos_beta', 'ave_pos_gamma']: 1121 return [0.0, 2*pi] 1122 1123 # Axis spherical coordinate theta. 1124 if param == 'axis_theta': 1125 return [0.0, pi] 1126 1127 # Axis spherical coordinate phi. 1128 if param == 'axis_phi': 1129 return [0.0, 2*pi] 1130 1131 # Cone angle. 1132 if param == 'cone_theta': 1133 return [0.0, pi]
1134 1135
1136 - 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):
1137 """Minimisation function. 1138 1139 @param min_algor: The minimisation algorithm to use. 1140 @type min_algor: str 1141 @param min_options: An array of options to be used by the minimisation algorithm. 1142 @type min_options: array of str 1143 @param func_tol: The function tolerance which, when reached, terminates optimisation. 1144 Setting this to None turns of the check. 1145 @type func_tol: None or float 1146 @param grad_tol: The gradient tolerance which, when reached, terminates optimisation. 1147 Setting this to None turns of the check. 1148 @type grad_tol: None or float 1149 @param max_iterations: The maximum number of iterations for the algorithm. 1150 @type max_iterations: int 1151 @param constraints: If True, constraints are used during optimisation. 1152 @type constraints: bool 1153 @param scaling: If True, diagonal scaling is enabled during optimisation to allow 1154 the problem to be better conditioned. 1155 @type scaling: bool 1156 @param verbosity: A flag specifying the amount of information to print. The higher 1157 the value, the greater the verbosity. 1158 @type verbosity: int 1159 @param sim_index: The index of the simulation to optimise. This should be None if 1160 normal optimisation is desired. 1161 @type sim_index: None or int 1162 @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. 1163 @type lower: array of numbers 1164 @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. 1165 @type upper: array of numbers 1166 @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. 1167 @type inc: array of int 1168 """ 1169 1170 # Constraints not implemented yet. 1171 if constraints: 1172 # Turn the constraints off. 1173 constraints = False 1174 1175 # Remove the method of multipliers arg. 1176 if not search('^[Gg]rid', min_algor): 1177 min_algor = min_options[0] 1178 min_options = min_options[1:] 1179 1180 # Throw a warning. 1181 warn(RelaxWarning("Constraints are as of yet not implemented - turning this option off.")) 1182 1183 # Simulated annealing constraints. 1184 lower, upper = self._assemble_limit_arrays() 1185 1186 # Assemble the parameter vector. 1187 param_vector = self._assemble_param_vector() 1188 1189 # Get the data structures for optimisation using the tensors as base data sets. 1190 full_tensors, red_tensors, red_tensor_err, full_in_ref_frame = self._minimise_setup_tensors(sim_index) 1191 1192 # Set up the optimisation function. 1193 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) 1194 1195 # Grid search. 1196 if search('^[Gg]rid', min_algor): 1197 results = grid_point_array(func=target.func, args=(), points=min_options, verbosity=verbosity) 1198 1199 # Minimisation. 1200 else: 1201 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) 1202 1203 # Unpack the results. 1204 self._unpack_opt_results(results, sim_index) 1205 1206 # Store the back-calculated tensors. 1207 self._store_bc_tensors(target)
1208 1209 1210
1211 - def model_loop(self):
1212 """Dummy generator method. 1213 1214 In this case only a single model per spin system is assumed. Hence the yielded data is the 1215 spin container object. 1216 1217 1218 @return: Information about the model which for this analysis is the spin container. 1219 @rtype: SpinContainer instance 1220 """ 1221 1222 # Don't return anything, just loop once. 1223 yield None
1224 1225
1226 - def model_statistics(self, model_info=None, spin_id=None, global_stats=None):
1227 """Return the k, n, and chi2 model statistics. 1228 1229 k - number of parameters. 1230 n - number of data points. 1231 chi2 - the chi-squared value. 1232 1233 1234 @keyword model_info: Unused. 1235 @type model_info: None 1236 @keyword spin_id: The spin identification string (unused). 1237 @type spin_id: None 1238 @keyword global_stats: Unused. 1239 @type global_stats: None 1240 @return: The optimisation statistics, in tuple format, of the number of 1241 parameters (k), the number of data points (n), and the chi-squared 1242 value (chi2). 1243 @rtype: tuple of (int, int, float) 1244 """ 1245 1246 # Count the number of parameters. 1247 param_vector = self._assemble_param_vector() 1248 k = len(param_vector) 1249 1250 # The number of data points. 1251 n = len(cdp.align_tensors.reduction) 1252 1253 # The chi2 value. 1254 if not hasattr(cdp, 'chi2'): 1255 raise RelaxError("Statistics are not available, most likely because the model has not been optimised.") 1256 chi2 = cdp.chi2 1257 1258 # Return the data. 1259 return k, n, chi2
1260 1261
1262 - def return_error(self, data_id):
1263 """Return the alignment tensor error structure. 1264 1265 @param data_id: The information yielded by the base_data_loop() generator method. 1266 @type data_id: None 1267 @return: The array of tensor error values. 1268 @rtype: list of float 1269 """ 1270 1271 # Get the tensor data structures. 1272 full_tensors, red_tensors, red_tensor_err, full_in_ref_frame = self._minimise_setup_tensors() 1273 1274 # Return the errors. 1275 return red_tensor_err
1276 1277
1278 - def return_units(self, param):
1279 """Return a string representing the parameters units. 1280 1281 @param param: The name of the parameter to return the units string for. 1282 @type param: str 1283 @return: The parameter units string. 1284 @rtype: str 1285 """ 1286 1287 # Average position Euler angles. 1288 if param in ['ave_pos_alpha', 'ave_pos_beta', 'ave_pos_gamma']: 1289 return 'rad' 1290 1291 # Eigenframe angles. 1292 if param in ['eigen_alpha', 'eigen_beta', 'eigen_gamma', 'axis_theta', 'axis_phi']: 1293 return 'rad' 1294 1295 # Cone angles. 1296 if param in ['cone_theta_x', 'cone_theta_y', 'cone_sigma_max', 'cone_theta']: 1297 return 'rad'
1298 1299
1300 - def set_error(self, model_info, index, error):
1301 """Set the parameter errors. 1302 1303 @param model_info: The model information originating from model_loop() (unused). 1304 @type model_info: None 1305 @param index: The index of the parameter to set the errors for. 1306 @type index: int 1307 @param error: The error value. 1308 @type error: float 1309 """ 1310 1311 # Parameter increment counter. 1312 inc = 0 1313 1314 # Loop over the residue specific parameters. 1315 for param in self.data_names(set='params'): 1316 # Return the parameter array. 1317 if index == inc: 1318 setattr(cdp, param + "_err", error) 1319 1320 # Increment. 1321 inc = inc + 1 1322 1323 # Add some additional parameters. 1324 if cdp.model == 'iso cone, free rotor' and inc == index: 1325 setattr(cdp, 'cone_theta_err', error)
1326 1327 1328
1329 - def set_selected_sim(self, model_info, select_sim):
1330 """Set the simulation selection flag for the spin. 1331 1332 @param model_info: The model information originating from model_loop() (unused). 1333 @type model_info: None 1334 @param select_sim: The selection flag for the simulations. 1335 @type select_sim: bool 1336 """ 1337 1338 # Set the array. 1339 cdp.select_sim = deepcopy(select_sim)
1340 1341
1342 - def sim_init_values(self):
1343 """Initialise the Monte Carlo parameter values.""" 1344 1345 # Get the parameter object names. 1346 param_names = self.data_names(set='params') 1347 1348 # Add some additional parameters. 1349 if cdp.model == 'iso cone, free rotor': 1350 param_names.append('cone_theta') 1351 1352 # Get the minimisation statistic object names. 1353 min_names = self.data_names(set='min') 1354 1355 1356 # Test if Monte Carlo parameter values have already been set. 1357 ############################################################# 1358 1359 # Loop over all the parameter names. 1360 for object_name in param_names: 1361 # Name for the simulation object. 1362 sim_object_name = object_name + '_sim' 1363 1364 # Test if the simulation object already exists. 1365 if hasattr(cdp, sim_object_name): 1366 raise RelaxError("Monte Carlo parameter values have already been set.") 1367 1368 1369 # Set the Monte Carlo parameter values. 1370 ####################################### 1371 1372 # Loop over all the data names. 1373 for object_name in param_names: 1374 # Skip non-existent objects. 1375 if not hasattr(cdp, object_name): 1376 continue 1377 1378 # Name for the simulation object. 1379 sim_object_name = object_name + '_sim' 1380 1381 # Create the simulation object. 1382 setattr(cdp, sim_object_name, []) 1383 1384 # Get the simulation object. 1385 sim_object = getattr(cdp, sim_object_name) 1386 1387 # Loop over the simulations. 1388 for j in range(cdp.sim_number): 1389 # Copy and append the data. 1390 sim_object.append(deepcopy(getattr(cdp, object_name))) 1391 1392 # Loop over all the minimisation object names. 1393 for object_name in min_names: 1394 # Name for the simulation object. 1395 sim_object_name = object_name + '_sim' 1396 1397 # Create the simulation object. 1398 setattr(cdp, sim_object_name, []) 1399 1400 # Get the simulation object. 1401 sim_object = getattr(cdp, sim_object_name) 1402 1403 # Loop over the simulations. 1404 for j in range(cdp.sim_number): 1405 # Copy and append the data. 1406 sim_object.append(deepcopy(getattr(cdp, object_name)))
1407 1408
1409 - def sim_pack_data(self, data_id, sim_data):
1410 """Pack the Monte Carlo simulation data. 1411 1412 @param data_id: The spin identification string, as yielded by the base_data_loop() generator method. 1413 @type data_id: None 1414 @param sim_data: The Monte Carlo simulation data. 1415 @type sim_data: list of float 1416 """ 1417 1418 # Transpose the data. 1419 sim_data = transpose(sim_data) 1420 1421 # Loop over the reduced tensors. 1422 for i, tensor in self._tensor_loop(red=True): 1423 # Set the simulation number. 1424 tensor.set_sim_num(cdp.sim_number) 1425 1426 # Loop over the simulations. 1427 for j in range(cdp.sim_number): 1428 # Set the reduced tensor simulation data. 1429 tensor.set(param='Axx', value=sim_data[5*i + 0][j], category='sim', sim_index=j) 1430 tensor.set(param='Ayy', value=sim_data[5*i + 1][j], category='sim', sim_index=j) 1431 tensor.set(param='Axy', value=sim_data[5*i + 2][j], category='sim', sim_index=j) 1432 tensor.set(param='Axz', value=sim_data[5*i + 3][j], category='sim', sim_index=j) 1433 tensor.set(param='Ayz', value=sim_data[5*i + 4][j], category='sim', sim_index=j)
1434 1435
1436 - def sim_return_param(self, model_info, index):
1437 """Return the array of simulation parameter values. 1438 1439 @param model_info: The model information originating from model_loop() (unused). 1440 @type model_info: None 1441 @param index: The index of the parameter to return the array of values for. 1442 @type index: int 1443 """ 1444 1445 # Parameter increment counter. 1446 inc = 0 1447 1448 # Get the parameter object names. 1449 param_names = self.data_names(set='params') 1450 1451 # Loop over the parameters. 1452 for param in param_names: 1453 # Skip non-existent objects. 1454 if not hasattr(cdp, param): 1455 continue 1456 1457 # Return the parameter array. 1458 if index == inc: 1459 return getattr(cdp, param + "_sim") 1460 1461 # Increment. 1462 inc = inc + 1 1463 1464 # Add some additional parameters. 1465 if cdp.model == 'iso cone, free rotor' and inc == index: 1466 return getattr(cdp, 'cone_theta_sim')
1467 1468
1469 - def sim_return_selected(self, model_info):
1470 """Return the array of selected simulation flags for the spin. 1471 1472 @param model_info: The model information originating from model_loop() (unused). 1473 @type model_info: None 1474 @return: The array of selected simulation flags. 1475 @rtype: list of int 1476 """ 1477 1478 # Return the array. 1479 return cdp.select_sim
1480