specific_analyses.frame_order.optimisation

1 ############################################################################### 2 # # 3 # Copyright (C) 2009-2014 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 numpy import arccos, array, dot, float64, ones, zeros 28 from numpy.linalg import inv, norm 29 import sys 30 from warnings import warn 31 32 # relax module imports. 33 from lib.float import isNaN, isInf 34 from lib.errors import RelaxError, RelaxInfError, RelaxNaNError, RelaxNoPCSError, RelaxNoRDCError 35 from lib.geometry.angles import wrap_angles 36 from lib.order import order_parameters 37 from lib.periodic_table import periodic_table 38 from lib.physical_constants import dipolar_constant 39 from lib.warnings import RelaxWarning 40 from pipe_control.interatomic import interatomic_loop 41 from pipe_control.mol_res_spin import return_spin, spin_loop 42 from pipe_control.rdc import check_rdcs 43 from pipe_control.structure.mass import pipe_centre_of_mass 44 from specific_analyses.frame_order.checks import check_ave_domain_setup 45 from specific_analyses.frame_order.data import base_data_types, domain_moving, pivot_fixed, tensor_loop, translation_fixed 46 from specific_analyses.frame_order.parameters import assemble_param_vector 47 from target_functions import frame_order 48 49

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

51 """Set up a row of the grid search for a given parameter. 52 53 @param incs: The number of increments. 54 @type incs: int 55 @param lower: The lower bounds. 56 @type lower: float 57 @param upper: The upper bounds. 58 @type upper: float 59 @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). 60 @type dist_type: None or str 61 @keyword end_point: A flag which if False will cause the end point to be removed. 62 @type end_point: bool 63 @return: The row of the grid. 64 @rtype: list of float 65 """ 66 67 # Init. 68 row = [] 69 70 # Linear grid. 71 if dist_type == None: 72 # Loop over the increments. 73 for i in range(incs): 74 # The row. 75 row.append(lower + i * (upper - lower) / (incs - 1.0)) 76 77 # Inverse cos distribution. 78 elif dist_type == 'acos': 79 # Generate the increment values of v from cos(upper) to cos(lower). 80 v = zeros(incs, float64) 81 val = (cos(lower) - cos(upper)) / (incs - 1.0) 82 for i in range(incs): 83 v[-i-1] = cos(upper) + float(i) * val 84 85 # Generate the distribution. 86 row = arccos(v) 87 88 # Remove the last point. 89 if not end_point: 90 row = row[:-1] 91 92 # Return the row (as a list). 93 return list(row)

94 95

96 -def minimise_setup_atomic_pos(sim_index=None):

97 """Set up the atomic position data structures for optimisation using PCSs and PREs as base data sets. 98 99 @keyword sim_index: The index of the simulation to optimise. This should be None if normal optimisation is desired. 100 @type sim_index: None or int 101 @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. 102 @rtype: numpy rank-3 array, numpy rank-1 array. 103 """ 104 105 # Initialise. 106 atomic_pos = [] 107 108 # Store the atomic positions. 109 for spin, spin_id in spin_loop(selection=domain_moving(), return_id=True): 110 # Skip deselected spins. 111 if not spin.select: 112 continue 113 114 # Only use spins with alignment/paramagnetic data. 115 if not hasattr(spin, 'pcs') and not hasattr(spin, 'pre'): 116 continue 117 118 # A single atomic position. 119 if spin.pos.shape == (3,): 120 atomic_pos.append(spin.pos) 121 122 # A single model (rank-2 array of a single position). 123 elif spin.pos.shape == (1, 3): 124 atomic_pos.append(spin.pos[0]) 125 126 # Average multiple atomic positions. 127 else: 128 # First throw a warning to tell the user what is happening. 129 if sim_index == None: 130 warn(RelaxWarning("Averaging the %s atomic positions for the PCS for the spin '%s'." % (len(spin.pos), spin_id))) 131 132 # The average position. 133 ave_pos = zeros(3, float64) 134 for i in range(len(spin.pos)): 135 ave_pos += spin.pos[i] 136 ave_pos = ave_pos / len(spin.pos) 137 138 # Store. 139 atomic_pos.append(ave_pos) 140 141 # Convert to numpy objects. 142 atomic_pos = array(atomic_pos, float64) 143 144 # The paramagnetic centre. 145 if not hasattr(cdp, 'paramagnetic_centre'): 146 paramag_centre = zeros(3, float64) 147 else: 148 paramag_centre = array(cdp.paramagnetic_centre) 149 150 # Return the data structures. 151 return atomic_pos, paramag_centre

152 153

154 -def minimise_setup_pcs(sim_index=None):

155 """Set up the data structures for optimisation using PCSs as base data sets. 156 157 @keyword sim_index: The index of the simulation to optimise. This should be None if normal optimisation is desired. 158 @type sim_index: None or int 159 @return: The assembled data structures for using PCSs as the base data for optimisation. These include: 160 - the PCS values. 161 - the unit vectors connecting the paramagnetic centre (the electron spin) to 162 - the PCS weight. 163 - the nuclear spin. 164 - the pseudocontact shift constants. 165 @rtype: tuple of (numpy rank-2 array, numpy rank-2 array, numpy rank-2 array, numpy rank-1 array, numpy rank-1 array) 166 """ 167 168 # Data setup tests. 169 if not hasattr(cdp, 'paramagnetic_centre'): 170 raise RelaxError("The paramagnetic centre has not yet been specified.") 171 if not hasattr(cdp, 'temperature'): 172 raise RelaxError("The experimental temperatures have not been set.") 173 if not hasattr(cdp, 'spectrometer_frq'): 174 raise RelaxError("The spectrometer frequencies of the experiments have not been set.") 175 176 # Initialise. 177 pcs = [] 178 pcs_err = [] 179 pcs_weight = [] 180 temp = [] 181 frq = [] 182 183 # The PCS data. 184 for i in range(len(cdp.align_ids)): 185 # Alias the ID. 186 align_id = cdp.align_ids[i] 187 188 # Skip non-optimised data. 189 if not opt_uses_align_data(align_id): 190 continue 191 192 # Append empty arrays to the PCS structures. 193 pcs.append([]) 194 pcs_err.append([]) 195 pcs_weight.append([]) 196 197 # Get the temperature for the PCS constant. 198 if align_id in cdp.temperature: 199 temp.append(cdp.temperature[align_id]) 200 201 # The temperature must be given! 202 else: 203 raise RelaxError("The experimental temperature for the alignment ID '%s' has not been set." % align_id) 204 205 # Get the spectrometer frequency in Tesla units for the PCS constant. 206 if align_id in cdp.spectrometer_frq: 207 frq.append(cdp.spectrometer_frq[align_id] * 2.0 * pi / periodic_table.gyromagnetic_ratio('1H')) 208 209 # The frequency must be given! 210 else: 211 raise RelaxError("The spectrometer frequency for the alignment ID '%s' has not been set." % align_id) 212 213 # Spin loop over the domain. 214 j = 0 215 for spin in spin_loop(selection=domain_moving()): 216 # Skip deselected spins. 217 if not spin.select: 218 continue 219 220 # Skip spins without PCS data. 221 if not hasattr(spin, 'pcs'): 222 continue 223 224 # Append the PCSs to the list. 225 if align_id in spin.pcs: 226 if sim_index != None: 227 pcs[-1].append(spin.pcs_sim[align_id][sim_index]) 228 else: 229 pcs[-1].append(spin.pcs[align_id]) 230 else: 231 pcs[-1].append(None) 232 233 # Append the PCS errors. 234 if hasattr(spin, 'pcs_err') and align_id in spin.pcs_err: 235 pcs_err[-1].append(spin.pcs_err[align_id]) 236 else: 237 pcs_err[-1].append(None) 238 239 # Append the weight. 240 if hasattr(spin, 'pcs_weight') and align_id in spin.pcs_weight: 241 pcs_weight[-1].append(spin.pcs_weight[align_id]) 242 else: 243 pcs_weight[-1].append(1.0) 244 245 # Spin index. 246 j = j + 1 247 248 # Convert to numpy objects. 249 pcs = array(pcs, float64) 250 pcs_err = array(pcs_err, float64) 251 pcs_weight = array(pcs_weight, float64) 252 253 # Convert the PCS from ppm to no units. 254 pcs = pcs * 1e-6 255 pcs_err = pcs_err * 1e-6 256 257 # Return the data structures. 258 return pcs, pcs_err, pcs_weight, temp, frq

259 260

261 -def minimise_setup_rdcs(sim_index=None):

262 """Set up the data structures for optimisation using RDCs as base data sets. 263 264 @keyword sim_index: The index of the simulation to optimise. This should be None if normal optimisation is desired. 265 @type sim_index: None or int 266 @return: The assembled data structures for using RDCs as the base data for optimisation. These include: 267 - rdc, the RDC values. 268 - rdc_err, the RDC errors. 269 - rdc_weight, the RDC weights. 270 - vectors, the interatomic vectors. 271 - rdc_const, the dipolar constants. 272 - absolute, the absolute value flags (as 1's and 0's). 273 @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) 274 """ 275 276 # Initialise. 277 rdc = [] 278 rdc_err = [] 279 rdc_weight = [] 280 unit_vect = [] 281 rdc_const = [] 282 absolute = [] 283 284 # The unit vectors and RDC constants. 285 for interatom in interatomic_loop(selection1=domain_moving()): 286 # RDC checks. 287 if not check_rdcs(interatom): 288 continue 289 290 # Get the spins. 291 spin1 = return_spin(interatom.spin_id1) 292 spin2 = return_spin(interatom.spin_id2) 293 294 # A single unit vector. 295 if interatom.vector.shape == (3,): 296 unit_vect.append(interatom.vector) 297 298 # Average multiple unit vectors. 299 else: 300 # First throw a warning to tell the user what is happening. 301 if sim_index == None: 302 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))) 303 304 # The average position. 305 ave_vector = zeros(3, float64) 306 for i in range(len(interatom.vector)): 307 ave_vector += interatom.vector[i] 308 309 # Store. 310 unit_vect.append(ave_vector) 311 312 # Normalise (to be sure). 313 unit_vect[-1] = unit_vect[-1] / norm(unit_vect[-1]) 314 315 # Gyromagnetic ratios. 316 g1 = periodic_table.gyromagnetic_ratio(spin1.isotope) 317 g2 = periodic_table.gyromagnetic_ratio(spin2.isotope) 318 319 # Calculate the RDC dipolar constant (in Hertz, and the 3 comes from the alignment tensor), and append it to the list. 320 rdc_const.append(3.0/(2.0*pi) * dipolar_constant(g1, g2, interatom.r)) 321 322 # Fix the unit vector data structure. 323 num = None 324 for rdc_index in range(len(unit_vect)): 325 # Number of vectors. 326 if num == None: 327 if unit_vect[rdc_index] != None: 328 num = len(unit_vect[rdc_index]) 329 continue 330 331 # Check. 332 if unit_vect[rdc_index] != None and len(unit_vect[rdc_index]) != num: 333 raise RelaxError("The number of interatomic vectors for all no match:\n%s" % unit_vect) 334 335 # Missing unit vectors. 336 if num == None: 337 raise RelaxError("No interatomic vectors could be found.") 338 339 # Update None entries. 340 for i in range(len(unit_vect)): 341 if unit_vect[i] == None: 342 unit_vect[i] = [[None, None, None]]*num 343 344 # The RDC data. 345 for i in range(len(cdp.align_ids)): 346 # Alias the ID. 347 align_id = cdp.align_ids[i] 348 349 # Skip non-optimised data. 350 if not opt_uses_align_data(align_id): 351 continue 352 353 # Append empty arrays to the RDC structures. 354 rdc.append([]) 355 rdc_err.append([]) 356 rdc_weight.append([]) 357 absolute.append([]) 358 359 # Interatom loop over the domain. 360 for interatom in interatomic_loop(domain_moving()): 361 # Get the spins. 362 spin1 = return_spin(interatom.spin_id1) 363 spin2 = return_spin(interatom.spin_id2) 364 365 # Skip deselected spins. 366 if not spin1.select or not spin2.select: 367 continue 368 369 # Only use interatomic data containers with RDC and vector data. 370 if not hasattr(interatom, 'rdc') or not hasattr(interatom, 'vector'): 371 continue 372 373 # Defaults of None. 374 value = None 375 error = None 376 377 # Pseudo-atom set up. 378 if (hasattr(spin1, 'members') or hasattr(spin2, 'members')) and align_id in interatom.rdc: 379 raise RelaxError("Psuedo-atoms are currently not supported for the frame order analysis.") 380 381 # Normal set up. 382 elif align_id in interatom.rdc: 383 # The RDC. 384 if sim_index != None: 385 value = interatom.rdc_sim[align_id][sim_index] 386 else: 387 value = interatom.rdc[align_id] 388 389 # The error. 390 if hasattr(interatom, 'rdc_err') and align_id in interatom.rdc_err: 391 error = interatom.rdc_err[align_id] 392 393 # Append the RDCs to the list. 394 rdc[-1].append(value) 395 396 # Append the RDC errors. 397 rdc_err[-1].append(error) 398 399 # Append the weight. 400 if hasattr(interatom, 'rdc_weight') and align_id in interatom.rdc_weight: 401 rdc_weight[-1].append(interatom.rdc_weight[align_id]) 402 else: 403 rdc_weight[-1].append(1.0) 404 405 # Append the absolute value flag. 406 if hasattr(interatom, 'absolute_rdc') and align_id in interatom.absolute_rdc: 407 absolute[-1].append(interatom.absolute_rdc[align_id]) 408 else: 409 absolute[-1].append(False) 410 411 # Convert to numpy objects. 412 rdc = array(rdc, float64) 413 rdc_err = array(rdc_err, float64) 414 rdc_weight = array(rdc_weight, float64) 415 unit_vect = array(unit_vect, float64) 416 rdc_const = array(rdc_const, float64) 417 absolute = array(absolute, float64) 418 419 # Return the data structures. 420 return rdc, rdc_err, rdc_weight, unit_vect, rdc_const, absolute

421 422

423 -def minimise_setup_tensors(sim_index=None):

424 """Set up the data structures for optimisation using alignment tensors as base data sets. 425 426 @keyword sim_index: The simulation index. This should be None if normal optimisation is desired. 427 @type sim_index: None or int 428 @return: The assembled data structures for using alignment tensors as the base data for optimisation. These include: 429 - full_tensors, the full tensors as concatenated arrays. 430 - full_err, the full tensor errors as concatenated arrays. 431 - full_in_ref_frame, the flags specifying if the tensor is the full or reduced tensor in the non-moving reference domain. 432 @rtype: tuple of 3 numpy nx5D, rank-1 arrays 433 """ 434 435 # Checks. 436 if not hasattr(cdp, 'ref_domain'): 437 raise RelaxError("The reference domain has not been set up.") 438 if not hasattr(cdp.align_tensors, 'reduction'): 439 raise RelaxError("The tensor reductions have not been specified.") 440 for i, tensor in tensor_loop(): 441 if not hasattr(tensor, 'domain'): 442 raise RelaxError("The domain that the '%s' tensor is attached to has not been set" % tensor.name) 443 444 # Initialise. 445 n = len(cdp.align_tensors.reduction) 446 full_tensors = zeros(n*5, float64) 447 full_err = ones(n*5, float64) * 1e-5 448 full_in_ref_frame = zeros(n, float64) 449 450 # Loop over the full tensors. 451 for i, tensor in tensor_loop(red=False): 452 # The full tensor. 453 full_tensors[5*i + 0] = tensor.Axx 454 full_tensors[5*i + 1] = tensor.Ayy 455 full_tensors[5*i + 2] = tensor.Axy 456 full_tensors[5*i + 3] = tensor.Axz 457 full_tensors[5*i + 4] = tensor.Ayz 458 459 # The full tensor corresponds to the frame of reference. 460 if cdp.ref_domain == tensor.domain: 461 full_in_ref_frame[i] = 1 462 463 # The full tensor errors. 464 if hasattr(tensor, 'Axx_err'): 465 full_err[5*i + 0] = tensor.Axx_err 466 full_err[5*i + 1] = tensor.Ayy_err 467 full_err[5*i + 2] = tensor.Axy_err 468 full_err[5*i + 3] = tensor.Axz_err 469 full_err[5*i + 4] = tensor.Ayz_err 470 471 # Return the data structures. 472 return full_tensors, full_err, full_in_ref_frame

473 474

475 -def opt_uses_align_data(align_id=None):

476 """Determine if the PCS or RDC data for the given alignment ID is needed for optimisation. 477 478 @keyword align_id: The optional alignment ID string. 479 @type align_id: str 480 @return: True if alignment data is to be used for optimisation, False otherwise. 481 @rtype: bool 482 """ 483 484 # No alignment IDs. 485 if not hasattr(cdp, 'align_ids'): 486 return False 487 488 # Convert the align IDs to an array, or take all IDs. 489 if align_id: 490 align_ids = [align_id] 491 else: 492 align_ids = cdp.align_ids 493 494 # Check the PCS and RDC. 495 for align_id in align_ids: 496 if opt_uses_pcs(align_id) or opt_uses_rdc(align_id): 497 return True 498 499 # No alignment data is used for optimisation. 500 return False

501 502

503 -def opt_uses_pcs(align_id):

504 """Determine if the PCS data for the given alignment ID is needed for optimisation. 505 506 @param align_id: The alignment ID string. 507 @type align_id: str 508 @return: True if the PCS data is to be used for optimisation, False otherwise. 509 @rtype: bool 510 """ 511 512 # No alignment IDs. 513 if not hasattr(cdp, 'pcs_ids'): 514 return False 515 516 # No PCS data for the alignment. 517 if align_id not in cdp.pcs_ids: 518 return False 519 520 # The PCS data is to be used for optimisation. 521 return True

522 523

524 -def opt_uses_rdc(align_id):

525 """Determine if the RDC data for the given alignment ID is needed for optimisation. 526 527 @param align_id: The alignment ID string. 528 @type align_id: str 529 @return: True if the RDC data is to be used for optimisation, False otherwise. 530 @rtype: bool 531 """ 532 533 # No alignment IDs. 534 if not hasattr(cdp, 'rdc_ids'): 535 return False 536 537 # No RDC data for the alignment. 538 if align_id not in cdp.rdc_ids: 539 return False 540 541 # The RDC data is to be used for optimisation. 542 return True

543 544

545 -def store_bc_data(target_fn):

546 """Store the back-calculated data. 547 548 @param target_fn: The frame-order target function class. 549 @type target_fn: class instance 550 """ 551 552 # Loop over the reduced tensors. 553 for i, tensor in tensor_loop(red=True): 554 # Store the values. 555 tensor.set(param='Axx', value=target_fn.A_5D_bc[5*i + 0]) 556 tensor.set(param='Ayy', value=target_fn.A_5D_bc[5*i + 1]) 557 tensor.set(param='Axy', value=target_fn.A_5D_bc[5*i + 2]) 558 tensor.set(param='Axz', value=target_fn.A_5D_bc[5*i + 3]) 559 tensor.set(param='Ayz', value=target_fn.A_5D_bc[5*i + 4]) 560 561 # The RDC data. 562 for i in range(len(cdp.align_ids)): 563 # The alignment ID. 564 align_id = cdp.align_ids[i] 565 566 # Data flags 567 rdc_flag = False 568 if hasattr(cdp, 'rdc_ids') and align_id in cdp.rdc_ids: 569 rdc_flag = True 570 pcs_flag = False 571 if hasattr(cdp, 'pcs_ids') and align_id in cdp.pcs_ids: 572 pcs_flag = True 573 574 # Spin loop over the domain. 575 pcs_index = 0 576 for spin in spin_loop(domain_moving()): 577 # Skip deselected spins. 578 if not spin.select: 579 continue 580 581 # Spins with PCS data. 582 if pcs_flag and hasattr(spin, 'pcs'): 583 # Initialise the data structure. 584 if not hasattr(spin, 'pcs_bc'): 585 spin.pcs_bc = {} 586 587 # Store the back-calculated value (in ppm). 588 spin.pcs_bc[align_id] = target_fn.pcs_theta[i, pcs_index] * 1e6 589 590 # Increment the index. 591 pcs_index += 1 592 593 # Interatomic data container loop. 594 rdc_index = 0 595 for interatom in interatomic_loop(domain_moving()): 596 # RDC checks. 597 if not check_rdcs(interatom): 598 continue 599 600 # Initialise the data structure. 601 if not hasattr(interatom, 'rdc_bc'): 602 interatom.rdc_bc = {} 603 604 # Store the back-calculated value. 605 interatom.rdc_bc[align_id] = target_fn.rdc_theta[i, rdc_index] 606 607 # Increment the index. 608 rdc_index += 1

609 610

611 -def target_fn_setup(sim_index=None, verbosity=1, scaling_matrix=None):

612 """Initialise the target function for optimisation or direct calculation. 613 614 @keyword sim_index: The index of the simulation to optimise. This should be None if normal optimisation is desired. 615 @type sim_index: None or int 616 @keyword verbosity: The amount of information to print. The higher the value, the greater the verbosity. 617 @type verbosity: int 618 @keyword scaling_matrix: The diagonal and square scaling matrices. 619 @type scaling_matrix: numpy rank-2, float64 array or None 620 """ 621 622 # Check for the average domain displacement information. 623 check_ave_domain_setup() 624 625 # Assemble the parameter vector. 626 param_vector = assemble_param_vector(sim_index=sim_index) 627 628 # Determine the base data types (RDCs and/or PCSs). 629 data_types = base_data_types() 630 631 # Diagonal scaling. 632 if scaling_matrix != None and len(param_vector): 633 param_vector = dot(inv(scaling_matrix), param_vector) 634 635 # Get the data structures for optimisation using the tensors as base data sets. 636 full_tensors, full_tensor_err, full_in_ref_frame = minimise_setup_tensors(sim_index) 637 638 # Get the data structures for optimisation using PCSs as base data sets. 639 pcs, pcs_err, pcs_weight, temp, frq = None, None, None, None, None 640 if 'pcs' in data_types: 641 pcs, pcs_err, pcs_weight, temp, frq = minimise_setup_pcs(sim_index=sim_index) 642 643 # Get the data structures for optimisation using RDCs as base data sets. 644 rdcs, rdc_err, rdc_weight, rdc_vect, rdc_const, absolute_rdc = None, None, None, None, None, None 645 if 'rdc' in data_types: 646 rdcs, rdc_err, rdc_weight, rdc_vect, rdc_const, absolute_rdc = minimise_setup_rdcs(sim_index=sim_index) 647 648 # Data checks. 649 if pcs != None and not len(pcs): 650 raise RelaxNoPCSError 651 if rdcs != None and not len(rdcs): 652 raise RelaxNoRDCError 653 654 # Get the atomic_positions. 655 atomic_pos, paramag_centre = None, None 656 if 'pcs' in data_types or 'pre' in data_types: 657 atomic_pos, paramag_centre = minimise_setup_atomic_pos(sim_index=sim_index) 658 659 # Average domain translation. 660 translation_opt = False 661 if not translation_fixed(): 662 translation_opt = True 663 664 # The fixed pivot point. 665 pivot = None 666 if hasattr(cdp, 'pivot_x'): 667 pivot = array([cdp.pivot_x, cdp.pivot_y, cdp.pivot_z]) 668 669 # The second pivot. 670 pivot2 = None 671 if hasattr(cdp, 'pivot_x_2'): 672 pivot2 = array([cdp.pivot_x_2, cdp.pivot_y_2, cdp.pivot_z_2]) 673 674 # Pivot optimisation. 675 pivot_opt = True 676 if pivot_fixed(): 677 pivot_opt = False 678 679 # The number of integration points. 680 if not hasattr(cdp, 'num_int_pts'): 681 cdp.num_int_pts = 200000 682 683 # The centre of mass, for use in the rotor models. 684 com = None 685 if cdp.model in ['rotor', 'double rotor']: 686 # The centre of mass of all objects in the data pipe. 687 com = pipe_centre_of_mass(verbosity=0) 688 com = array(com, float64) 689 690 # Printout. 691 if verbosity: 692 print("The centre of mass reference coordinate for the rotor models is at:\n %s" % list(com)) 693 694 # The centre of mass of the moving domain - to use as the centroid for the average domain position rotation. 695 if cdp.ave_pos_pivot == 'com': 696 ave_pos_pivot = pipe_centre_of_mass(atom_id=domain_moving(), verbosity=0) 697 ave_pos_piv_sync = False 698 699 # The centre of mass of the moving domain - to use as the centroid for the average domain position rotation. 700 if cdp.ave_pos_pivot == 'motional': 701 ave_pos_pivot = pivot 702 ave_pos_piv_sync = True 703 704 # Print outs. 705 if sim_index == None: 706 if cdp.model != 'rigid': 707 if cdp.quad_int: 708 sys.stdout.write("Numerical integration via Scipy quadratic integration.\n") 709 else: 710 sys.stdout.write("Numerical integration via the quasi-random Sobol' sequence.\n") 711 sys.stdout.write("Number of integration points: %s\n" % cdp.num_int_pts) 712 base_data = [] 713 if rdcs != None and len(rdcs): 714 base_data.append("RDCs") 715 if pcs != None and len(pcs): 716 base_data.append("PCSs") 717 sys.stdout.write("Base data: %s\n" % repr(base_data)) 718 719 # Set up the optimisation function. 720 target = frame_order.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=scaling_matrix, com=com, ave_pos_pivot=ave_pos_pivot, ave_pos_piv_sync=ave_pos_piv_sync, translation_opt=translation_opt, pivot=pivot, pivot2=pivot2, pivot_opt=pivot_opt, num_int_pts=cdp.num_int_pts, quad_int=cdp.quad_int) 721 722 # Return the data. 723 return target, param_vector

724 725

726 -def unpack_opt_results(results, scaling_matrix=None, sim_index=None):

727 """Unpack and store the Frame Order optimisation results. 728 729 @param results: The results tuple returned by the minfx generic_minimise() function. 730 @type results: tuple 731 @keyword scaling_matrix: The diagonal and square scaling matrices. 732 @type scaling_matrix: numpy rank-2, float64 array or None 733 @keyword sim_index: The index of the simulation to optimise. This should be None for normal optimisation. 734 @type sim_index: None or int 735 """ 736 737 # Disassemble the results. 738 if len(results) == 4: # Grid search. 739 param_vector, func, iter_count, warning = results 740 f_count = iter_count 741 g_count = 0.0 742 h_count = 0.0 743 else: 744 param_vector, func, iter_count, f_count, g_count, h_count, warning = results 745 746 # Catch infinite chi-squared values. 747 if isInf(func): 748 raise RelaxInfError('chi-squared') 749 750 # Catch chi-squared values of NaN. 751 if isNaN(func): 752 raise RelaxNaNError('chi-squared') 753 754 # Scaling. 755 if scaling_matrix != None: 756 param_vector = dot(scaling_matrix, param_vector) 757 758 # The parameters to wrap. 759 wrap = [ 760 'ave_pos_alpha', 761 'ave_pos_beta', 762 'ave_pos_gamma', 763 'eigen_alpha', 764 'eigen_beta', 765 'eigen_gamma', 766 'axis_theta', 767 'axis_phi' 768 ] 769 770 # Monte Carlo simulation data structures. 771 if sim_index != None: 772 # Loop over the parameters. 773 for i in range(len(cdp.params)): 774 # Angle wrapping around the real value. 775 if cdp.params[i] in wrap or cdp.params[i] == 'axis_alpha': 776 val = getattr(cdp, cdp.params[i]) 777 param_vector[i] = wrap_angles(param_vector[i], val-pi, val+pi) 778 779 # FIXME: Implement linear constraints via the log-barrier algorithm, then delete this. 780 # Handle negative values of the cone_sigma_max parameter. 781 if cdp.params[i] == 'cone_sigma_max': 782 param_vector[i] = abs(param_vector[i]) 783 784 # Store the value. 785 obj = getattr(cdp, cdp.params[i]+'_sim') 786 obj[sim_index] = param_vector[i] 787 788 # Order parameter to angle conversion. 789 if cdp.params[i] == 'cone_s1': 790 cdp.cone_theta_sim[sim_index] = order_parameters.iso_cone_S_to_theta(param_vector[i]) 791 792 # Optimisation stats. 793 cdp.chi2_sim[sim_index] = func 794 cdp.iter_sim[sim_index] = iter_count 795 cdp.f_count_sim[sim_index] = f_count 796 cdp.g_count_sim[sim_index] = g_count 797 cdp.h_count_sim[sim_index] = h_count 798 cdp.warning_sim[sim_index] = warning 799 800 # Normal data structures. 801 else: 802 # Loop over the parameters. 803 for i in range(len(cdp.params)): 804 # Angle wrapping. 805 if cdp.params[i] in wrap: 806 param_vector[i] = wrap_angles(param_vector[i], 0.0, 2.0*pi) 807 if cdp.params[i] == 'axis_alpha': 808 param_vector[i] = wrap_angles(param_vector[i], -pi, pi) 809 810 # FIXME: Implement linear constraints via the log-barrier algorithm, then delete this. 811 # Handle negative values of the cone_sigma_max parameter. 812 if cdp.params[i] == 'cone_sigma_max': 813 param_vector[i] = abs(param_vector[i]) 814 815 # Store the value. 816 setattr(cdp, cdp.params[i], param_vector[i]) 817 818 # Order parameter to angle conversion. 819 if cdp.params[i] == 'cone_s1': 820 cdp.cone_theta = order_parameters.iso_cone_S_to_theta(param_vector[i]) 821 822 # Optimisation stats. 823 cdp.chi2 = func 824 cdp.iter = iter_count 825 cdp.f_count = f_count 826 cdp.g_count = g_count 827 cdp.h_count = h_count 828 cdp.warning = warning

829

Source Code for Module specific_analyses.frame_order.optimisation