auto_analyses.relax

49 """The relaxation dispersion auto-analysis.""" 50 51 # Some class variables. 52 opt_func_tol = 1e-25 53 opt_max_iterations = int(1e7) 54

55 - 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):

56 """Perform a full relaxation dispersion analysis for the given list of models. 57 58 @keyword pipe_name: The name of the data pipe containing all of the data for the analysis. 59 @type pipe_name: str 60 @keyword pipe_bundle: The data pipe bundle to associate all spawned data pipes with. 61 @type pipe_bundle: str 62 @keyword results_dir: The directory where results files are saved. 63 @type results_dir: str 64 @keyword models: The list of relaxation dispersion models to optimise. 65 @type models: list of str 66 @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. 67 @type grid_inc: int or None 68 @keyword mc_sim_num: The number of Monte Carlo simulations to be used for error analysis at the end of the analysis. 69 @type mc_sim_num: int 70 @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. 71 @type exp_mc_sim_num: int or None 72 @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'. 73 @type modsel: str 74 @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. 75 @type pre_run_dir: None or str 76 @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. 77 @type optimise_r2eff: bool 78 @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. 79 @type insignificance: float 80 @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. 81 @type numeric_only: bool 82 @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. 83 @type mc_sim_all_models: bool 84 @keyword eliminate: A flag which if True will enable the elimination of failed models and failed Monte Carlo simulations through the eliminate user function. 85 @type eliminate: bool 86 @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. 87 @type set_grid_r20: bool 88 @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. 89 """ 90 91 # Printout. 92 title(file=sys.stdout, text="Relaxation dispersion auto-analysis", prespace=4) 93 94 # Execution lock. 95 status.exec_lock.acquire(pipe_bundle, mode='auto-analysis') 96 97 # Set up the analysis status object. 98 status.init_auto_analysis(pipe_bundle, type='relax_disp') 99 status.current_analysis = pipe_bundle 100 101 # Store the args. 102 self.pipe_name = pipe_name 103 self.pipe_bundle = pipe_bundle 104 self.results_dir = results_dir 105 self.grid_inc = grid_inc 106 self.mc_sim_num = mc_sim_num 107 self.exp_mc_sim_num = exp_mc_sim_num 108 self.models = models 109 self.modsel = modsel 110 self.pre_run_dir = pre_run_dir 111 self.optimise_r2eff = optimise_r2eff 112 self.insignificance = insignificance 113 self.set_grid_r20 = set_grid_r20 114 self.numeric_only = numeric_only 115 self.mc_sim_all_models = mc_sim_all_models 116 self.eliminate = eliminate 117 self.r1_fit = r1_fit 118 119 # No results directory, so default to the current directory. 120 if not self.results_dir: 121 self.results_dir = getcwd() 122 123 # Data checks. 124 self.check_vars() 125 126 # Check for numerical model using numpy version under 1.8. 127 # This will result in slow "for loop" calculation through data, making the analysis 5-6 times slower. 128 self.check_numpy_less_1_8_and_numerical_model() 129 130 # Load the interpreter. 131 self.interpreter = Interpreter(show_script=False, raise_relax_error=True) 132 self.interpreter.populate_self() 133 self.interpreter.on(verbose=False) 134 135 # Execute. 136 try: 137 self.run() 138 139 # Finish and unlock execution. 140 finally: 141 status.auto_analysis[self.pipe_bundle].fin = True 142 status.current_analysis = None 143 status.exec_lock.release()

144 145

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

147 """Determine if the model should be used for model selection. 148 149 @keyword model: The model to check. 150 @type model: str 151 @return: True if the model should be included in the model selection list, False if not. 152 @rtype: bool 153 """ 154 155 # Skip the 'R2eff' base model. 156 if model == MODEL_R2EFF: 157 return False 158 159 # Do not use the analytic models. 160 if self.numeric_only and model in MODEL_LIST_ANALYTIC: 161 return False 162 163 # All models allowed. 164 return True

165 166

167 - def check_vars(self):

168 """Check that the user has set the variables correctly.""" 169 170 # Printout. 171 section(file=sys.stdout, text="Variable checking", prespace=2) 172 173 # The pipe name. 174 if not has_pipe(self.pipe_name): 175 raise RelaxNoPipeError(self.pipe_name) 176 177 # Check the model selection. 178 allowed = ['AIC', 'AICc', 'BIC'] 179 if self.modsel not in allowed: 180 raise RelaxError("The model selection technique '%s' is not in the allowed list of %s." % (self.modsel, allowed)) 181 182 # Some warning for the user if the pure numeric solution is selected. 183 if self.numeric_only: 184 # Loop over all models. 185 for model in self.models: 186 # Skip the models used for nesting. 187 if model in MODEL_LIST_NEST: 188 continue 189 190 # Warnings for all other analytic models. 191 if model in MODEL_LIST_ANALYTIC: 192 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)) 193 194 # Printout. 195 print("The dispersion auto-analysis variables are OK.")

196 197

198 - def check_numpy_less_1_8_and_numerical_model(self):

199 """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.""" 200 201 # Some warning for the user if the pure numeric solution is selected. 202 if float(version.version[:3]) < 1.8: 203 # Store which models are in numeric. 204 models = [] 205 206 # Loop through models. 207 for model in self.models: 208 if model in MODEL_LIST_NUMERIC: 209 models.append(model) 210 211 # Write system message if numerical models is present and numpy version is below 1.8. 212 if len(models) > 0: 213 # Printout. 214 section(file=sys.stdout, text="Numpy version checking for numerical models.", prespace=2) 215 warn(RelaxWarning("Your version of numpy is %s, and below the recommended version of 1.8 for numerical models." % (version.version))) 216 warn(RelaxWarning("Please consider upgrading your numpy version to 1.8.")) 217 218 # Loop over models. 219 for model in models: 220 warn(RelaxWarning("This could make the numerical analysis with model '%s', 5 to 6 times slower." % (model)))

221 222

223 - def error_analysis(self):

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

250 251

252 - def name_pipe(self, prefix):

253 """Generate a unique name for the data pipe. 254 255 @param prefix: The prefix of the data pipe name. 256 @type prefix: str 257 """ 258 259 # The unique pipe name. 260 name = "%s - %s" % (prefix, self.pipe_bundle) 261 262 # Return the name. 263 return name

264 265

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

267 """Support for model nesting. 268 269 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. 270 271 272 @keyword model: The model to be optimised. 273 @type model: str 274 @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. 275 @rtype: bool 276 """ 277 278 # Printout. 279 subsection(file=sys.stdout, text="Nesting and model equivalence checks", prespace=1) 280 281 # The simpler model. 282 model_info, comparable_model_info = nesting_model(self_models=self.models, model=model) 283 if comparable_model_info != None: 284 nested_pipe = self.name_pipe(comparable_model_info.model) 285 else: 286 nested_pipe = None 287 288 # No nesting. 289 if not nested_pipe: 290 print("No model nesting or model equivalence detected.") 291 return False 292 293 # Copying the parameters to a numerical model from an analytic solution. 294 if model_info.eq in [EQ_NUMERIC, EQ_SILICO] and comparable_model_info.eq == EQ_ANALYTIC: 295 analytic = True 296 else: 297 analytic = False 298 299 # Determine if model is equivalent or nested. 300 if model_info.params == comparable_model_info.params: 301 equivalent = True 302 else: 303 equivalent = False 304 305 # Printout. 306 if equivalent: 307 print("Model equivalence detected, copying the optimised parameters from the '%s' model rather than performing a grid search." % comparable_model_info.model) 308 else: 309 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)) 310 if analytic: 311 print("The parameters are copied from a %s model to a %s model." % (comparable_model_info.eq, model_info.eq)) 312 313 # Get the dictionary of how the model parameters of the current model can be copied. 314 par_dic = nesting_param(model_params=model_info.params, nested_model_params=comparable_model_info.params) 315 316 # Loop over the parameters in current model. 317 for param in model_info.params: 318 # Extract how parameter is translated. 319 param_conv = par_dic[param] 320 321 # If the param_conv is None, then continue. 322 if param_conv == None: 323 continue 324 325 print("Copying from parameter '%s' to '%s'." % (param_conv, param)) 326 327 # Loop over the spins to copy the parameters. 328 for spin, spin_id in spin_loop(return_id=True, skip_desel=True): 329 # Get the nested spin. 330 nested_spin = return_spin(spin_id=spin_id, pipe=nested_pipe) 331 332 # Set value. 333 # Some special conversions. 334 if param_conv == '1 - pA': 335 val = 1.0 - getattr(nested_spin, 'pA') 336 337 elif param_conv == '0.0': 338 val = 0.0 339 340 else: 341 val = deepcopy(getattr(nested_spin, param_conv)) 342 343 # Set the attribute. 344 setattr(spin, param, val) 345 346 # 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. 347 return equivalent

348 349

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

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

515 516

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

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

549 550

551 - def run(self):

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

648 649

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

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

774 775

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

777 """Create a set of results, text and Grace files for the current data pipe. 778 779 @keyword path: The directory to place the files into. 780 @type path: str 781 @keyword model: The model tested. 782 @type model: None or str 783 @keyword model_tested: List of models tested, if the pipe is final. 784 @type model_tested: None or list of str. 785 @keyword param: The param to write out. 786 @type param: None or list of str. 787 @keyword file_name_ini: The initial part of the file name for the grace and text files. 788 @type file_name_ini: None or str. 789 """ 790 791 # If not set, use the name of the parameter. 792 if file_name_ini == None: 793 file_name_ini = param 794 795 # If the model is in the list of models which support the parameter. 796 write_result = False 797 if model != None: 798 # Get the model params. 799 model_params = MODEL_PARAMS[model] 800 801 if param in model_params: 802 write_result = True 803 804 # If this is the final pipe, then check if the model has been tested at any time. 805 elif model == None: 806 # Loop through all tested models. 807 for model_tested in models_tested: 808 # If one of the models tested has a parameter which belong in the list of models which support the parameter, then write it out. 809 model_params = MODEL_PARAMS[model_tested] 810 811 if param in model_params: 812 write_result = True 813 break 814 815 # Write results if some of the models supports the parameter. 816 if write_result: 817 self.interpreter.value.write(param=param, file='%s.out'%file_name_ini, dir=path, force=True) 818 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