auto_analyses.dauvergne

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

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

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

243 244

245 - def check_vars(self):

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

327 328

329 - def convergence(self):

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

480 481

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

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

539 540

541 - def execute(self):

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

737 738

739 - def load_tensor(self):

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

754 755

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

757 """Model selection function.""" 758 759 # Model selection (delete the model selection pipe if it already exists). 760 if has_pipe(modsel_pipe): 761 self.interpreter.pipe.delete(modsel_pipe) 762 self.interpreter.model_selection(method='AIC', modsel_pipe=modsel_pipe, bundle=self.pipe_bundle, pipes=self.pipes) 763 764 # Write the results. 765 if write_flag: 766 self.interpreter.results.write(file='results', dir=dir, force=True)

767 768

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

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

814 815

816 - def name_pipe(self, prefix):

817 """Generate a unique name for the data pipe. 818 819 @param prefix: The prefix of the data pipe name. 820 @type prefix: str 821 """ 822 823 # The unique pipe name. 824 name = "%s - %s" % (prefix, self.pipe_bundle) 825 826 # Return the name. 827 return name

828 829

830 - def status_setup(self):

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

857 858

859 - def write_results(self):

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

Source Code for Module auto_analyses.dauvergne_protocol