auto_analyses.relax

53 """The relaxation dispersion auto-analysis.""" 54 55 # Some class variables. 56 opt_func_tol = 1e-25 57 opt_max_iterations = int(1e7) 58

59 - def __init__(self, pipe_name=None, pipe_bundle=None, results_dir=None, models=[MODEL_R2EFF], grid_inc=11, mc_sim_num=500, exp_mc_sim_num=None, modsel='AIC', pre_run_dir=None, optimise_r2eff=False, insignificance=0.0, numeric_only=False, mc_sim_all_models=False, eliminate=True, set_grid_r20=False, r1_fit=False):

60 """Perform a full relaxation dispersion analysis for the given list of models. 61 62 @keyword pipe_name: The name of the data pipe containing all of the data for the analysis. 63 @type pipe_name: str 64 @keyword pipe_bundle: The data pipe bundle to associate all spawned data pipes with. 65 @type pipe_bundle: str 66 @keyword results_dir: The directory where results files are saved. 67 @type results_dir: str 68 @keyword models: The list of relaxation dispersion models to optimise. 69 @type models: list of str 70 @keyword grid_inc: Number of grid search increments. If set to None, then the grid search will be turned off and the default parameter values will be used instead. 71 @type grid_inc: int or None 72 @keyword mc_sim_num: The number of Monte Carlo simulations to be used for error analysis at the end of the analysis. 73 @type mc_sim_num: int 74 @keyword exp_mc_sim_num: The number of Monte Carlo simulations for the error analysis in the 'R2eff' model when exponential curves are fitted. This defaults to the value of the mc_sim_num argument when not given. When set to '-1', the R2eff errors are estimated from the Covariance matrix. For the 2-point fixed-time calculation for the 'R2eff' model, this argument is ignored. 75 @type exp_mc_sim_num: int or None 76 @keyword modsel: The model selection technique to use in the analysis to determine which model is the best for each spin cluster. This can currently be one of 'AIC', 'AICc', and 'BIC'. 77 @type modsel: str 78 @keyword pre_run_dir: The optional directory containing the dispersion auto-analysis results from a previous run. The optimised parameters from these previous results will be used as the starting point for optimisation rather than performing a grid search. This is essential for when large spin clusters are specified, as a grid search becomes prohibitively expensive with clusters of three or more spins. At some point a RelaxError will occur because the grid search is impossibly large. For the cluster specific parameters, i.e. the populations of the states and the exchange parameters, an average value will be used as the starting point. For all other parameters, the R20 values for each spin and magnetic field, as well as the parameters related to the chemical shift difference dw, the optimised values of the previous run will be directly copied. 79 @type pre_run_dir: None or str 80 @keyword optimise_r2eff: Flag to specify if the read previous R2eff results should be optimised. For R1rho models where the error of R2eff values are determined by Monte-Carlo simulations, it can be valuable to make an initial R2eff run with a high number of Monte-Carlo simulations. Any subsequent model analysis can then be based on these R2eff values, without optimising the R2eff values. 81 @type optimise_r2eff: bool 82 @keyword insignificance: The R2eff/R1rho value in rad/s by which to judge insignificance. If the maximum difference between two points on all dispersion curves for a spin is less than this value, that spin will be deselected. This does not affect the 'No Rex' model. Set this value to 0.0 to use all data. The value will be passed on to the relax_disp.insignificance user function. 83 @type insignificance: float 84 @keyword numeric_only: The class of models to use in the model selection. The default of False allows all dispersion models to be used in the analysis (no exchange, the analytic models and the numeric models). The value of True will activate a pure numeric solution - the analytic models will be optimised, as they are very useful for replacing the grid search for the numeric models, but the final model selection will not include them. 85 @type numeric_only: bool 86 @keyword mc_sim_all_models: A flag which if True will cause Monte Carlo simulations to be performed for each individual model. Otherwise Monte Carlo simulations will be reserved for the final model. 87 @type mc_sim_all_models: bool 88 @keyword eliminate: A flag which if True will enable the elimination of failed models and failed Monte Carlo simulations through the eliminate user function. 89 @type eliminate: bool 90 @keyword set_grid_r20: A flag which if True will set the grid R20 values from the minimum R2eff values through the r20_from_min_r2eff user function. This will speed up the grid search with a factor GRID_INC^(Nr_spec_freq). For a CPMG experiment with two fields and standard GRID_INC=21, the speed-up is a factor 441. 91 @type set_grid_r20: bool 92 @keyword r1_fit: A flag which if True will activate R1 parameter fitting via relax_disp.r1_fit for the models that support it. If False, then the relax_disp.r1_fit user function will not be called. 93 """ 94 95 # Initial printout. 96 title(file=sys.stdout, text="Relaxation dispersion auto-analysis", prespace=4) 97 98 # Safely execute the full protocol. 99 try: 100 # Execution lock. 101 status.exec_lock.acquire(pipe_bundle, mode='auto-analysis') 102 103 # Set up the analysis status object. 104 status.init_auto_analysis(pipe_bundle, type='relax_disp') 105 status.current_analysis = pipe_bundle 106 107 # Store the args. 108 self.pipe_name = pipe_name 109 self.pipe_bundle = pipe_bundle 110 self.results_dir = results_dir 111 self.grid_inc = grid_inc 112 self.mc_sim_num = mc_sim_num 113 self.exp_mc_sim_num = exp_mc_sim_num 114 self.models = models 115 self.modsel = modsel 116 self.pre_run_dir = pre_run_dir 117 self.optimise_r2eff = optimise_r2eff 118 self.insignificance = insignificance 119 self.set_grid_r20 = set_grid_r20 120 self.numeric_only = numeric_only 121 self.mc_sim_all_models = mc_sim_all_models 122 self.eliminate = eliminate 123 self.r1_fit = r1_fit 124 125 # No results directory, so default to the current directory. 126 if not self.results_dir: 127 self.results_dir = getcwd() 128 129 # Data checks. 130 self.check_vars() 131 132 # Check for numerical model using numpy version under 1.8. 133 # This will result in slow "for loop" calculation through data, making the analysis 5-6 times slower. 134 self.check_numpy_less_1_8_and_numerical_model() 135 136 # Load the interpreter. 137 self.interpreter = Interpreter(show_script=False, raise_relax_error=True) 138 self.interpreter.populate_self() 139 self.interpreter.on(verbose=False) 140 141 # Execute. 142 self.run() 143 144 # Clean up. 145 finally: 146 # Final printout. 147 title(file=sys.stdout, text="Completion of the relaxation dispersion auto-analysis", prespace=4) 148 print_elapsed_time(time() - status.start_time) 149 150 # Finish and unlock execution. 151 status.auto_analysis[self.pipe_bundle].fin = True 152 status.current_analysis = None 153 status.exec_lock.release()

154 155

156 - def is_model_for_selection(self, model=None):

157 """Determine if the model should be used for model selection. 158 159 @keyword model: The model to check. 160 @type model: str 161 @return: True if the model should be included in the model selection list, False if not. 162 @rtype: bool 163 """ 164 165 # Skip the 'R2eff' base model. 166 if model == MODEL_R2EFF: 167 return False 168 169 # Do not use the analytic models. 170 if self.numeric_only and model in MODEL_LIST_ANALYTIC: 171 return False 172 173 # All models allowed. 174 return True

175 176

177 - def check_vars(self):

178 """Check that the user has set the variables correctly.""" 179 180 # Printout. 181 section(file=sys.stdout, text="Variable checking", prespace=2) 182 183 # The pipe name. 184 if not has_pipe(self.pipe_name): 185 raise RelaxNoPipeError(self.pipe_name) 186 187 # Check the model selection. 188 allowed = ['AIC', 'AICc', 'BIC'] 189 if self.modsel not in allowed: 190 raise RelaxError("The model selection technique '%s' is not in the allowed list of %s." % (self.modsel, allowed)) 191 192 # Some warning for the user if the pure numeric solution is selected. 193 if self.numeric_only: 194 # Loop over all models. 195 for model in self.models: 196 # Skip the models used for nesting. 197 if model in MODEL_LIST_NEST: 198 continue 199 200 # Warnings for all other analytic models. 201 if model in MODEL_LIST_ANALYTIC: 202 warn(RelaxWarning("The analytic model '%s' will be optimised but will not be used in any way in this numeric model only auto-analysis." % model)) 203 204 # Printout. 205 print("The dispersion auto-analysis variables are OK.")

206 207

208 - def check_numpy_less_1_8_and_numerical_model(self):

209 """Check for numerical model using numpy version under 1.8. This will result in slow "for loop" calculation through data, making the analysis 5-6 times slower.""" 210 211 # Some warning for the user if the pure numeric solution is selected. 212 if version_comparison(version.version, '1.8') < 0: 213 # Store which models are in numeric. 214 models = [] 215 216 # Loop through models. 217 for model in self.models: 218 if model in MODEL_LIST_NUMERIC: 219 models.append(model) 220 221 # Write system message if numerical models is present and numpy version is below 1.8. 222 if len(models) > 0: 223 # Printout. 224 section(file=sys.stdout, text="Numpy version checking for numerical models.", prespace=2) 225 warn(RelaxWarning("Your version of numpy is %s, and below the recommended version of 1.8 for numerical models." % (version.version))) 226 warn(RelaxWarning("Please consider upgrading your numpy version to 1.8.")) 227 228 # Loop over models. 229 for model in models: 230 warn(RelaxWarning("This could make the numerical analysis with model '%s', 5 to 6 times slower." % (model)))

231 232

233 - def error_analysis(self):

234 """Perform an error analysis of the peak intensities for each field strength separately.""" 235 236 # Printout. 237 section(file=sys.stdout, text="Error analysis", prespace=2) 238 239 # Check if intensity errors have already been calculated by the user. 240 precalc = True 241 for spin in spin_loop(skip_desel=True): 242 # No structure. 243 if not hasattr(spin, 'peak_intensity_err'): 244 precalc = False 245 break 246 247 # Determine if a spectrum ID is missing from the list. 248 for id in cdp.spectrum_ids: 249 if id not in spin.peak_intensity_err: 250 precalc = False 251 break 252 253 # Skip. 254 if precalc: 255 print("Skipping the error analysis as it has already been performed.") 256 return 257 258 # Perform the error analysis. 259 self.interpreter.spectrum.error_analysis_per_field()

260 261

262 - def name_pipe(self, prefix):

263 """Generate a unique name for the data pipe. 264 265 @param prefix: The prefix of the data pipe name. 266 @type prefix: str 267 """ 268 269 # The unique pipe name. 270 name = "%s - %s" % (prefix, self.pipe_bundle) 271 272 # Return the name. 273 return name

274 275

276 - def nesting(self, model=None):

277 """Support for model nesting. 278 279 If model nesting is detected, the optimised parameters from the simpler model will be used for the more complex model. The method will then signal if the nesting condition is met for the model, allowing the grid search to be skipped. 280 281 282 @keyword model: The model to be optimised. 283 @type model: str 284 @return: True if the model parameters is equivalent to the nested model, and all parameters are copied. False if none or some of the parameters have been translated from the nested model. Here the Grid search should still be performed. 285 @rtype: bool 286 """ 287 288 # Printout. 289 subsection(file=sys.stdout, text="Nesting and model equivalence checks", prespace=1) 290 291 # The simpler model. 292 model_info, comparable_model_info = nesting_model(self_models=self.models, model=model) 293 if comparable_model_info != None: 294 nested_pipe = self.name_pipe(comparable_model_info.model) 295 else: 296 nested_pipe = None 297 298 # No nesting. 299 if not nested_pipe: 300 print("No model nesting or model equivalence detected.") 301 return False 302 303 # Copying the parameters to a numerical model from an analytic solution. 304 if model_info.eq in [EQ_NUMERIC, EQ_SILICO] and comparable_model_info.eq == EQ_ANALYTIC: 305 analytic = True 306 else: 307 analytic = False 308 309 # Determine if model is equivalent or nested. 310 if model_info.params == comparable_model_info.params: 311 equivalent = True 312 else: 313 equivalent = False 314 315 # Printout. 316 if equivalent: 317 print("Model equivalence detected, copying the optimised parameters from the '%s' model rather than performing a grid search." % comparable_model_info.model) 318 else: 319 print("Model nesting detected, translating the optimised parameters %s from the '%s' model to the parameters %s of model '%s'. A grid search is issued for the remaining parameters." % (comparable_model_info.params, comparable_model_info.model, model_info.params, model)) 320 if analytic: 321 print("The parameters are copied from a %s model to a %s model." % (comparable_model_info.eq, model_info.eq)) 322 323 # Get the dictionary of how the model parameters of the current model can be copied. 324 par_dic = nesting_param(model_params=model_info.params, nested_model_params=comparable_model_info.params) 325 326 # Loop over the parameters in current model. 327 for param in model_info.params: 328 # Extract how parameter is translated. 329 param_conv = par_dic[param] 330 331 # If the param_conv is None, then continue. 332 if param_conv == None: 333 continue 334 335 print("Copying from parameter '%s' to '%s'." % (param_conv, param)) 336 337 # Loop over the spins to copy the parameters. 338 for spin, spin_id in spin_loop(return_id=True, skip_desel=True): 339 # Get the nested spin. 340 nested_spin = return_spin(spin_id=spin_id, pipe=nested_pipe) 341 342 # Set value. 343 # Some special conversions. 344 if param_conv == '1 - pA': 345 val = 1.0 - getattr(nested_spin, 'pA') 346 347 elif param_conv == '0.0': 348 val = 0.0 349 350 else: 351 val = deepcopy(getattr(nested_spin, param_conv)) 352 353 # Set the attribute. 354 setattr(spin, param, val) 355 356 # Determine if model is equivalent, and should not be Grid searched, or if nested, and some parameters are pre-set. Here Grid search should still be issued. 357 return equivalent

358 359

360 - def optimise(self, model=None, model_path=None):

361 """Optimise the model, taking model nesting into account. 362 363 @keyword model: The model to be optimised. 364 @type model: str 365 @keyword model_path: The folder name for the model, where possible spaces has been replaced with underscore. 366 @type model_path: str 367 """ 368 369 # Printout. 370 section(file=sys.stdout, text="Optimisation", prespace=2) 371 372 # Deselect insignificant spins. 373 if model not in [MODEL_R2EFF, MODEL_NOREX]: 374 self.interpreter.relax_disp.insignificance(level=self.insignificance) 375 376 # Speed-up grid-search by using minium R2eff value. 377 if self.set_grid_r20 and model != MODEL_R2EFF: 378 self.interpreter.relax_disp.r20_from_min_r2eff(force=True) 379 380 # Use pre-run results as the optimisation starting point. 381 # Test if file exists. 382 if self.pre_run_dir: 383 path = self.pre_run_dir + sep + model_path 384 # File path. 385 file_path = get_file_path('results', path) 386 387 # Test if the file exists and determine the compression type. 388 try: 389 compress_type, file_path = determine_compression(file_path) 390 res_file_exists = True 391 392 except RelaxFileError: 393 res_file_exists = False 394 395 if self.pre_run_dir and res_file_exists: 396 self.pre_run_parameters(model=model, model_path=model_path) 397 398 # Otherwise use the normal nesting check and grid search if not nested. 399 else: 400 # Nested model simplification. 401 nested = self.nesting(model=model) 402 403 # Otherwise use a grid search of default values to start optimisation with. 404 if not nested: 405 # Grid search. 406 if self.grid_inc: 407 self.interpreter.minimise.grid_search(inc=self.grid_inc) 408 409 # Default values. 410 else: 411 # The standard parameters. 412 for param in MODEL_PARAMS[model]: 413 self.interpreter.value.set(param=param, index=None) 414 415 # The optional R1 parameter. 416 if is_r1_optimised(model=model): 417 self.interpreter.value.set(param='r1', index=None) 418 419 # 'R2eff' model minimisation flags. 420 do_minimise = False 421 if model == MODEL_R2EFF: 422 # The constraints flag. 423 constraints = False 424 425 # The minimisation algorithm to use. 426 # Both the Jacobian and Hessian matrix has been specified for exponential curve-fitting, allowing for the much faster algorithms to be used. 427 min_algor = 'Newton' 428 429 # Check if all spins contains 'r2eff and it associated error. 430 has_r2eff = False 431 432 # Loop over all spins. 433 for cur_spin, spin_id in spin_loop(return_id=True, skip_desel=True): 434 # Check 'r2eff' 435 if hasattr(cur_spin, 'r2eff') and hasattr(cur_spin, 'r2eff_err'): 436 has_r2eff = True 437 else: 438 has_r2eff = False 439 break 440 441 # Skip optimisation, if 'r2eff' + 'r2eff_err' is present and flag for forcing optimisation is not raised. 442 if has_r2eff and not self.optimise_r2eff: 443 pass 444 445 # Do optimisation, if 'r2eff' + 'r2eff_err' is present and flag for forcing optimisation is raised. 446 elif has_r2eff and self.optimise_r2eff: 447 do_minimise = True 448 449 # Optimise, if no R2eff and error is present. 450 elif not has_r2eff: 451 do_minimise = True 452 453 # Dispersion model minimisation flags. 454 else: 455 do_minimise = True 456 constraints = True 457 # The minimisation algorithm to use. If the Jacobian and Hessian matrix have not been specified for fitting, 'simplex' should be used. 458 min_algor = 'simplex' 459 460 # Do the minimisation. 461 if do_minimise: 462 self.interpreter.minimise.execute(min_algor=min_algor, func_tol=self.opt_func_tol, max_iter=self.opt_max_iterations, constraints=constraints) 463 464 # Model elimination. 465 if self.eliminate: 466 self.interpreter.eliminate() 467 468 # Monte Carlo simulations. 469 do_monte_carlo = False 470 if model == MODEL_R2EFF: 471 # The constraints flag. 472 constraints = False 473 474 # Both the Jacobian and Hessian matrix has been specified for exponential curve-fitting, allowing for the much faster algorithms to be used. 475 min_algor = 'Newton' 476 477 # Skip optimisation, if 'r2eff' + 'r2eff_err' is present and flag for forcing optimisation is not raised. 478 if has_r2eff and not self.optimise_r2eff: 479 pass 480 481 # Do optimisation, if 'r2eff' + 'r2eff_err' is present and flag for forcing optimisation is raised. 482 elif has_r2eff and self.optimise_r2eff: 483 do_monte_carlo = True 484 485 # Optimise, if no R2eff and error is present. 486 elif not has_r2eff: 487 do_monte_carlo = True 488 489 elif self.mc_sim_all_models or len(self.models) < 2: 490 do_monte_carlo = True 491 # The constraints flag. 492 constraints = True 493 # The minimisation algorithm to use. If the Jacobian and Hessian matrix have not been specified for fitting, 'simplex' should be used. 494 min_algor = 'simplex' 495 496 # Error estimation by Monte Carlo simulations. 497 if do_monte_carlo: 498 # Set the number of Monte-Carlo simulations. 499 monte_carlo_sim = self.mc_sim_num 500 501 # If the number for exponential curve fitting has been set. 502 if model == MODEL_R2EFF and self.exp_mc_sim_num != None: 503 monte_carlo_sim = self.exp_mc_sim_num 504 505 # When set to minus 1, estimation of the errors will be extracted from the covariance matrix. 506 # This is HIGHLY likely to be wrong, but can be used in an initial test fase. 507 if model == MODEL_R2EFF and self.exp_mc_sim_num == -1: 508 # Print 509 subsection(file=sys.stdout, text="Estimating errors from Covariance matrix", prespace=1) 510 511 # Raise warning. 512 text = 'Estimating errors from the Covariance matrix is highly likely to be "quite" wrong. Use only with extreme care, and for initial rapid testing of your data.' 513 warn(RelaxWarning(text)) 514 515 # Estimate errors 516 self.interpreter.relax_disp.r2eff_err_estimate() 517 else: 518 self.interpreter.monte_carlo.setup(number=monte_carlo_sim) 519 self.interpreter.monte_carlo.create_data() 520 self.interpreter.monte_carlo.initial_values() 521 self.interpreter.minimise.execute(min_algor=min_algor, func_tol=self.opt_func_tol, max_iter=self.opt_max_iterations, constraints=constraints) 522 if self.eliminate: 523 self.interpreter.eliminate() 524 self.interpreter.monte_carlo.error_analysis()

525 526

527 - def pre_run_parameters(self, model=None, model_path=None):

528 """Copy parameters from an earlier analysis. 529 530 @keyword model: The model to be optimised. 531 @type model: str 532 @keyword model_path: The folder name for the model, where possible spaces has been replaced with underscore. 533 @type model_path: str 534 """ 535 536 # Printout. 537 subsection(file=sys.stdout, text="Pre-run parameters", prespace=1) 538 539 # The data pipe name. 540 pipe_name = self.name_pipe('pre') 541 542 # Create a temporary data pipe for the previous run. 543 self.interpreter.pipe.create(pipe_name=pipe_name, pipe_type='relax_disp') 544 545 # Load the previous results. 546 path = self.pre_run_dir + sep + model_path 547 self.interpreter.results.read(file='results', dir=path) 548 549 # Force copy of the R2eff values. 550 if model == MODEL_R2EFF: 551 self.interpreter.value.copy(pipe_from=pipe_name, pipe_to=self.name_pipe(model), param='r2eff', force=True) 552 553 # Copy the parameters. 554 self.interpreter.relax_disp.parameter_copy(pipe_from=pipe_name, pipe_to=self.name_pipe(model)) 555 556 # Finally, switch back to the original data pipe and delete the temporary one. 557 self.interpreter.pipe.switch(pipe_name=self.name_pipe(model)) 558 self.interpreter.pipe.delete(pipe_name=pipe_name)

559 560

561 - def run(self):

562 """Execute the auto-analysis.""" 563 564 # Peak intensity error analysis. 565 if MODEL_R2EFF in self.models: 566 self.error_analysis() 567 568 # R1 parameter fitting. 569 if self.r1_fit: 570 subtitle(file=sys.stdout, text="R1 parameter optimisation activation", prespace=3) 571 self.interpreter.relax_disp.r1_fit(fit=self.r1_fit) 572 else: 573 # No print out. 574 self.interpreter.relax_disp.r1_fit(fit=self.r1_fit) 575 576 # Loop over the models. 577 self.model_pipes = [] 578 for model in self.models: 579 # Printout. 580 subtitle(file=sys.stdout, text="The '%s' model" % model, prespace=3) 581 582 # The results directory path. 583 model_path = model.replace(" ", "_") 584 path = self.results_dir+sep+model_path 585 586 # The name of the data pipe for the model. 587 model_pipe = self.name_pipe(model) 588 if self.is_model_for_selection(model): 589 self.model_pipes.append(model_pipe) 590 591 # Check that results do not already exist - i.e. a previous run was interrupted. 592 path1 = path + sep + 'results' 593 path2 = path1 + '.bz2' 594 path3 = path1 + '.gz' 595 if access(path1, F_OK) or access(path2, F_OK) or access(path2, F_OK): 596 # Printout. 597 print("Detected the presence of results files for the '%s' model - loading these instead of performing optimisation for a second time." % model) 598 599 # Create a data pipe and switch to it. 600 self.interpreter.pipe.create(pipe_name=model_pipe, pipe_type='relax_disp', bundle=self.pipe_bundle) 601 self.interpreter.pipe.switch(model_pipe) 602 603 # Load the results. 604 self.interpreter.results.read(file='results', dir=path) 605 606 # Jump to the next model. 607 continue 608 609 # Create the data pipe by copying the base pipe, then switching to it. 610 self.interpreter.pipe.copy(pipe_from=self.pipe_name, pipe_to=model_pipe, bundle_to=self.pipe_bundle) 611 self.interpreter.pipe.switch(model_pipe) 612 613 # Select the model. 614 self.interpreter.relax_disp.select_model(model) 615 616 # Copy the R2eff values from the R2eff model data pipe. 617 if model != MODEL_R2EFF and MODEL_R2EFF in self.models: 618 self.interpreter.value.copy(pipe_from=self.name_pipe(MODEL_R2EFF), pipe_to=model_pipe, param='r2eff') 619 620 # Calculate the R2eff values for the fixed relaxation time period data types. 621 if model == MODEL_R2EFF and not has_exponential_exp_type(): 622 self.interpreter.minimise.calculate() 623 624 # Optimise the model. 625 else: 626 self.optimise(model=model, model_path=model_path) 627 628 # Write out the results. 629 self.write_results(path=path, model=model) 630 631 # The final model selection data pipe. 632 if len(self.models) >= 2: 633 # Printout. 634 section(file=sys.stdout, text="Final results", prespace=2) 635 636 # Perform model selection. 637 self.interpreter.model_selection(method=self.modsel, modsel_pipe=self.name_pipe('final'), bundle=self.pipe_bundle, pipes=self.model_pipes) 638 639 # Final Monte Carlo simulations only. 640 if not self.mc_sim_all_models: 641 self.interpreter.monte_carlo.setup(number=self.mc_sim_num) 642 self.interpreter.monte_carlo.create_data() 643 self.interpreter.monte_carlo.initial_values() 644 self.interpreter.minimise.execute('simplex', func_tol=self.opt_func_tol, max_iter=self.opt_max_iterations, constraints=True) 645 if self.eliminate: 646 self.interpreter.eliminate() 647 self.interpreter.monte_carlo.error_analysis() 648 649 # Writing out the final results. 650 self.write_results(path=self.results_dir+sep+'final') 651 652 # No model selection. 653 else: 654 warn(RelaxWarning("Model selection in the dispersion auto-analysis has been skipped as only %s models have been optimised." % len(self.model_pipes))) 655 656 # Finally save the program state. 657 self.interpreter.state.save(state='final_state', dir=self.results_dir, force=True)

658 659

660 - def write_results(self, path=None, model=None):

661 """Create a set of results, text and Grace files for the current data pipe. 662 663 @keyword path: The directory to place the files into. 664 @type path: str 665 """ 666 667 # Printout. 668 section(file=sys.stdout, text="Results writing", prespace=2) 669 670 # If this is the final model selection round, check which models have been tested. 671 if model == None: 672 models_tested = [] 673 for spin, spin_id in spin_loop(return_id=True, skip_desel=True): 674 spin_model = spin.model 675 676 # Add to list, if not in already. 677 if spin_model not in models_tested: 678 models_tested.append(spin_model) 679 else: 680 models_tested = None 681 682 # Special for R2eff model. 683 if model == MODEL_R2EFF: 684 # The R2eff parameter. 685 self.interpreter.value.write(param='r2eff', file='r2eff.out', dir=path, force=True) 686 self.interpreter.grace.write(x_data_type='res_num', y_data_type='r2eff', file='r2eff.agr', dir=path, force=True) 687 688 # Exponential curves. 689 if has_exponential_exp_type(): 690 self.interpreter.relax_disp.plot_exp_curves(file='intensities.agr', dir=path, force=True) # Average peak intensities. 691 self.interpreter.relax_disp.plot_exp_curves(file='intensities_norm.agr', dir=path, force=True, norm=True) # Average peak intensities (normalised). 692 693 # The I0 parameter. 694 self.interpreter.value.write(param='i0', file='i0.out', dir=path, force=True) 695 self.interpreter.grace.write(x_data_type='res_num', y_data_type='i0', file='i0.agr', dir=path, force=True) 696 697 # Dispersion curves. 698 self.interpreter.relax_disp.plot_disp_curves(dir=path, force=True) 699 self.interpreter.relax_disp.write_disp_curves(dir=path, force=True) 700 701 # The selected models for the final run. 702 if model == None: 703 self.interpreter.value.write(param='model', file='model.out', dir=path, force=True) 704 705 # For CPMG models. 706 if has_cpmg_exp_type(): 707 # The R20 parameter. 708 self.write_results_test(path=path, model=model, models_tested=models_tested, param='r2', file_name_ini='r20') 709 710 # The R20A and R20B parameters. 711 self.write_results_test(path=path, model=model, models_tested=models_tested, param='r2a', file_name_ini='r20a') 712 self.write_results_test(path=path, model=model, models_tested=models_tested, param='r2b', file_name_ini='r20b') 713 714 # For R1ho models. 715 if has_r1rho_exp_type(): 716 # The R1 parameter. 717 self.write_results_test(path=path, model=model, models_tested=models_tested, param='r1') 718 719 # The R1rho prime parameter. 720 self.write_results_test(path=path, model=model, models_tested=models_tested, param='r2', file_name_ini='r1rho_prime') 721 722 # Plot specific R1rho graphs. 723 if model in [None] + MODEL_LIST_R1RHO: 724 self.interpreter.relax_disp.plot_disp_curves(dir=path, x_axis=X_AXIS_THETA, force=True) 725 self.interpreter.relax_disp.plot_disp_curves(dir=path, y_axis=Y_AXIS_R2_R1RHO, x_axis=X_AXIS_W_EFF, force=True) 726 self.interpreter.relax_disp.plot_disp_curves(dir=path, y_axis=Y_AXIS_R2_EFF, x_axis=X_AXIS_THETA, interpolate=INTERPOLATE_OFFSET, force=True) 727 728 # The calculation of theta and w_eff parameter in R1rho experiments. 729 if model in MODEL_LIST_R1RHO_FULL: 730 self.interpreter.value.write(param='theta', file='theta.out', dir=path, force=True) 731 self.interpreter.value.write(param='w_eff', file='w_eff.out', dir=path, force=True) 732 733 # The pA and pB parameters. 734 self.write_results_test(path=path, model=model, models_tested=models_tested, param='pA') 735 self.write_results_test(path=path, model=model, models_tested=models_tested, param='pB') 736 737 # The pC parameter. 738 self.write_results_test(path=path, model=model, models_tested=models_tested, param='pC') 739 740 # The phi_ex parameter. 741 self.write_results_test(path=path, model=model, models_tested=models_tested, param='phi_ex') 742 743 # The phi_ex_B nd phi_ex_C parameters. 744 self.write_results_test(path=path, model=model, models_tested=models_tested, param='phi_ex_B') 745 self.write_results_test(path=path, model=model, models_tested=models_tested, param='phi_ex_C') 746 747 # The dw parameter. 748 self.write_results_test(path=path, model=model, models_tested=models_tested, param='dw') 749 750 # The dw_AB, dw_BC and dw_AC parameter. 751 self.write_results_test(path=path, model=model, models_tested=models_tested, param='dw_AB') 752 self.write_results_test(path=path, model=model, models_tested=models_tested, param='dw_BC') 753 self.write_results_test(path=path, model=model, models_tested=models_tested, param='dw_AC') 754 755 # The dwH parameter. 756 self.write_results_test(path=path, model=model, models_tested=models_tested, param='dwH') 757 758 # The dwH_AB, dwH_BC and dwH_AC parameter. 759 self.write_results_test(path=path, model=model, models_tested=models_tested, param='dwH_AB') 760 self.write_results_test(path=path, model=model, models_tested=models_tested, param='dwH_BC') 761 self.write_results_test(path=path, model=model, models_tested=models_tested, param='dwH_AC') 762 763 # The k_AB, kex and tex parameters. 764 self.write_results_test(path=path, model=model, models_tested=models_tested, param='k_AB') 765 self.write_results_test(path=path, model=model, models_tested=models_tested, param='kex') 766 self.write_results_test(path=path, model=model, models_tested=models_tested, param='tex') 767 768 # The kex_AB, kex_BC, kex_AC parameters. 769 self.write_results_test(path=path, model=model, models_tested=models_tested, param='kex_AB') 770 self.write_results_test(path=path, model=model, models_tested=models_tested, param='kex_BC') 771 self.write_results_test(path=path, model=model, models_tested=models_tested, param='kex_AC') 772 773 # The kB and kC parameters. 774 self.write_results_test(path=path, model=model, models_tested=models_tested, param='kB') 775 self.write_results_test(path=path, model=model, models_tested=models_tested, param='kC') 776 777 # Minimisation statistics. 778 if not (model == MODEL_R2EFF and has_fixed_time_exp_type()): 779 self.interpreter.value.write(param='chi2', file='chi2.out', dir=path, force=True) 780 self.interpreter.grace.write(y_data_type='chi2', file='chi2.agr', dir=path, force=True) 781 782 # Finally save the results. This is last to allow the continuation of an interrupted analysis while ensuring that all results files have been created. 783 self.interpreter.results.write(file='results', dir=path, force=True)

784 785

786 - def write_results_test(self, path=None, model=None, models_tested=None, param=None, file_name_ini=None):

787 """Create a set of results, text and Grace files for the current data pipe. 788 789 @keyword path: The directory to place the files into. 790 @type path: str 791 @keyword model: The model tested. 792 @type model: None or str 793 @keyword model_tested: List of models tested, if the pipe is final. 794 @type model_tested: None or list of str. 795 @keyword param: The param to write out. 796 @type param: None or list of str. 797 @keyword file_name_ini: The initial part of the file name for the grace and text files. 798 @type file_name_ini: None or str. 799 """ 800 801 # If not set, use the name of the parameter. 802 if file_name_ini == None: 803 file_name_ini = param 804 805 # Normal model. 806 write_result = False 807 if model != None: 808 if param in Model_class(model=model).params: 809 write_result = True 810 811 # If this is the final pipe, then check if the model has been tested at any time. 812 elif model == None: 813 for model_tested in models_tested: 814 if param in Model_class(model=model_tested).params: 815 write_result = True 816 break 817 818 # Write results if some of the models supports the parameter. 819 if write_result: 820 self.interpreter.value.write(param=param, file='%s.out'%file_name_ini, dir=path, force=True) 821 self.interpreter.grace.write(x_data_type='res_num', y_data_type=param, file='%s.agr'%file_name_ini, dir=path, force=True)

Source Code for Module auto_analyses.relax_disp