auto_analyses.dauvergne

116 """The model-free auto-analysis.""" 117 118 # Some class variables. 119 opt_func_tol = 1e-25 120 opt_max_iterations = int(1e7) 121

122 - def __init__(self, pipe_name=None, pipe_bundle=None, 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):

123 """Perform the full model-free analysis protocol of d'Auvergne and Gooley, 2008b. 124 125 @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. 126 @type pipe_name: str 127 @keyword pipe_bundle: The data pipe bundle to associate all spawned data pipes with. 128 @type pipe_bundle: str 129 @keyword results_dir: The directory, where files are saved in. 130 @type results_dir: str 131 @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. 132 @type diff_model: str or list of str 133 @keyword mf_models: The model-free models. 134 @type mf_models: list of str 135 @keyword local_tm_models: The model-free models. 136 @type local_tm_models: list of str 137 @keyword grid_inc: The grid search size (the number of increments per dimension). 138 @type grid_inc: int 139 @keyword diff_tensor_grid_inc: A list of grid search sizes for the optimisation of the sphere, prolate spheroid, oblate spheroid, and ellipsoid. 140 @type diff_tensor_grid_inc: list of int 141 @keyword min_algor: The minimisation algorithm (in most cases this should not be changed). 142 @type min_algor: str 143 @keyword mc_sim_num: The number of Monte Carlo simulations to be used for error analysis at the end of the analysis. 144 @type mc_sim_num: int 145 @keyword max_iter: The maximum number of iterations for the global iteration. Set to None, then the algorithm iterates until convergence. 146 @type max_iter: int or None. 147 @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. 148 @type user_fns: dict 149 @keyword conv_loop: Automatic looping over all rounds until convergence. 150 @type conv_loop: bool 151 """ 152 153 # Execution lock. 154 status.exec_lock.acquire(pipe_bundle, mode='auto-analysis') 155 156 # Store the args. 157 self.pipe_name = pipe_name 158 self.pipe_bundle = pipe_bundle 159 self.mf_models = mf_models 160 self.local_tm_models = local_tm_models 161 self.grid_inc = grid_inc 162 self.diff_tensor_grid_inc = diff_tensor_grid_inc 163 self.min_algor = min_algor 164 self.mc_sim_num = mc_sim_num 165 self.max_iter = max_iter 166 self.conv_loop = conv_loop 167 168 # The model-free data pipe names. 169 self.mf_model_pipes = [] 170 for i in range(len(self.mf_models)): 171 self.mf_model_pipes.append(self.name_pipe(self.mf_models[i])) 172 self.local_tm_model_pipes = [] 173 for i in range(len(self.local_tm_models)): 174 self.local_tm_model_pipes.append(self.name_pipe(self.local_tm_models[i])) 175 176 # The diffusion models. 177 if isinstance(diff_model, list): 178 self.diff_model_list = diff_model 179 else: 180 self.diff_model_list = [diff_model] 181 182 # Project directory (i.e. directory containing the model-free model results and the newly generated files) 183 if results_dir: 184 self.results_dir = results_dir + sep 185 else: 186 self.results_dir = getcwd() + sep 187 188 # Data checks. 189 self.check_vars() 190 191 # Set the data pipe to the current data pipe. 192 if self.pipe_name != cdp_name(): 193 switch(self.pipe_name) 194 195 # Some info for the status. 196 self.status_setup() 197 198 # Load the interpreter. 199 self.interpreter = Interpreter(show_script=False, quit=False, raise_relax_error=True) 200 self.interpreter.populate_self() 201 self.interpreter.on(verbose=False) 202 203 # Replacement user functions. 204 if user_fns: 205 for name in user_fns: 206 setattr(self.interpreter, name, user_fns[name]) 207 208 # Execute the protocol. 209 try: 210 # Loop over the models. 211 for self.diff_model in self.diff_model_list: 212 # Wait a little while between diffusion models. 213 sleep(1) 214 215 # Set the global model name. 216 status.auto_analysis[self.pipe_bundle].diff_model = self.diff_model 217 218 # Initialise the convergence data structures. 219 self.conv_data = Container() 220 self.conv_data.chi2 = [] 221 self.conv_data.models = [] 222 self.conv_data.diff_vals = [] 223 if self.diff_model == 'sphere': 224 self.conv_data.diff_params = ['tm'] 225 elif self.diff_model == 'oblate' or self.diff_model == 'prolate': 226 self.conv_data.diff_params = ['tm', 'Da', 'theta', 'phi'] 227 elif self.diff_model == 'ellipsoid': 228 self.conv_data.diff_params = ['tm', 'Da', 'Dr', 'alpha', 'beta', 'gamma'] 229 self.conv_data.spin_ids = [] 230 self.conv_data.mf_params = [] 231 self.conv_data.mf_vals = [] 232 233 # Execute the analysis for each diffusion model. 234 self.execute() 235 236 # Clean up. 237 finally: 238 # Finish and unlock execution. 239 status.auto_analysis[self.pipe_bundle].fin = True 240 status.current_analysis = None 241 status.exec_lock.release()

242 243

244 - def check_vars(self):

245 """Check that the user has set the variables correctly.""" 246 247 # The pipe bundle. 248 if not isinstance(self.pipe_bundle, str): 249 raise RelaxError("The pipe bundle name '%s' is invalid." % self.pipe_bundle) 250 251 # The diff model. 252 valid_models = ['local_tm', 'sphere', 'oblate', 'prolate', 'ellipsoid', 'final'] 253 for i in range(len(self.diff_model_list)): 254 if self.diff_model_list[i] not in valid_models: 255 raise RelaxError("The diff_model value '%s' is incorrectly set. It must be one of %s." % (self.diff_model_list[i], valid_models)) 256 257 # Model-free models. 258 mf_models = ['m0', 'm1', 'm2', 'm3', 'm4', 'm5', 'm6', 'm7', 'm8', 'm9'] 259 local_tm_models = ['tm0', 'tm1', 'tm2', 'tm3', 'tm4', 'tm5', 'tm6', 'tm7', 'tm8', 'tm9'] 260 if not isinstance(self.mf_models, list): 261 raise RelaxError("The self.mf_models user variable must be a list.") 262 if not isinstance(self.local_tm_models, list): 263 raise RelaxError("The self.local_tm_models user variable must be a list.") 264 for i in range(len(self.mf_models)): 265 if self.mf_models[i] not in mf_models: 266 raise RelaxError("The self.mf_models user variable '%s' is incorrectly set. It must be one of %s." % (self.mf_models, mf_models)) 267 for i in range(len(self.local_tm_models)): 268 if self.local_tm_models[i] not in local_tm_models: 269 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)) 270 271 # Sequence data. 272 if not exists_mol_res_spin_data(): 273 raise RelaxNoSequenceError(self.pipe_name) 274 275 # Relaxation data. 276 if not hasattr(cdp, 'ri_ids') or len(cdp.ri_ids) == 0: 277 raise RelaxNoRiError(ri_id) 278 279 # Insufficient data. 280 if len(cdp.ri_ids) <= 3: 281 raise RelaxError("Insufficient relaxation data, 4 or more data sets are essential for the execution of this script.") 282 283 # Spin vars. 284 for spin, spin_id in spin_loop(return_id=True): 285 # Skip deselected spins. 286 if not spin.select: 287 continue 288 289 # Test if the isotope type has been set. 290 if not hasattr(spin, 'isotope') or spin.isotope == None: 291 raise RelaxNoValueError("nuclear isotope type", spin_id=spin_id) 292 293 # Skip spins with no relaxation data. 294 if not hasattr(spin, 'ri_data') or spin.ri_data == None: 295 continue 296 297 # Test if the CSA value has been set. 298 if not hasattr(spin, 'csa') or spin.csa == None: 299 raise RelaxNoValueError("CSA", spin_id=spin_id) 300 301 # Interatomic vars. 302 for interatom in interatomic_loop(): 303 # Get the corresponding spins. 304 spin1 = return_spin(interatom.spin_id1) 305 spin2 = return_spin(interatom.spin_id2) 306 307 # Skip deselected spins. 308 if not spin1.select or not spin2.select: 309 continue 310 311 # Test if the interatomic distance has been set. 312 if not hasattr(interatom, 'r') or interatom.r == None: 313 raise RelaxNoValueError("interatomic distance", spin_id=interatom.spin_id1, spin_id2=interatom.spin_id2) 314 315 # Min vars. 316 if not isinstance(self.grid_inc, int): 317 raise RelaxError("The grid_inc user variable '%s' is incorrectly set. It should be an integer." % self.grid_inc) 318 if not isinstance(self.diff_tensor_grid_inc, dict): 319 raise RelaxError("The diff_tensor_grid_inc user variable '%s' is incorrectly set. It should be a dictionary." % self.diff_tensor_grid_inc) 320 for tensor in ['sphere', 'prolate', 'oblate', 'ellipsoid']: 321 if not tensor in self.diff_tensor_grid_inc: 322 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)) 323 if not isinstance(self.diff_tensor_grid_inc[tensor], int): 324 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)) 325 if not isinstance(self.min_algor, str): 326 raise RelaxError("The min_algor user variable '%s' is incorrectly set. It should be a string." % self.min_algor) 327 if not isinstance(self.mc_sim_num, int): 328 raise RelaxError("The mc_sim_num user variable '%s' is incorrectly set. It should be an integer." % self.mc_sim_num) 329 330 # Looping. 331 if not isinstance(self.conv_loop, bool): 332 raise RelaxError("The conv_loop user variable '%s' is incorrectly set. It should be one of the booleans True or False." % self.conv_loop)

333 334

335 - def convergence(self):

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

486 487

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

489 """Function for returning the name of next round of optimisation.""" 490 491 # Get a list of all files in the directory model. If no directory exists, set the round to 'init' or 0. 492 try: 493 dir_list = listdir(self.results_dir+sep+model) 494 except: 495 return 0 496 497 # Set the round to 'init' or 0 if there is no directory called 'init'. 498 if 'init' not in dir_list: 499 return 0 500 501 # Create a list of all files which begin with 'round_'. 502 rnd_dirs = [] 503 for file in dir_list: 504 if search('^round_', file): 505 rnd_dirs.append(file) 506 507 # Create a sorted list of integer round numbers. 508 numbers = [] 509 for dir in rnd_dirs: 510 try: 511 numbers.append(int(dir[6:])) 512 except: 513 pass 514 numbers.sort() 515 516 # No directories beginning with 'round_' exist, set the round to 1. 517 if not len(numbers): 518 return 1 519 520 # The highest number. 521 max_round = numbers[-1] 522 523 # Check that the opt/results file exists for the round (working backwards). 524 for i in range(max_round, -1, -1): 525 # Assume the round is complete. 526 complete_round = i 527 528 # The file root. 529 file_root = self.results_dir + sep + model + sep + "round_%i" % i + sep + 'opt' + sep + 'results' 530 531 # Stop looping when the opt/results file is found. 532 if access(file_root + '.bz2', F_OK): 533 break 534 if access(file_root + '.gz', F_OK): 535 break 536 if access(file_root, F_OK): 537 break 538 539 # No round, so assume the initial state. 540 if complete_round == 0: 541 return 0 542 543 # Determine the number for the next round (add 1 to the highest completed round). 544 return complete_round + 1

545 546

547 - def execute(self):

548 """Execute the protocol.""" 549 550 # MI - Local tm. 551 ################ 552 553 if self.diff_model == 'local_tm': 554 # Base directory to place files into. 555 self.base_dir = self.results_dir+'local_tm'+sep 556 557 # Sequential optimisation of all model-free models (function must be modified to suit). 558 self.multi_model(local_tm=True) 559 560 # Model selection. 561 self.model_selection(modsel_pipe=self.name_pipe('aic'), dir=self.base_dir + 'aic') 562 563 564 # Diffusion models MII to MV. 565 ############################# 566 567 elif self.diff_model == 'sphere' or self.diff_model == 'prolate' or self.diff_model == 'oblate' or self.diff_model == 'ellipsoid': 568 # No local_tm directory! 569 dir_list = listdir(self.results_dir) 570 if 'local_tm' not in dir_list: 571 raise RelaxError("The local_tm model must be optimised first.") 572 573 # The initial round of optimisation - not zero if calculations were interrupted. 574 self.start_round = self.determine_rnd(model=self.diff_model) 575 576 # Loop until convergence if conv_loop is set, otherwise just loop once. 577 # This looping could be made much cleaner by removing the dependence on the determine_rnd() function. 578 while True: 579 # Determine which round of optimisation to do (init, round_1, round_2, etc). 580 self.round = self.determine_rnd(model=self.diff_model) 581 status.auto_analysis[self.pipe_bundle].round = self.round 582 583 # Inital round of optimisation for diffusion models MII to MV. 584 if self.round == 0: 585 # Base directory to place files into. 586 self.base_dir = self.results_dir+self.diff_model+sep+'init'+sep 587 588 # Run name. 589 name = self.name_pipe(self.diff_model) 590 591 # Create the data pipe (deleting the old one if it exists). 592 if has_pipe(name): 593 self.interpreter.pipe.delete(name) 594 self.interpreter.pipe.create(name, 'mf', bundle=self.pipe_bundle) 595 596 # Load the local tm diffusion model MI results. 597 self.interpreter.results.read(file='results', dir=self.results_dir+'local_tm'+sep+'aic') 598 599 # Remove the tm parameter. 600 self.interpreter.model_free.remove_tm() 601 602 # Add an arbitrary diffusion tensor which will be optimised. 603 if self.diff_model == 'sphere': 604 self.interpreter.diffusion_tensor.init(10e-9, fixed=False) 605 inc = self.diff_tensor_grid_inc['sphere'] 606 elif self.diff_model == 'prolate': 607 self.interpreter.diffusion_tensor.init((10e-9, 0, 0, 0), spheroid_type='prolate', fixed=False) 608 inc = self.diff_tensor_grid_inc['prolate'] 609 elif self.diff_model == 'oblate': 610 self.interpreter.diffusion_tensor.init((10e-9, 0, 0, 0), spheroid_type='oblate', fixed=False) 611 inc = self.diff_tensor_grid_inc['oblate'] 612 elif self.diff_model == 'ellipsoid': 613 self.interpreter.diffusion_tensor.init((10e-09, 0, 0, 0, 0, 0), fixed=False) 614 inc = self.diff_tensor_grid_inc['ellipsoid'] 615 616 # Minimise just the diffusion tensor. 617 self.interpreter.fix('all_spins') 618 self.interpreter.grid_search(inc=inc) 619 self.interpreter.minimise(self.min_algor, func_tol=self.opt_func_tol, max_iter=self.opt_max_iterations) 620 621 # Write the results. 622 self.interpreter.results.write(file='results', dir=self.base_dir, force=True) 623 624 625 # Normal round of optimisation for diffusion models MII to MV. 626 else: 627 # Base directory to place files into. 628 self.base_dir = self.results_dir+self.diff_model + sep+'round_'+repr(self.round)+sep 629 630 # Load the optimised diffusion tensor from either the previous round. 631 self.load_tensor() 632 633 # Sequential optimisation of all model-free models (function must be modified to suit). 634 self.multi_model() 635 636 # Model selection. 637 self.model_selection(modsel_pipe=self.name_pipe('aic'), dir=self.base_dir + 'aic') 638 639 # Final optimisation of all diffusion and model-free parameters. 640 self.interpreter.fix('all', fixed=False) 641 642 # Minimise all parameters. 643 self.interpreter.minimise(self.min_algor, func_tol=self.opt_func_tol, max_iter=self.opt_max_iterations) 644 645 # Write the results. 646 dir = self.base_dir + 'opt' 647 self.interpreter.results.write(file='results', dir=dir, force=True) 648 649 # Test for convergence. 650 converged = self.convergence() 651 652 # Break out of the infinite while loop if automatic looping is not activated or if convergence has occurred. 653 if converged or not self.conv_loop: 654 break 655 656 # Unset the status. 657 status.auto_analysis[self.pipe_bundle].round = None 658 659 660 # Final run. 661 ############ 662 663 elif self.diff_model == 'final': 664 # Diffusion model selection. 665 ############################ 666 667 # All the global diffusion models to be used in the model selection. 668 models = ['local_tm', 'sphere', 'prolate', 'oblate', 'ellipsoid'] 669 self.pipes = [] 670 for name in models: 671 self.pipes.append(self.name_pipe(name)) 672 673 # Remove all temporary pipes used in this auto-analysis. 674 for name in pipe_names(bundle=self.pipe_bundle): 675 if name in self.pipes + self.mf_model_pipes + self.local_tm_model_pipes + [self.name_pipe('aic'), self.name_pipe('previous')]: 676 self.interpreter.pipe.delete(name) 677 678 # Missing optimised model. 679 dir_list = listdir(self.results_dir) 680 for name in models: 681 if name not in dir_list: 682 raise RelaxError("The %s model must be optimised first." % name) 683 684 # Create the local_tm data pipe. 685 self.interpreter.pipe.create(self.name_pipe('local_tm'), 'mf', bundle=self.pipe_bundle) 686 687 # Load the local tm diffusion model MI results. 688 self.interpreter.results.read(file='results', dir=self.results_dir+'local_tm'+sep+'aic') 689 690 # Loop over models MII to MV. 691 for model in ['sphere', 'prolate', 'oblate', 'ellipsoid']: 692 # Determine which was the last round of optimisation for each of the models. 693 self.round = self.determine_rnd(model=model) - 1 694 695 # If no directories begining with 'round_' exist, the script has not been properly utilised! 696 if self.round < 1: 697 # Construct the name of the diffusion tensor. 698 name = model 699 if model == 'prolate' or model == 'oblate': 700 name = name + ' spheroid' 701 702 # Throw an error to prevent misuse of the script. 703 raise RelaxError("Multiple rounds of optimisation of the " + name + " (between 8 to 15) are required for the proper execution of this script.") 704 705 # Create the data pipe. 706 self.interpreter.pipe.create(self.name_pipe(model), 'mf', bundle=self.pipe_bundle) 707 708 # Load the diffusion model results. 709 self.interpreter.results.read(file='results', dir=self.results_dir+model + sep+'round_'+repr(self.round)+sep+'opt') 710 711 # Model selection between MI to MV. 712 self.model_selection(modsel_pipe=self.name_pipe('final'), write_flag=False) 713 714 715 # Monte Carlo simulations. 716 ########################## 717 718 # Fix the diffusion tensor, if it exists. 719 if hasattr(get_pipe(self.name_pipe('final')), 'diff_tensor'): 720 self.interpreter.fix('diff') 721 722 # Simulations. 723 self.interpreter.monte_carlo.setup(number=self.mc_sim_num) 724 self.interpreter.monte_carlo.create_data() 725 self.interpreter.monte_carlo.initial_values() 726 self.interpreter.minimise(self.min_algor, func_tol=self.opt_func_tol, max_iter=self.opt_max_iterations) 727 self.interpreter.eliminate() 728 self.interpreter.monte_carlo.error_analysis() 729 730 731 # Write the final results. 732 ########################## 733 734 # Create results files and plots of the data. 735 self.write_results() 736 737 738 # Unknown script behaviour. 739 ########################### 740 741 else: 742 raise RelaxError("Unknown diffusion model, change the value of 'self.diff_model'")

743 744

745 - def load_tensor(self):

746 """Function for loading the optimised diffusion tensor.""" 747 748 # Create the data pipe for the previous data (deleting the old data pipe first if necessary). 749 if has_pipe(self.name_pipe('previous')): 750 self.interpreter.pipe.delete(self.name_pipe('previous')) 751 self.interpreter.pipe.create(self.name_pipe('previous'), 'mf', bundle=self.pipe_bundle) 752 753 # Load the optimised diffusion tensor from the initial round. 754 if self.round == 1: 755 self.interpreter.results.read('results', self.results_dir+self.diff_model + sep+'init') 756 757 # Load the optimised diffusion tensor from the previous round. 758 else: 759 self.interpreter.results.read('results', self.results_dir+self.diff_model + sep+'round_'+repr(self.round-1)+sep+'opt')

760 761

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

763 """Model selection function.""" 764 765 # Model selection (delete the model selection pipe if it already exists). 766 if has_pipe(modsel_pipe): 767 self.interpreter.pipe.delete(modsel_pipe) 768 self.interpreter.model_selection(method='AIC', modsel_pipe=modsel_pipe, bundle=self.pipe_bundle, pipes=self.pipes) 769 770 # Write the results. 771 if write_flag: 772 self.interpreter.results.write(file='results', dir=dir, force=True)

773 774

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

776 """Function for optimisation of all model-free models.""" 777 778 # Set the data pipe names (also the names of preset model-free models). 779 if local_tm: 780 models = self.local_tm_models 781 self.pipes = self.local_tm_models 782 else: 783 models = self.mf_models 784 self.pipes = [] 785 for i in range(len(models)): 786 self.pipes.append(self.name_pipe(models[i])) 787 788 # Loop over the data pipes. 789 for i in range(len(models)): 790 # Place the model name into the status container. 791 status.auto_analysis[self.pipe_bundle].current_model = models[i] 792 793 # Create the data pipe (by copying). 794 if has_pipe(self.pipes[i]): 795 self.interpreter.pipe.delete(self.pipes[i]) 796 self.interpreter.pipe.copy(self.pipe_name, self.pipes[i], bundle_to=self.pipe_bundle) 797 self.interpreter.pipe.switch(self.pipes[i]) 798 799 # Copy the diffusion tensor from the 'opt' data pipe and prevent it from being minimised. 800 if not local_tm: 801 self.interpreter.diffusion_tensor.copy(self.name_pipe('previous')) 802 self.interpreter.fix('diff') 803 804 # Select the model-free model. 805 self.interpreter.model_free.select_model(model=models[i]) 806 807 # Minimise. 808 self.interpreter.grid_search(inc=self.grid_inc) 809 self.interpreter.minimise(self.min_algor, func_tol=self.opt_func_tol, max_iter=self.opt_max_iterations) 810 811 # Model elimination. 812 self.interpreter.eliminate() 813 814 # Write the results. 815 dir = self.base_dir + models[i] 816 self.interpreter.results.write(file='results', dir=dir, force=True) 817 818 # Unset the status. 819 status.auto_analysis[self.pipe_bundle].current_model = None

820 821

822 - def name_pipe(self, prefix):

823 """Generate a unique name for the data pipe. 824 825 @param prefix: The prefix of the data pipe name. 826 @type prefix: str 827 """ 828 829 # The unique pipe name. 830 name = "%s - %s" % (prefix, self.pipe_bundle) 831 832 # Return the name. 833 return name

834 835

836 - def status_setup(self):

837 """Initialise the status object.""" 838 839 # Initialise the status object for this auto-analysis. 840 status.init_auto_analysis(self.pipe_bundle, type='dauvergne_protocol') 841 status.current_analysis = self.pipe_bundle 842 843 # The global diffusion model. 844 status.auto_analysis[self.pipe_bundle].diff_model = None 845 846 # The round of optimisation, i.e. the global iteration. 847 status.auto_analysis[self.pipe_bundle].round = None 848 849 # The list of model-free local tm models for optimisation, i.e. the global iteration. 850 status.auto_analysis[self.pipe_bundle].local_tm_models = self.local_tm_models 851 852 # The list of model-free models for optimisation, i.e. the global iteration. 853 status.auto_analysis[self.pipe_bundle].mf_models = self.mf_models 854 855 # The current model-free model. 856 status.auto_analysis[self.pipe_bundle].current_model = None 857 858 # The maximum number of iterations of the global model. 859 status.auto_analysis[self.pipe_bundle].max_iter = self.max_iter 860 861 # The convergence of the global model. 862 status.auto_analysis[self.pipe_bundle].convergence = False

863 864

865 - def write_results(self):

866 """Create Grace plots of the final model-free results.""" 867 868 # Save the results file. 869 dir = self.results_dir + 'final' 870 self.interpreter.results.write(file='results', dir=dir, force=True) 871 872 # The Grace plots. 873 dir = self.results_dir + 'final' + sep + 'grace' 874 self.interpreter.grace.write(x_data_type='spin', y_data_type='s2', file='s2.agr', dir=dir, force=True) 875 self.interpreter.grace.write(x_data_type='spin', y_data_type='s2f', file='s2f.agr', dir=dir, force=True) 876 self.interpreter.grace.write(x_data_type='spin', y_data_type='s2s', file='s2s.agr', dir=dir, force=True) 877 self.interpreter.grace.write(x_data_type='spin', y_data_type='te', file='te.agr', dir=dir, force=True) 878 self.interpreter.grace.write(x_data_type='spin', y_data_type='tf', file='tf.agr', dir=dir, force=True) 879 self.interpreter.grace.write(x_data_type='spin', y_data_type='ts', file='ts.agr', dir=dir, force=True) 880 self.interpreter.grace.write(x_data_type='spin', y_data_type='rex', file='rex.agr', dir=dir, force=True) 881 self.interpreter.grace.write(x_data_type='s2', y_data_type='te', file='s2_vs_te.agr', dir=dir, force=True) 882 self.interpreter.grace.write(x_data_type='s2', y_data_type='rex', file='s2_vs_rex.agr', dir=dir, force=True) 883 self.interpreter.grace.write(x_data_type='te', y_data_type='rex', file='te_vs_rex.agr', dir=dir, force=True) 884 885 # Write the values to text files. 886 dir = self.results_dir + 'final' 887 self.interpreter.value.write(param='s2', file='s2.txt', dir=dir, force=True) 888 self.interpreter.value.write(param='s2f', file='s2f.txt', dir=dir, force=True) 889 self.interpreter.value.write(param='s2s', file='s2s.txt', dir=dir, force=True) 890 self.interpreter.value.write(param='te', file='te.txt', dir=dir, force=True) 891 self.interpreter.value.write(param='tf', file='tf.txt', dir=dir, force=True) 892 self.interpreter.value.write(param='ts', file='ts.txt', dir=dir, force=True) 893 self.interpreter.value.write(param='rex', file='rex.txt', dir=dir, force=True) 894 self.interpreter.value.write(param='local_tm', file='local_tm.txt', dir=dir, force=True) 895 896 # Create the PyMOL macros. 897 dir = self.results_dir + 'final' + sep + 'pymol' 898 self.interpreter.pymol.macro_write(data_type='s2', dir=dir, force=True) 899 self.interpreter.pymol.macro_write(data_type='s2f', dir=dir, force=True) 900 self.interpreter.pymol.macro_write(data_type='s2s', dir=dir, force=True) 901 self.interpreter.pymol.macro_write(data_type='amp_fast', dir=dir, force=True) 902 self.interpreter.pymol.macro_write(data_type='amp_slow', dir=dir, force=True) 903 self.interpreter.pymol.macro_write(data_type='te', dir=dir, force=True) 904 self.interpreter.pymol.macro_write(data_type='tf', dir=dir, force=True) 905 self.interpreter.pymol.macro_write(data_type='ts', dir=dir, force=True) 906 self.interpreter.pymol.macro_write(data_type='time_fast', dir=dir, force=True) 907 self.interpreter.pymol.macro_write(data_type='time_slow', dir=dir, force=True) 908 self.interpreter.pymol.macro_write(data_type='rex', dir=dir, force=True) 909 910 # Create the Molmol macros. 911 dir = self.results_dir + 'final' + sep + 'molmol' 912 self.interpreter.molmol.macro_write(data_type='s2', dir=dir, force=True) 913 self.interpreter.molmol.macro_write(data_type='s2f', dir=dir, force=True) 914 self.interpreter.molmol.macro_write(data_type='s2s', dir=dir, force=True) 915 self.interpreter.molmol.macro_write(data_type='amp_fast', dir=dir, force=True) 916 self.interpreter.molmol.macro_write(data_type='amp_slow', dir=dir, force=True) 917 self.interpreter.molmol.macro_write(data_type='te', dir=dir, force=True) 918 self.interpreter.molmol.macro_write(data_type='tf', dir=dir, force=True) 919 self.interpreter.molmol.macro_write(data_type='ts', dir=dir, force=True) 920 self.interpreter.molmol.macro_write(data_type='time_fast', dir=dir, force=True) 921 self.interpreter.molmol.macro_write(data_type='time_slow', dir=dir, force=True) 922 self.interpreter.molmol.macro_write(data_type='rex', dir=dir, force=True) 923 924 # Create a diffusion tensor representation of the tensor. 925 dir = self.results_dir + 'final' 926 self.interpreter.structure.create_diff_tensor_pdb(file="tensor.pdb", dir=dir, force=True)

Source Code for Module auto_analyses.dauvergne_protocol