specific_analyses.frame_order.optimisation

1 ############################################################################### 2 # # 3 # Copyright (C) 2009-2011,2013-2015 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 """Module for the optimisation of the frame order models.""" 24 25 # Python module imports. 26 from math import cos, pi 27 from minfx.generic import generic_minimise 28 from minfx.grid import grid_point_array 29 from numpy import arccos, array, dot, float64, ndarray, ones, zeros 30 from numpy.linalg import inv, norm 31 import sys 32 from warnings import warn 33 34 # relax module imports. 35 from lib.check_types import is_float 36 from lib.float import isNaN, isInf 37 from lib.errors import RelaxError, RelaxInfError, RelaxNaNError, RelaxNoPCSError, RelaxNoRDCError 38 from lib.frame_order.pseudo_ellipse import tmax_pseudo_ellipse_array 39 from lib.frame_order.variables import MODEL_DOUBLE_ROTOR, MODEL_FREE_ROTOR, MODEL_ISO_CONE, MODEL_ISO_CONE_FREE_ROTOR, MODEL_ISO_CONE_TORSIONLESS, MODEL_LIST_FREE_ROTORS, MODEL_LIST_PSEUDO_ELLIPSE, MODEL_PSEUDO_ELLIPSE, MODEL_PSEUDO_ELLIPSE_FREE_ROTOR, MODEL_PSEUDO_ELLIPSE_TORSIONLESS, MODEL_RIGID, MODEL_ROTOR 40 from lib.geometry.angles import wrap_angles 41 from lib.periodic_table import periodic_table 42 from lib.physical_constants import dipolar_constant 43 from lib.warnings import RelaxWarning 44 from multi import Memo, Result_command, Slave_command 45 from pipe_control.interatomic import interatomic_loop 46 from pipe_control.mol_res_spin import return_spin, spin_loop 47 from pipe_control.structure.mass import pipe_centre_of_mass 48 from specific_analyses.frame_order.checks import check_domain, check_model, check_parameters 49 from specific_analyses.frame_order.data import base_data_types, domain_moving, pivot_fixed, tensor_loop 50 from specific_analyses.frame_order.parameters import assemble_param_vector, linear_constraints 51 from target_functions.frame_order import Frame_order, sobol_data 52 53

54 -def count_sobol_points(target_fn=None, verbosity=1):

55 """Count the number of Sobol' points for the current parameter values of the model. 56 57 The count will be stored in the current data pipe and printed out. 58 59 60 @keyword target_fn: The pre-initialised frame order target function class. 61 @type target_fn: target_functions.frame_order.Frame_order instance 62 @keyword verbosity: If greater than 0, lots of information will be printed out. 63 @type verbosity: int 64 """ 65 66 # Printout. 67 if verbosity: 68 print("\nSobol' quasi-random integration point counting for the current parameter values.") 69 70 # Checks. 71 if not check_model(escalate=1): 72 return 73 if not check_parameters(escalate=1): 74 return 75 if not check_domain(escalate=1): 76 return 77 78 # Handle the rigid model. 79 if cdp.model == MODEL_RIGID: 80 if verbosity: 81 print("Sobol' quasi-random integration points are not used for the rigid frame order model.\n") 82 return 83 84 # Set up the target function, if required. 85 if target_fn == None: 86 # Set up the data structures for the target function. 87 param_vector, full_tensors, full_in_ref_frame, rdcs, rdc_err, rdc_weight, rdc_vect, rdc_const, pcs, pcs_err, pcs_weight, atomic_pos, temp, frq, paramag_centre, com, ave_pos_pivot, pivot, pivot_opt = target_fn_data_setup(verbosity=0, unset_fail=True) 88 89 # The numeric integration information. 90 if not hasattr(cdp, 'quad_int'): 91 cdp.quad_int = False 92 sobol_max_points, sobol_oversample = None, None 93 if hasattr(cdp, 'sobol_max_points'): 94 sobol_max_points = cdp.sobol_max_points 95 sobol_oversample = cdp.sobol_oversample 96 97 # Set up the optimisation target function class. 98 target_fn = Frame_order(model=cdp.model, init_params=param_vector, full_tensors=full_tensors, full_in_ref_frame=full_in_ref_frame, rdcs=rdcs, rdc_errors=rdc_err, rdc_weights=rdc_weight, rdc_vect=rdc_vect, dip_const=rdc_const, pcs=pcs, pcs_errors=pcs_err, pcs_weights=pcs_weight, atomic_pos=atomic_pos, temp=temp, frq=frq, paramag_centre=paramag_centre, scaling_matrix=None, com=com, ave_pos_pivot=ave_pos_pivot, pivot=pivot, pivot_opt=pivot_opt, sobol_max_points=sobol_max_points, sobol_oversample=sobol_oversample, quad_int=cdp.quad_int) 99 100 # The Sobol' sequence dimensions. 101 if cdp.model in [MODEL_ISO_CONE, MODEL_ISO_CONE_FREE_ROTOR, MODEL_PSEUDO_ELLIPSE, MODEL_PSEUDO_ELLIPSE_FREE_ROTOR]: 102 dims = ['theta', 'phi', 'sigma'] 103 elif cdp.model in [MODEL_ISO_CONE_TORSIONLESS, MODEL_PSEUDO_ELLIPSE_TORSIONLESS]: 104 dims = ['theta', 'phi'] 105 elif cdp.model in [MODEL_ROTOR, MODEL_FREE_ROTOR]: 106 dims = ['sigma'] 107 elif cdp.model in [MODEL_DOUBLE_ROTOR]: 108 dims = ['sigma', 'sigma2'] 109 110 # Unpack the points. 111 theta, phi, sigma, sigma2 = None, None, None, None 112 if dims == ['theta', 'phi', 'sigma']: 113 theta, phi, sigma = sobol_data.sobol_angles 114 elif dims == ['theta', 'phi']: 115 theta, phi = sobol_data.sobol_angles 116 elif dims == ['sigma']: 117 sigma = sobol_data.sobol_angles[0] 118 elif dims == ['sigma', 'sigma2']: 119 sigma, sigma2 = sobol_data.sobol_angles 120 121 # Pseudo-ellipse. 122 pe = False 123 if cdp.model in MODEL_LIST_PSEUDO_ELLIPSE: 124 pe = True 125 126 # Calculate theta_max. 127 theta_max = tmax_pseudo_ellipse_array(phi, cdp.cone_theta_x, cdp.cone_theta_y) 128 129 # The torsion angle. 130 if cdp.model in MODEL_LIST_FREE_ROTORS: 131 cone_sigma_max = pi 132 elif 'sigma' in dims: 133 cone_sigma_max = cdp.cone_sigma_max 134 135 # The isotropic cone angle. 136 if not pe and 'theta' in dims: 137 cone_theta = cdp.cone_theta 138 139 # Loop over the Sobol' points to count them. 140 total_num = len(sobol_data.sobol_angles[0]) 141 count = 0 142 for i in range(total_num): 143 # Pseudo-elliptic cone opening angle. 144 if pe and theta[i] > theta_max[i]: 145 continue 146 147 # Isotropic cones. 148 if not pe and 'theta' in dims and theta[i] > cone_theta: 149 continue 150 151 # 1st torsion angle. 152 if 'sigma' in dims and abs(sigma[i]) > cone_sigma_max: 153 continue 154 155 # 2nd torsion angle. 156 if 'sigma2' in dims and abs(sigma2[i]) > cdp.cone_sigma_max_2: 157 continue 158 159 # Increment the point count. 160 count += 1 161 162 # Maximum reached. 163 if count == cdp.sobol_max_points: 164 break 165 166 # Store the count. 167 cdp.sobol_points_used = count 168 169 # Printout. 170 if verbosity: 171 format = " %-30s %20s\n" 172 percent = "%s" % (float(count)/float(cdp.sobol_max_points)*100) + '%' 173 sys.stdout.write(format % ("Maximum number of points:", cdp.sobol_max_points)) 174 sys.stdout.write(format % ("Oversampling factor:", cdp.sobol_oversample)) 175 sys.stdout.write(format % ("Total points:", total_num)) 176 sys.stdout.write(format % ("Used points:", count)) 177 sys.stdout.write(format % ("Percentage:", percent)) 178 sys.stdout.write('\n')

179 180

181 -def grid_row(incs, lower, upper, dist_type=None, end_point=True):

182 """Set up a row of the grid search for a given parameter. 183 184 @param incs: The number of increments. 185 @type incs: int 186 @param lower: The lower bounds. 187 @type lower: float 188 @param upper: The upper bounds. 189 @type upper: float 190 @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). 191 @type dist_type: None or str 192 @keyword end_point: A flag which if False will cause the end point to be removed. 193 @type end_point: bool 194 @return: The row of the grid. 195 @rtype: list of float 196 """ 197 198 # Init. 199 row = [] 200 201 # Linear grid. 202 if dist_type == None: 203 # Handle a single increment. 204 if incs == 1: 205 row.append((lower + upper) / 2.0) 206 207 # Loop over the increments. 208 else: 209 for i in range(incs): 210 row.append(lower + i * (upper - lower) / (incs - 1.0)) 211 212 # Inverse cos distribution. 213 elif dist_type == 'acos': 214 # Handle a single increment. 215 if incs == 1: 216 row.append((lower + upper) / 2.0) 217 218 # Generate the increment values of v from cos(upper) to cos(lower). 219 else: 220 v = zeros(incs, float64) 221 val = (cos(lower) - cos(upper)) / (incs - 1.0) 222 for i in range(incs): 223 v[-i-1] = cos(upper) + float(i) * val 224 225 # Generate the distribution. 226 row = arccos(v) 227 228 # Remove the last point. 229 if incs != 1 and not end_point: 230 row = row[:-1] 231 232 # Return the row (as a list). 233 return list(row)

234 235

236 -def minimise_setup_atomic_pos(sim_index=None, verbosity=1):

237 """Set up the atomic position data structures for optimisation using PCSs and PREs as base data sets. 238 239 @keyword sim_index: The index of the simulation to optimise. This should be None if normal optimisation is desired. 240 @type sim_index: None or int 241 @keyword verbosity: If set to 1 or higher, then printouts and warnings will be active. 242 @type verbosity: int 243 @return: The atomic positions (the first index is the spins, the second is the structures, and the third is the atomic coordinates) and the paramagnetic centre. 244 @rtype: numpy rank-3 array, numpy rank-1 array. 245 """ 246 247 # Initialise. 248 atomic_pos = [] 249 250 # Store the atomic positions. 251 ave_warning_spin_ids = [] 252 ave_warning_num = None 253 for spin, spin_id in spin_loop(selection=domain_moving(), return_id=True): 254 # Skip deselected spins. 255 if not spin.select: 256 continue 257 258 # Only use spins with alignment/paramagnetic data. 259 if not hasattr(spin, 'pcs') and not hasattr(spin, 'pre'): 260 continue 261 262 # A single atomic position. 263 if (isinstance(spin.pos, list) or isinstance(spin.pos, ndarray)) and len(spin.pos) == 3 and is_float(spin.pos[0]): 264 atomic_pos.append(spin.pos) 265 266 # Average multiple atomic positions. 267 else: 268 # First throw a warning to tell the user what is happening. 269 if sim_index == None: 270 ave_warning_spin_ids.append(spin_id) 271 if ave_warning_num == None: 272 ave_warning_num = len(spin.pos) 273 elif ave_warning_num != len(spin.pos): 274 ave_warning_num = 'multiple' 275 276 # The average position. 277 ave_pos = zeros(3, float64) 278 count = 0 279 for i in range(len(spin.pos)): 280 if spin.pos[i] is None: 281 continue 282 ave_pos += spin.pos[i] 283 count += 1 284 ave_pos = ave_pos / float(count) 285 286 # Store. 287 atomic_pos.append(ave_pos) 288 289 # Give a warning about the atomic position averaging. 290 if verbosity and ave_warning_num != 1 and len(ave_warning_spin_ids): 291 warn(RelaxWarning("Averaging the %s atomic positions for the PCS for the spins %s." % (ave_warning_num, ave_warning_spin_ids))) 292 293 # Convert to numpy objects. 294 atomic_pos = array(atomic_pos, float64) 295 296 # The paramagnetic centre. 297 if not hasattr(cdp, 'paramagnetic_centre'): 298 paramag_centre = zeros(3, float64) 299 else: 300 paramag_centre = array(cdp.paramagnetic_centre) 301 302 # Return the data structures. 303 return atomic_pos, paramag_centre

304 305

306 -def minimise_setup_pcs(sim_index=None):

307 """Set up the data structures for optimisation using PCSs as base data sets. 308 309 @keyword sim_index: The index of the simulation to optimise. This should be None if normal optimisation is desired. 310 @type sim_index: None or int 311 @return: The assembled data structures for using PCSs as the base data for optimisation. These include: 312 - the PCS values. 313 - the unit vectors connecting the paramagnetic centre (the electron spin) to 314 - the PCS weight. 315 - the nuclear spin. 316 - the pseudocontact shift constants. 317 @rtype: tuple of (numpy rank-2 array, numpy rank-2 array, numpy rank-2 array, numpy rank-1 array, numpy rank-1 array) 318 """ 319 320 # Data setup tests. 321 if not hasattr(cdp, 'paramagnetic_centre'): 322 raise RelaxError("The paramagnetic centre has not yet been specified.") 323 if not hasattr(cdp, 'temperature'): 324 raise RelaxError("The experimental temperatures have not been set.") 325 if not hasattr(cdp, 'spectrometer_frq'): 326 raise RelaxError("The spectrometer frequencies of the experiments have not been set.") 327 328 # Initialise. 329 pcs = [] 330 pcs_err = [] 331 pcs_weight = [] 332 temp = [] 333 frq = [] 334 335 # The PCS data. 336 for i in range(len(cdp.align_ids)): 337 # Alias the ID. 338 align_id = cdp.align_ids[i] 339 340 # Skip non-optimised data. 341 if not opt_uses_align_data(align_id): 342 continue 343 344 # Append empty arrays to the PCS structures. 345 pcs.append([]) 346 pcs_err.append([]) 347 pcs_weight.append([]) 348 349 # Get the temperature for the PCS constant. 350 if align_id in cdp.temperature: 351 temp.append(cdp.temperature[align_id]) 352 353 # The temperature must be given! 354 else: 355 raise RelaxError("The experimental temperature for the alignment ID '%s' has not been set." % align_id) 356 357 # Get the spectrometer frequency in Tesla units for the PCS constant. 358 if align_id in cdp.spectrometer_frq: 359 frq.append(cdp.spectrometer_frq[align_id] * 2.0 * pi / periodic_table.gyromagnetic_ratio('1H')) 360 361 # The frequency must be given! 362 else: 363 raise RelaxError("The spectrometer frequency for the alignment ID '%s' has not been set." % align_id) 364 365 # Spin loop over the domain. 366 j = 0 367 for spin in spin_loop(selection=domain_moving()): 368 # Skip deselected spins. 369 if not spin.select: 370 continue 371 372 # Skip spins without PCS data. 373 if not hasattr(spin, 'pcs'): 374 continue 375 376 # Append the PCSs to the list. 377 if align_id in spin.pcs: 378 if sim_index != None: 379 pcs[-1].append(spin.pcs_sim[align_id][sim_index]) 380 else: 381 pcs[-1].append(spin.pcs[align_id]) 382 else: 383 pcs[-1].append(None) 384 385 # Append the PCS errors. 386 if hasattr(spin, 'pcs_err') and align_id in spin.pcs_err: 387 pcs_err[-1].append(spin.pcs_err[align_id]) 388 else: 389 pcs_err[-1].append(None) 390 391 # Append the weight. 392 if hasattr(spin, 'pcs_weight') and align_id in spin.pcs_weight: 393 pcs_weight[-1].append(spin.pcs_weight[align_id]) 394 else: 395 pcs_weight[-1].append(1.0) 396 397 # Spin index. 398 j = j + 1 399 400 # Convert to numpy objects. 401 pcs = array(pcs, float64) 402 pcs_err = array(pcs_err, float64) 403 pcs_weight = array(pcs_weight, float64) 404 405 # Convert the PCS from ppm to no units. 406 pcs = pcs * 1e-6 407 pcs_err = pcs_err * 1e-6 408 409 # Return the data structures. 410 return pcs, pcs_err, pcs_weight, temp, frq

411 412

413 -def minimise_setup_rdcs(sim_index=None):

414 """Set up the data structures for optimisation using RDCs as base data sets. 415 416 @keyword sim_index: The index of the simulation to optimise. This should be None if normal optimisation is desired. 417 @type sim_index: None or int 418 @return: The assembled data structures for using RDCs as the base data for optimisation. These include: 419 - rdc, the RDC values. 420 - rdc_err, the RDC errors. 421 - rdc_weight, the RDC weights. 422 - vectors, the interatomic vectors. 423 - rdc_const, the dipolar constants. 424 - absolute, the absolute value flags (as 1's and 0's). 425 @rtype: tuple of (numpy rank-2 array, numpy rank-2 array, numpy rank-2 array, numpy rank-3 array, numpy rank-2 array, numpy rank-2 array) 426 """ 427 428 # Initialise. 429 rdc = [] 430 rdc_err = [] 431 rdc_weight = [] 432 unit_vect = [] 433 rdc_const = [] 434 absolute = [] 435 436 # The unit vectors and RDC constants. 437 for interatom in interatomic_loop(selection1=domain_moving()): 438 # Get the spins. 439 spin1 = return_spin(spin_hash=interatom._spin_hash1) 440 spin2 = return_spin(spin_hash=interatom._spin_hash2) 441 442 # A single unit vector. 443 if interatom.vector.shape == (3,): 444 unit_vect.append(interatom.vector) 445 446 # Average multiple unit vectors. 447 else: 448 # First throw a warning to tell the user what is happening. 449 if sim_index == None: 450 warn(RelaxWarning("Averaging the %s unit vectors for the RDC for the spin pair '%s' and '%s'." % (len(interatom.vector), interatom.spin_id1, interatom.spin_id2))) 451 452 # The average position. 453 ave_vector = zeros(3, float64) 454 for i in range(len(interatom.vector)): 455 ave_vector += interatom.vector[i] 456 457 # Store. 458 unit_vect.append(ave_vector) 459 460 # Normalise (to be sure). 461 unit_vect[-1] = unit_vect[-1] / norm(unit_vect[-1]) 462 463 # Gyromagnetic ratios. 464 g1 = periodic_table.gyromagnetic_ratio(spin1.isotope) 465 g2 = periodic_table.gyromagnetic_ratio(spin2.isotope) 466 467 # Calculate the RDC dipolar constant (in Hertz, and the 3 comes from the alignment tensor), and append it to the list. 468 rdc_const.append(3.0/(2.0*pi) * dipolar_constant(g1, g2, interatom.r)) 469 470 # Fix the unit vector data structure. 471 num = None 472 for rdc_index in range(len(unit_vect)): 473 # Number of vectors. 474 if num == None: 475 if unit_vect[rdc_index] is not None: 476 num = len(unit_vect[rdc_index]) 477 continue 478 479 # Check. 480 if unit_vect[rdc_index] is not None and len(unit_vect[rdc_index]) != num: 481 raise RelaxError("The number of interatomic vectors for all no match:\n%s" % unit_vect) 482 483 # Missing unit vectors. 484 if num == None: 485 raise RelaxError("No interatomic vectors could be found.") 486 487 # Update None entries. 488 for i in range(len(unit_vect)): 489 if unit_vect[i] is None: 490 unit_vect[i] = [[None, None, None]]*num 491 492 # The RDC data. 493 for i in range(len(cdp.align_ids)): 494 # Alias the ID. 495 align_id = cdp.align_ids[i] 496 497 # Skip non-optimised data. 498 if not opt_uses_align_data(align_id): 499 continue 500 501 # Append empty arrays to the RDC structures. 502 rdc.append([]) 503 rdc_err.append([]) 504 rdc_weight.append([]) 505 absolute.append([]) 506 507 # Interatom loop over the domain. 508 for interatom in interatomic_loop(domain_moving(), skip_desel=True): 509 # Get the spins. 510 spin1 = return_spin(spin_hash=interatom._spin_hash1) 511 spin2 = return_spin(spin_hash=interatom._spin_hash2) 512 513 # Only use interatomic data containers with RDC and vector data. 514 if not hasattr(interatom, 'rdc') or not hasattr(interatom, 'vector'): 515 continue 516 517 # Defaults of None. 518 value = None 519 error = None 520 521 # Pseudo-atom set up. 522 if (hasattr(spin1, 'members') or hasattr(spin2, 'members')) and align_id in interatom.rdc: 523 raise RelaxError("Psuedo-atoms are currently not supported for the frame order analysis.") 524 525 # Normal set up. 526 elif align_id in interatom.rdc: 527 # The RDC. 528 if sim_index != None: 529 value = interatom.rdc_sim[align_id][sim_index] 530 else: 531 value = interatom.rdc[align_id] 532 533 # The error. 534 if hasattr(interatom, 'rdc_err') and align_id in interatom.rdc_err: 535 error = interatom.rdc_err[align_id] 536 537 # Append the RDCs to the list. 538 rdc[-1].append(value) 539 540 # Append the RDC errors. 541 rdc_err[-1].append(error) 542 543 # Append the weight. 544 if hasattr(interatom, 'rdc_weight') and align_id in interatom.rdc_weight: 545 rdc_weight[-1].append(interatom.rdc_weight[align_id]) 546 else: 547 rdc_weight[-1].append(1.0) 548 549 # Append the absolute value flag. 550 if hasattr(interatom, 'absolute_rdc') and align_id in interatom.absolute_rdc: 551 absolute[-1].append(interatom.absolute_rdc[align_id]) 552 else: 553 absolute[-1].append(False) 554 555 # Convert to numpy objects. 556 rdc = array(rdc, float64) 557 rdc_err = array(rdc_err, float64) 558 rdc_weight = array(rdc_weight, float64) 559 unit_vect = array(unit_vect, float64) 560 rdc_const = array(rdc_const, float64) 561 absolute = array(absolute, float64) 562 563 # Return the data structures. 564 return rdc, rdc_err, rdc_weight, unit_vect, rdc_const, absolute

565 566

567 -def minimise_setup_tensors(sim_index=None):

568 """Set up the data structures for optimisation using alignment tensors as base data sets. 569 570 @keyword sim_index: The simulation index. This should be None if normal optimisation is desired. 571 @type sim_index: None or int 572 @return: The assembled data structures for using alignment tensors as the base data for optimisation. These include: 573 - full_tensors, the full tensors as concatenated arrays. 574 - full_err, the full tensor errors as concatenated arrays. 575 - full_in_ref_frame, the flags specifying if the tensor is the full or reduced tensor in the non-moving reference domain. 576 @rtype: tuple of 3 numpy nx5D, rank-1 arrays 577 """ 578 579 # Checks. 580 if not hasattr(cdp, 'ref_domain'): 581 raise RelaxError("The reference non-moving domain has not been specified.") 582 if not hasattr(cdp.align_tensors, 'reduction'): 583 raise RelaxError("The tensor reductions have not been specified.") 584 for i, tensor in tensor_loop(): 585 if not hasattr(tensor, 'domain'): 586 raise RelaxError("The domain that the '%s' tensor is attached to has not been set" % tensor.name) 587 588 # Initialise. 589 n = len(cdp.align_tensors.reduction) 590 full_tensors = zeros(n*5, float64) 591 full_err = ones(n*5, float64) * 1e-5 592 full_in_ref_frame = zeros(n, float64) 593 594 # Loop over the full tensors. 595 for i, tensor in tensor_loop(red=False): 596 # The full tensor. 597 full_tensors[5*i + 0] = tensor.Axx 598 full_tensors[5*i + 1] = tensor.Ayy 599 full_tensors[5*i + 2] = tensor.Axy 600 full_tensors[5*i + 3] = tensor.Axz 601 full_tensors[5*i + 4] = tensor.Ayz 602 603 # The full tensor corresponds to the frame of reference. 604 if cdp.ref_domain == tensor.domain: 605 full_in_ref_frame[i] = 1 606 607 # The full tensor errors. 608 if hasattr(tensor, 'Axx_err'): 609 full_err[5*i + 0] = tensor.Axx_err 610 full_err[5*i + 1] = tensor.Ayy_err 611 full_err[5*i + 2] = tensor.Axy_err 612 full_err[5*i + 3] = tensor.Axz_err 613 full_err[5*i + 4] = tensor.Ayz_err 614 615 # Return the data structures. 616 return full_tensors, full_err, full_in_ref_frame

617 618

619 -def opt_uses_align_data(align_id=None):

620 """Determine if the PCS or RDC data for the given alignment ID is needed for optimisation. 621 622 @keyword align_id: The optional alignment ID string. 623 @type align_id: str 624 @return: True if alignment data is to be used for optimisation, False otherwise. 625 @rtype: bool 626 """ 627 628 # No alignment IDs. 629 if not hasattr(cdp, 'align_ids'): 630 return False 631 632 # Convert the align IDs to an array, or take all IDs. 633 if align_id: 634 align_ids = [align_id] 635 else: 636 align_ids = cdp.align_ids 637 638 # Check the PCS and RDC. 639 for align_id in align_ids: 640 if opt_uses_pcs(align_id) or opt_uses_rdc(align_id): 641 return True 642 643 # No alignment data is used for optimisation. 644 return False

645 646

647 -def opt_uses_pcs(align_id):

648 """Determine if the PCS data for the given alignment ID is needed for optimisation. 649 650 @param align_id: The alignment ID string. 651 @type align_id: str 652 @return: True if the PCS data is to be used for optimisation, False otherwise. 653 @rtype: bool 654 """ 655 656 # No alignment IDs. 657 if not hasattr(cdp, 'pcs_ids'): 658 return False 659 660 # No PCS data for the alignment. 661 if align_id not in cdp.pcs_ids: 662 return False 663 664 # The PCS data is to be used for optimisation. 665 return True

666 667

668 -def opt_uses_rdc(align_id):

669 """Determine if the RDC data for the given alignment ID is needed for optimisation. 670 671 @param align_id: The alignment ID string. 672 @type align_id: str 673 @return: True if the RDC data is to be used for optimisation, False otherwise. 674 @rtype: bool 675 """ 676 677 # No alignment IDs. 678 if not hasattr(cdp, 'rdc_ids'): 679 return False 680 681 # No RDC data for the alignment. 682 if align_id not in cdp.rdc_ids: 683 return False 684 685 # The RDC data is to be used for optimisation. 686 return True

687 688

689 -def store_bc_data(A_5D_bc=None, pcs_theta=None, rdc_theta=None):

690 """Store the back-calculated data. 691 692 @keyword A_5D_bc: The reduced back-calculated alignment tensors from the target function. 693 @type A_5D_bc: numpy float64 array 694 @keyword pcs_theta: The back calculated PCS values from the target function. 695 @type pcs_theta: numpy float64 array 696 @keyword rdc_theta: The back calculated RDC values from the target function. 697 @type rdc_theta: numpy float64 array 698 """ 699 700 # Loop over the reduced tensors. 701 for i, tensor in tensor_loop(red=True): 702 # Store the values. 703 tensor.set(param='Axx', value=A_5D_bc[5*i + 0]) 704 tensor.set(param='Ayy', value=A_5D_bc[5*i + 1]) 705 tensor.set(param='Axy', value=A_5D_bc[5*i + 2]) 706 tensor.set(param='Axz', value=A_5D_bc[5*i + 3]) 707 tensor.set(param='Ayz', value=A_5D_bc[5*i + 4]) 708 709 # The RDC data. 710 for i in range(len(cdp.align_ids)): 711 # The alignment ID. 712 align_id = cdp.align_ids[i] 713 714 # Data flags 715 rdc_flag = False 716 if hasattr(cdp, 'rdc_ids') and align_id in cdp.rdc_ids: 717 rdc_flag = True 718 pcs_flag = False 719 if hasattr(cdp, 'pcs_ids') and align_id in cdp.pcs_ids: 720 pcs_flag = True 721 722 # Spin loop over the domain. 723 pcs_index = 0 724 for spin in spin_loop(domain_moving()): 725 # Skip deselected spins. 726 if not spin.select: 727 continue 728 729 # Spins with PCS data. 730 if pcs_flag and hasattr(spin, 'pcs'): 731 # Initialise the data structure. 732 if not hasattr(spin, 'pcs_bc'): 733 spin.pcs_bc = {} 734 735 # Store the back-calculated value (in ppm). 736 spin.pcs_bc[align_id] = pcs_theta[i, pcs_index] * 1e6 737 738 # Increment the index. 739 pcs_index += 1 740 741 # Interatomic data container loop. 742 rdc_index = 0 743 for interatom in interatomic_loop(domain_moving(), skip_desel=True): 744 # Initialise the data structure. 745 if not hasattr(interatom, 'rdc_bc'): 746 interatom.rdc_bc = {} 747 748 # Store the back-calculated value. 749 interatom.rdc_bc[align_id] = rdc_theta[i, rdc_index] 750 751 # Increment the index. 752 rdc_index += 1

753 754

755 -def target_fn_data_setup(sim_index=None, verbosity=1, scaling_matrix=None, unset_fail=False):

756 """Initialise the target function for optimisation or direct calculation. 757 758 @keyword sim_index: The index of the simulation to optimise. This should be None if normal optimisation is desired. 759 @type sim_index: None or int 760 @keyword verbosity: The amount of information to print. The higher the value, the greater the verbosity. 761 @type verbosity: int 762 @keyword scaling_matrix: The diagonal and square scaling matrices. 763 @type scaling_matrix: numpy rank-2, float64 array or None 764 @keyword unset_fail: A flag which if True will cause a RelaxError to be raised if the parameter is not set yet. 765 @type unset_fail: bool 766 """ 767 768 # Model printout. 769 if verbosity: 770 print("Frame order model: %s" % cdp.model) 771 772 # Assemble the parameter vector. 773 param_vector = assemble_param_vector(sim_index=sim_index, unset_fail=unset_fail) 774 775 # Determine the base data types (RDCs and/or PCSs). 776 data_types = base_data_types() 777 778 # Diagonal scaling. 779 if scaling_matrix != None and len(param_vector): 780 param_vector = dot(inv(scaling_matrix), param_vector) 781 782 # Get the data structures for optimisation using the tensors as base data sets. 783 full_tensors, full_tensor_err, full_in_ref_frame = minimise_setup_tensors(sim_index) 784 785 # Get the data structures for optimisation using PCSs as base data sets. 786 pcs, pcs_err, pcs_weight, temp, frq = None, None, None, None, None 787 if 'pcs' in data_types: 788 pcs, pcs_err, pcs_weight, temp, frq = minimise_setup_pcs(sim_index=sim_index) 789 790 # Get the data structures for optimisation using RDCs as base data sets. 791 rdcs, rdc_err, rdc_weight, rdc_vect, rdc_const, absolute_rdc = None, None, None, None, None, None 792 if 'rdc' in data_types: 793 rdcs, rdc_err, rdc_weight, rdc_vect, rdc_const, absolute_rdc = minimise_setup_rdcs(sim_index=sim_index) 794 795 # Data checks. 796 if pcs is not None and not len(pcs): 797 raise RelaxNoPCSError 798 if rdcs is not None and not len(rdcs): 799 raise RelaxNoRDCError 800 801 # Get the atomic_positions. 802 atomic_pos, paramag_centre = None, None 803 if 'pcs' in data_types or 'pre' in data_types: 804 atomic_pos, paramag_centre = minimise_setup_atomic_pos(sim_index=sim_index, verbosity=verbosity) 805 806 # The fixed pivot point. 807 pivot = None 808 if hasattr(cdp, 'pivot_x'): 809 pivot = array([cdp.pivot_x, cdp.pivot_y, cdp.pivot_z]) 810 811 # Pivot optimisation. 812 pivot_opt = True 813 if pivot_fixed(): 814 pivot_opt = False 815 816 # The centre of mass of the moving domain - to use as the centroid for the average domain position rotation. 817 ave_pos_pivot = pipe_centre_of_mass(atom_id=domain_moving(), verbosity=0) 818 819 # The centre of mass, for use in the rotor models. 820 com = None 821 if cdp.model in [MODEL_ROTOR, MODEL_FREE_ROTOR, MODEL_DOUBLE_ROTOR]: 822 # The centre of mass of all objects in the data pipe. 823 com = pipe_centre_of_mass(verbosity=0) 824 com = array(com, float64) 825 826 # Information printout. 827 if verbosity and sim_index == None: 828 sys.stdout.write("\nThe average domain rotation centroid, taken as the CoM of the atoms defined as the moving domain, is:\n %s\n" % list(ave_pos_pivot)) 829 if com is not None: 830 sys.stdout.write("The centre of mass reference coordinate for the rotor models is:\n %s\n" % list(com)) 831 if cdp.model != MODEL_RIGID: 832 if hasattr(cdp, 'quad_int') and cdp.quad_int: 833 sys.stdout.write("Numerical PCS integration: SciPy quadratic integration.\n") 834 else: 835 sys.stdout.write("Numerical PCS integration: Quasi-random Sobol' sequence.\n") 836 base_data = [] 837 if rdcs is not None and len(rdcs): 838 base_data.append("RDCs") 839 if pcs is not None and len(pcs): 840 base_data.append("PCSs") 841 sys.stdout.write("Base data: %s\n" % repr(base_data)) 842 sys.stdout.write("\n") 843 844 # Return the data. 845 return param_vector, full_tensors, full_in_ref_frame, rdcs, rdc_err, rdc_weight, rdc_vect, rdc_const, pcs, pcs_err, pcs_weight, atomic_pos, temp, frq, paramag_centre, com, ave_pos_pivot, pivot, pivot_opt

846 847

848 -def unpack_opt_results(param_vector=None, func=None, iter_count=None, f_count=None, g_count=None, h_count=None, warning=None, scaling_matrix=None, sim_index=None):

849 """Unpack and store the Frame Order optimisation results. 850 851 @keyword param_vector: The model-free parameter vector. 852 @type param_vector: numpy array 853 @keyword func: The optimised chi-squared value. 854 @type func: float 855 @keyword iter_count: The number of optimisation steps required to find the minimum. 856 @type iter_count: int 857 @keyword f_count: The function count. 858 @type f_count: int 859 @keyword g_count: The gradient count. 860 @type g_count: int 861 @keyword h_count: The Hessian count. 862 @type h_count: int 863 @keyword warning: Any optimisation warnings. 864 @type warning: str or None 865 @keyword scaling_matrix: The diagonal and square scaling matrices. 866 @type scaling_matrix: numpy rank-2, float64 array or None 867 @keyword sim_index: The index of the simulation to optimise. This should be None for normal optimisation. 868 @type sim_index: None or int 869 """ 870 871 # Catch infinite chi-squared values. 872 if isInf(func): 873 raise RelaxInfError('chi-squared') 874 875 # Catch chi-squared values of NaN. 876 if isNaN(func): 877 raise RelaxNaNError('chi-squared') 878 879 # Scaling. 880 if scaling_matrix is not None: 881 param_vector = dot(scaling_matrix, param_vector) 882 883 # The parameters to wrap. 884 wrap = [ 885 'ave_pos_alpha', 886 'ave_pos_beta', 887 'ave_pos_gamma', 888 'eigen_alpha', 889 'eigen_beta', 890 'eigen_gamma', 891 'axis_theta', 892 'axis_phi' 893 ] 894 895 # Monte Carlo simulation data structures. 896 if sim_index != None: 897 # Loop over the parameters. 898 for i in range(len(cdp.params)): 899 # Angle wrapping around the real value. 900 if cdp.params[i] in wrap or cdp.params[i] == 'axis_alpha': 901 val = getattr(cdp, cdp.params[i]) 902 param_vector[i] = wrap_angles(param_vector[i], val-pi, val+pi) 903 904 # Store the value. 905 obj = getattr(cdp, cdp.params[i]+'_sim') 906 obj[sim_index] = param_vector[i] 907 908 # Optimisation stats. 909 cdp.chi2_sim[sim_index] = func 910 cdp.iter_sim[sim_index] = iter_count 911 cdp.f_count_sim[sim_index] = f_count 912 cdp.g_count_sim[sim_index] = g_count 913 cdp.h_count_sim[sim_index] = h_count 914 cdp.warning_sim[sim_index] = warning 915 916 # Normal data structures. 917 else: 918 # Loop over the parameters. 919 for i in range(len(cdp.params)): 920 # Angle wrapping. 921 if cdp.params[i] in wrap: 922 param_vector[i] = wrap_angles(param_vector[i], 0.0, 2.0*pi) 923 if cdp.params[i] == 'axis_alpha': 924 param_vector[i] = wrap_angles(param_vector[i], -pi, pi) 925 926 # Store the value. 927 setattr(cdp, cdp.params[i], param_vector[i]) 928 929 # Optimisation stats. 930 cdp.chi2 = func 931 cdp.iter = iter_count 932 cdp.f_count = f_count 933 cdp.g_count = g_count 934 cdp.h_count = h_count 935 cdp.warning = warning

936 937 938

939 -class Frame_order_grid_command(Slave_command):

940 """Command class for relaxation dispersion optimisation on the slave processor.""" 941

942 - def __init__(self, points=None, scaling_matrix=None, sim_index=None, model=None, param_vector=None, full_tensors=None, full_in_ref_frame=None, rdcs=None, rdc_err=None, rdc_weight=None, rdc_vect=None, rdc_const=None, pcs=None, pcs_err=None, pcs_weight=None, atomic_pos=None, temp=None, frq=None, paramag_centre=None, com=None, ave_pos_pivot=None, pivot=None, pivot_opt=None, sobol_max_points=None, sobol_oversample=None, verbosity=None, quad_int=False):

943 """Initialise the base class, storing all the master data to be sent to the slave processor. 944 945 This method is run on the master processor whereas the run() method is run on the slave processor. 946 947 @keyword points: The points of the grid search subdivision to search over. 948 @type points: numpy rank-2 array 949 @keyword scaling_matrix: The diagonal, square scaling matrix. 950 @type scaling_matrix: numpy diagonal matrix 951 @keyword sim_index: The index of the simulation to optimise. This should be None if normal optimisation is desired. 952 @type sim_index: None or int 953 @keyword model: The name of the Frame Order model. 954 @type model: str 955 @keyword param_vector: The initial parameter values. 956 @type param_vector: numpy float64 array 957 @keyword full_tensors: An array of the {Axx, Ayy, Axy, Axz, Ayz} values for all full alignment tensors. The format is [Axx1, Ayy1, Axy1, Axz1, Ayz1, Axx2, Ayy2, Axy2, Axz2, Ayz2, ..., Axxn, Ayyn, Axyn, Axzn, Ayzn]. 958 @type full_tensors: numpy nx5D, rank-1 float64 array 959 @keyword full_in_ref_frame: An array of flags specifying if the tensor in the reference frame is the full or reduced tensor. 960 @type full_in_ref_frame: numpy rank-1 array 961 @keyword rdcs: The RDC lists. The first index must correspond to the different alignment media i and the second index to the spin systems j. 962 @type rdcs: numpy rank-2 array 963 @keyword rdc_err: The RDC error lists. The dimensions of this argument are the same as for 'rdcs'. 964 @type rdc_err: numpy rank-2 array 965 @keyword rdc_weight: The RDC weight lists. The dimensions of this argument are the same as for 'rdcs'. 966 @type rdc_weight: numpy rank-2 array 967 @keyword rdc_vect: The unit XH vector lists corresponding to the RDC values. The first index must correspond to the spin systems and the second index to the x, y, z elements. 968 @type rdc_vect: numpy rank-2 array 969 @keyword rdc_const: The dipolar constants for each RDC. The indices correspond to the spin systems j. 970 @type rdc_const: numpy rank-1 array 971 @keyword pcs: The PCS lists. The first index must correspond to the different alignment media i and the second index to the spin systems j. 972 @type pcs: numpy rank-2 array 973 @keyword pcs_err: The PCS error lists. The dimensions of this argument are the same as for 'pcs'. 974 @type pcs_err: numpy rank-2 array 975 @keyword pcs_weight: The PCS weight lists. The dimensions of this argument are the same as for 'pcs'. 976 @type pcs_weight: numpy rank-2 array 977 @keyword atomic_pos: The atomic positions of all spins for the PCS and PRE data. The first index is the spin systems j and the second is the structure or state c. 978 @type atomic_pos: numpy rank-3 array 979 @keyword temp: The temperature of each PCS data set. 980 @type temp: numpy rank-1 array 981 @keyword frq: The frequency of each PCS data set. 982 @type frq: numpy rank-1 array 983 @keyword paramag_centre: The paramagnetic centre position (or positions). 984 @type paramag_centre: numpy rank-1, 3D array or rank-2, Nx3 array 985 @keyword com: The centre of mass of the system. This is used for defining the rotor model systems. 986 @type com: numpy 3D rank-1 array 987 @keyword ave_pos_pivot: The pivot point to rotate all atoms about to the average domain position. In most cases this will be the centre of mass of the moving domain. This pivot is shifted by the translation vector. 988 @type ave_pos_pivot: numpy 3D rank-1 array 989 @keyword pivot: The pivot point for the ball-and-socket joint motion. This is needed if PCS or PRE values are used. 990 @type pivot: numpy rank-1, 3D array or None 991 @keyword pivot_opt: A flag which if True will allow the pivot point of the motion to be optimised. 992 @type pivot_opt: bool 993 @keyword sobol_max_points: The maximum number of Sobol' points to use for the numerical PCS integration technique. 994 @type sobol_max_points: int 995 @keyword sobol_oversample: The oversampling factor Ov used for the total number of points N * Ov * 10**M, where N is the maximum number of Sobol' points and M is the number of dimensions or torsion-tilt angles for the system. 996 @type sobol_oversample: int 997 @keyword verbosity: The verbosity level. This is used by the result command returned to the master for printouts. 998 @type verbosity: int 999 @keyword quad_int: A flag which if True will perform high precision numerical integration via the scipy.integrate quad(), dblquad() and tplquad() integration methods rather than the rough quasi-random numerical integration. 1000 @type quad_int: bool 1001 """ 1002 1003 # Store the arguments. 1004 self.points = points 1005 self.scaling_matrix = scaling_matrix 1006 self.model = model 1007 self.param_vector = param_vector 1008 self.full_tensors = full_tensors 1009 self.full_in_ref_frame = full_in_ref_frame 1010 self.rdcs = rdcs 1011 self.rdc_err = rdc_err 1012 self.rdc_weight = rdc_weight 1013 self.rdc_vect = rdc_vect 1014 self.rdc_const = rdc_const 1015 self.pcs = pcs 1016 self.pcs_err = pcs_err 1017 self.pcs_weight = pcs_weight 1018 self.atomic_pos = atomic_pos 1019 self.temp = temp 1020 self.frq = frq 1021 self.paramag_centre = paramag_centre 1022 self.com = com 1023 self.ave_pos_pivot = ave_pos_pivot 1024 self.pivot = pivot 1025 self.pivot_opt = pivot_opt 1026 self.sobol_max_points = sobol_max_points 1027 self.sobol_oversample = sobol_oversample 1028 self.verbosity = verbosity 1029 self.quad_int = quad_int

1030 1031

1032 - def run(self, processor, completed):

1033 """Set up and perform the optimisation.""" 1034 1035 # Set up the optimisation target function class. 1036 target_fn = Frame_order(model=self.model, init_params=self.param_vector, full_tensors=self.full_tensors, full_in_ref_frame=self.full_in_ref_frame, rdcs=self.rdcs, rdc_errors=self.rdc_err, rdc_weights=self.rdc_weight, rdc_vect=self.rdc_vect, dip_const=self.rdc_const, pcs=self.pcs, pcs_errors=self.pcs_err, pcs_weights=self.pcs_weight, atomic_pos=self.atomic_pos, temp=self.temp, frq=self.frq, paramag_centre=self.paramag_centre, scaling_matrix=self.scaling_matrix, com=self.com, ave_pos_pivot=self.ave_pos_pivot, pivot=self.pivot, pivot_opt=self.pivot_opt, sobol_max_points=self.sobol_max_points, sobol_oversample=self.sobol_oversample, quad_int=self.quad_int) 1037 1038 # Grid search. 1039 results = grid_point_array(func=target_fn.func, args=(), points=self.points, verbosity=self.verbosity) 1040 1041 # Create the result command object on the slave to send back to the master. 1042 processor.return_object(Frame_order_result_command(processor=processor, memo_id=self.memo_id, results=results, A_5D_bc=target_fn.A_5D_bc, pcs_theta=target_fn.pcs_theta, rdc_theta=target_fn.rdc_theta, completed=completed))

1043 1044 1045

1046 -class Frame_order_memo(Memo):

1047 """The frame order memo class.""" 1048

1049 - def __init__(self, spins=None, spin_ids=None, sim_index=None, scaling_matrix=None, verbosity=None):

1050 """Initialise the relaxation dispersion memo class. 1051 1052 This is used for handling the optimisation results returned from a slave processor. It runs on the master processor and is used to store data which is passed to the slave processor and then passed back to the master via the results command. 1053 1054 1055 @keyword spins: The list of spin data container for the cluster. If this argument is supplied, then the spin_id argument will be ignored. 1056 @type spins: list of SpinContainer instances 1057 @keyword spin_ids: The spin ID strings for the cluster. 1058 @type spin_ids: list of str 1059 @keyword sim_index: The optional MC simulation index. 1060 @type sim_index: int 1061 @keyword scaling_matrix: The diagonal, square scaling matrix. 1062 @type scaling_matrix: numpy diagonal matrix 1063 @keyword verbosity: The verbosity level. This is used by the result command returned to the master for printouts. 1064 @type verbosity: int 1065 """ 1066 1067 # Execute the base class __init__() method. 1068 super(Frame_order_memo, self).__init__() 1069 1070 # Store the arguments. 1071 self.spins = spins 1072 self.spin_ids = spin_ids 1073 self.sim_index = sim_index 1074 self.scaling_matrix = scaling_matrix

1075 1076 1077

1078 -class Frame_order_minimise_command(Slave_command):

1079 """Command class for relaxation dispersion optimisation on the slave processor.""" 1080

1081 - def __init__(self, min_algor=None, min_options=None, func_tol=None, grad_tol=None, max_iterations=None, scaling_matrix=None, constraints=False, sim_index=None, model=None, param_vector=None, full_tensors=None, full_in_ref_frame=None, rdcs=None, rdc_err=None, rdc_weight=None, rdc_vect=None, rdc_const=None, pcs=None, pcs_err=None, pcs_weight=None, atomic_pos=None, temp=None, frq=None, paramag_centre=None, com=None, ave_pos_pivot=None, pivot=None, pivot_opt=None, sobol_max_points=None, sobol_oversample=None, verbosity=None, quad_int=False):

1082 """Initialise the base class, storing all the master data to be sent to the slave processor. 1083 1084 This method is run on the master processor whereas the run() method is run on the slave processor. 1085 1086 @keyword min_algor: The minimisation algorithm to use. 1087 @type min_algor: str 1088 @keyword min_options: An array of options to be used by the minimisation algorithm. 1089 @type min_options: array of str 1090 @keyword func_tol: The function tolerance which, when reached, terminates optimisation. Setting this to None turns of the check. 1091 @type func_tol: None or float 1092 @keyword grad_tol: The gradient tolerance which, when reached, terminates optimisation. Setting this to None turns of the check. 1093 @type grad_tol: None or float 1094 @keyword max_iterations: The maximum number of iterations for the algorithm. 1095 @type max_iterations: int 1096 @keyword constraints: If True, constraints are used during optimisation. 1097 @type constraints: bool 1098 @keyword sim_index: The index of the simulation to optimise. This should be None if normal optimisation is desired. 1099 @type sim_index: None or int 1100 @keyword model: The name of the Frame Order model. 1101 @type model: str 1102 @keyword param_vector: The initial parameter values. 1103 @type param_vector: numpy float64 array 1104 @keyword full_tensors: An array of the {Axx, Ayy, Axy, Axz, Ayz} values for all full alignment tensors. The format is [Axx1, Ayy1, Axy1, Axz1, Ayz1, Axx2, Ayy2, Axy2, Axz2, Ayz2, ..., Axxn, Ayyn, Axyn, Axzn, Ayzn]. 1105 @type full_tensors: numpy nx5D, rank-1 float64 array 1106 @keyword full_in_ref_frame: An array of flags specifying if the tensor in the reference frame is the full or reduced tensor. 1107 @type full_in_ref_frame: numpy rank-1 array 1108 @keyword rdcs: The RDC lists. The first index must correspond to the different alignment media i and the second index to the spin systems j. 1109 @type rdcs: numpy rank-2 array 1110 @keyword rdc_err: The RDC error lists. The dimensions of this argument are the same as for 'rdcs'. 1111 @type rdc_err: numpy rank-2 array 1112 @keyword rdc_weight: The RDC weight lists. The dimensions of this argument are the same as for 'rdcs'. 1113 @type rdc_weight: numpy rank-2 array 1114 @keyword rdc_vect: The unit XH vector lists corresponding to the RDC values. The first index must correspond to the spin systems and the second index to the x, y, z elements. 1115 @type rdc_vect: numpy rank-2 array 1116 @keyword rdc_const: The dipolar constants for each RDC. The indices correspond to the spin systems j. 1117 @type rdc_const: numpy rank-1 array 1118 @keyword pcs: The PCS lists. The first index must correspond to the different alignment media i and the second index to the spin systems j. 1119 @type pcs: numpy rank-2 array 1120 @keyword pcs_err: The PCS error lists. The dimensions of this argument are the same as for 'pcs'. 1121 @type pcs_err: numpy rank-2 array 1122 @keyword pcs_weight: The PCS weight lists. The dimensions of this argument are the same as for 'pcs'. 1123 @type pcs_weight: numpy rank-2 array 1124 @keyword atomic_pos: The atomic positions of all spins for the PCS and PRE data. The first index is the spin systems j and the second is the structure or state c. 1125 @type atomic_pos: numpy rank-3 array 1126 @keyword temp: The temperature of each PCS data set. 1127 @type temp: numpy rank-1 array 1128 @keyword frq: The frequency of each PCS data set. 1129 @type frq: numpy rank-1 array 1130 @keyword paramag_centre: The paramagnetic centre position (or positions). 1131 @type paramag_centre: numpy rank-1, 3D array or rank-2, Nx3 array 1132 @keyword com: The centre of mass of the system. This is used for defining the rotor model systems. 1133 @type com: numpy 3D rank-1 array 1134 @keyword ave_pos_pivot: The pivot point to rotate all atoms about to the average domain position. In most cases this will be the centre of mass of the moving domain. This pivot is shifted by the translation vector. 1135 @type ave_pos_pivot: numpy 3D rank-1 array 1136 @keyword pivot: The pivot point for the ball-and-socket joint motion. This is needed if PCS or PRE values are used. 1137 @type pivot: numpy rank-1, 3D array or None 1138 @keyword pivot_opt: A flag which if True will allow the pivot point of the motion to be optimised. 1139 @type pivot_opt: bool 1140 @keyword sobol_max_points: The maximum number of Sobol' points to use for the numerical PCS integration technique. 1141 @type sobol_max_points: int 1142 @keyword sobol_oversample: The oversampling factor Ov used for the total number of points N * Ov * 10**M, where N is the maximum number of Sobol' points and M is the number of dimensions or torsion-tilt angles for the system. 1143 @type sobol_oversample: int 1144 @keyword scaling_matrix: The diagonal, square scaling matrix. 1145 @type scaling_matrix: numpy diagonal matrix 1146 @keyword quad_int: A flag which if True will perform high precision numerical integration via the scipy.integrate quad(), dblquad() and tplquad() integration methods rather than the rough quasi-random numerical integration. 1147 @type quad_int: bool 1148 """ 1149 1150 # Store some arguments. 1151 self.min_algor = min_algor 1152 self.min_options = min_options 1153 self.func_tol = func_tol 1154 self.grad_tol = grad_tol 1155 self.max_iterations = max_iterations 1156 self.scaling_matrix = scaling_matrix 1157 self.model = model 1158 self.param_vector = param_vector 1159 self.full_tensors = full_tensors 1160 self.full_in_ref_frame = full_in_ref_frame 1161 self.rdcs = rdcs 1162 self.rdc_err = rdc_err 1163 self.rdc_weight = rdc_weight 1164 self.rdc_vect = rdc_vect 1165 self.rdc_const = rdc_const 1166 self.pcs = pcs 1167 self.pcs_err = pcs_err 1168 self.pcs_weight = pcs_weight 1169 self.atomic_pos = atomic_pos 1170 self.temp = temp 1171 self.frq = frq 1172 self.paramag_centre = paramag_centre 1173 self.com = com 1174 self.ave_pos_pivot = ave_pos_pivot 1175 self.pivot = pivot 1176 self.pivot_opt = pivot_opt 1177 self.sobol_max_points = sobol_max_points 1178 self.sobol_oversample = sobol_oversample 1179 self.verbosity = verbosity 1180 self.quad_int = quad_int 1181 1182 # Feedback on the number of integration points used (target function setup required). This must be run here on the master and not in run() on the slave. 1183 target_fn = Frame_order(model=self.model, init_params=self.param_vector, full_tensors=self.full_tensors, full_in_ref_frame=self.full_in_ref_frame, rdcs=self.rdcs, rdc_errors=self.rdc_err, rdc_weights=self.rdc_weight, rdc_vect=self.rdc_vect, dip_const=self.rdc_const, pcs=self.pcs, pcs_errors=self.pcs_err, pcs_weights=self.pcs_weight, atomic_pos=self.atomic_pos, temp=self.temp, frq=self.frq, paramag_centre=self.paramag_centre, scaling_matrix=self.scaling_matrix, com=self.com, ave_pos_pivot=self.ave_pos_pivot, pivot=self.pivot, pivot_opt=self.pivot_opt, sobol_max_points=self.sobol_max_points, sobol_oversample=self.sobol_oversample, quad_int=self.quad_int) 1184 if not self.quad_int: 1185 count_sobol_points(target_fn=target_fn, verbosity=self.verbosity) 1186 1187 # Linear constraints. 1188 self.A, self.b = None, None 1189 if constraints: 1190 # Obtain the constraints. 1191 self.A, self.b = linear_constraints(scaling_matrix=scaling_matrix) 1192 1193 # Constraint flag set but no constraints present. 1194 if self.A is None: 1195 if verbosity: 1196 warn(RelaxWarning("The '%s' model parameters are not constrained, turning the linear constraint algorithm off." % cdp.model)) 1197 1198 # Pop out the log barrier algorithm. 1199 if self.min_algor == 'Log barrier': 1200 self.min_algor = self.min_options[0] 1201 self.min_options = self.min_options[1:]

1202 1203

1204 - def run(self, processor, completed):

1205 """Set up and perform the optimisation.""" 1206 1207 # Set up the optimisation target function class. 1208 target_fn = Frame_order(model=self.model, init_params=self.param_vector, full_tensors=self.full_tensors, full_in_ref_frame=self.full_in_ref_frame, rdcs=self.rdcs, rdc_errors=self.rdc_err, rdc_weights=self.rdc_weight, rdc_vect=self.rdc_vect, dip_const=self.rdc_const, pcs=self.pcs, pcs_errors=self.pcs_err, pcs_weights=self.pcs_weight, atomic_pos=self.atomic_pos, temp=self.temp, frq=self.frq, paramag_centre=self.paramag_centre, scaling_matrix=self.scaling_matrix, com=self.com, ave_pos_pivot=self.ave_pos_pivot, pivot=self.pivot, pivot_opt=self.pivot_opt, sobol_max_points=self.sobol_max_points, sobol_oversample=self.sobol_oversample, quad_int=self.quad_int) 1209 1210 # Minimisation. 1211 results = generic_minimise(func=target_fn.func, args=(), x0=self.param_vector, min_algor=self.min_algor, min_options=self.min_options, func_tol=self.func_tol, grad_tol=self.grad_tol, maxiter=self.max_iterations, A=self.A, b=self.b, full_output=True, print_flag=self.verbosity) 1212 1213 # Create the result command object on the slave to send back to the master. 1214 processor.return_object(Frame_order_result_command(processor=processor, memo_id=self.memo_id, results=results, A_5D_bc=target_fn.A_5D_bc, pcs_theta=target_fn.pcs_theta, rdc_theta=target_fn.rdc_theta, completed=completed))

1215 1216 1217

1218 -class Frame_order_result_command(Result_command):

1219 """Class for processing the frame order results. 1220 1221 This object will be sent from the slave back to the master to have its run() method executed. 1222 """ 1223

1224 - def __init__(self, processor=None, memo_id=None, results=None, A_5D_bc=None, pcs_theta=None, rdc_theta=None, completed=False):

1225 """Set up the class, placing the minimisation results here. 1226 1227 @keyword processor: The processor object. 1228 @type processor: multi.processor.Processor instance 1229 @keyword memo_id: The memo identification string. 1230 @type memo_id: str 1231 @keyword results: The results as returned by minfx. 1232 @type results: tuple 1233 @keyword A_5D_bc: The reduced back-calculated alignment tensors from the target function. 1234 @type A_5D_bc: numpy float64 array 1235 @keyword pcs_theta: The back calculated PCS values from the target function. 1236 @type pcs_theta: numpy float64 array 1237 @keyword rdc_theta: The back calculated RDC values from the target function. 1238 @type rdc_theta: numpy float64 array 1239 @keyword model: The target function class instance which has been optimised. 1240 @type model: class instance 1241 @keyword completed: A flag which if True signals that the optimisation successfully completed. 1242 @type completed: bool 1243 """ 1244 1245 # Execute the base class __init__() method. 1246 super(Frame_order_result_command, self).__init__(processor=processor, completed=completed) 1247 1248 # Store the arguments. 1249 self.memo_id = memo_id 1250 self.results = results 1251 self.A_5D_bc = A_5D_bc 1252 self.pcs_theta = pcs_theta 1253 self.rdc_theta = rdc_theta

1254 1255

1256 - def run(self, processor, memo):

1257 """Disassemble the frame order optimisation results. 1258 1259 @param processor: Unused! 1260 @type processor: None 1261 @param memo: The model-free memo. 1262 @type memo: memo 1263 """ 1264 1265 # Printout. 1266 if memo.sim_index != None: 1267 print("Simulation %i" % (memo.sim_index+1)) 1268 1269 # Disassemble the results. 1270 if len(self.results) == 4: # Grid search. 1271 param_vector, func, iter_count, warning = self.results 1272 f_count = iter_count 1273 g_count = 0.0 1274 h_count = 0.0 1275 else: 1276 param_vector, func, iter_count, f_count, g_count, h_count, warning = self.results 1277 1278 # Check if the chi-squared value is lower - required for a parallelised grid search to work. 1279 if memo.sim_index == None: 1280 if hasattr(cdp, 'chi2') and cdp.chi2 != None: 1281 # No improvement, so exit the function. 1282 if func > cdp.chi2: 1283 print("Discarding the optimisation results, the optimised chi-squared value is higher than the current value (%s > %s)." % (func, cdp.chi2)) 1284 return 1285 1286 # New minimum. 1287 print("Storing the optimisation results, the optimised chi-squared value is lower than the current value (%s < %s)." % (func, cdp.chi2)) 1288 1289 # Just print something to avoid user confusion in parallelised grid searches. 1290 else: 1291 print("Storing the optimisation results, no optimised values currently exist.") 1292 1293 # Unpack the results. 1294 unpack_opt_results(param_vector, func, iter_count, f_count, g_count, h_count, warning, memo.scaling_matrix, memo.sim_index) 1295 1296 # Store the back-calculated data. 1297 store_bc_data(A_5D_bc=self.A_5D_bc, pcs_theta=self.pcs_theta, rdc_theta=self.rdc_theta)

1298

Source Code for Module specific_analyses.frame_order.optimisation