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