Package specific_analyses :: Package n_state_model :: Module api
[hide private]
[frames] | no frames]

Source Code for Module specific_analyses.n_state_model.api

  1  ############################################################################### 
  2  #                                                                             # 
  3  # Copyright (C) 2007-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  """The N-state model or structural ensemble analysis.""" 
 24   
 25  # Python module imports. 
 26  from copy import deepcopy 
 27  from minfx.generic import generic_minimise 
 28  from minfx.grid import grid 
 29  from numpy import dot, float64, zeros 
 30  from re import search 
 31  from warnings import warn 
 32   
 33  # relax module imports. 
 34  import lib.arg_check 
 35  from lib.errors import RelaxError, RelaxInfError, RelaxNaNError, RelaxNoModelError 
 36  from lib.float import isNaN, isInf 
 37  from lib.warnings import RelaxWarning 
 38  from pipe_control import align_tensor, pcs, rdc 
 39  from pipe_control.align_tensor import opt_uses_align_data, opt_uses_tensor 
 40  from pipe_control.interatomic import interatomic_loop 
 41  from pipe_control.mol_res_spin import spin_loop 
 42  from specific_analyses.api_base import API_base 
 43  from specific_analyses.api_common import API_common 
 44  from specific_analyses.n_state_model.data import base_data_types, calc_ave_dist, num_data_points 
 45  from specific_analyses.n_state_model.optimisation import minimise_bc_data, target_fn_setup 
 46  from specific_analyses.n_state_model.parameter_object import N_state_params 
 47  from specific_analyses.n_state_model.parameters import assemble_param_vector, disassemble_param_vector, linear_constraints, param_num 
 48  from target_functions.potential import quad_pot 
 49   
 50   
51 -class N_state_model(API_base, API_common):
52 """Class containing functions for the N-state model.""" 53 54 # Class variable for storing the class instance (for the singleton design pattern). 55 instance = None 56
57 - def __init__(self):
58 """Initialise the class by placing API_common methods into the API.""" 59 60 # Place methods into the API. 61 self.get_model_container = self._get_model_container_cdp 62 self.model_loop = self._model_loop_single_global 63 self.overfit_deselect = self._overfit_deselect_dummy 64 self.print_model_title = self._print_model_title_global 65 self.set_selected_sim = self._set_selected_sim_global 66 self.sim_return_selected = self._sim_return_selected_global 67 68 # Place a copy of the parameter list object in the instance namespace. 69 self._PARAMS = N_state_params()
70 71
72 - def base_data_loop(self):
73 """Loop over the base data of the spins - RDCs, PCSs, and NOESY data. 74 75 This loop iterates for each data point (RDC, PCS, NOESY) for each spin or interatomic data container, returning the identification information. 76 77 @return: A list of the spin or interatomic data container, the data type ('rdc', 'pcs', 'noesy'), and the alignment ID if required. 78 @rtype: list of [SpinContainer instance, str, str] or [InteratomContainer instance, str, str] 79 """ 80 81 # Loop over the interatomic data containers. 82 for interatom in interatomic_loop(): 83 # Skip deselected data. 84 if not interatom.select: 85 continue 86 87 # Re-initialise the data structure. 88 data = [interatom, None, None] 89 90 # RDC data. 91 if hasattr(interatom, 'rdc'): 92 data[1] = 'rdc' 93 94 # Loop over the alignment IDs. 95 for id in cdp.rdc_ids: 96 # Add the ID. 97 data[2] = id 98 99 # Yield the set. 100 yield data 101 102 # NOESY data. 103 if hasattr(interatom, 'noesy'): 104 data[1] = 'noesy' 105 106 # Loop over the alignment IDs. 107 for id in cdp.noesy_ids: 108 # Add the ID. 109 data[2] = id 110 111 # Yield the set. 112 yield data 113 114 # Loop over the spins. 115 for spin in spin_loop(): 116 # Skip deselected data. 117 if not spin.select: 118 continue 119 120 # Re-initialise the data structure. 121 data = [spin, None, None] 122 123 # PCS data. 124 if hasattr(spin, 'pcs'): 125 data[1] = 'pcs' 126 127 # Loop over the alignment IDs. 128 for id in cdp.pcs_ids: 129 # Add the ID. 130 data[2] = id 131 132 # Yield the set. 133 yield data
134 135
136 - def calculate(self, spin_id=None, scaling_matrix=None, verbosity=1, sim_index=None):
137 """Calculation function. 138 139 Currently this function simply calculates the NOESY flat-bottom quadratic energy potential, 140 if NOE restraints are available. 141 142 @keyword spin_id: The spin identification string (unused). 143 @type spin_id: None or str 144 @keyword scaling_matrix: The per-model list of diagonal and square scaling matrices. 145 @type scaling_matrix: list of numpy rank-2, float64 array or list of None 146 @keyword verbosity: The amount of information to print. The higher the value, the greater the verbosity. 147 @type verbosity: int 148 @keyword sim_index: The MC simulation index (unused). 149 @type sim_index: None 150 """ 151 152 # Set up the target function for direct calculation. 153 if scaling_matrix != None: 154 scaling_matrix = scaling_matrix[0] 155 model, param_vector, data_types = target_fn_setup(scaling_matrix=scaling_matrix, verbosity=verbosity) 156 157 # Calculate the chi-squared value. 158 if model: 159 # Make a function call. 160 chi2 = model.func(param_vector) 161 162 # Store the global chi-squared value. 163 cdp.chi2 = chi2 164 165 # Store the back-calculated data. 166 minimise_bc_data(model, sim_index=sim_index) 167 168 # Calculate the RDC Q factors. 169 if 'rdc' in data_types: 170 rdc.q_factors(sim_index=sim_index, verbosity=verbosity) 171 172 # Calculate the PCS Q factors. 173 if 'pcs' in data_types: 174 pcs.q_factors(sim_index=sim_index, verbosity=verbosity) 175 176 # NOE potential. 177 if hasattr(cdp, 'noe_restraints'): 178 # Init some numpy arrays. 179 num_restraints = len(cdp.noe_restraints) 180 dist = zeros(num_restraints, float64) 181 pot = zeros(num_restraints, float64) 182 lower = zeros(num_restraints, float64) 183 upper = zeros(num_restraints, float64) 184 185 # Loop over the NOEs. 186 for i in range(num_restraints): 187 # Create arrays of the NOEs. 188 lower[i] = cdp.noe_restraints[i][2] 189 upper[i] = cdp.noe_restraints[i][3] 190 191 # Calculate the average distances, using -6 power averaging. 192 dist[i] = calc_ave_dist(cdp.noe_restraints[i][0], cdp.noe_restraints[i][1], exp=-6) 193 194 # Calculate the quadratic potential. 195 quad_pot(dist, pot, lower, upper) 196 197 # Store the distance and potential information. 198 cdp.ave_dist = [] 199 cdp.quad_pot = [] 200 for i in range(num_restraints): 201 cdp.ave_dist.append([cdp.noe_restraints[i][0], cdp.noe_restraints[i][1], dist[i]]) 202 cdp.quad_pot.append([cdp.noe_restraints[i][0], cdp.noe_restraints[i][1], pot[i]])
203 204
205 - def create_mc_data(self, data_id=None):
206 """Create the Monte Carlo data by back-calculation. 207 208 @keyword data_id: The list of spin ID, data type, and alignment ID, as yielded by the base_data_loop() generator method. 209 @type data_id: str 210 @return: The Monte Carlo Ri data. 211 @rtype: list of floats 212 """ 213 214 # Initialise the MC data structure. 215 mc_data = [] 216 217 # Alias the spin or interatomic data container. 218 container = data_id[0] 219 220 # RDC data. 221 if data_id[1] == 'rdc' and hasattr(container, 'rdc'): 222 # Does back-calculated data exist? 223 if not hasattr(container, 'rdc_bc'): 224 self.calculate() 225 226 # The data. 227 if not hasattr(container, 'rdc_bc') or not data_id[2] in container.rdc_bc: 228 data = None 229 else: 230 data = container.rdc_bc[data_id[2]] 231 232 # Append the data. 233 mc_data.append(data) 234 235 # NOESY data. 236 elif data_id[1] == 'noesy' and hasattr(container, 'noesy'): 237 # Does back-calculated data exist? 238 if not hasattr(container, 'noesy_bc'): 239 self.calculate() 240 241 # Append the data. 242 mc_data.append(container.noesy_bc) 243 244 # PCS data. 245 elif data_id[1] == 'pcs' and hasattr(container, 'pcs'): 246 # Does back-calculated data exist? 247 if not hasattr(container, 'pcs_bc'): 248 self.calculate() 249 250 # The data. 251 if not hasattr(container, 'pcs_bc') or not data_id[2] in container.pcs_bc: 252 data = None 253 else: 254 data = container.pcs_bc[data_id[2]] 255 256 # Append the data. 257 mc_data.append(data) 258 259 # Return the data. 260 return mc_data
261 262
263 - def get_param_names(self, model_info=None):
264 """Return a vector of parameter names. 265 266 @keyword model_info: The model information from model_loop(). This is unused. 267 @type model_info: None 268 @return: The vector of parameter names. 269 @rtype: list of str 270 """ 271 272 # Init. 273 param_names = [] 274 275 # A RDC or PCS data type requires the alignment tensors to be at the start of the parameter vector (unless the tensors are fixed). 276 if opt_uses_align_data(): 277 for i in range(len(cdp.align_tensors)): 278 # Skip non-optimised tensors. 279 if not opt_uses_tensor(cdp.align_tensors[i]): 280 continue 281 282 # Add the parameters. 283 param_names += ['Axx', 'Ayy', 'Axy', 'Axz', 'Ayz'] 284 285 # Populations. 286 if cdp.model in ['2-domain', 'population']: 287 for i in range(cdp.N - 1): 288 param_names.append('probs') 289 290 # The Euler angles. 291 if cdp.model == '2-domain': 292 for i in range(cdp.N): 293 param_names.append('alpha') 294 param_names.append('beta') 295 param_names.append('gamma') 296 297 # The paramagnetic centre. 298 if hasattr(cdp, 'paramag_centre_fixed') and not cdp.paramag_centre_fixed: 299 for i in range(3): 300 param_names.append('paramagnetic_centre') 301 302 # Return the list. 303 return param_names
304 305
306 - def get_param_values(self, model_info=None, sim_index=None):
307 """Return a vector of parameter values. 308 309 @keyword model_info: The model information from model_loop(). This is unused. 310 @type model_info: None 311 @keyword sim_index: The optional Monte Carlo simulation index. 312 @type sim_index: int 313 @return: The vector of parameter values. 314 @rtype: list of str 315 """ 316 317 # Return the vector. 318 return assemble_param_vector(sim_index=sim_index)
319 320
321 - def grid_search(self, lower=None, upper=None, inc=None, scaling_matrix=None, constraints=False, verbosity=0, sim_index=None):
322 """The grid search function. 323 324 @param lower: The lower bounds of the grid search which must be equal to the number of parameters in the model. 325 @type lower: list of lists of floats 326 @param upper: The upper bounds of the grid search which must be equal to the number of parameters in the model. 327 @type upper: list of lists of floats 328 @param 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. 329 @type inc: list of lists of int 330 @keyword scaling_matrix: The per-model list of diagonal and square scaling matrices. 331 @type scaling_matrix: list of numpy rank-2, float64 array or list of None 332 @param constraints: If True, constraints are applied during the grid search (elinating parts of the grid). If False, no constraints are used. 333 @type constraints: bool 334 @param verbosity: A flag specifying the amount of information to print. The higher the value, the greater the verbosity. 335 @type verbosity: int 336 """ 337 338 # Test if the N-state model has been set up. 339 if not hasattr(cdp, 'model'): 340 raise RelaxNoModelError('N-state') 341 342 # The number of parameters. 343 n = param_num() 344 345 # Determine the data type. 346 data_types = base_data_types() 347 348 # The number of tensors to optimise. 349 tensor_num = align_tensor.num_tensors(skip_fixed=True) 350 351 # Custom sub-grid search for when only tensors are optimised (as each tensor is independent, the number of points collapses from inc**(5*N) to N*inc**5). 352 if cdp.model == 'fixed' and tensor_num > 1 and ('rdc' in data_types or 'pcs' in data_types) and not align_tensor.all_tensors_fixed() and hasattr(cdp, 'paramag_centre_fixed') and cdp.paramag_centre_fixed: 353 # Print out. 354 print("Optimising each alignment tensor separately.") 355 356 # Store the alignment tensor fixed flags. 357 fixed_flags = [] 358 for i in range(len(cdp.align_ids)): 359 # Get the tensor object. 360 tensor = align_tensor.return_tensor(index=i, skip_fixed=False) 361 362 # Store the flag. 363 fixed_flags.append(tensor.fixed) 364 365 # Fix the tensor. 366 tensor.set('fixed', True) 367 368 # Loop over each sub-grid. 369 for i in range(len(cdp.align_ids)): 370 # Skip the tensor if originally fixed. 371 if fixed_flags[i]: 372 continue 373 374 # Get the tensor object. 375 tensor = align_tensor.return_tensor(index=i, skip_fixed=False) 376 377 # Unfix the current tensor. 378 tensor.set('fixed', False) 379 380 # Grid search parameter subsets. 381 lower_sub = lower[0][i*5:i*5+5] 382 upper_sub = upper[0][i*5:i*5+5] 383 inc_sub = inc[0][i*5:i*5+5] 384 385 # Minimisation of the sub-grid. 386 self.minimise(min_algor='grid', lower=[lower_sub], upper=[upper_sub], inc=[inc_sub], scaling_matrix=[None], constraints=constraints, verbosity=verbosity, sim_index=sim_index) 387 388 # Fix the tensor again. 389 tensor.set('fixed', True) 390 391 # Reset the state of the tensors. 392 for i in range(len(cdp.align_ids)): 393 # Get the tensor object. 394 tensor = align_tensor.return_tensor(index=i, skip_fixed=False) 395 396 # Fix the tensor. 397 tensor.set('fixed', fixed_flags[i]) 398 399 # All other minimisation. 400 else: 401 self.minimise(min_algor='grid', lower=lower, upper=upper, inc=inc, scaling_matrix=scaling_matrix, constraints=constraints, verbosity=verbosity, sim_index=sim_index)
402 403
404 - def is_spin_param(self, name):
405 """Determine whether the given parameter is spin specific. 406 407 @param name: The name of the parameter. 408 @type name: str 409 @return: False 410 @rtype: bool 411 """ 412 413 # Spin specific parameters. 414 if name in ['r', 'heteronuc_type', 'proton_type']: 415 return True 416 417 # All other parameters are global. 418 return False
419 420
421 - def map_bounds(self, param, spin_id=None):
422 """Create bounds for the OpenDX mapping function. 423 424 @param param: The name of the parameter to return the lower and upper bounds of. 425 @type param: str 426 @param spin_id: The spin identification string (unused). 427 @type spin_id: None 428 @return: The upper and lower bounds of the parameter. 429 @rtype: list of float 430 """ 431 432 # Paramagnetic centre. 433 if search('^paramag_[xyz]$', param): 434 return [-100.0, 100.0]
435 436
437 - def minimise(self, min_algor=None, min_options=None, func_tol=None, grad_tol=None, max_iterations=None, constraints=False, scaling_matrix=None, verbosity=0, sim_index=None, lower=None, upper=None, inc=None):
438 """Minimisation function. 439 440 @param min_algor: The minimisation algorithm to use. 441 @type min_algor: str 442 @param min_options: An array of options to be used by the minimisation algorithm. 443 @type min_options: array of str 444 @param func_tol: The function tolerance which, when reached, terminates optimisation. Setting this to None turns of the check. 445 @type func_tol: None or float 446 @param grad_tol: The gradient tolerance which, when reached, terminates optimisation. Setting this to None turns of the check. 447 @type grad_tol: None or float 448 @param max_iterations: The maximum number of iterations for the algorithm. 449 @type max_iterations: int 450 @param constraints: If True, constraints are used during optimisation. 451 @type constraints: bool 452 @keyword scaling_matrix: The per-model list of diagonal and square scaling matrices. 453 @type scaling_matrix: list of numpy rank-2, float64 array or list of None 454 @param verbosity: A flag specifying the amount of information to print. The higher the value, the greater the verbosity. 455 @type verbosity: int 456 @param sim_index: The index of the simulation to optimise. This should be None if normal optimisation is desired. 457 @type sim_index: None or int 458 @keyword lower: The per-model 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. 459 @type lower: list of lists of numbers 460 @keyword upper: The per-model 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. 461 @type upper: list of lists of numbers 462 @keyword inc: The per-model 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. 463 @type inc: list of lists of int 464 """ 465 466 # Set up the target function for direct calculation. 467 model, param_vector, data_types = target_fn_setup(sim_index=sim_index, scaling_matrix=scaling_matrix[0], verbosity=verbosity) 468 469 # Nothing to do! 470 if not len(param_vector): 471 warn(RelaxWarning("The model has no parameters, minimisation cannot be performed.")) 472 return 473 474 # Right, constraints cannot be used for the 'fixed' model. 475 if constraints and cdp.model == 'fixed': 476 if verbosity: 477 warn(RelaxWarning("Turning constraints off. These cannot be used for the 'fixed' model.")) 478 constraints = False 479 480 # Pop out the Method of Multipliers algorithm. 481 if min_algor == 'Method of Multipliers': 482 min_algor = min_options[0] 483 min_options = min_options[1:] 484 485 # And constraints absolutely must be used for the 'population' model. 486 if not constraints and cdp.model == 'population': 487 warn(RelaxWarning("Turning constraints on. These absolutely must be used for the 'population' model.")) 488 constraints = True 489 490 # Add the Method of Multipliers algorithm. 491 min_options = (min_algor,) + min_options 492 min_algor = 'Method of Multipliers' 493 494 # Disallow Newton optimisation and other Hessian optimisers for the paramagnetic centre position optimisation (the PCS Hessian is not yet implemented). 495 if hasattr(cdp, 'paramag_centre_fixed') and not cdp.paramag_centre_fixed: 496 if min_algor in ['newton']: 497 raise RelaxError("For the paramagnetic centre position, as the Hessians are not yet implemented Newton optimisation cannot be performed.") 498 499 # Linear constraints. 500 A, b = None, None 501 if constraints: 502 A, b = linear_constraints(data_types=data_types, scaling_matrix=scaling_matrix[0]) 503 504 # Grid search. 505 if search('^[Gg]rid', min_algor): 506 # The search. 507 results = grid(func=model.func, args=(), num_incs=inc[0], lower=lower[0], upper=upper[0], A=A, b=b, verbosity=verbosity) 508 509 # Unpack the results. 510 param_vector, func, iter_count, warning = results 511 f_count = iter_count 512 g_count = 0.0 513 h_count = 0.0 514 515 # Minimisation. 516 else: 517 results = generic_minimise(func=model.func, dfunc=model.dfunc, d2func=model.d2func, args=(), x0=param_vector, min_algor=min_algor, min_options=min_options, func_tol=func_tol, grad_tol=grad_tol, maxiter=max_iterations, A=A, b=b, full_output=1, print_flag=verbosity) 518 519 # Unpack the results. 520 if results == None: 521 return 522 param_vector, func, iter_count, f_count, g_count, h_count, warning = results 523 524 # Catch infinite chi-squared values. 525 if isInf(func): 526 raise RelaxInfError('chi-squared') 527 528 # Catch chi-squared values of NaN. 529 if isNaN(func): 530 raise RelaxNaNError('chi-squared') 531 532 # Make a last function call to update the back-calculated RDC and PCS structures to the optimal values. 533 chi2 = model.func(param_vector) 534 535 # Scaling. 536 if scaling_matrix[0] is not None: 537 param_vector = dot(scaling_matrix[0], param_vector) 538 539 # Disassemble the parameter vector. 540 disassemble_param_vector(param_vector=param_vector, data_types=data_types, sim_index=sim_index) 541 542 # Monte Carlo minimisation statistics. 543 if sim_index != None: 544 # Chi-squared statistic. 545 cdp.chi2_sim[sim_index] = func 546 547 # Iterations. 548 cdp.iter_sim[sim_index] = iter_count 549 550 # Function evaluations. 551 cdp.f_count_sim[sim_index] = f_count 552 553 # Gradient evaluations. 554 cdp.g_count_sim[sim_index] = g_count 555 556 # Hessian evaluations. 557 cdp.h_count_sim[sim_index] = h_count 558 559 # Warning. 560 cdp.warning_sim[sim_index] = warning 561 562 # Normal statistics. 563 else: 564 # Chi-squared statistic. 565 cdp.chi2 = func 566 567 # Iterations. 568 cdp.iter = iter_count 569 570 # Function evaluations. 571 cdp.f_count = f_count 572 573 # Gradient evaluations. 574 cdp.g_count = g_count 575 576 # Hessian evaluations. 577 cdp.h_count = h_count 578 579 # Warning. 580 cdp.warning = warning 581 582 # Statistical analysis. 583 if 'rdc' in data_types or 'pcs' in data_types: 584 # Get the final back calculated data (for the Q factor and 585 minimise_bc_data(model, sim_index=sim_index) 586 587 # Calculate the RDC Q factors. 588 if 'rdc' in data_types: 589 rdc.q_factors(sim_index=sim_index, verbosity=verbosity) 590 591 # Calculate the PCS Q factors. 592 if 'pcs' in data_types: 593 pcs.q_factors(sim_index=sim_index, verbosity=verbosity)
594 595
596 - def model_statistics(self, model_info=None, spin_id=None, global_stats=None):
597 """Return the k, n, and chi2 model statistics. 598 599 k - number of parameters. 600 n - number of data points. 601 chi2 - the chi-squared value. 602 603 604 @keyword model_info: The model information from model_loop(). This is unused. 605 @type model_info: None 606 @keyword spin_id: The spin identification string. This is ignored in the N-state model. 607 @type spin_id: None or str 608 @keyword global_stats: A parameter which determines if global or local statistics are returned. For the N-state model, this argument is ignored. 609 @type global_stats: None or bool 610 @return: The optimisation statistics, in tuple format, of the number of parameters (k), the number of data points (n), and the chi-squared value (chi2). 611 @rtype: tuple of (int, int, float) 612 """ 613 614 # Return the values. 615 return param_num(), num_data_points(), cdp.chi2
616 617
618 - def return_data(self, data_id):
619 """Return the base data for the given data ID. 620 621 @keyword data_id: The list of spin ID, data type, and alignment ID, as yielded by the base_data_loop() generator method. 622 @type data_id: list of str 623 @return: The base data. 624 @rtype: list of (float or None) 625 """ 626 627 # Alias the spin or interatomic data container, data type and alignment ID. 628 container = data_id[0] 629 data_type = data_id[1] 630 align_id = data_id[2] 631 632 # The data structure to return. 633 data = [] 634 635 # Skip deselected spins. 636 if data_id[1] == 'pcs' and not container.select: 637 return 638 639 # Return the RDC data. 640 if data_type == 'rdc' and hasattr(container, 'rdc'): 641 if align_id not in container.rdc: 642 data.append(None) 643 else: 644 data.append(container.rdc[align_id]) 645 646 # Return the NOESY data. 647 elif data_type == 'noesy' and hasattr(container, 'noesy'): 648 data.append(container.noesy) 649 650 # Return the PCS data. 651 elif data_id[1] == 'pcs' and hasattr(container, 'pcs'): 652 if align_id not in container.pcs: 653 data.append(None) 654 else: 655 data.append(container.pcs[align_id]) 656 657 # Return the data. 658 return data
659 660
661 - def return_error(self, data_id=None):
662 """Create and return the spin specific Monte Carlo Ri error structure. 663 664 @keyword data_id: The list of spin ID, data type, and alignment ID, as yielded by the base_data_loop() generator method. 665 @type data_id: str 666 @return: The Monte Carlo simulation data errors. 667 @rtype: list of float 668 """ 669 670 # Initialise the MC data structure. 671 mc_errors = [] 672 673 # Alias the spin or interatomic data container. 674 container = data_id[0] 675 676 # Skip deselected spins. 677 if data_id[1] == 'pcs' and not container.select: 678 return 679 680 # RDC data. 681 if data_id[1] == 'rdc' and hasattr(container, 'rdc'): 682 # Do errors exist? 683 if not hasattr(container, 'rdc_err'): 684 raise RelaxError("The RDC errors are missing for the spin pair '%s' and '%s'." % (container.spin_id1, container.spin_id2)) 685 686 # The error. 687 if data_id[2] not in container.rdc_err: 688 err = None 689 else: 690 err = container.rdc_err[data_id[2]] 691 692 # Append the data. 693 mc_errors.append(err) 694 695 # NOESY data. 696 elif data_id[1] == 'noesy' and hasattr(container, 'noesy'): 697 # Do errors exist? 698 if not hasattr(container, 'noesy_err'): 699 raise RelaxError("The NOESY errors are missing for the spin pair '%s' and '%s'." % (container.spin_id1, container.spin_id2)) 700 701 # Append the data. 702 mc_errors.append(container.noesy_err) 703 704 # PCS data. 705 elif data_id[1] == 'pcs' and hasattr(container, 'pcs'): 706 # Do errors exist? 707 if not hasattr(container, 'pcs_err'): 708 raise RelaxError("The PCS errors are missing for spin '%s'." % data_id[0]) 709 710 # The error. 711 if data_id[2] not in container.pcs_err: 712 err = None 713 else: 714 err = container.pcs_err[data_id[2]] 715 716 # Append the data. 717 mc_errors.append(err) 718 719 # Return the errors. 720 return mc_errors
721 722
723 - def set_error(self, index, error, model_info=None):
724 """Set the parameter errors. 725 726 @param index: The index of the parameter to set the errors for. 727 @type index: int 728 @param error: The error value. 729 @type error: float 730 @keyword model_info: The model information from model_loop(). This is unused. 731 @type model_info: None 732 """ 733 734 # Align parameters. 735 names = ['Axx', 'Ayy', 'Axy', 'Axz', 'Ayz'] 736 737 # Alignment tensor parameters. 738 if index < len(cdp.align_ids)*5: 739 # The tensor and parameter index. 740 param_index = index % 5 741 tensor_index = (index - index % 5) / 5 742 743 # Set the error. 744 tensor = align_tensor.return_tensor(index=tensor_index, skip_fixed=True) 745 tensor.set(param=names[param_index], value=error, category='err') 746 747 # Return the object. 748 return getattr(tensor, names[param_index]+'_err')
749 750
751 - def set_param_values(self, param=None, value=None, index=None, spin_id=None, error=False, force=True):
752 """Set the N-state model parameter values. 753 754 @keyword param: The parameter name list. 755 @type param: list of str 756 @keyword value: The parameter value list. 757 @type value: list 758 @keyword index: The index for parameters which are of the list-type (probs, alpha, beta, and gamma). This is ignored for all other types. 759 @type index: None or int 760 @keyword spin_id: The spin identification string (unused). 761 @type spin_id: None 762 @keyword error: A flag which if True will allow the parameter errors to be set instead of the values. 763 @type error: bool 764 @keyword force: A flag which if True will cause current values to be overwritten. If False, a RelaxError will raised if the parameter value is already set. 765 @type force: bool 766 """ 767 768 # Checks. 769 lib.arg_check.is_str_list(param, 'parameter name') 770 lib.arg_check.is_list(value, 'parameter value') 771 772 # Loop over the parameters. 773 for i in range(len(param)): 774 # Is the parameter is valid? 775 if not param[i]: 776 raise RelaxError("The parameter '%s' is not valid for this data pipe type." % param[i]) 777 778 # Error object. 779 if error: 780 param[i] += '_err' 781 782 # Set the indexed parameter. 783 if param[i] in ['probs', 'alpha', 'beta', 'gamma']: 784 obj = getattr(cdp, param[i]) 785 obj[index] = value[i] 786 787 # Set the spin parameters. 788 else: 789 for spin in spin_loop(spin_id): 790 setattr(spin, param[i], value[i])
791 792
793 - def sim_init_values(self):
794 """Initialise the Monte Carlo parameter values.""" 795 796 # Get the minimisation statistic object names. 797 sim_names = self.data_names(set='min') 798 799 # Add the paramagnetic centre, if optimised. 800 if hasattr(cdp, 'paramag_centre_fixed') and not cdp.paramag_centre_fixed: 801 sim_names += ['paramagnetic_centre'] 802 803 # Alignments. 804 if hasattr(cdp, 'align_tensors'): 805 # The parameter names. 806 names = ['Axx', 'Ayy', 'Axy', 'Axz', 'Ayz'] 807 808 # Loop over the alignments, adding the alignment tensor parameters to the tensor data container. 809 for i in range(len(cdp.align_tensors)): 810 # Skip non-optimised tensors. 811 if not opt_uses_tensor(cdp.align_tensors[i]): 812 continue 813 814 # Set up the number of simulations. 815 cdp.align_tensors[i].set_sim_num(cdp.sim_number) 816 817 # Loop over all the parameter names, setting the initial simulation values to those of the parameter value. 818 for object_name in names: 819 for j in range(cdp.sim_number): 820 cdp.align_tensors[i].set(param=object_name, value=deepcopy(getattr(cdp.align_tensors[i], object_name)), category='sim', sim_index=j) 821 822 # Create all other simulation objects. 823 for object_name in sim_names: 824 # Name for the simulation object. 825 sim_object_name = object_name + '_sim' 826 827 # Create the simulation object. 828 setattr(cdp, sim_object_name, []) 829 830 # Get the simulation object. 831 sim_object = getattr(cdp, sim_object_name) 832 833 # Loop over the simulations. 834 for j in range(cdp.sim_number): 835 # Append None to fill the structure. 836 sim_object.append(None) 837 838 # Set the simulation paramagnetic centre positions to the optimised values. 839 if hasattr(cdp, 'paramag_centre_fixed') and not cdp.paramag_centre_fixed: 840 for j in range(cdp.sim_number): 841 cdp.paramagnetic_centre_sim[j] = deepcopy(cdp.paramagnetic_centre)
842 843
844 - def sim_pack_data(self, data_id, sim_data):
845 """Pack the Monte Carlo simulation data. 846 847 @keyword data_id: The list of spin ID, data type, and alignment ID, as yielded by the base_data_loop() generator method. 848 @type data_id: list of str 849 @param sim_data: The Monte Carlo simulation data. 850 @type sim_data: list of float 851 """ 852 853 # Alias the spin or interatomic data container. 854 container = data_id[0] 855 856 # RDC data. 857 if data_id[1] == 'rdc' and hasattr(container, 'rdc'): 858 # Initialise. 859 if not hasattr(container, 'rdc_sim'): 860 container.rdc_sim = {} 861 862 # Store the data structure. 863 container.rdc_sim[data_id[2]] = [] 864 for i in range(cdp.sim_number): 865 container.rdc_sim[data_id[2]].append(sim_data[i][0]) 866 867 # NOESY data. 868 elif data_id[1] == 'noesy' and hasattr(container, 'noesy'): 869 # Store the data structure. 870 container.noesy_sim = [] 871 for i in range(cdp.sim_number): 872 container.noesy_sim[data_id[2]].append(sim_data[i][0]) 873 874 # PCS data. 875 elif data_id[1] == 'pcs' and hasattr(container, 'pcs'): 876 # Initialise. 877 if not hasattr(container, 'pcs_sim'): 878 container.pcs_sim = {} 879 880 # Store the data structure. 881 container.pcs_sim[data_id[2]] = [] 882 for i in range(cdp.sim_number): 883 container.pcs_sim[data_id[2]].append(sim_data[i][0])
884 885
886 - def sim_return_param(self, index, model_info=None):
887 """Return the array of simulation parameter values. 888 889 @param index: The index of the parameter to return the array of values for. 890 @type index: int 891 @keyword model_info: The model information from model_loop(). This is unused. 892 @type model_info: None 893 @return: The array of simulation parameter values. 894 @rtype: list of float 895 """ 896 897 # Align parameters. 898 names = ['Axx', 'Ayy', 'Axy', 'Axz', 'Ayz'] 899 900 # Alignment tensor parameters. 901 if index < align_tensor.num_tensors(skip_fixed=True)*5: 902 # The tensor and parameter index. 903 param_index = index % 5 904 tensor_index = (index - index % 5) / 5 905 906 # Return the simulation parameter array. 907 tensor = align_tensor.return_tensor(index=tensor_index, skip_fixed=True) 908 return getattr(tensor, names[param_index]+'_sim')
909