auto_analyses.dauvergne

122 """The model-free auto-analysis.""" 123 124 # Some class variables. 125 opt_func_tol = 1e-25 126 opt_max_iterations = int(1e7) 127

128 - def __init__(self, pipe_name=None, pipe_bundle=None, results_dir=None, write_results_dir=None, diff_model=None, mf_models=['m0', 'm1', 'm2', 'm3', 'm4', 'm5', 'm6', 'm7', 'm8', 'm9'], local_tm_models=['tm0', 'tm1', 'tm2', 'tm3', 'tm4', 'tm5', 'tm6', 'tm7', 'tm8', 'tm9'], grid_inc=11, diff_tensor_grid_inc={'sphere': 11, 'prolate': 11, 'oblate': 11, 'ellipsoid': 6}, min_algor='newton', mc_sim_num=500, max_iter=None, user_fns=None, conv_loop=True):

129 """Perform the full model-free analysis protocol of d'Auvergne and Gooley, 2008b. 130 131 @keyword pipe_name: The name of the data pipe containing the sequence info. This data pipe should have all values set including the CSA value, the bond length, the heteronucleus name and proton name. It should also have all relaxation data loaded. 132 @type pipe_name: str 133 @keyword pipe_bundle: The data pipe bundle to associate all spawned data pipes with. 134 @type pipe_bundle: str 135 @keyword results_dir: The directory where optimisation results will read from. Results will also be saved to this directory if the write_results_dir argument is not given. 136 @type results_dir: str 137 @keyword write_results_dir: The directory where optimisation results will be saved in. If None, it will default to the value of the results_dir argument. This is mainly used for debugging. 138 @type write_results_dir: str or None 139 @keyword diff_model: The global diffusion model to optimise. This can be one of 'local_tm', 'sphere', 'oblate', 'prolate', 'ellipsoid', or 'final'. If all or a subset of these are supplied as a list, then these will be automatically looped over and calculated. 140 @type diff_model: str or list of str 141 @keyword mf_models: The model-free models. 142 @type mf_models: list of str 143 @keyword local_tm_models: The model-free models. 144 @type local_tm_models: list of str 145 @keyword grid_inc: The grid search size (the number of increments per dimension). 146 @type grid_inc: int 147 @keyword diff_tensor_grid_inc: A list of grid search sizes for the optimisation of the sphere, prolate spheroid, oblate spheroid, and ellipsoid. 148 @type diff_tensor_grid_inc: list of int 149 @keyword min_algor: The minimisation algorithm (in most cases this should not be changed). 150 @type min_algor: str 151 @keyword mc_sim_num: The number of Monte Carlo simulations to be used for error analysis at the end of the analysis. 152 @type mc_sim_num: int 153 @keyword max_iter: The maximum number of iterations for the global iteration. Set to None, then the algorithm iterates until convergence. 154 @type max_iter: int or None. 155 @keyword user_fns: A dictionary of replacement user functions. These will overwrite the standard user functions. The key should be the name of the user function or user function class and the value should be the function or class instance. 156 @type user_fns: dict 157 @keyword conv_loop: Automatic looping over all rounds until convergence. 158 @type conv_loop: bool 159 """ 160 161 # Initial printout. 162 title(file=sys.stdout, text="d'Auvergne protocol model-free auto-analysis") 163 164 # Safely execute the full protocol. 165 try: 166 # Execution lock. 167 status.exec_lock.acquire(pipe_bundle, mode='auto-analysis') 168 169 # Store the args. 170 self.pipe_name = pipe_name 171 self.pipe_bundle = pipe_bundle 172 self.mf_models = mf_models 173 self.local_tm_models = local_tm_models 174 self.grid_inc = grid_inc 175 self.diff_tensor_grid_inc = diff_tensor_grid_inc 176 self.min_algor = min_algor 177 self.mc_sim_num = mc_sim_num 178 self.max_iter = max_iter 179 self.conv_loop = conv_loop 180 181 # The model-free data pipe names. 182 self.mf_model_pipes = [] 183 for i in range(len(self.mf_models)): 184 self.mf_model_pipes.append(self.name_pipe(self.mf_models[i])) 185 self.local_tm_model_pipes = [] 186 for i in range(len(self.local_tm_models)): 187 self.local_tm_model_pipes.append(self.name_pipe(self.local_tm_models[i])) 188 189 # The diffusion models. 190 if isinstance(diff_model, list): 191 self.diff_model_list = diff_model 192 else: 193 self.diff_model_list = [diff_model] 194 195 # Project directory (i.e. directory containing the model-free model results and the newly generated files) 196 if results_dir: 197 self.results_dir = results_dir + sep 198 else: 199 self.results_dir = getcwd() + sep 200 if write_results_dir: 201 self.write_results_dir = write_results_dir + sep 202 else: 203 self.write_results_dir = self.results_dir 204 205 # Data checks. 206 self.check_vars() 207 208 # Set the data pipe to the current data pipe. 209 if self.pipe_name != cdp_name(): 210 switch(self.pipe_name) 211 212 # Some info for the status. 213 self.status_setup() 214 215 # Load the interpreter. 216 self.interpreter = Interpreter(show_script=False, raise_relax_error=True) 217 self.interpreter.populate_self() 218 self.interpreter.on(verbose=False) 219 220 # Replacement user functions. 221 if user_fns: 222 for name in user_fns: 223 setattr(self.interpreter, name, user_fns[name]) 224 225 # Execute the protocol. 226 for self.diff_model in self.diff_model_list: 227 # Wait a little while between diffusion models. 228 sleep(1) 229 230 # Set the global model name. 231 status.auto_analysis[self.pipe_bundle].diff_model = self.diff_model 232 233 # Initialise the convergence data structures. 234 self.conv_data = Container() 235 self.conv_data.chi2 = [] 236 self.conv_data.models = [] 237 self.conv_data.diff_vals = [] 238 if self.diff_model == 'sphere': 239 self.conv_data.diff_params = ['tm'] 240 elif self.diff_model == 'oblate' or self.diff_model == 'prolate': 241 self.conv_data.diff_params = ['tm', 'Da', 'theta', 'phi'] 242 elif self.diff_model == 'ellipsoid': 243 self.conv_data.diff_params = ['tm', 'Da', 'Dr', 'alpha', 'beta', 'gamma'] 244 self.conv_data.spin_ids = [] 245 self.conv_data.mf_params = [] 246 self.conv_data.mf_vals = [] 247 248 # Execute the analysis for each diffusion model. 249 self.execute() 250 251 # Clean up. 252 finally: 253 # Final printout. 254 title(file=sys.stdout, text="Completion of the d'Auvergne protocol model-free auto-analysis") 255 print_elapsed_time(time() - status.start_time) 256 257 # Finish and unlock execution. 258 status.auto_analysis[self.pipe_bundle].fin = True 259 status.current_analysis = None 260 status.exec_lock.release()

261 262

263 - def check_vars(self):

264 """Check that the user has set the variables correctly.""" 265 266 # Printout. 267 subtitle(file=sys.stdout, text="Auto-analysis variable checking") 268 269 # The pipe bundle. 270 if not isinstance(self.pipe_bundle, str): 271 raise RelaxError("The pipe bundle name '%s' is invalid." % self.pipe_bundle) 272 273 # The diff model. 274 valid_models = ['local_tm', 'sphere', 'oblate', 'prolate', 'ellipsoid', 'final'] 275 for i in range(len(self.diff_model_list)): 276 if self.diff_model_list[i] not in valid_models: 277 raise RelaxError("The diff_model value '%s' is incorrectly set. It must be one of %s." % (self.diff_model_list[i], valid_models)) 278 279 # Model-free models. 280 mf_models = ['m0', 'm1', 'm2', 'm3', 'm4', 'm5', 'm6', 'm7', 'm8', 'm9'] 281 local_tm_models = ['tm0', 'tm1', 'tm2', 'tm3', 'tm4', 'tm5', 'tm6', 'tm7', 'tm8', 'tm9'] 282 if not isinstance(self.mf_models, list): 283 raise RelaxError("The self.mf_models user variable must be a list.") 284 if not isinstance(self.local_tm_models, list): 285 raise RelaxError("The self.local_tm_models user variable must be a list.") 286 for i in range(len(self.mf_models)): 287 if self.mf_models[i] not in mf_models: 288 raise RelaxError("The self.mf_models user variable '%s' is incorrectly set. It must be one of %s." % (self.mf_models, mf_models)) 289 for i in range(len(self.local_tm_models)): 290 if self.local_tm_models[i] not in local_tm_models: 291 raise RelaxError("The self.local_tm_models user variable '%s' is incorrectly set. It must be one of %s." % (self.local_tm_models, local_tm_models)) 292 293 # Sequence data. 294 if not exists_mol_res_spin_data(): 295 raise RelaxNoSequenceError(self.pipe_name) 296 297 # Relaxation data. 298 if not hasattr(cdp, 'ri_ids') or len(cdp.ri_ids) == 0: 299 raise RelaxNoRiError(ri_id) 300 301 # Insufficient data. 302 if len(cdp.ri_ids) <= 3: 303 raise RelaxError("Insufficient relaxation data, 4 or more data sets are essential for the execution of this script.") 304 305 # Spin vars. 306 for spin, spin_id in spin_loop(return_id=True): 307 # Skip deselected spins. 308 if not spin.select: 309 continue 310 311 # Test if the isotope type has been set. 312 if not hasattr(spin, 'isotope') or spin.isotope == None: 313 raise RelaxNoValueError("nuclear isotope type", spin_id=spin_id) 314 315 # Skip spins with no relaxation data. 316 if not hasattr(spin, 'ri_data') or spin.ri_data == None: 317 continue 318 319 # Test if the CSA value has been set. 320 if not hasattr(spin, 'csa') or spin.csa == None: 321 raise RelaxNoValueError("CSA", spin_id=spin_id) 322 323 # Interatomic vars. 324 for interatom in interatomic_loop(): 325 # Get the corresponding spins. 326 spin1 = return_spin(spin_hash=interatom._spin_hash1) 327 spin2 = return_spin(spin_hash=interatom._spin_hash2) 328 329 # Skip deselected spins. 330 if not spin1.select or not spin2.select: 331 continue 332 333 # Test if the interatomic distance has been set. 334 if not hasattr(interatom, 'r') or interatom.r == None: 335 raise RelaxNoValueError("interatomic distance", spin_id=interatom.spin_id1, spin_id2=interatom.spin_id2) 336 337 # Min vars. 338 if not isinstance(self.grid_inc, int): 339 raise RelaxError("The grid_inc user variable '%s' is incorrectly set. It should be an integer." % self.grid_inc) 340 if not isinstance(self.diff_tensor_grid_inc, dict): 341 raise RelaxError("The diff_tensor_grid_inc user variable '%s' is incorrectly set. It should be a dictionary." % self.diff_tensor_grid_inc) 342 for tensor in ['sphere', 'prolate', 'oblate', 'ellipsoid']: 343 if not tensor in self.diff_tensor_grid_inc: 344 raise RelaxError("The diff_tensor_grid_inc user variable '%s' is incorrectly set. It should contain the '%s' key." % (self.diff_tensor_grid_inc, tensor)) 345 if not isinstance(self.diff_tensor_grid_inc[tensor], int): 346 raise RelaxError("The diff_tensor_grid_inc user variable '%s' is incorrectly set. The value corresponding to the key '%s' should be an integer." % (self.diff_tensor_grid_inc, tensor)) 347 if not isinstance(self.min_algor, str): 348 raise RelaxError("The min_algor user variable '%s' is incorrectly set. It should be a string." % self.min_algor) 349 if not isinstance(self.mc_sim_num, int): 350 raise RelaxError("The mc_sim_num user variable '%s' is incorrectly set. It should be an integer." % self.mc_sim_num) 351 352 # Looping. 353 if not isinstance(self.conv_loop, bool): 354 raise RelaxError("The conv_loop user variable '%s' is incorrectly set. It should be one of the booleans True or False." % self.conv_loop)

355 356

357 - def convergence(self):

358 """Test for the convergence of the global model.""" 359 360 # Printout. 361 title(file=sys.stdout, text="Convergence testing") 362 363 # Maximum number of iterations reached. 364 if self.max_iter and self.round > self.max_iter: 365 print("Maximum number of global iterations reached. Terminating the protocol before convergence has been reached.") 366 return True 367 368 # Store the data of the current data pipe. 369 self.conv_data.chi2.append(cdp.chi2) 370 371 # Create a string representation of the model-free models of the current data pipe. 372 curr_models = '' 373 for spin in spin_loop(): 374 if hasattr(spin, 'model'): 375 if not spin.model == 'None': 376 curr_models = curr_models + spin.model 377 self.conv_data.models.append(curr_models) 378 379 # Store the diffusion tensor parameters. 380 self.conv_data.diff_vals.append([]) 381 for param in self.conv_data.diff_params: 382 # Get the parameter values. 383 self.conv_data.diff_vals[-1].append(getattr(cdp.diff_tensor, param)) 384 385 # Store the model-free parameters. 386 self.conv_data.mf_vals.append([]) 387 self.conv_data.mf_params.append([]) 388 self.conv_data.spin_ids.append([]) 389 for spin, spin_id in spin_loop(return_id=True): 390 # Skip spin systems with no 'params' object. 391 if not hasattr(spin, 'params'): 392 continue 393 394 # Add the spin ID, parameters, and empty value list. 395 self.conv_data.spin_ids[-1].append(spin_id) 396 self.conv_data.mf_params[-1].append([]) 397 self.conv_data.mf_vals[-1].append([]) 398 399 # Loop over the parameters. 400 for j in range(len(spin.params)): 401 # Get the parameters and values. 402 self.conv_data.mf_params[-1][-1].append(spin.params[j]) 403 self.conv_data.mf_vals[-1][-1].append(getattr(spin, spin.params[j].lower())) 404 405 # No need for tests. 406 if self.round == 1: 407 print("First round of optimisation, skipping the convergence tests.\n\n\n") 408 return False 409 410 # Loop over the iterations. 411 converged = False 412 for i in range(self.start_round, self.round - 1): 413 # Print out. 414 print("\n\n\n# Comparing the current iteration to iteration %i.\n" % (i+1)) 415 416 # Index. 417 index = i - self.start_round 418 419 # Chi-squared test. 420 print("Chi-squared test:") 421 print(" chi2 (iter %i): %s" % (i+1, self.conv_data.chi2[index])) 422 print(" (as an IEEE-754 byte array: %s)" % floatAsByteArray(self.conv_data.chi2[index])) 423 print(" chi2 (iter %i): %s" % (self.round, self.conv_data.chi2[-1])) 424 print(" (as an IEEE-754 byte array: %s)" % floatAsByteArray(self.conv_data.chi2[-1])) 425 print(" chi2 (difference): %s" % (self.conv_data.chi2[index] - self.conv_data.chi2[-1])) 426 if self.conv_data.chi2[index] == self.conv_data.chi2[-1]: 427 print(" The chi-squared value has converged.\n") 428 else: 429 print(" The chi-squared value has not converged.\n") 430 continue 431 432 # Identical model-free model test. 433 print("Identical model-free models test:") 434 if self.conv_data.models[index] == self.conv_data.models[-1]: 435 print(" The model-free models have converged.\n") 436 else: 437 print(" The model-free models have not converged.\n") 438 continue 439 440 # Identical diffusion tensor parameter value test. 441 print("Identical diffusion tensor parameter test:") 442 params_converged = True 443 for k in range(len(self.conv_data.diff_params)): 444 # Test if not identical. 445 if self.conv_data.diff_vals[index][k] != self.conv_data.diff_vals[-1][k]: 446 print(" Parameter: %s" % param) 447 print(" Value (iter %i): %s" % (i+1, self.conv_data.diff_vals[index][k])) 448 print(" (as an IEEE-754 byte array: %s)" % floatAsByteArray(self.conv_data.diff_vals[index][k])) 449 print(" Value (iter %i): %s" % (self.round, self.conv_data.diff_vals[-1][k])) 450 print(" (as an IEEE-754 byte array: %s)" % floatAsByteArray(self.conv_data.diff_vals[-1][k])) 451 print(" The diffusion parameters have not converged.\n") 452 params_converged = False 453 break 454 if not params_converged: 455 continue 456 print(" The diffusion tensor parameters have converged.\n") 457 458 # Identical model-free parameter value test. 459 print("\nIdentical model-free parameter test:") 460 if len(self.conv_data.spin_ids[index]) != len(self.conv_data.spin_ids[-1]): 461 print(" Different number of spins.") 462 continue 463 for j in range(len(self.conv_data.spin_ids[-1])): 464 # Loop over the parameters. 465 for k in range(len(self.conv_data.mf_params[-1][j])): 466 # Test if not identical. 467 if self.conv_data.mf_vals[index][j][k] != self.conv_data.mf_vals[-1][j][k]: 468 print(" Spin ID: %s" % self.conv_data.spin_ids[-1][j]) 469 print(" Parameter: %s" % self.conv_data.mf_params[-1][j][k]) 470 print(" Value (iter %i): %s" % (i+1, self.conv_data.mf_vals[index][j][k])) 471 print(" (as an IEEE-754 byte array: %s)" % floatAsByteArray(self.conv_data.mf_vals[index][j][k])) 472 print(" Value (iter %i): %s" % (self.round, self.conv_data.mf_vals[-1][j][k])) 473 print(" (as an IEEE-754 byte array: %s)" % floatAsByteArray(self.conv_data.mf_vals[index][j][k])) 474 print(" The model-free parameters have not converged.\n") 475 params_converged = False 476 break 477 if not params_converged: 478 continue 479 print(" The model-free parameters have converged.\n") 480 481 # Convergence. 482 converged = True 483 break 484 485 486 # Final printout. 487 ################## 488 489 print("\nConvergence:") 490 if converged: 491 # Update the status. 492 status.auto_analysis[self.pipe_bundle].convergence = True 493 494 # Print out. 495 print(" [ Yes ]") 496 497 # Return the termination condition. 498 return True 499 else: 500 # Print out. 501 print(" [ No ]") 502 503 # Return False to not terminate. 504 return False

505 506

507 - def determine_rnd(self, model=None):

508 """Function for returning the name of next round of optimisation.""" 509 510 # The base model directory. 511 base_dir = self.results_dir+sep+model 512 513 # Printout. 514 subtitle(file=sys.stdout, text="Determining the next round of optimisation") 515 print("%-30s %s" % ("Base model directory:", base_dir)) 516 517 # Catch if a file exists with the name of the directory. 518 if not isdir(base_dir) and access(base_dir, F_OK): 519 raise RelaxError("The base model directory '%s' is not usable as a file with the same name already exists." % base_dir) 520 521 # If no directory exists, set the round to 'init' or 0. 522 if not isdir(base_dir): 523 print("%-30s %i" % ("Round:", 0)) 524 print("The base directory does not exist.") 525 return 0 526 527 # Is the directory readable, writable, and executable. 528 if not access(base_dir, R_OK): 529 raise RelaxError("The base model directory '%s' is not readable." % base_dir) 530 if not access(base_dir, X_OK): 531 raise RelaxError("The base model directory '%s' is not executable." % base_dir) 532 533 # Get a list of all files in the directory model. 534 dir_list = listdir(base_dir) 535 536 # Set the round to 'init' or 0 if there is no directory called 'init'. 537 if 'init' not in dir_list: 538 print("%-30s %i" % ("Round:", 0)) 539 print("No 'init' directory present in the base directory.") 540 return 0 541 542 # Create a list of all files which begin with 'round_'. 543 rnd_dirs = [] 544 for file in dir_list: 545 if search('^round_', file): 546 rnd_dirs.append(file) 547 548 # Create a sorted list of integer round numbers. 549 numbers = [] 550 for dir in rnd_dirs: 551 try: 552 numbers.append(int(dir[6:])) 553 except: 554 pass 555 numbers.sort() 556 557 # No directories beginning with 'round_' exist, set the round to 1. 558 if not len(numbers): 559 print("%-30s %i" % ("Round:", 1)) 560 print("No directories beginning with 'round_' exist.") 561 return 1 562 563 # The highest number. 564 max_round = numbers[-1] 565 566 # Check that the opt/results file exists for the round (working backwards). 567 for i in range(max_round, -1, -1): 568 # Assume the round is complete. 569 complete_round = i 570 571 # The file root. 572 file_root = base_dir + sep + "round_%i" % i + sep + 'opt' + sep + 'results' 573 574 # Stop looping when the opt/results file is found. 575 if access(file_root + '.bz2', F_OK): 576 break 577 if access(file_root + '.gz', F_OK): 578 break 579 if access(file_root, F_OK): 580 break 581 582 # No round, so assume the initial state. 583 if complete_round == 0: 584 print("%-30s %i" % ("Round:", 0)) 585 print("No opt/results files can be found.") 586 return 0 587 588 # Determine the number for the next round (add 1 to the highest completed round). 589 print("%-30s %i" % ("Round:", complete_round + 1)) 590 return complete_round + 1

591 592

593 - def execute(self):

594 """Execute the protocol.""" 595 596 # MI - Local tm. 597 ################ 598 599 if self.diff_model == 'local_tm': 600 # Printout. 601 title(file=sys.stdout, text="Model MI - Local tm") 602 603 # Base directory to place files into. 604 self.base_dir = self.results_dir+'local_tm'+sep 605 606 # Sequential optimisation of all model-free models (function must be modified to suit). 607 self.multi_model(local_tm=True) 608 609 # Model selection. 610 self.model_selection(modsel_pipe=self.name_pipe('aic'), dir=self.base_dir + 'aic') 611 612 613 # Diffusion models MII to MV. 614 ############################# 615 616 elif self.diff_model == 'sphere' or self.diff_model == 'prolate' or self.diff_model == 'oblate' or self.diff_model == 'ellipsoid': 617 # Printout. 618 if self.diff_model == 'sphere': 619 title(file=sys.stdout, text="Model MII - Spherical diffusion") 620 elif self.diff_model == 'prolate': 621 title(file=sys.stdout, text="Model MIII - Prolate spheroidal diffusion") 622 elif self.diff_model == 'oblate': 623 title(file=sys.stdout, text="Model MIV - Oblate spheroidal diffusion") 624 elif self.diff_model == 'ellipsoid': 625 title(file=sys.stdout, text="Model MV - Ellipsoidal diffusion") 626 627 # No local_tm directory! 628 dir_list = listdir(self.results_dir) 629 if 'local_tm' not in dir_list: 630 raise RelaxError("The local_tm model must be optimised first.") 631 632 # The initial round of optimisation - not zero if calculations were interrupted. 633 self.start_round = self.determine_rnd(model=self.diff_model) 634 635 # Loop until convergence if conv_loop is set, otherwise just loop once. 636 # This looping could be made much cleaner by removing the dependence on the determine_rnd() function. 637 while True: 638 # Determine which round of optimisation to do (init, round_1, round_2, etc). 639 self.round = self.determine_rnd(model=self.diff_model) 640 status.auto_analysis[self.pipe_bundle].round = self.round 641 642 # Inital round of optimisation for diffusion models MII to MV. 643 if self.round == 0: 644 # Printout. 645 subtitle(file=sys.stdout, text="Initial round of optimisation") 646 647 # Base directory to place files into. 648 self.base_dir = self.results_dir+self.diff_model+sep+'init'+sep 649 650 # Run name. 651 name = self.name_pipe(self.diff_model) 652 653 # Create the data pipe (deleting the old one if it exists). 654 if has_pipe(name): 655 self.interpreter.pipe.delete(name) 656 self.interpreter.pipe.create(name, 'mf', bundle=self.pipe_bundle) 657 658 # Load the local tm diffusion model MI results. 659 self.interpreter.results.read(file='results', dir=self.results_dir+'local_tm'+sep+'aic') 660 661 # Remove the tm parameter. 662 self.interpreter.model_free.remove_tm() 663 664 # Initialise the diffusion tensor. 665 if self.diff_model == 'sphere': 666 self.interpreter.diffusion_tensor.init(None, fixed=False) 667 inc = self.diff_tensor_grid_inc['sphere'] 668 elif self.diff_model == 'prolate': 669 self.interpreter.diffusion_tensor.init((None, None, None, None), spheroid_type='prolate', fixed=False) 670 inc = self.diff_tensor_grid_inc['prolate'] 671 elif self.diff_model == 'oblate': 672 self.interpreter.diffusion_tensor.init((None, None, None, None), spheroid_type='oblate', fixed=False) 673 inc = self.diff_tensor_grid_inc['oblate'] 674 elif self.diff_model == 'ellipsoid': 675 self.interpreter.diffusion_tensor.init((None, None, None, None, None, None), fixed=False) 676 inc = self.diff_tensor_grid_inc['ellipsoid'] 677 678 # Minimise just the diffusion tensor. 679 self.interpreter.fix('all_spins') 680 self.interpreter.minimise.grid_search(inc=inc) 681 self.interpreter.minimise.execute(self.min_algor, func_tol=self.opt_func_tol, max_iter=self.opt_max_iterations) 682 683 # Write the results. 684 self.interpreter.results.write(file='results', dir=self.base_dir, force=True) 685 686 687 # Normal round of optimisation for diffusion models MII to MV. 688 else: 689 # Printout. 690 subtitle(file=sys.stdout, text="Round %i of optimisation" % self.round) 691 692 # Base directory to place files into. 693 self.base_dir = self.results_dir+self.diff_model + sep+'round_'+repr(self.round)+sep 694 695 # Load the optimised diffusion tensor from either the previous round. 696 self.load_tensor() 697 698 # Sequential optimisation of all model-free models (function must be modified to suit). 699 self.multi_model() 700 701 # Model selection. 702 self.model_selection(modsel_pipe=self.name_pipe('aic'), dir=self.base_dir + 'aic') 703 704 # Final optimisation of all diffusion and model-free parameters. 705 self.interpreter.fix('all', fixed=False) 706 707 # Minimise all parameters. 708 self.interpreter.minimise.execute(self.min_algor, func_tol=self.opt_func_tol, max_iter=self.opt_max_iterations) 709 710 # Write the results. 711 dir = self.base_dir + 'opt' 712 self.interpreter.results.write(file='results', dir=dir, force=True) 713 714 # Test for convergence. 715 converged = self.convergence() 716 717 # Break out of the infinite while loop if automatic looping is not activated or if convergence has occurred. 718 if converged or not self.conv_loop: 719 break 720 721 # Unset the status. 722 status.auto_analysis[self.pipe_bundle].round = None 723 724 725 # Final run. 726 ############ 727 728 elif self.diff_model == 'final': 729 # Printout. 730 title(file=sys.stdout, text="Final run") 731 732 # Diffusion model selection. 733 ############################ 734 735 # Printout. 736 subtitle(file=sys.stdout, text="Diffusion model selection") 737 738 # The contents of the results directory. 739 dir_list = listdir(self.results_dir) 740 741 # Check that the minimal set of global diffusion models required for the protocol has been optimised. 742 min_models = ['local_tm', 'sphere'] 743 for model in min_models: 744 if model not in dir_list: 745 raise RelaxError("The minimum set of global diffusion models required for the protocol have not been optimised, the '%s' model results cannot be found." % model) 746 747 # Build a list of all global diffusion models optimised. 748 all_models = ['local_tm', 'sphere', 'prolate', 'oblate', 'ellipsoid'] 749 self.opt_models = [] 750 self.pipes = [] 751 for model in all_models: 752 if model in dir_list: 753 self.opt_models.append(model) 754 self.pipes.append(self.name_pipe(model)) 755 756 # Remove all temporary pipes used in this auto-analysis. 757 for name in pipe_names(bundle=self.pipe_bundle): 758 if name in self.pipes + self.mf_model_pipes + self.local_tm_model_pipes + [self.name_pipe('aic'), self.name_pipe('previous')]: 759 self.interpreter.pipe.delete(name) 760 761 # Create the local_tm data pipe. 762 self.interpreter.pipe.create(self.name_pipe('local_tm'), 'mf', bundle=self.pipe_bundle) 763 764 # Load the local tm diffusion model MI results. 765 self.interpreter.results.read(file='results', dir=self.results_dir+'local_tm'+sep+'aic') 766 767 # Loop over models MII to MV. 768 for model in ['sphere', 'prolate', 'oblate', 'ellipsoid']: 769 # Skip missing models. 770 if model not in self.opt_models: 771 continue 772 773 # Determine which was the last round of optimisation for each of the models. 774 self.round = self.determine_rnd(model=model) - 1 775 776 # If no directories begining with 'round_' exist, the script has not been properly utilised! 777 if self.round < 1: 778 # Construct the name of the diffusion tensor. 779 name = model 780 if model == 'prolate' or model == 'oblate': 781 name = name + ' spheroid' 782 783 # Throw an error to prevent misuse of the script. 784 raise RelaxError("Multiple rounds of optimisation of the " + name + " (between 8 to 15) are required for the proper execution of this script.") 785 786 # Create the data pipe. 787 self.interpreter.pipe.create(self.name_pipe(model), 'mf', bundle=self.pipe_bundle) 788 789 # Load the diffusion model results. 790 self.interpreter.results.read(file='results', dir=self.results_dir+model + sep+'round_'+repr(self.round)+sep+'opt') 791 792 # Model selection between MI to MV. 793 self.model_selection(modsel_pipe=self.name_pipe('final'), write_flag=False) 794 795 796 # Monte Carlo simulations. 797 ########################## 798 799 # Printout. 800 subtitle(file=sys.stdout, text="Monte Carlo simulations") 801 802 # Fix the diffusion tensor, if it exists. 803 if hasattr(get_pipe(self.name_pipe('final')), 'diff_tensor'): 804 self.interpreter.fix('diff') 805 806 # Simulations. 807 self.interpreter.monte_carlo.setup(number=self.mc_sim_num) 808 self.interpreter.monte_carlo.create_data() 809 self.interpreter.monte_carlo.initial_values() 810 self.interpreter.minimise.execute(self.min_algor, func_tol=self.opt_func_tol, max_iter=self.opt_max_iterations) 811 self.interpreter.eliminate() 812 self.interpreter.monte_carlo.error_analysis() 813 814 815 # Write the final results. 816 ########################## 817 818 # Printout. 819 subtitle(file=sys.stdout, text="Writing the final results") 820 821 # Create results files and plots of the data. 822 self.write_results() 823 824 825 # Unknown script behaviour. 826 ########################### 827 828 else: 829 raise RelaxError("Unknown diffusion model, change the value of 'self.diff_model'")

830 831

832 - def load_tensor(self):

833 """Function for loading the optimised diffusion tensor.""" 834 835 # Create the data pipe for the previous data (deleting the old data pipe first if necessary). 836 if has_pipe(self.name_pipe('previous')): 837 self.interpreter.pipe.delete(self.name_pipe('previous')) 838 self.interpreter.pipe.create(self.name_pipe('previous'), 'mf', bundle=self.pipe_bundle) 839 840 # Load the optimised diffusion tensor from the initial round. 841 if self.round == 1: 842 self.interpreter.results.read('results', self.results_dir+self.diff_model + sep+'init') 843 844 # Load the optimised diffusion tensor from the previous round. 845 else: 846 self.interpreter.results.read('results', self.results_dir+self.diff_model + sep+'round_'+repr(self.round-1)+sep+'opt')

847 848

849 - def model_selection(self, modsel_pipe=None, dir=None, write_flag=True):

850 """Model selection function.""" 851 852 # Model selection (delete the model selection pipe if it already exists). 853 if has_pipe(modsel_pipe): 854 self.interpreter.pipe.delete(modsel_pipe) 855 self.interpreter.model_selection(method='AIC', modsel_pipe=modsel_pipe, bundle=self.pipe_bundle, pipes=self.pipes) 856 857 # Write the results. 858 if write_flag: 859 self.interpreter.results.write(file='results', dir=dir, force=True)

860 861

862 - def multi_model(self, local_tm=False):

863 """Function for optimisation of all model-free models.""" 864 865 # Set the data pipe names (also the names of preset model-free models). 866 if local_tm: 867 models = self.local_tm_models 868 self.pipes = self.local_tm_models 869 else: 870 models = self.mf_models 871 self.pipes = [] 872 for i in range(len(models)): 873 self.pipes.append(self.name_pipe(models[i])) 874 875 # Loop over the data pipes. 876 for i in range(len(models)): 877 # Place the model name into the status container. 878 status.auto_analysis[self.pipe_bundle].current_model = models[i] 879 880 # Create the data pipe (by copying). 881 if has_pipe(self.pipes[i]): 882 self.interpreter.pipe.delete(self.pipes[i]) 883 self.interpreter.pipe.copy(self.pipe_name, self.pipes[i], bundle_to=self.pipe_bundle) 884 self.interpreter.pipe.switch(self.pipes[i]) 885 886 # Copy the diffusion tensor from the 'opt' data pipe and prevent it from being minimised. 887 if not local_tm: 888 self.interpreter.diffusion_tensor.copy(self.name_pipe('previous')) 889 self.interpreter.fix('diff') 890 891 # Select the model-free model. 892 self.interpreter.model_free.select_model(model=models[i]) 893 894 # Minimise. 895 self.interpreter.minimise.grid_search(inc=self.grid_inc) 896 self.interpreter.minimise.execute(self.min_algor, func_tol=self.opt_func_tol, max_iter=self.opt_max_iterations) 897 898 # Model elimination. 899 self.interpreter.eliminate() 900 901 # Write the results. 902 dir = self.base_dir + models[i] 903 self.interpreter.results.write(file='results', dir=dir, force=True) 904 905 # Unset the status. 906 status.auto_analysis[self.pipe_bundle].current_model = None

907 908

909 - def name_pipe(self, prefix):

910 """Generate a unique name for the data pipe. 911 912 @param prefix: The prefix of the data pipe name. 913 @type prefix: str 914 """ 915 916 # The unique pipe name. 917 name = "%s - %s" % (prefix, self.pipe_bundle) 918 919 # Return the name. 920 return name

921 922

923 - def status_setup(self):

924 """Initialise the status object.""" 925 926 # Initialise the status object for this auto-analysis. 927 status.init_auto_analysis(self.pipe_bundle, type='dauvergne_protocol') 928 status.current_analysis = self.pipe_bundle 929 930 # The global diffusion model. 931 status.auto_analysis[self.pipe_bundle].diff_model = None 932 933 # The round of optimisation, i.e. the global iteration. 934 status.auto_analysis[self.pipe_bundle].round = None 935 936 # The list of model-free local tm models for optimisation, i.e. the global iteration. 937 status.auto_analysis[self.pipe_bundle].local_tm_models = self.local_tm_models 938 939 # The list of model-free models for optimisation, i.e. the global iteration. 940 status.auto_analysis[self.pipe_bundle].mf_models = self.mf_models 941 942 # The current model-free model. 943 status.auto_analysis[self.pipe_bundle].current_model = None 944 945 # The maximum number of iterations of the global model. 946 status.auto_analysis[self.pipe_bundle].max_iter = self.max_iter 947 948 # The convergence of the global model. 949 status.auto_analysis[self.pipe_bundle].convergence = False

950 951

952 - def write_results(self):

953 """Create Grace plots of the final model-free results.""" 954 955 # Save the results file. 956 dir = self.write_results_dir + 'final' 957 self.interpreter.results.write(file='results', dir=dir, force=True) 958 959 # The Grace plots. 960 dir = self.write_results_dir + 'final' + sep + 'grace' 961 self.interpreter.grace.write(x_data_type='res_num', y_data_type='s2', file='s2.agr', dir=dir, force=True) 962 self.interpreter.grace.write(x_data_type='res_num', y_data_type='s2f', file='s2f.agr', dir=dir, force=True) 963 self.interpreter.grace.write(x_data_type='res_num', y_data_type='s2s', file='s2s.agr', dir=dir, force=True) 964 self.interpreter.grace.write(x_data_type='res_num', y_data_type='te', file='te.agr', dir=dir, force=True) 965 self.interpreter.grace.write(x_data_type='res_num', y_data_type='tf', file='tf.agr', dir=dir, force=True) 966 self.interpreter.grace.write(x_data_type='res_num', y_data_type='ts', file='ts.agr', dir=dir, force=True) 967 self.interpreter.grace.write(x_data_type='res_num', y_data_type='rex', file='rex.agr', dir=dir, force=True) 968 self.interpreter.grace.write(x_data_type='s2', y_data_type='te', file='s2_vs_te.agr', dir=dir, force=True) 969 self.interpreter.grace.write(x_data_type='s2', y_data_type='rex', file='s2_vs_rex.agr', dir=dir, force=True) 970 self.interpreter.grace.write(x_data_type='te', y_data_type='rex', file='te_vs_rex.agr', dir=dir, force=True) 971 972 # Write the values to text files. 973 dir = self.write_results_dir + 'final' 974 self.interpreter.value.write(param='s2', file='s2.txt', dir=dir, force=True) 975 self.interpreter.value.write(param='s2f', file='s2f.txt', dir=dir, force=True) 976 self.interpreter.value.write(param='s2s', file='s2s.txt', dir=dir, force=True) 977 self.interpreter.value.write(param='te', file='te.txt', dir=dir, force=True) 978 self.interpreter.value.write(param='tf', file='tf.txt', dir=dir, force=True) 979 self.interpreter.value.write(param='ts', file='ts.txt', dir=dir, force=True) 980 self.interpreter.value.write(param='rex', file='rex.txt', dir=dir, force=True) 981 self.interpreter.value.write(param='local_tm', file='local_tm.txt', dir=dir, force=True) 982 frqs = get_frequencies() 983 for i in range(len(frqs)): 984 comment = "This is the Rex value with units rad.s^-1 scaled to a magnetic field strength of %s MHz." % (frqs[i]/1e6) 985 self.interpreter.value.write(param='rex', file='rex_%s.txt'%int(frqs[i]/1e6), dir=dir, scaling=(2.0*pi*frqs[i])**2, comment=comment, force=True) 986 987 # Create the PyMOL macros. 988 dir = self.write_results_dir + 'final' + sep + 'pymol' 989 self.interpreter.pymol.macro_write(data_type='s2', dir=dir, force=True) 990 self.interpreter.pymol.macro_write(data_type='s2f', dir=dir, force=True) 991 self.interpreter.pymol.macro_write(data_type='s2s', dir=dir, force=True) 992 self.interpreter.pymol.macro_write(data_type='amp_fast', dir=dir, force=True) 993 self.interpreter.pymol.macro_write(data_type='amp_slow', dir=dir, force=True) 994 self.interpreter.pymol.macro_write(data_type='te', dir=dir, force=True) 995 self.interpreter.pymol.macro_write(data_type='tf', dir=dir, force=True) 996 self.interpreter.pymol.macro_write(data_type='ts', dir=dir, force=True) 997 self.interpreter.pymol.macro_write(data_type='time_fast', dir=dir, force=True) 998 self.interpreter.pymol.macro_write(data_type='time_slow', dir=dir, force=True) 999 self.interpreter.pymol.macro_write(data_type='rex', dir=dir, force=True) 1000 1001 # Create the Molmol macros. 1002 dir = self.write_results_dir + 'final' + sep + 'molmol' 1003 self.interpreter.molmol.macro_write(data_type='s2', dir=dir, force=True) 1004 self.interpreter.molmol.macro_write(data_type='s2f', dir=dir, force=True) 1005 self.interpreter.molmol.macro_write(data_type='s2s', dir=dir, force=True) 1006 self.interpreter.molmol.macro_write(data_type='amp_fast', dir=dir, force=True) 1007 self.interpreter.molmol.macro_write(data_type='amp_slow', dir=dir, force=True) 1008 self.interpreter.molmol.macro_write(data_type='te', dir=dir, force=True) 1009 self.interpreter.molmol.macro_write(data_type='tf', dir=dir, force=True) 1010 self.interpreter.molmol.macro_write(data_type='ts', dir=dir, force=True) 1011 self.interpreter.molmol.macro_write(data_type='time_fast', dir=dir, force=True) 1012 self.interpreter.molmol.macro_write(data_type='time_slow', dir=dir, force=True) 1013 self.interpreter.molmol.macro_write(data_type='rex', dir=dir, force=True) 1014 1015 # Create a diffusion tensor representation of the tensor, if a PDB file is present and the local tm global model has not been selected. 1016 if hasattr(cdp, 'structure') and hasattr(cdp, 'diff_tensor'): 1017 dir = self.write_results_dir + 'final' 1018 self.interpreter.structure.create_diff_tensor_pdb(file="tensor.pdb", dir=dir, force=True)

Source Code for Module auto_analyses.dauvergne_protocol