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

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

248 249

250 - def check_vars(self):

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

339 340

341 - def convergence(self):

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

492 493

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

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

551 552

553 - def execute(self):

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

759 760

761 - def load_tensor(self):

762 """Function for loading the optimised diffusion tensor.""" 763 764 # Create the data pipe for the previous data (deleting the old data pipe first if necessary). 765 if has_pipe(self.name_pipe('previous')): 766 self.interpreter.pipe.delete(self.name_pipe('previous')) 767 self.interpreter.pipe.create(self.name_pipe('previous'), 'mf', bundle=self.pipe_bundle) 768 769 # Load the optimised diffusion tensor from the initial round. 770 if self.round == 1: 771 self.interpreter.results.read('results', self.results_dir+self.diff_model + sep+'init') 772 773 # Load the optimised diffusion tensor from the previous round. 774 else: 775 self.interpreter.results.read('results', self.results_dir+self.diff_model + sep+'round_'+repr(self.round-1)+sep+'opt')

776 777

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

779 """Model selection function.""" 780 781 # Model selection (delete the model selection pipe if it already exists). 782 if has_pipe(modsel_pipe): 783 self.interpreter.pipe.delete(modsel_pipe) 784 self.interpreter.model_selection(method='AIC', modsel_pipe=modsel_pipe, bundle=self.pipe_bundle, pipes=self.pipes) 785 786 # Write the results. 787 if write_flag: 788 self.interpreter.results.write(file='results', dir=dir, force=True)

789 790

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

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

836 837

838 - def name_pipe(self, prefix):

839 """Generate a unique name for the data pipe. 840 841 @param prefix: The prefix of the data pipe name. 842 @type prefix: str 843 """ 844 845 # The unique pipe name. 846 name = "%s - %s" % (prefix, self.pipe_bundle) 847 848 # Return the name. 849 return name

850 851

852 - def status_setup(self):

853 """Initialise the status object.""" 854 855 # Initialise the status object for this auto-analysis. 856 status.init_auto_analysis(self.pipe_bundle, type='dauvergne_protocol') 857 status.current_analysis = self.pipe_bundle 858 859 # The global diffusion model. 860 status.auto_analysis[self.pipe_bundle].diff_model = None 861 862 # The round of optimisation, i.e. the global iteration. 863 status.auto_analysis[self.pipe_bundle].round = None 864 865 # The list of model-free local tm models for optimisation, i.e. the global iteration. 866 status.auto_analysis[self.pipe_bundle].local_tm_models = self.local_tm_models 867 868 # The list of model-free models for optimisation, i.e. the global iteration. 869 status.auto_analysis[self.pipe_bundle].mf_models = self.mf_models 870 871 # The current model-free model. 872 status.auto_analysis[self.pipe_bundle].current_model = None 873 874 # The maximum number of iterations of the global model. 875 status.auto_analysis[self.pipe_bundle].max_iter = self.max_iter 876 877 # The convergence of the global model. 878 status.auto_analysis[self.pipe_bundle].convergence = False

879 880

881 - def write_results(self):

882 """Create Grace plots of the final model-free results.""" 883 884 # Save the results file. 885 dir = self.write_results_dir + 'final' 886 self.interpreter.results.write(file='results', dir=dir, force=True) 887 888 # The Grace plots. 889 dir = self.write_results_dir + 'final' + sep + 'grace' 890 self.interpreter.grace.write(x_data_type='res_num', y_data_type='s2', file='s2.agr', dir=dir, force=True) 891 self.interpreter.grace.write(x_data_type='res_num', y_data_type='s2f', file='s2f.agr', dir=dir, force=True) 892 self.interpreter.grace.write(x_data_type='res_num', y_data_type='s2s', file='s2s.agr', dir=dir, force=True) 893 self.interpreter.grace.write(x_data_type='res_num', y_data_type='te', file='te.agr', dir=dir, force=True) 894 self.interpreter.grace.write(x_data_type='res_num', y_data_type='tf', file='tf.agr', dir=dir, force=True) 895 self.interpreter.grace.write(x_data_type='res_num', y_data_type='ts', file='ts.agr', dir=dir, force=True) 896 self.interpreter.grace.write(x_data_type='res_num', y_data_type='rex', file='rex.agr', dir=dir, force=True) 897 self.interpreter.grace.write(x_data_type='s2', y_data_type='te', file='s2_vs_te.agr', dir=dir, force=True) 898 self.interpreter.grace.write(x_data_type='s2', y_data_type='rex', file='s2_vs_rex.agr', dir=dir, force=True) 899 self.interpreter.grace.write(x_data_type='te', y_data_type='rex', file='te_vs_rex.agr', dir=dir, force=True) 900 901 # Write the values to text files. 902 dir = self.write_results_dir + 'final' 903 self.interpreter.value.write(param='s2', file='s2.txt', dir=dir, force=True) 904 self.interpreter.value.write(param='s2f', file='s2f.txt', dir=dir, force=True) 905 self.interpreter.value.write(param='s2s', file='s2s.txt', dir=dir, force=True) 906 self.interpreter.value.write(param='te', file='te.txt', dir=dir, force=True) 907 self.interpreter.value.write(param='tf', file='tf.txt', dir=dir, force=True) 908 self.interpreter.value.write(param='ts', file='ts.txt', dir=dir, force=True) 909 self.interpreter.value.write(param='rex', file='rex.txt', dir=dir, force=True) 910 self.interpreter.value.write(param='local_tm', file='local_tm.txt', dir=dir, force=True) 911 frqs = spectrometer.get_frequencies() 912 for i in range(len(frqs)): 913 comment = "This is the Rex value with units rad.s^-1 scaled to a magnetic field strength of %s MHz." % (frqs[i]/1e6) 914 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) 915 916 # Create the PyMOL macros. 917 dir = self.write_results_dir + 'final' + sep + 'pymol' 918 self.interpreter.pymol.macro_write(data_type='s2', dir=dir, force=True) 919 self.interpreter.pymol.macro_write(data_type='s2f', dir=dir, force=True) 920 self.interpreter.pymol.macro_write(data_type='s2s', dir=dir, force=True) 921 self.interpreter.pymol.macro_write(data_type='amp_fast', dir=dir, force=True) 922 self.interpreter.pymol.macro_write(data_type='amp_slow', dir=dir, force=True) 923 self.interpreter.pymol.macro_write(data_type='te', dir=dir, force=True) 924 self.interpreter.pymol.macro_write(data_type='tf', dir=dir, force=True) 925 self.interpreter.pymol.macro_write(data_type='ts', dir=dir, force=True) 926 self.interpreter.pymol.macro_write(data_type='time_fast', dir=dir, force=True) 927 self.interpreter.pymol.macro_write(data_type='time_slow', dir=dir, force=True) 928 self.interpreter.pymol.macro_write(data_type='rex', dir=dir, force=True) 929 930 # Create the Molmol macros. 931 dir = self.write_results_dir + 'final' + sep + 'molmol' 932 self.interpreter.molmol.macro_write(data_type='s2', dir=dir, force=True) 933 self.interpreter.molmol.macro_write(data_type='s2f', dir=dir, force=True) 934 self.interpreter.molmol.macro_write(data_type='s2s', dir=dir, force=True) 935 self.interpreter.molmol.macro_write(data_type='amp_fast', dir=dir, force=True) 936 self.interpreter.molmol.macro_write(data_type='amp_slow', dir=dir, force=True) 937 self.interpreter.molmol.macro_write(data_type='te', dir=dir, force=True) 938 self.interpreter.molmol.macro_write(data_type='tf', dir=dir, force=True) 939 self.interpreter.molmol.macro_write(data_type='ts', dir=dir, force=True) 940 self.interpreter.molmol.macro_write(data_type='time_fast', dir=dir, force=True) 941 self.interpreter.molmol.macro_write(data_type='time_slow', dir=dir, force=True) 942 self.interpreter.molmol.macro_write(data_type='rex', dir=dir, force=True) 943 944 # Create a diffusion tensor representation of the tensor, if a PDB file is present and the local tm global model has not been selected. 945 if hasattr(cdp, 'structure') and hasattr(cdp, 'diff_tensor'): 946 dir = self.write_results_dir + 'final' 947 self.interpreter.structure.create_diff_tensor_pdb(file="tensor.pdb", dir=dir, force=True)

Source Code for Module auto_analyses.dauvergne_protocol