auto_analyses.relax_disp_repeat

106 107 """The relaxation dispersion analysis for repeated data.""" 108 109 # Some class variables. 110 opt_func_tol = 1e-25 111 opt_max_iterations = int(1e7) 112 113

114 - def __init__(self, settings):

115 """Perform a repeated dispersion analysis for settings given.""" 116 117 # Store settings. 118 self.settings = settings 119 120 # Unpack settings from dictionary to self. 121 for setting in self.settings: 122 setattr(self, setting, self.settings[setting]) 123 124 if 'pipe_type' not in self.settings: 125 self.set_self(key='pipe_type', value='relax_disp') 126 127 if 'time' not in self.settings: 128 self.set_self(key='time', value=datetime.now().strftime('%Y_%m_%d_%H_%M')) 129 130 # No results directory, so default to the current directory. 131 if 'results_dir' not in self.settings: 132 self.set_self(key='results_dir', value=getcwd() + sep + 'results' + sep + self.time ) 133 134 # Standard Monte-Carlo simulations. 135 if 'mc_sim_num' not in self.settings: 136 self.set_self(key='mc_sim_num', value=40) 137 138 # Standard Monte-Carlo simulations for exponential fit. '-1' is getting R2eff err from Co-variance. 139 if 'exp_mc_sim_num' not in self.settings: 140 self.set_self(key='exp_mc_sim_num', value=-1) 141 142 # Standard Monte-Carlo simulations for exponential fit. '-1' is getting R2eff err from Co-variance. 143 if 'modsel' not in self.settings: 144 self.set_self(key='modsel', value='AIC') 145 146 # The R2eff/R1rho value in rad/s by which to judge insignificance. If the maximum difference between two points on all dispersion curves for a spin is less than this value, that spin will be deselected. 147 if 'insignificance' not in self.settings: 148 self.set_self(key='insignificance', value=0.0) 149 150 # A flag which if True will activate R1 parameter fitting via relax_disp.r1_fit for the models that support it. 151 # If False, then the relax_disp.r1_fit user function will not be called. 152 if 'r1_fit' not in self.settings: 153 self.set_self(key='r1_fit', value=False) 154 155 # The minimisation algorithm. 156 if 'min_algor' not in self.settings: 157 self.set_self(key='min_algor', value='simplex') 158 159 # The constraints settings. 160 if 'constraints' not in self.settings: 161 self.set_self(key='constraints', value=True) 162 163 # The base setup. 164 if 'base_setup_pipe_name' not in self.settings: 165 base_setup_pipe_name = self.name_pipe(method='setup', model='setup', analysis='setup', glob_ini='setup') 166 self.set_self(key='base_setup_pipe_name', value=base_setup_pipe_name) 167 168 # Start interpreter. 169 self.interpreter_start()

170 171

172 - def set_base_cpmg(self, method=None, glob_ini=None, force=False):

173 """ Setup base information, but do not load intensity. """ 174 175 # Define model 176 model = 'setup' 177 analysis = 'setup' 178 179 # Check previous, and get the pipe name. 180 found, pipe_name, resfile, path = self.check_previous_result(method='setup', model=model, analysis=analysis, glob_ini='setup', bundle='setup') 181 182 # If found, then pass, else calculate it. 183 if found: 184 pass 185 else: 186 # Create the data pipe. 187 self.interpreter.pipe.create(pipe_name=pipe_name, pipe_type=self.pipe_type, bundle=None) 188 189 # Loop over frequency, store spectrum ids. 190 dic_spectrum_ids = {} 191 dic_spectrum_ids_replicates = {} 192 for i, sfrq in enumerate(self.sfrqs): 193 # Access the key in self. 194 key = DIC_KEY_FORMAT % (sfrq) 195 196 # Loop over cpmg_frqs. 197 cpmg_frqs = getattr(self, key)['cpmg_frqs'] 198 199 # Get the folder for peak files. 200 peaks_folder = getattr(self, key)['peaks_folder'] + sep + method 201 202 # Define glop pattern for peak files. 203 peaks_glob_pat = '%s_%s.ser' % (glob_ini, method) 204 205 # Get the file list. 206 peaks_file_list = glob(peaks_folder + sep + peaks_glob_pat) 207 208 # Sort the file list Alphanumeric. 209 peaks_file_list = sort_filenames(filenames=peaks_file_list) 210 211 # Create the spins. 212 for peaks_file in peaks_file_list: 213 self.interpreter.spectrum.read_spins(file=peaks_file, dir=None) 214 215 # Collect data keys. 216 dic_spectrum_ids[key] = [] 217 for j, cpmg_frq in enumerate(cpmg_frqs): 218 # Define the key. 219 data_key = return_param_key_from_data(exp_type=self.exp_type, frq=sfrq, point=cpmg_frq) 220 spectrum_id = data_key + '_%i'%j 221 222 # Store data key 223 dic_spectrum_ids[key].append(spectrum_id) 224 225 # Set the current experiment type. 226 self.interpreter.relax_disp.exp_type(spectrum_id=spectrum_id, exp_type=self.exp_type) 227 228 # Set the relaxation dispersion CPMG frequencies. 229 if cpmg_frq == 0.0: 230 cpmg_frq = None 231 self.interpreter.relax_disp.cpmg_setup(spectrum_id=spectrum_id, cpmg_frq=cpmg_frq) 232 233 # Relaxation dispersion CPMG constant time delay T (in s). 234 time_T2 = getattr(self, key)['time_T2'] 235 self.interpreter.relax_disp.relax_time(spectrum_id=spectrum_id, time=time_T2) 236 237 # Set the NMR field strength of the spectrum. 238 self.interpreter.spectrometer.frequency(id=spectrum_id, frq=sfrq, units=self.sfrq_unit) 239 240 # Get the list of duplications 241 list_dub = self.get_dublicates(dic_spectrum_ids[key], cpmg_frqs) 242 243 # Store to dic 244 dic_spectrum_ids_replicates[key] = list_dub 245 246 # Store to current data pipe. 247 cdp.dic_spectrum_ids = dic_spectrum_ids 248 cdp.dic_spectrum_ids_replicates = dic_spectrum_ids_replicates 249 250 # Name the isotope for field strength scaling. 251 self.interpreter.spin.isotope(isotope=self.isotope) 252 253 # Finally store the pipe name and save the setup pipe. 254 self.interpreter.results.write(file=resfile, dir=path, force=force)

255 256

257 - def set_intensity_and_error(self, pipe_name, glob_ini=None, set_rmsd=None):

258 # Read the intensity per spectrum id and set the RMSD error. 259 260 # Switch to the pipe. 261 if pipes.cdp_name() != pipe_name: 262 self.interpreter.pipe.switch(pipe_name) 263 264 # Loop over spectrometer frequencies. 265 finished = len(self.sfrqs) * [False] 266 for i, sfrq in enumerate(self.sfrqs): 267 # Access the key in self. 268 key = DIC_KEY_FORMAT % (sfrq) 269 270 # Get the spectrum ids. 271 spectrum_ids = cdp.dic_spectrum_ids[key] 272 273 # Get the folder for peak files. 274 peaks_folder = getattr(self, key)['peaks_folder'] + sep + self.method 275 276 # Define glop pattern for peak files. 277 peaks_glob_pat = '%s_%s.ser' % (glob_ini, self.method) 278 279 # Get the file list. 280 peaks_file_list = glob(peaks_folder + sep + peaks_glob_pat) 281 282 # Sort the file list Alphanumeric. 283 peaks_file_list = sort_filenames(filenames=peaks_file_list) 284 285 # If there is no peak list, then continue. 286 if len(peaks_file_list) == 0: 287 finished[i] = False 288 continue 289 290 # There should only be one peak file. 291 for peaks_file in peaks_file_list: 292 self.interpreter.spectrum.read_intensities(file=peaks_file, spectrum_id=spectrum_ids, int_method=self.int_method, int_col=list(range(len(spectrum_ids)))) 293 294 if set_rmsd: 295 # Get the folder for rmsd files. 296 rmsd_folder = getattr(self, key)['rmsd_folder'] 297 298 # Define glop pattern for rmsd files. 299 rmsd_glob_pat = '%s_*_%s.rmsd' % (glob_ini, self.method) 300 301 # Get the file list. 302 rmsd_file_list = glob(rmsd_folder + sep + rmsd_glob_pat) 303 304 # Sort the file list Alphanumeric. 305 rmsd_file_list = sort_filenames(filenames=rmsd_file_list) 306 307 # Loop over spectrum ids 308 for j, spectrum_id in enumerate(spectrum_ids): 309 # Set the peak intensity errors, as defined as the baseplane RMSD. 310 rmsd_file = rmsd_file_list[j] 311 312 # Extract rmsd from line 0, and column 0. 313 rmsd = float(extract_data(file=rmsd_file)[0][0]) 314 self.interpreter.spectrum.baseplane_rmsd(error=rmsd, spectrum_id=spectrum_id) 315 316 finished[i] = True 317 318 return all(finished)

319 320

321 - def do_spectrum_error_analysis(self, pipe_name, set_rep=None):

322 """Do spectrum error analysis, where both replicates per spectrometer frequency and subset is taken into consideration.""" 323 324 325 # Switch to the pipe. 326 if pipes.cdp_name() != pipe_name: 327 self.interpreter.pipe.switch(pipe_name) 328 329 # Loop over spectrometer frequencies. 330 for i, sfrq in enumerate(self.sfrqs): 331 # Access the key in self. 332 key = DIC_KEY_FORMAT % (sfrq) 333 334 # Printout. 335 section(file=sys.stdout, text="Error analysis for pipe='%s' and sfr:%3.2f"%(pipe_name, sfrq), prespace=2) 336 337 # Get the spectrum ids. 338 spectrum_ids = cdp.dic_spectrum_ids[key] 339 340 if set_rep: 341 # Get the spectrum ids replicates. 342 spectrum_ids_replicates = cdp.dic_spectrum_ids_replicates[key] 343 344 # Check if there are any replicates. 345 for replicate in spectrum_ids_replicates: 346 spectrum_id, rep_list = replicate 347 348 # If there is a replicated list, specify it. 349 if len(rep_list) > 0: 350 # Define the replicates. 351 self.interpreter.spectrum.replicated(spectrum_ids=rep_list) 352 353 # Run the error analysis on the subset. 354 self.interpreter.spectrum.error_analysis(subset=spectrum_ids)

355 356

357 - def set_int(self, methods=None, list_glob_ini=None, set_rmsd=True, set_rep=False, force=False):

358 """Call both the setup of data and the error analysis""" 359 360 # Define model 361 model = 'setup' 362 analysis = 'int' 363 364 # Loop over the methods. 365 finished = len(methods) * [False] 366 for i, method in enumerate(methods): 367 # Change the self key. 368 self.set_self(key='method', value=method) 369 370 # Loop over the glob ini: 371 for glob_ini in list_glob_ini: 372 # Check previous, and get the pipe name. 373 found, pipe_name, resfile, path = self.check_previous_result(method=self.method, model=model, analysis=analysis, glob_ini=glob_ini, bundle=self.method) 374 375 if not found: 376 calculate = True 377 finished[i] = False 378 elif found: 379 calculate = False 380 finished[i] = True 381 382 if calculate: 383 # Create the data pipe, by copying setup pipe. 384 self.interpreter.pipe.copy(pipe_from=self.base_setup_pipe_name, pipe_to=pipe_name, bundle_to=self.method) 385 self.interpreter.pipe.switch(pipe_name) 386 387 # Call set intensity. 388 finished_int = self.set_intensity_and_error(pipe_name=pipe_name, glob_ini=glob_ini, set_rmsd=set_rmsd) 389 390 if finished_int: 391 # Call error analysis. 392 self.do_spectrum_error_analysis(pipe_name=pipe_name, set_rep=set_rep) 393 394 # Save results, and store the current settings dic to pipe. 395 cdp.settings = self.settings 396 self.interpreter.results.write(file=resfile, dir=path, force=force) 397 398 finished[i] = True 399 400 else: 401 pipe_name = pipes.cdp_name() 402 self.interpreter.pipe.delete(pipe_name=pipe_name) 403 404 return all(finished)

405 406

407 - def calc_r2eff(self, methods=None, list_glob_ini=None, force=False):

408 """Method to calculate R2eff or read previous results.""" 409 410 model = MODEL_R2EFF 411 analysis = 'int' 412 413 # Loop over the methods. 414 for method in methods: 415 # Change the self key. 416 self.set_self(key='method', value=method) 417 418 # Loop over the glob ini: 419 for glob_ini in list_glob_ini: 420 # Check previous, and get the pipe name. 421 found, pipe_name, resfile, path = self.check_previous_result(method=self.method, model=model, analysis=analysis, glob_ini=glob_ini, bundle=self.method) 422 423 if not found: 424 calculate = True 425 elif found: 426 calculate = False 427 428 if calculate: 429 # Create the data pipe by copying the intensity pipe, then switching to it. 430 # If not intensity pipe name pipe exists, then calculate it. 431 intensity_pipe_name = self.name_pipe(method=self.method, model='setup', analysis='int', glob_ini=glob_ini) 432 433 if not pipes.has_pipe(intensity_pipe_name): 434 finished = self.set_int(methods=[method], list_glob_ini=[glob_ini]) 435 if not finished: 436 continue 437 438 self.interpreter.pipe.copy(pipe_from=intensity_pipe_name, pipe_to=pipe_name, bundle_to=self.method) 439 self.interpreter.pipe.switch(pipe_name) 440 441 # Select the model. 442 self.interpreter.relax_disp.select_model(model) 443 444 # Print 445 subtitle(file=sys.stdout, text="The '%s' model for pipe='%s'" % (model, pipe_name), prespace=3) 446 447 # Calculate the R2eff values for the fixed relaxation time period data types. 448 if model == MODEL_R2EFF and not has_exponential_exp_type(): 449 self.interpreter.minimise.calculate() 450 451 # Save results, and store the current settings dic to pipe. 452 cdp.settings = self.settings 453 self.interpreter.results.write(file=resfile, dir=path, force=force)

454 455

456 - def deselect_all(self, methods=None, model=None, model_from=None, analysis=None, analysis_from=None, list_glob_ini=None, force=False):

457 """Method to deselect all spins for a pipe.""" 458 459 # Set default 460 if model_from == None: 461 model_from = model 462 if analysis_from == None: 463 analysis_from = analysis 464 465 # Loop over the methods. 466 for method in methods: 467 # Change the self key. 468 self.set_self(key='method', value=method) 469 470 # Loop over the glob ini: 471 for glob_ini in list_glob_ini: 472 # Check previous, and get the pipe name. 473 found, pipe_name, resfile, path = self.check_previous_result(method=self.method, model=model, analysis=analysis, glob_ini=glob_ini, bundle=self.method) 474 475 if not found: 476 # If previous pipe not found, then create it. 477 model_from_pipe_name = self.name_pipe(method=self.method, model=model_from, analysis=analysis_from, glob_ini=glob_ini) 478 479 # Copy pipe and switch. 480 self.interpreter.pipe.copy(pipe_from=model_from_pipe_name, pipe_to=pipe_name, bundle_to=self.method) 481 self.interpreter.pipe.switch(pipe_name) 482 483 # Print 484 subtitle(file=sys.stdout, text="Deselect all spins for pipe='%s'" % (pipe_name), prespace=3) 485 486 # Deselect spins. 487 self.interpreter.deselect.all() 488 489 # Save results, and store the current settings dic to pipe. 490 cdp.settings = self.settings 491 492 if found and not force: 493 file_path = get_file_path(file_name=resfile, dir=path) 494 text = "The file '%s' already exists. Set the force flag to True to overwrite." % (file_path) 495 warn(RelaxWarning(text)) 496 else: 497 self.interpreter.results.write(file=resfile, dir=path, force=force) 498 499 # Show selected spins 500 self.interpreter.spin.display()

501 502

503 - def select_spin(self, spin_id=None, methods=None, model=None, model_from=None, analysis=None, analysis_from=None, list_glob_ini=None, force=False):

504 """Method to select spins for a pipe.""" 505 506 # Set default 507 if model_from == None: 508 model_from = model 509 if analysis_from == None: 510 analysis_from = analysis 511 512 # Loop over the methods. 513 for method in methods: 514 # Change the self key. 515 self.set_self(key='method', value=method) 516 517 # Loop over the glob ini: 518 for glob_ini in list_glob_ini: 519 # Check previous, and get the pipe name. 520 found, pipe_name, resfile, path = self.check_previous_result(method=self.method, model=model, analysis=analysis, glob_ini=glob_ini, bundle=self.method) 521 522 if not found: 523 # If previous pipe not found, then create it. 524 model_from_pipe_name = self.name_pipe(method=self.method, model=model_from, analysis=analysis_from, glob_ini=glob_ini) 525 526 # Copy pipe and switch. 527 self.interpreter.pipe.copy(pipe_from=model_from_pipe_name, pipe_to=pipe_name, bundle_to=self.method) 528 self.interpreter.pipe.switch(pipe_name) 529 530 # Print 531 subtitle(file=sys.stdout, text="Select spins '%s' for pipe='%s'" % (spin_id, pipe_name), prespace=3) 532 533 # Select spins. 534 self.interpreter.select.spin(spin_id=spin_id) 535 536 # Save results, and store the current settings dic to pipe. 537 cdp.settings = self.settings 538 539 if found and not force: 540 file_path = get_file_path(file_name=resfile, dir=path) 541 text = "The file '%s' already exists. Set the force flag to True to overwrite." % (file_path) 542 warn(RelaxWarning(text)) 543 else: 544 self.interpreter.results.write(file=resfile, dir=path, force=force) 545 546 # Show selected spins 547 self.interpreter.spin.display()

548 549

550 - def value_set(self, spin_id=None, val=None, param=None, methods=None, model=None, model_from=None, analysis=None, analysis_from=None, list_glob_ini=None, force=False):

551 """Use value.set on all pipes.""" 552 553 # Set default 554 if model_from == None: 555 model_from = model 556 if analysis_from == None: 557 analysis_from = analysis 558 559 # Loop over the methods. 560 for method in methods: 561 # Change the self key. 562 self.set_self(key='method', value=method) 563 564 # Loop over the glob ini: 565 for glob_ini in list_glob_ini: 566 # Check previous, and get the pipe name. 567 found, pipe_name, resfile, path = self.check_previous_result(method=self.method, model=model, analysis=analysis, glob_ini=glob_ini, bundle=self.method) 568 569 if not found: 570 # If previous pipe not found, then create it. 571 model_from_pipe_name = self.name_pipe(method=self.method, model=model_from, analysis=analysis_from, glob_ini=glob_ini) 572 573 # Copy pipe and switch. 574 self.interpreter.pipe.copy(pipe_from=model_from_pipe_name, pipe_to=pipe_name, bundle_to=self.method) 575 self.interpreter.pipe.switch(pipe_name) 576 577 # Print 578 subtitle(file=sys.stdout, text="For param '%s' set value '%3.2f' for pipe='%s'" % (param, val, pipe_name), prespace=3) 579 580 # Select the model. 581 self.interpreter.relax_disp.select_model(model) 582 583 # Set value 584 self.interpreter.value.set(val=val, param=param, spin_id=spin_id) 585 586 # Save results, and store the current settings dic to pipe. 587 cdp.settings = self.settings 588 589 if found and not force: 590 file_path = get_file_path(file_name=resfile, dir=path) 591 text = "The file '%s' already exists. Set the force flag to True to overwrite." % (file_path) 592 warn(RelaxWarning(text)) 593 else: 594 self.interpreter.results.write(file=resfile, dir=path, force=force) 595 596 # Print for pipe name 597 self.spin_display_params(pipe_name=pipe_name)

598 599

600 - def r20_from_min_r2eff(self, spin_id=None, methods=None, model=None, model_from=None, analysis=None, analysis_from=None, list_glob_ini=None, force=False):

601 """Will set the grid R20 values from the minimum R2eff values through the r20_from_min_r2eff user function. 602 This will speed up the grid search with a factor GRID_INC^(Nr_spec_freq). For a CPMG experiment with two fields and standard GRID_INC = 21, the speed-up is a factor 441. 603 """ 604 605 # Set default 606 if model_from == None: 607 model_from = model 608 if analysis_from == None: 609 analysis_from = analysis 610 611 # Loop over the methods. 612 for method in methods: 613 # Change the self key. 614 self.set_self(key='method', value=method) 615 616 # Loop over the glob ini: 617 for glob_ini in list_glob_ini: 618 # Check previous, and get the pipe name. 619 found, pipe_name, resfile, path = self.check_previous_result(method=self.method, model=model, analysis=analysis, glob_ini=glob_ini, bundle=self.method) 620 621 if not found: 622 # If previous pipe not found, then create it. 623 model_from_pipe_name = self.name_pipe(method=self.method, model=model_from, analysis=analysis_from, glob_ini=glob_ini) 624 625 # Copy pipe and switch. 626 self.interpreter.pipe.copy(pipe_from=model_from_pipe_name, pipe_to=pipe_name, bundle_to=self.method) 627 self.interpreter.pipe.switch(pipe_name) 628 629 # Print 630 subtitle(file=sys.stdout, text="Set grid r20 for pipe='%s'" % (pipe_name), prespace=3) 631 632 # Select the model. 633 self.interpreter.relax_disp.select_model(model) 634 635 # Set r20 from min r2eff. 636 self.interpreter.relax_disp.r20_from_min_r2eff(force=True) 637 638 # Save results, and store the current settings dic to pipe. 639 cdp.settings = self.settings 640 641 if found and not force: 642 file_path = get_file_path(file_name=resfile, dir=path) 643 text = "The file '%s' already exists. Set the force flag to True to overwrite." % (file_path) 644 warn(RelaxWarning(text)) 645 else: 646 self.interpreter.results.write(file=resfile, dir=path, force=force) 647 648 # Print for pipe name 649 self.spin_display_params(pipe_name=pipe_name)

650 651

652 - def minimise_grid_search(self, inc=11, verbosity=0, methods=None, model=None, model_from=None, analysis=None, analysis_from=None, list_glob_ini=None, force=False):

653 """Use value.set on all pipes.""" 654 655 # Set default 656 if model_from == None: 657 model_from = model 658 if analysis_from == None: 659 analysis_from = analysis 660 661 # Loop over the methods. 662 for method in methods: 663 # Change the self key. 664 self.set_self(key='method', value=method) 665 666 # Loop over the glob ini: 667 for glob_ini in list_glob_ini: 668 # Check previous, and get the pipe name. 669 found, pipe_name, resfile, path = self.check_previous_result(method=self.method, model=model, analysis=analysis, glob_ini=glob_ini, bundle=self.method) 670 671 # Try from analysis 672 if not found: 673 # Check previous, and get the pipe name. 674 found, pipe_name, resfile, path = self.check_previous_result(method=self.method, model=model, analysis=analysis_from, glob_ini=glob_ini, bundle=self.method) 675 676 if not found: 677 # If previous pipe not found, then create it. 678 model_from_pipe_name = self.name_pipe(method=self.method, model=model_from, analysis=analysis_from, glob_ini=glob_ini) 679 680 # Check if pipe exists. If not, try the R2eff pipe. 681 if not pipes.has_pipe(model_from_pipe_name): 682 model_from_pipe_name = self.name_pipe(method=self.method, model=MODEL_R2EFF, analysis='int', glob_ini=glob_ini) 683 684 # If not the R2eff pipe exist, then calculate it. 685 if not pipes.has_pipe(model_from_pipe_name): 686 self.calc_r2eff(methods=[self.method], list_glob_ini=[glob_ini]) 687 688 # Copy pipe and switch. 689 self.interpreter.pipe.copy(pipe_from=model_from_pipe_name, pipe_to=pipe_name, bundle_to=self.method) 690 self.interpreter.pipe.switch(pipe_name) 691 692 # Select the model. 693 self.interpreter.relax_disp.select_model(model) 694 695 # Deselect insignificant spins. 696 if model not in [MODEL_R2EFF, MODEL_NOREX]: 697 self.interpreter.relax_disp.insignificance(level=self.insignificance) 698 699 # Print 700 subtitle(file=sys.stdout, text="Grid search for pipe='%s'" % (pipe_name), prespace=3) 701 702 # Grid search. 703 if inc: 704 self.interpreter.minimise.grid_search(inc=inc, verbosity=verbosity, constraints=self.constraints, skip_preset=True) 705 706 # Default values. 707 else: 708 # The standard parameters. 709 for param in MODEL_PARAMS[model]: 710 self.interpreter.value.set(param=param, index=None, force=False) 711 712 # The optional R1 parameter. 713 if is_r1_optimised(model=model): 714 self.interpreter.value.set(param='r1', index=None) 715 716 717 # Save results, and store the current settings dic to pipe. 718 cdp.settings = self.settings 719 720 # Define new pipe names. 721 pipe_name = self.name_pipe(method=self.method, model=model, analysis=analysis, glob_ini=glob_ini) 722 resfile = pipe_name.replace(" ", "_") 723 model_path = model.replace(" ", "_") 724 path = self.results_dir+sep+model_path 725 726 if found and not force: 727 file_path = get_file_path(file_name=resfile, dir=path) 728 text = "The file '%s' already exists. Set the force flag to True to overwrite." % (file_path) 729 warn(RelaxWarning(text)) 730 else: 731 self.interpreter.results.write(file=resfile, dir=path, force=force) 732 733 # Print for pipe name 734 self.spin_display_params(pipe_name=pipe_name)

735 736

737 - def cluster_spins(self, spin_id=None, methods=None, model=None, model_from=None, analysis=None, analysis_from=None, list_glob_ini=None, force=False):

738 """Method to select spins for a pipe.""" 739 740 # Set default 741 if model_from == None: 742 model_from = model 743 if analysis_from == None: 744 analysis_from = analysis 745 746 # Loop over the methods. 747 for method in methods: 748 # Change the self key. 749 self.set_self(key='method', value=method) 750 751 # Loop over the glob ini: 752 for glob_ini in list_glob_ini: 753 # Check previous, and get the pipe name. 754 found, pipe_name, resfile, path = self.check_previous_result(method=self.method, model=model, analysis=analysis, glob_ini=glob_ini, bundle=self.method) 755 756 if not found: 757 # If previous pipe not found, then create it. 758 model_from_pipe_name = self.name_pipe(method=self.method, model=model_from, analysis=analysis_from, glob_ini=glob_ini) 759 760 # Copy pipe and switch. 761 self.interpreter.pipe.copy(pipe_from=model_from_pipe_name, pipe_to=pipe_name, bundle_to=self.method) 762 self.interpreter.pipe.switch(pipe_name) 763 764 # Print 765 subtitle(file=sys.stdout, text="Cluster spins '%s' for pipe='%s'" % (spin_id, pipe_name), prespace=3) 766 767 # Select spins. 768 self.interpreter.relax_disp.cluster(cluster_id='sel', spin_id=spin_id) 769 770 # Save results, and store the current settings dic to pipe. 771 cdp.settings = self.settings 772 773 if found and not force: 774 file_path = get_file_path(file_name=resfile, dir=path) 775 text = "The file '%s' already exists. Set the force flag to True to overwrite." % (file_path) 776 warn(RelaxWarning(text)) 777 else: 778 self.interpreter.results.write(file=resfile, dir=path, force=force) 779 780 # print clustered spins 781 print("Clustered spins are:", cdp.clustering)

782 783

784 - def minimise_execute(self, verbosity=1, methods=None, model=None, model_from=None, analysis=None, analysis_from=None, list_glob_ini=None, force=False, mc_err_analysis=False):

785 """Use value.set on all pipes.""" 786 787 # Set default 788 if model_from == None: 789 model_from = model 790 if analysis_from == None: 791 analysis_from = analysis 792 793 # Loop over the methods. 794 for method in methods: 795 # Change the self key. 796 self.set_self(key='method', value=method) 797 798 # Loop over the glob ini: 799 for glob_ini in list_glob_ini: 800 # Check previous, and get the pipe name. 801 found_pipe, pipe_name, resfile, path = self.check_previous_result(method=self.method, model=model, analysis=analysis, glob_ini=glob_ini, bundle=self.method) 802 803 # Try from analysis 804 if not found_pipe: 805 # Check previous, and get the pipe name. 806 found_analysis, pipe_name, resfile, path = self.check_previous_result(method=self.method, model=model, analysis=analysis_from, glob_ini=glob_ini, bundle=self.method) 807 808 if not (found_pipe or found_analysis): 809 # If previous pipe not found, then create it. 810 model_from_pipe_name = self.name_pipe(method=self.method, model=model_from, analysis=analysis_from, glob_ini=glob_ini) 811 812 # Check if pipe exists. If not, try grid pipe. 813 if not pipes.has_pipe(model_from_pipe_name): 814 model_from_pipe_name = self.name_pipe(method=self.method, model=model_from, analysis='grid', glob_ini=glob_ini) 815 816 # Copy pipe and switch. 817 self.interpreter.pipe.copy(pipe_from=model_from_pipe_name, pipe_to=pipe_name, bundle_to=self.method) 818 self.interpreter.pipe.switch(pipe_name) 819 820 # Select the model. 821 self.interpreter.relax_disp.select_model(model) 822 823 # Print 824 subtitle(file=sys.stdout, text="Minimise for pipe='%s'" % (pipe_name), prespace=3) 825 if hasattr(cdp, "sim_number"): 826 subsection(file=sys.stdout, text="Performing Monte-Carlo minimisations on %i simulations"%(getattr(cdp, "sim_number")), prespace=0) 827 828 # Do the minimisation. 829 self.interpreter.minimise.execute(min_algor=self.min_algor, func_tol=self.opt_func_tol, max_iter=self.opt_max_iterations, constraints=self.constraints, scaling=True, verbosity=verbosity) 830 831 # Do Monte-Carlo error analysis 832 if mc_err_analysis: 833 self.interpreter.monte_carlo.error_analysis() 834 835 # Save results, and store the current settings dic to pipe. 836 cdp.settings = self.settings 837 838 # Define new pipe names. 839 pipe_name = self.name_pipe(method=self.method, model=model, analysis=analysis, glob_ini=glob_ini) 840 resfile = pipe_name.replace(" ", "_") 841 model_path = model.replace(" ", "_") 842 path = self.results_dir+sep+model_path 843 844 if found_pipe and not force: 845 file_path = get_file_path(file_name=resfile, dir=path) 846 text = "The file '%s' already exists. Set the force flag to True to overwrite." % (file_path) 847 warn(RelaxWarning(text)) 848 else: 849 self.interpreter.results.write(file=resfile, dir=path, force=force) 850 851 # Print for pipe name 852 self.spin_display_params(pipe_name=pipe_name)

853 854

855 - def name_pipe(self, method, model, analysis, glob_ini, clusterid=None):

856 """Generate a unique name for the data pipe.""" 857 858 # Cluster group is none, set to standard free spins. 859 # cdp.clustering['free spins'] 860 if clusterid == None: 861 clusterid = 'free spins' 862 863 # The unique pipe name. 864 name = "%s - %s - %s - %s - %s" % (method, model, analysis, glob_ini, clusterid) 865 866 # Replace name with underscore. 867 name = name.replace(" ", "_") 868 869 # Return the name. 870 return name

871 872

873 - def check_previous_result(self, method=None, model=None, analysis=None, glob_ini=None, clusterid=None, bundle=None):

874 875 # Define if found and loaded 876 found = False 877 if bundle == None: 878 bundle = self.pipe_bundle 879 880 # Define the pipe name. 881 pipe_name = self.name_pipe(method=method, model=model, analysis=analysis, glob_ini=glob_ini, clusterid=clusterid) 882 883 # The results directory path. 884 model_path = model.replace(" ", "_") 885 path = self.results_dir+sep+model_path 886 887 # The result file. 888 resfile = pipe_name.replace(" ", "_") 889 890 # First check if the pipe already exists. Then switch to it. 891 if pipes.has_pipe(pipe_name): 892 print("Detected the presence of previous '%s' pipe - switching to it." % pipe_name) 893 self.interpreter.pipe.switch(pipe_name) 894 found = True 895 896 else: 897 # Check that results do not already exist - i.e. a previous run was interrupted. 898 path1 = get_file_path(file_name=resfile, dir=path) 899 path2 = path1 + '.bz2' 900 path3 = path1 + '.gz' 901 902 #print("Path to R2eff file is: %s"%path1) 903 if access(path1, F_OK) or access(path2, F_OK) or access(path2, F_OK): 904 # Printout. 905 print("Detected the presence of results files for the '%s' pipe - loading these instead of performing optimisation for a second time." % pipe_name) 906 907 # Create a data new pipe and switch to it. 908 self.interpreter.pipe.create(pipe_name=pipe_name, pipe_type=self.pipe_type, bundle=bundle) 909 self.interpreter.pipe.switch(pipe_name) 910 911 # Load the results. 912 self.interpreter.results.read(file=resfile, dir=path) 913 914 # Set found to True 915 found = True 916 917 return found, pipe_name, resfile, path

918 919

920 - def spin_display_params(self, spin_id=None, pipe_name=None):

921 """Display parameters for model in pipe.""" 922 923 924 # First check if the pipe already exists. Then switch to it. 925 if pipes.has_pipe(pipe_name): 926 # Switch to the pipe. 927 if pipes.cdp_name() != pipe_name: 928 print("Detected the presence of previous '%s' pipe - switching to it." % pipe_name) 929 self.interpreter.pipe.switch(pipe_name) 930 931 else: 932 # The result file. 933 pipe_name_split = pipe_name.split("_-_") 934 method = pipe_name_split[0] 935 model = pipe_name_split[1] 936 analysis = pipe_name_split[2] 937 bundle = method 938 939 model_path = model.replace(" ", "_") 940 path = self.results_dir+sep+model_path 941 # The result file. 942 resfile = pipe_name.replace(" ", "_") 943 944 # Check that results do not already exist - i.e. a previous run was interrupted. 945 path1 = get_file_path(file_name=resfile, dir=path) 946 path2 = path1 + '.bz2' 947 path3 = path1 + '.gz' 948 949 #print("Path to R2eff file is: %s"%path1) 950 if access(path1, F_OK) or access(path2, F_OK) or access(path2, F_OK): 951 # Printout. 952 print("Detected the presence of results files for the '%s' pipe - loading these instead of performing optimisation for a second time." % pipe_name) 953 954 # Create a data new pipe and switch to it. 955 self.interpreter.pipe.create(pipe_name=pipe_name, pipe_type=self.pipe_type, bundle=bundle) 956 self.interpreter.pipe.switch(pipe_name) 957 958 # Load the results. 959 self.interpreter.results.read(file=resfile, dir=path) 960 961 # Set found to True 962 found = True 963 964 965 # Start dic. 966 my_dic = {} 967 spin_id_list = [] 968 969 # Define data list. 970 data = [] 971 972 for cur_spin, mol_name, resi, resn, spin_id in spin_loop(selection=spin_id, full_info=True, return_id=True, skip_desel=True): 973 # Add key to dic. 974 my_dic[spin_id] = {} 975 # Store id, for ordering. 976 spin_id_list.append(spin_id) 977 978 # Add list to data. 979 cur_data_list = [repr(mol_name), repr(resi), repr(resn), repr(cur_spin.num), repr(cur_spin.name), spin_id] 980 981 # Get the parameters fitted in the model. 982 params = cur_spin.params 983 my_dic[spin_id]['params'] = params 984 my_dic[spin_id]['params_err'] = [] 985 986 # Loop over params. 987 for i, param in enumerate(params): 988 # Set the param error name 989 param_err = param + '_err' 990 my_dic[spin_id]['params_err'].append(param_err) 991 992 # If param in PARAMS_R20, values are stored in with parameter key. 993 param_key_list = [] 994 if param in PARAMS_R20: 995 # Loop over frq key. 996 for exp_type, frq, offset, ei, mi, oi, in loop_exp_frq_offset(return_indices=True): 997 # Get the parameter key. 998 param_key = generate_r20_key(exp_type=exp_type, frq=frq) 999 param_key_list.append(param_key) 1000 my_dic[spin_id]['param_key_list'] = param_key_list 1001 my_dic[spin_id][param] = {} 1002 1003 # Get the Value. 1004 if len(getattr(cur_spin, param)) == 0: 1005 param_val = None 1006 else: 1007 param_val = deepcopy(getattr(cur_spin, param)[param_key]) 1008 my_dic[spin_id][param][param_key] = param_val 1009 1010 # Add information to data. 1011 if param_val == None: 1012 cur_data_list.append("%s"%param_val) 1013 else: 1014 cur_data_list.append("%3.3f"%param_val) 1015 1016 else: 1017 # Get the value. 1018 param_val = deepcopy(getattr(cur_spin, param)) 1019 my_dic[spin_id][param] = param_val 1020 1021 # Add information to data. 1022 if param_val == None: 1023 cur_data_list.append("%s"%param_val) 1024 else: 1025 cur_data_list.append("%.3f"%param_val) 1026 1027 # Add data. 1028 data.append(cur_data_list) 1029 1030 # Collect header from spin 0. 1031 cur_spin_id = spin_id_list[0] 1032 cur_spin_params = my_dic[cur_spin_id]['params'] 1033 cur_param_keys = my_dic[cur_spin_id]['param_key_list'] 1034 1035 # Define header. 1036 param_header = ["Molecule", "Res number", "Res name", "Spin number", "Spin name", "Spin id"] 1037 1038 # Loop over params, add to header. 1039 for param in cur_spin_params: 1040 if param in PARAMS_R20: 1041 for param_key in cur_param_keys: 1042 # Take the second last part of key. 1043 cur_key = "%3.1f" % float(param_key.split()[-2]) 1044 hstring = "%s_%s" % (param, cur_key) 1045 param_header.append(hstring) 1046 else: 1047 hstring = "%s" % (param) 1048 param_header.append(hstring) 1049 1050 write_data(out=sys.stdout, headings=param_header, data=data)

1051 1052

1053 - def get_dublicates(self, spectrum_ids, cpmg_frqs):

1054 """Method which return a list of tubles, where each tuble is a spectrum id and a list of spectrum ids which are replicated""" 1055 1056 # Get the dublicates. 1057 dublicates = [(val, [i for i in range(len(cpmg_frqs)) if cpmg_frqs[i] == val]) for val in cpmg_frqs] 1058 1059 # Loop over the list of the mapping of cpmg frequency and duplications. 1060 list_dub_mapping = [] 1061 for i, dub in enumerate(dublicates): 1062 # Get current spectum id. 1063 cur_spectrum_id = spectrum_ids[i] 1064 1065 # Get the tuple of cpmg frequency and indexes of duplications. 1066 cpmg_frq, list_index_occur = dub 1067 1068 # Collect mapping of index to id. 1069 id_list = [] 1070 if len(list_index_occur) > 1: 1071 for list_index in list_index_occur: 1072 id_list.append(spectrum_ids[list_index]) 1073 1074 # Store to list 1075 list_dub_mapping.append((cur_spectrum_id, id_list)) 1076 1077 return list_dub_mapping

1078 1079

1080 - def col_int(self, method=None, list_glob_ini=None, selection=None):

1081 1082 # Loop over the glob ini: 1083 res_dic = {} 1084 res_dic['method'] = method 1085 res_dic['selection'] = selection 1086 for glob_ini in list_glob_ini: 1087 # Get the pipe name for peak_intensity values. 1088 pipe_name = self.name_pipe(method=method, model='setup', analysis='int', glob_ini=glob_ini) 1089 1090 # Check if pipe exists, or else calculate. 1091 if not pipes.has_pipe(pipe_name): 1092 self.set_int(methods=[method], list_glob_ini=[glob_ini]) 1093 1094 if pipes.cdp_name() != pipe_name and pipes.has_pipe(pipe_name): 1095 self.interpreter.pipe.switch(pipe_name) 1096 1097 elif pipes.has_pipe(pipe_name) == False: 1098 continue 1099 1100 # Results dictionary. 1101 res_dic[str(glob_ini)] = {} 1102 res_dic[str(glob_ini)]['peak_intensity'] = {} 1103 res_dic[str(glob_ini)]['peak_intensity_err'] = {} 1104 spin_point_peak_intensity_list = [] 1105 spin_point_peak_intensity_err_list = [] 1106 1107 # Loop over the spins. 1108 for cur_spin, mol_name, resi, resn, spin_id in spin_loop(selection=selection, full_info=True, return_id=True, skip_desel=True): 1109 # Make spin dic. 1110 res_dic[str(glob_ini)]['peak_intensity'][spin_id] = {} 1111 res_dic[str(glob_ini)]['peak_intensity_err'][spin_id] = {} 1112 1113 # Loop over spectrum_ids. 1114 for s_id in cdp.spectrum_ids: 1115 # Check for bad data has skipped peak_intensity points 1116 if s_id in cur_spin.peak_intensity: 1117 peak_intensity_point = cur_spin.peak_intensity[s_id] 1118 peak_intensity_err_point = cur_spin.peak_intensity_err[s_id] 1119 1120 res_dic[str(glob_ini)]['peak_intensity'][spin_id][s_id] = peak_intensity_point 1121 res_dic[str(glob_ini)]['peak_intensity_err'][spin_id][s_id] = peak_intensity_err_point 1122 spin_point_peak_intensity_list.append(peak_intensity_point) 1123 spin_point_peak_intensity_err_list.append(peak_intensity_err_point) 1124 1125 res_dic[str(glob_ini)]['peak_intensity_arr'] = asarray(spin_point_peak_intensity_list) 1126 res_dic[str(glob_ini)]['peak_intensity_err_arr'] = asarray(spin_point_peak_intensity_err_list) 1127 1128 return res_dic

1129 1130

1131 - def plot_int_corr(self, corr_data, show=False, write_stats=False):

1132 1133 # Define figure. 1134 # Nr of columns is number of datasets. 1135 nr_cols = len(corr_data) 1136 # Nr of rows, is 2. With and without scaling. 1137 nr_rows = 4 1138 1139 # Define figure 1140 fig, axises = plt.subplots(nrows=nr_rows, ncols=nr_cols) 1141 fig.suptitle('Correlation plot') 1142 1143 # axises is a tuple with number of elements corresponding to number of rows. 1144 # Each sub-tuple contains axis for each column. 1145 1146 # For writing out stats. 1147 data_dic = {} 1148 1149 # Loop over the rows. 1150 for i, row_axises in enumerate(axises): 1151 # Loop over the columns. 1152 for j, ax in enumerate(row_axises) : 1153 # Extract from lists. 1154 data, methods, glob_inis = corr_data[j] 1155 data_x, data_y = data 1156 method_x, method_y = methods 1157 glob_ini_x, glob_ini_y = glob_inis 1158 1159 x = data_x[str(glob_ini_x)]['peak_intensity_arr'] 1160 x_err = data_x[str(glob_ini_x)]['peak_intensity_err_arr'] 1161 np = len(x) 1162 1163 y = data_y[str(glob_ini_y)]['peak_intensity_arr'] 1164 y_err = data_y[str(glob_ini_y)]['peak_intensity_err_arr'] 1165 1166 # If row 1. 1167 if i == 0: 1168 # Add to data dic. 1169 method_xy_NI = "int_%s%s_%s%s" % (method_x, glob_ini_x, method_y, glob_ini_y) 1170 data_dic[method_xy_NI] = [] 1171 1172 ax.plot(x, x, 'o', label='%s vs. %s' % (method_x, method_x)) 1173 ax.plot(x, y, '.', label='%s vs. %s' % (method_y, method_x) ) 1174 1175 np = len(y) 1176 ax.set_title(r'$I$' + ' for %s %i vs. %s %i. np=%i' % (method_y, glob_ini_y, method_x, glob_ini_x, np), fontsize=10) 1177 ax.legend(loc='upper center', shadow=True, prop=fontP) 1178 ax.ticklabel_format(style='sci', axis='x', scilimits=(0, 0)) 1179 ax.ticklabel_format(style='sci', axis='y', scilimits=(0, 0)) 1180 ax.set_xlabel(r'$I$') 1181 ax.set_ylabel(r'$I$') 1182 1183 # Calculate straight line. 1184 # Linear a, with no intercept. 1185 a = sum(x * y) / sum(x**2) 1186 min_x = min(x) 1187 max_x = max(x) 1188 step_x = (max_x - min_x) / np 1189 x_arange = arange(min_x, max_x, step_x) 1190 y_arange = a * x_arange 1191 1192 # Add to data. 1193 for k, x_k in enumerate(x): 1194 y_k = y[k] 1195 x_arange_k = x_arange[k] 1196 y_arange_k = y_arange[k] 1197 data_dic[method_xy_NI].append(["%3.5f"%x_k, "%3.5f"%y_k, "%3.5f"%x_arange_k, "%3.5f"%y_arange_k]) 1198 1199 # Scale intensity 1200 if i == 1: 1201 1202 x_norm = x / x.max() 1203 y_norm = y / y.max() 1204 1205 ax.plot(x_norm, x_norm, 'o', label='%s vs. %s' % (method_x, method_x)) 1206 ax.plot(x_norm, y_norm, '.', label='%s vs. %s' % (method_y, method_x)) 1207 1208 np = len(y_norm) 1209 ax.set_title('Normalised intensity for %s %i vs. %s %i. np=%i' % (method_y, glob_ini_y, method_x, glob_ini_x, np), fontsize=10) 1210 ax.legend(loc='upper center', shadow=True, prop = fontP) 1211 ax.set_xlabel(r'$\mathrm{Norm.} I$') 1212 ax.set_ylabel(r'$\mathrm{Norm.} I$') 1213 1214 1215 # Error. 1216 if i == 2: 1217 # Add to data dic. 1218 method_xy_NI = "int_err_%s%s_%s%s" % (method_x, glob_ini_x, method_y, glob_ini_y) 1219 data_dic[method_xy_NI] = [] 1220 1221 ax.plot(x_err, x_err, 'o', label='%s vs. %s' % (method_x, method_x)) 1222 ax.plot(x_err, y_err, '.', label='%s vs. %s' % (method_y, method_x)) 1223 1224 np = len(y_err) 1225 ax.set_title(r'$\sigma(I)$' + ' for %s %i vs. %s %i. np=%i' % (method_y, glob_ini_y, method_x, glob_ini_x, np), fontsize=10) 1226 ax.legend(loc='upper center', shadow=True, prop = fontP) 1227 ax.set_xlabel(r'$\sigma(I)$') 1228 ax.set_ylabel(r'$\sigma(I)$') 1229 1230 # Calculate straight line. 1231 # Linear a, with no intercept. 1232 a = sum(x_err * y_err) / sum(x_err**2) 1233 min_x = min(x_err) 1234 max_x = max(x_err) 1235 step_x = (max_x - min_x) / np 1236 x_err_arange = arange(min_x, max_x, step_x) 1237 y_err_arange = a * x_err_arange 1238 1239 # Add to data. 1240 for k, x_err_k in enumerate(x_err): 1241 y_err_k = y_err[k] 1242 x_err_arange_k = x_err_arange[k] 1243 y_err_arange_k = y_err_arange[k] 1244 data_dic[method_xy_NI].append(["%3.5f"%x_err_k, "%3.5f"%y_err_k, "%3.5f"%x_err_arange_k, "%3.5f"%y_err_arange_k]) 1245 1246 1247 # Intensity to error. 1248 if i == 3: 1249 1250 x_to_x_err = x / x_err 1251 y_to_y_err = y / y_err 1252 1253 ax.plot(x_to_x_err, x_to_x_err, 'o', label='%s vs. %s' % (method_x, method_x)) 1254 ax.plot(x_to_x_err, y_to_y_err, '.', label='%s vs. %s' % (method_y, method_x)) 1255 1256 np = len(y_to_y_err) 1257 ax.set_title(r'$I/\sigma(I)$' + ' for %s %i vs. %s %i. np=%i' % (method_y, glob_ini_y, method_x, glob_ini_x, np), fontsize=10) 1258 ax.legend(loc='upper center', shadow=True, prop = fontP) 1259 ax.set_xlabel(r'$I/\sigma(I)$') 1260 ax.set_ylabel(r'$I/\sigma(I)$') 1261 1262 plt.tight_layout() 1263 1264 # Loop over columns for writing data. 1265 # Write to file. 1266 if write_stats: 1267 # Re-order the data. 1268 headings_all = [] 1269 method_xy_NI_all = [] 1270 # Loop over the columns. 1271 for j in range(nr_cols): 1272 headings_j = [] 1273 method_xy_NI_j = [] 1274 # Loop over rows 1275 for i in range(nr_rows): 1276 # Extract from lists. 1277 data, methods, glob_inis = corr_data[j] 1278 method_x, method_y = methods 1279 glob_ini_x, glob_ini_y = glob_inis 1280 1281 # If row 1. 1282 if i == 0: 1283 # Add to headings. 1284 method_x_NI = "int_%s%s" % (method_x, glob_ini_x) 1285 method_y_NI = "int_%s%s" % (method_y, glob_ini_y) 1286 method_x_NI_lin = "int_lin_%s%s" % (method_x, glob_ini_x) 1287 method_y_NI_lin = "int_lin_%s%s" % (method_y, glob_ini_y) 1288 headings_j = headings_j + [method_x_NI, method_y_NI, method_x_NI_lin, method_y_NI_lin] 1289 1290 method_xy_NI = "int_%s%s_%s%s" % (method_x, glob_ini_x, method_y, glob_ini_y) 1291 method_xy_NI_j.append(method_xy_NI) 1292 1293 # Error. 1294 if i == 2: 1295 # Add to headings 1296 method_x_NI = "int_err_%s%s" % (method_x, glob_ini_x) 1297 method_y_NI = "int_err_%s%s" % (method_y, glob_ini_y) 1298 method_x_NI_lin = "int_err_lin_%s%s" % (method_x, glob_ini_x) 1299 method_y_NI_lin = "int_err_lin_%s%s" % (method_y, glob_ini_y) 1300 headings_j = headings_j + [method_x_NI, method_y_NI, method_x_NI_lin, method_y_NI_lin] 1301 1302 method_xy_NI = "int_err_%s%s_%s%s" % (method_x, glob_ini_x, method_y, glob_ini_y) 1303 method_xy_NI_j.append(method_xy_NI) 1304 1305 headings_all.append(headings_j) 1306 method_xy_NI_all.append(method_xy_NI_j) 1307 1308 # Loop over the columns. 1309 for j, headings_j in enumerate(headings_all): 1310 method_xy_NI_j = method_xy_NI_all[j] 1311 1312 data_w = [] 1313 data_int = data_dic[method_xy_NI_j[0]] 1314 data_int_err = data_dic[method_xy_NI_j[1]] 1315 1316 for k, data_int_k in enumerate(data_int): 1317 data_int_err_k = data_int_err[k] 1318 data_w.append(data_int_k + data_int_err_k) 1319 1320 # Define file name. 1321 data, methods, glob_inis = corr_data[j] 1322 data_x, data_y = data 1323 method_x, method_y = methods 1324 glob_ini_x, glob_ini_y = glob_inis 1325 np = len(data_int) 1326 1327 # Get the spin selection for correlation. 1328 selection = data_x['selection'] 1329 1330 file_name_ini = 'int_corr_%s_%s_%s_%s_NP_%i' % (method_x, glob_ini_x, method_y, glob_ini_y, np) 1331 if selection == None: 1332 file_name_ini = file_name_ini + '_all' 1333 else: 1334 file_name_ini = file_name_ini + '_sel' 1335 1336 file_name = file_name_ini + '.txt' 1337 path = self.results_dir 1338 1339 # save figure 1340 # Write png. 1341 png_file_name = file_name_ini + '.png' 1342 png_file_path = get_file_path(file_name=png_file_name, dir=path) 1343 plt.savefig(png_file_path, bbox_inches='tight') 1344 1345 # Write file 1346 file_obj, file_path = open_write_file(file_name=file_name, dir=path, force=True, compress_type=0, verbosity=1, return_path=True) 1347 1348 # Write data. 1349 write_data(out=file_obj, headings=headings_j, data=data_w) 1350 1351 # Close file. 1352 file_obj.close() 1353 1354 if show: 1355 plt.show()

1356 1357

1358 - def get_int_stat_dic(self, list_int_dics=None, list_glob_ini=None):

1359 1360 # Loop over the result dictionaries: 1361 res_dic = {} 1362 for i, int_dic in enumerate(list_int_dics): 1363 # Let the reference dic be initial dic 1364 int_dic_ref = list_int_dics[0] 1365 method_ref = int_dic_ref['method'] 1366 res_dic['method_ref'] = method_ref 1367 glob_ini_ref = list_glob_ini[0] 1368 res_dic['glob_ini_ref'] = str(glob_ini_ref) 1369 selection = int_dic_ref['selection'] 1370 res_dic['selection'] = selection 1371 1372 # Let the reference int array be the initial glob. 1373 int_arr_ref = int_dic_ref[str(glob_ini_ref)]['peak_intensity_arr'] 1374 res_dic['int_arr_ref'] = int_arr_ref 1375 int_err_arr_ref = int_dic_ref[str(glob_ini_ref)]['peak_intensity_err_arr'] 1376 res_dic['int_err_arr_ref'] = int_err_arr_ref 1377 1378 # Get the current method 1379 method_cur = int_dic['method'] 1380 res_dic[method_cur] = {} 1381 res_dic[method_cur]['method'] = method_cur 1382 res_dic[method_cur]['sampling_sparseness'] = [] 1383 res_dic[method_cur]['glob_ini'] = [] 1384 res_dic[method_cur]['int_norm_std'] = [] 1385 1386 # Other stats. 1387 res_dic[method_cur]['r_xy_int'] = [] 1388 res_dic[method_cur]['a_int'] = [] 1389 res_dic[method_cur]['r_xy_int_err'] = [] 1390 res_dic[method_cur]['a_int_err'] = [] 1391 1392 # Now loop over glob_ini: 1393 for glob_ini in list_glob_ini: 1394 # Get the array, if it exists. 1395 if str(glob_ini) not in int_dic: 1396 continue 1397 1398 # Get the data. 1399 int_arr = int_dic[str(glob_ini)]['peak_intensity_arr'] 1400 int_err_arr = int_dic[str(glob_ini)]['peak_intensity_err_arr'] 1401 1402 # This require that all number of points are equal. 1403 # If they are not of same length, then dont even bother to continue. 1404 if len(int_arr) != len(int_arr_ref): 1405 continue 1406 1407 # Store x 1408 sampling_sparseness = float(glob_ini) / float(glob_ini_ref) * 100. 1409 res_dic[method_cur]['sampling_sparseness'].append(sampling_sparseness) 1410 res_dic[method_cur]['glob_ini'].append(glob_ini) 1411 1412 # Store to result dic. 1413 res_dic[method_cur][str(glob_ini)] = {} 1414 res_dic[method_cur][str(glob_ini)]['sampling_sparseness'] = sampling_sparseness 1415 res_dic[method_cur][str(glob_ini)]['int_arr'] = int_arr 1416 res_dic[method_cur][str(glob_ini)]['int_err_arr'] = int_err_arr 1417 1418 # Calculate sample correlation coefficient, measure of goodness-of-fit of linear regression 1419 # Without intercept. 1420 x = int_arr_ref 1421 y = int_arr 1422 1423 a_int = sum(x*y) / sum(x**2) 1424 r_xy_int = sum(x*y) / sqrt(sum(x**2) * sum(y**2)) 1425 1426 x = int_err_arr_ref 1427 y = int_err_arr 1428 a_int_err = sum(x*y) / sum(x**2) 1429 r_xy_int_err = sum(x*y) / sqrt(sum(x**2) * sum(y**2)) 1430 1431 print(method_ref, method_cur, sampling_sparseness, glob_ini, r_xy_int**2, a_int, r_xy_int_err**2, a_int_err) 1432 1433 # Store to result dic. 1434 res_dic[method_cur][str(glob_ini)]['r_xy_int'] = r_xy_int 1435 res_dic[method_cur]['r_xy_int'].append(r_xy_int) 1436 res_dic[method_cur][str(glob_ini)]['a_int'] = a_int 1437 res_dic[method_cur]['a_int'].append(a_int) 1438 1439 res_dic[method_cur][str(glob_ini)]['r_xy_int_err'] = r_xy_int_err 1440 res_dic[method_cur]['r_xy_int_err'].append(r_xy_int_err) 1441 res_dic[method_cur][str(glob_ini)]['a_int_err'] = a_int_err 1442 res_dic[method_cur]['a_int_err'].append(a_int_err) 1443 1444 res_dic[method_cur]['sampling_sparseness'] = asarray(res_dic[method_cur]['sampling_sparseness']) 1445 res_dic[method_cur]['glob_ini'] = asarray(res_dic[method_cur]['glob_ini']) 1446 1447 res_dic[method_cur]['r_xy_int'] = asarray(res_dic[method_cur]['r_xy_int']) 1448 res_dic[method_cur]['a_int'] = asarray(res_dic[method_cur]['a_int']) 1449 res_dic[method_cur]['r_xy_int_err'] = asarray(res_dic[method_cur]['r_xy_int_err']) 1450 res_dic[method_cur]['a_int_err'] = asarray(res_dic[method_cur]['a_int_err']) 1451 1452 return res_dic

1453 1454

1455 - def plot_int_stat(self, int_stat_dic=None, methods=[], list_glob_ini=[], show=False, write_stats=False):

1456 1457 # Define figure 1458 fig, axises = plt.subplots(nrows=2, ncols=1) 1459 fig.suptitle('Stats per NI') 1460 ax1, ax2 = axises 1461 1462 # Catch min and max values for all methods. 1463 min_a = 1.0 1464 max_a = 0.0 1465 1466 min_r_xy2 = 1.0 1467 max_r_xy2 = 0.0 1468 1469 # Prepare header for writing. 1470 selection = int_stat_dic['selection'] 1471 1472 # For writing out stats. 1473 headings = [] 1474 data_dic = {} 1475 data_dic_methods = [] 1476 i_max = 0 1477 1478 for method in methods: 1479 if method not in int_stat_dic: 1480 continue 1481 1482 # Use NI as x. 1483 NI = int_stat_dic[method]['glob_ini'] 1484 # Use sampling_sparseness as x. 1485 SS = int_stat_dic[method]['sampling_sparseness'] 1486 1487 # Add to headings. 1488 headings = headings + ['method', 'SS', 'NI', 'slope_int', 'rxy2_int', 'slope_int_err', 'rxy2_int_err'] 1489 1490 # Get stats. 1491 # Linear regression slope, without intercept 1492 a_int = int_stat_dic[method]['a_int'] 1493 1494 if max(a_int) > max_a: 1495 max_a = max(a_int) 1496 if min(a_int) < min_a: 1497 min_a = min(a_int) 1498 1499 # sample correlation coefficient, without intercept 1500 r_xy_int = int_stat_dic[method]['r_xy_int'] 1501 r_xy_int2 = r_xy_int**2 1502 1503 if max(r_xy_int2) > max_r_xy2: 1504 max_r_xy2 = max(r_xy_int2) 1505 if min(r_xy_int2) < min_r_xy2: 1506 min_r_xy2 = min(r_xy_int2) 1507 1508 # For just the int values 1509 a_int_err = int_stat_dic[method]['a_int_err'] 1510 r_xy_int_err = int_stat_dic[method]['r_xy_int_err'] 1511 r_xy_int_err2 = r_xy_int_err**2 1512 1513 # Add to data. 1514 data_dic[method] = {} 1515 data_dic_methods.append(method) 1516 for i, NI_i in enumerate(NI): 1517 SS_i = SS[i] 1518 a_int_i = a_int[i] 1519 r_xy_int2_i = r_xy_int2[i] 1520 a_int_err_i = a_int_err[i] 1521 r_xy_int_err2_i = r_xy_int_err2[i] 1522 data_dic[method][str(i)] = ["%3.5f"%SS_i, "%i"%NI_i, "%3.5f"%a_int_i, "%3.5f"%r_xy_int2_i, "%3.5f"%a_int_err_i, "%3.5f"%r_xy_int_err2_i] 1523 if i > i_max: 1524 i_max = i 1525 1526 t = ax1.plot(SS, a_int, ".--", label='%s slope int'%method) 1527 color = t[0].get_color() 1528 ax1.plot(SS, a_int_err, ".-", label='%s slope int_err'%method, color=color) 1529 1530 t = ax2.plot(SS, r_xy_int2, "o--", label='%s r2 int'%method) 1531 color = t[0].get_color() 1532 ax2.plot(SS, r_xy_int_err2, "o-", label='%s r2 int_err'%method, color=color) 1533 1534 # Loop over methods for writing data. 1535 data = [] 1536 1537 for i in range(0, i_max+1): 1538 data_i = [] 1539 for method in data_dic_methods: 1540 data_dic_m = data_dic[method] 1541 # Loop over all possible data points. 1542 if str(i) in data_dic_m: 1543 data_i = data_i + [method] + data_dic_m[str(i)] 1544 else: 1545 data_i = data_i + [method] + ["0", "0", "0", "0", "0", "0"] 1546 1547 data.append(data_i) 1548 1549 # Set legends. 1550 ax1.legend(loc='lower left', shadow=True, prop = fontP) 1551 #ax1.set_xlabel('NI') 1552 ax1.set_xlabel('SS') 1553 #ax1.set_ylabel(r'$\sigma ( R_{2,\mathrm{eff}} )$') 1554 ax1.set_ylabel('Linear regression slope, without intercept') 1555 #ax1.set_xticks(NI) 1556 #ax1.set_xticks(SS) 1557 ax1.set_ylim(min_a*0.95, max_a*1.05) 1558 ax1.invert_xaxis() 1559 1560 ax2.legend(loc='lower right', shadow=True, prop = fontP) 1561 ax2.set_ylabel('Sample correlation ' + r'$r_{xy}^2$') 1562 #ax2.set_xticks(NI) 1563 #ax2.set_xticks(SS) 1564 ax2.set_ylim(min_r_xy2*0.95, max_r_xy2*1.05) 1565 ax2.invert_xaxis() 1566 1567 # Determine filename. 1568 if selection == None: 1569 file_name_ini = 'int_stat_all' 1570 else: 1571 file_name_ini = 'int_stat_sel' 1572 1573 # Write png. 1574 png_file_name = file_name_ini + '.png' 1575 png_file_path = get_file_path(file_name=png_file_name, dir=self.results_dir) 1576 1577 # Write to file. 1578 if write_stats: 1579 # save figure 1580 plt.savefig(png_file_path, bbox_inches='tight') 1581 1582 file_name = file_name_ini + '.txt' 1583 path = self.results_dir 1584 file_obj, file_path = open_write_file(file_name=file_name, dir=path, force=True, compress_type=0, verbosity=1, return_path=True) 1585 1586 # Write data. 1587 write_data(out=file_obj, headings=headings, data=data) 1588 1589 # Close file. 1590 file_obj.close() 1591 1592 # Plot data. 1593 if show: 1594 plt.show()

1595 1596

1597 - def col_r2eff(self, method=None, list_glob_ini=None, selection=None):

1598 1599 # Loop over the glob ini: 1600 res_dic = {} 1601 res_dic['method'] = method 1602 res_dic['selection'] = selection 1603 for glob_ini in list_glob_ini: 1604 # Get the pipe name for R2eff values. 1605 pipe_name = self.name_pipe(method=method, model=MODEL_R2EFF, analysis='int', glob_ini=glob_ini) 1606 1607 # Check if pipe exists, or else calculate. 1608 if not pipes.has_pipe(pipe_name): 1609 self.calc_r2eff(methods=[method], list_glob_ini=[glob_ini]) 1610 1611 if pipes.cdp_name() != pipe_name and pipes.has_pipe(pipe_name): 1612 self.interpreter.pipe.switch(pipe_name) 1613 1614 elif pipes.has_pipe(pipe_name) == False: 1615 continue 1616 1617 # Results dictionary. 1618 res_dic[str(glob_ini)] = {} 1619 res_dic[str(glob_ini)]['r2eff'] = {} 1620 res_dic[str(glob_ini)]['r2eff_err'] = {} 1621 spin_point_r2eff_list = [] 1622 spin_point_r2eff_err_list = [] 1623 1624 # Loop over the spins. 1625 for cur_spin, mol_name, resi, resn, spin_id in spin_loop(selection=selection, full_info=True, return_id=True, skip_desel=True): 1626 # Make spin dic. 1627 res_dic[str(glob_ini)]['r2eff'][spin_id] = {} 1628 res_dic[str(glob_ini)]['r2eff_err'][spin_id] = {} 1629 1630 # Loop over the R2eff points 1631 for exp_type, frq, offset, point, ei, mi, oi, di in loop_exp_frq_offset_point(return_indices=True): 1632 # Define the key. 1633 data_key = return_param_key_from_data(exp_type=self.exp_type, frq=frq, offset=offset, point=point) 1634 1635 # Check for bad data has skipped r2eff points 1636 if data_key in cur_spin.r2eff: 1637 r2eff_point = cur_spin.r2eff[data_key] 1638 r2eff_err_point = cur_spin.r2eff_err[data_key] 1639 1640 res_dic[str(glob_ini)]['r2eff'][spin_id][data_key] = r2eff_point 1641 res_dic[str(glob_ini)]['r2eff_err'][spin_id][data_key] = r2eff_err_point 1642 spin_point_r2eff_list.append(r2eff_point) 1643 spin_point_r2eff_err_list.append(r2eff_err_point) 1644 1645 res_dic[str(glob_ini)]['r2eff_arr'] = asarray(spin_point_r2eff_list) 1646 res_dic[str(glob_ini)]['r2eff_err_arr'] = asarray(spin_point_r2eff_err_list) 1647 1648 return res_dic

1649 1650

1651 - def plot_r2eff_corr(self, corr_data, show=False, write_stats=False):

1652 1653 # Define figure. 1654 # Nr of columns is number of datasets. 1655 nr_cols = len(corr_data) 1656 # Nr of rows, is 2. With and without scaling. 1657 nr_rows = 2 1658 1659 # Define figure 1660 fig, axises = plt.subplots(nrows=nr_rows, ncols=nr_cols) 1661 fig.suptitle('Correlation plot') 1662 1663 # axises is a tuple with number of elements corresponding to number of rows. 1664 # Each sub-tuple contains axis for each column. 1665 1666 # For writing out stats. 1667 data_dic = {} 1668 1669 # Loop over the rows. 1670 for i, row_axises in enumerate(axises): 1671 # Loop over the columns. 1672 for j, ax in enumerate(row_axises): 1673 # Extract from lists. 1674 data, methods, glob_inis = corr_data[j] 1675 data_x, data_y = data 1676 method_x, method_y = methods 1677 glob_ini_x, glob_ini_y = glob_inis 1678 1679 x = data_x[str(glob_ini_x)]['r2eff_arr'] 1680 x_err = data_x[str(glob_ini_x)]['r2eff_err_arr'] 1681 np = len(x) 1682 1683 y = data_y[str(glob_ini_y)]['r2eff_arr'] 1684 y_err = data_y[str(glob_ini_y)]['r2eff_err_arr'] 1685 1686 # If row 1. 1687 if i == 0: 1688 # Add to data dic. 1689 method_xy_NI = "r2eff_%s%s_%s%s" % (method_x, glob_ini_x, method_y, glob_ini_y) 1690 data_dic[method_xy_NI] = [] 1691 1692 # Calculate straight line. 1693 # Linear a, with no intercept. 1694 a = sum(x * y) / sum(x**2) 1695 min_x = min(x) 1696 max_x = max(x) 1697 step_x = (max_x - min_x) / np 1698 x_arange = arange(min_x, max_x, step_x) 1699 y_arange = a * x_arange 1700 1701 ax.plot(x, x, 'o', label='%s vs. %s' % (method_x, method_x)) 1702 ax.plot(x_arange, x_arange, 'b--') 1703 if len(x) != len(y): 1704 print(len(x), len(y)) 1705 ax.plot(x, y, '.', label='%s vs. %s' % (method_y, method_x) ) 1706 ax.plot(x_arange, y_arange, 'g--') 1707 1708 ax.set_title(r'$R_{2,\mathrm{eff}}$' + ' for %s %i vs. %s %i. np=%i' % (method_y, glob_ini_y, method_x, glob_ini_x, np), fontsize=10) 1709 ax.legend(loc='upper left', shadow=True, prop = fontP) 1710 ax.ticklabel_format(style='sci', axis='x', scilimits=(0, 0)) 1711 ax.ticklabel_format(style='sci', axis='y', scilimits=(0, 0)) 1712 ax.set_xlabel(r'$R_{2,\mathrm{eff}}$') 1713 ax.set_ylabel(r'$R_{2,\mathrm{eff}}$') 1714 1715 # Add to data. 1716 for k, x_k in enumerate(x): 1717 y_k = y[k] 1718 x_arange_k = x_arange[k] 1719 y_arange_k = y_arange[k] 1720 data_dic[method_xy_NI].append(["%3.5f"%x_k, "%3.5f"%y_k, "%3.5f"%x_arange_k, "%3.5f"%y_arange_k]) 1721 1722 # R2eff to error. 1723 if i == 1: 1724 # Add to data dic. 1725 method_xy_NI = "r2eff_to_err_%s%s_%s%s" % (method_x, glob_ini_x, method_y, glob_ini_y) 1726 data_dic[method_xy_NI] = [] 1727 1728 x_to_x_err = x / x_err 1729 y_to_y_err = y / y_err 1730 1731 # Calculate straight line. 1732 # Linear a, with no intercept. 1733 a = sum(x_to_x_err * y_to_y_err) / sum(x_to_x_err**2) 1734 min_x = min(x_to_x_err) 1735 max_x = max(x_to_x_err) 1736 step_x = (max_x - min_x) / np 1737 x_to_x_err_arange = arange(min_x, max_x, step_x) 1738 y_to_x_err_arange = a * x_to_x_err_arange 1739 1740 ax.plot(x_to_x_err, x_to_x_err, 'o', label='%s vs. %s' % (method_x, method_x)) 1741 ax.plot(x_to_x_err_arange, x_to_x_err_arange, 'b--') 1742 ax.plot(x_to_x_err, y_to_y_err, '.', label='%s vs. %s' % (method_y, method_x) ) 1743 ax.plot(x_to_x_err_arange, y_to_x_err_arange, 'g--') 1744 1745 ax.set_title(r'$R_{2,\mathrm{eff}}/\sigma(R_{2,\mathrm{eff}})$' + ' for %s %i vs. %s %i. np=%i' % (method_y, glob_ini_y, method_x, glob_ini_x, np), fontsize=10) 1746 ax.legend(loc='upper left', shadow=True, prop = fontP) 1747 ax.set_xlabel(r'$R_{2,\mathrm{eff}}/\sigma(R_{2,\mathrm{eff}})$') 1748 ax.set_ylabel(r'$R_{2,\mathrm{eff}}/\sigma(R_{2,\mathrm{eff}})$') 1749 1750 # Add to data. 1751 for k, x_to_x_err_k in enumerate(x_to_x_err): 1752 y_to_y_err_k = y_to_y_err[k] 1753 x_to_x_err_arange_k = x_to_x_err_arange[k] 1754 y_to_x_err_arange_k = y_to_x_err_arange[k] 1755 data_dic[method_xy_NI].append(["%3.5f"%x_to_x_err_k, "%3.5f"%y_to_y_err_k, "%3.5f"%x_to_x_err_arange_k, "%3.5f"%y_to_x_err_arange_k]) 1756 1757 1758 plt.tight_layout() 1759 1760 # Loop over columns for writing data. 1761 # Write to file. 1762 if write_stats: 1763 # Re-order the data. 1764 headings_all = [] 1765 method_xy_NI_all = [] 1766 # Loop over the columns. 1767 for j in range(nr_cols): 1768 headings_j = [] 1769 method_xy_NI_j = [] 1770 # Loop over rows 1771 for i in range(nr_rows): 1772 # Extract from lists. 1773 data, methods, glob_inis = corr_data[j] 1774 method_x, method_y = methods 1775 glob_ini_x, glob_ini_y = glob_inis 1776 1777 # If row 1. 1778 if i == 0: 1779 # Add to headings. 1780 method_x_NI = "r2eff_%s%s" % (method_x, glob_ini_x) 1781 method_y_NI = "r2eff_%s%s" % (method_y, glob_ini_y) 1782 method_x_NI_lin = "r2eff_lin_%s%s" % (method_x, glob_ini_x) 1783 method_y_NI_lin = "r2eff_lin_%s%s" % (method_y, glob_ini_y) 1784 headings_j = headings_j + [method_x_NI, method_y_NI, method_x_NI_lin, method_y_NI_lin] 1785 1786 method_xy_NI = "r2eff_%s%s_%s%s" % (method_x, glob_ini_x, method_y, glob_ini_y) 1787 method_xy_NI_j.append(method_xy_NI) 1788 1789 # R2eff to error. 1790 if i == 1: 1791 # Add to headings 1792 method_x_NI = "r2eff_to_err_%s%s" % (method_x, glob_ini_x) 1793 method_y_NI = "r2eff_to_err_%s%s" % (method_y, glob_ini_y) 1794 method_x_NI_lin = "r2eff_to_err_lin_%s%s" % (method_x, glob_ini_x) 1795 method_y_NI_lin = "r2eff_to_err_lin_%s%s" % (method_y, glob_ini_y) 1796 headings_j = headings_j + [method_x_NI, method_y_NI, method_x_NI_lin, method_y_NI_lin] 1797 1798 method_xy_NI = "r2eff_to_err_%s%s_%s%s" % (method_x, glob_ini_x, method_y, glob_ini_y) 1799 method_xy_NI_j.append(method_xy_NI) 1800 1801 headings_all.append(headings_j) 1802 method_xy_NI_all.append(method_xy_NI_j) 1803 1804 # Loop over the columns. 1805 for j, headings_j in enumerate(headings_all): 1806 method_xy_NI_j = method_xy_NI_all[j] 1807 1808 data_w = [] 1809 data_r2eff = data_dic[method_xy_NI_j[0]] 1810 data_r2eff_to_err = data_dic[method_xy_NI_j[1]] 1811 1812 for k, data_r2eff_k in enumerate(data_r2eff): 1813 data_r2eff_to_err_k = data_r2eff_to_err[k] 1814 data_w.append(data_r2eff_k + data_r2eff_to_err_k) 1815 1816 # Define file name. 1817 data, methods, glob_inis = corr_data[j] 1818 data_x, data_y = data 1819 method_x, method_y = methods 1820 glob_ini_x, glob_ini_y = glob_inis 1821 np = len(data_r2eff) 1822 1823 # Get the spin selection for correlation. 1824 selection = data_x['selection'] 1825 1826 file_name_ini = 'r2eff_corr_%s_%s_%s_%s_NP_%i' % (method_x, glob_ini_x, method_y, glob_ini_y, np) 1827 if selection == None: 1828 file_name_ini = file_name_ini + '_all' 1829 else: 1830 file_name_ini = file_name_ini + '_sel' 1831 1832 file_name = file_name_ini + '.txt' 1833 path = self.results_dir 1834 1835 # save figure 1836 # Write png. 1837 png_file_name = file_name_ini + '.png' 1838 png_file_path = get_file_path(file_name=png_file_name, dir=path) 1839 plt.savefig(png_file_path, bbox_inches='tight') 1840 1841 # Write file 1842 file_obj, file_path = open_write_file(file_name=file_name, dir=path, force=True, compress_type=0, verbosity=1, return_path=True) 1843 1844 # Write data. 1845 write_data(out=file_obj, headings=headings_j, data=data_w) 1846 1847 # Close file. 1848 file_obj.close() 1849 1850 if show: 1851 plt.show()

1852 1853

1854 - def get_r2eff_stat_dic(self, list_r2eff_dics=None, list_glob_ini=None):

1855 1856 # Loop over the result dictionaries: 1857 res_dic = {} 1858 for i, r2eff_dic in enumerate(list_r2eff_dics): 1859 # Let the reference dic be initial dic 1860 r2eff_dic_ref = list_r2eff_dics[0] 1861 method_ref = r2eff_dic_ref['method'] 1862 res_dic['method_ref'] = method_ref 1863 glob_ini_ref = list_glob_ini[0] 1864 res_dic['glob_ini_ref'] = str(glob_ini_ref) 1865 selection = r2eff_dic_ref['selection'] 1866 res_dic['selection'] = selection 1867 1868 # Let the reference R2eff array be the initial glob. 1869 r2eff_arr_ref = r2eff_dic_ref[str(glob_ini_ref)]['r2eff_arr'] 1870 res_dic['r2eff_arr_ref'] = r2eff_arr_ref 1871 r2eff_err_arr_ref = r2eff_dic_ref[str(glob_ini_ref)]['r2eff_err_arr'] 1872 res_dic['r2eff_err_arr_ref'] = r2eff_err_arr_ref 1873 1874 # Get the current method 1875 method_cur = r2eff_dic['method'] 1876 res_dic[method_cur] = {} 1877 res_dic[method_cur]['method'] = method_cur 1878 res_dic[method_cur]['sampling_sparseness'] = [] 1879 res_dic[method_cur]['glob_ini'] = [] 1880 res_dic[method_cur]['r2eff_norm_std'] = [] 1881 1882 # Other stats. 1883 res_dic[method_cur]['pearsons_correlation_coefficient'] = [] 1884 res_dic[method_cur]['two_tailed_p_value'] = [] 1885 res_dic[method_cur]['r_xy'] = [] 1886 res_dic[method_cur]['a'] = [] 1887 res_dic[method_cur]['r_xy_2'] = [] 1888 res_dic[method_cur]['a_2'] = [] 1889 res_dic[method_cur]['r_xy_int'] = [] 1890 res_dic[method_cur]['a_int'] = [] 1891 res_dic[method_cur]['b_int'] = [] 1892 1893 # Now loop over glob_ini: 1894 for glob_ini in list_glob_ini: 1895 # Get the array, if it exists. 1896 if str(glob_ini) not in r2eff_dic: 1897 continue 1898 1899 # Get the data. 1900 r2eff_arr = r2eff_dic[str(glob_ini)]['r2eff_arr'] 1901 r2eff_err_arr = r2eff_dic[str(glob_ini)]['r2eff_err_arr'] 1902 1903 # Make a normalised R2eff array according to reference. 1904 # This require that all number of points are equal. 1905 # If they are not of same length, then dont even bother to continue. 1906 if len(r2eff_arr) != len(r2eff_arr_ref): 1907 continue 1908 1909 # Get the normalised array. 1910 r2eff_norm_arr = r2eff_arr/r2eff_arr_ref 1911 1912 # Calculate the standard deviation 1913 r2eff_norm_std = std(r2eff_norm_arr, ddof=1) 1914 1915 # Get the diff, then norm 1916 r2eff_diff_norm_arr = (r2eff_arr - r2eff_arr_ref) / r2eff_arr_ref 1917 r2eff_diff_norm_std = std(r2eff_diff_norm_arr, ddof=1) 1918 1919 # Store x 1920 sampling_sparseness = float(glob_ini) / float(glob_ini_ref) * 100. 1921 res_dic[method_cur]['sampling_sparseness'].append(sampling_sparseness) 1922 res_dic[method_cur]['glob_ini'].append(glob_ini) 1923 1924 # Store to result dic. 1925 res_dic[method_cur][str(glob_ini)] = {} 1926 res_dic[method_cur][str(glob_ini)]['sampling_sparseness'] = sampling_sparseness 1927 res_dic[method_cur][str(glob_ini)]['r2eff_arr'] = r2eff_arr 1928 res_dic[method_cur][str(glob_ini)]['r2eff_norm_arr'] = r2eff_norm_arr 1929 res_dic[method_cur][str(glob_ini)]['r2eff_norm_std'] = r2eff_norm_std 1930 res_dic[method_cur]['r2eff_norm_std'].append(r2eff_norm_std) 1931 1932 res_dic[method_cur][str(glob_ini)]['r2eff_diff_norm_arr'] = r2eff_diff_norm_arr 1933 res_dic[method_cur][str(glob_ini)]['r2eff_diff_norm_std'] = r2eff_diff_norm_std 1934 1935 ### Calculate for value over error. 1936 1937 # Calculate the R2eff versus R2eff error. 1938 r2eff_vs_err_ref = r2eff_arr_ref / r2eff_err_arr_ref 1939 r2eff_vs_err = r2eff_arr / r2eff_err_arr 1940 1941 # Get the statistics from scipy. 1942 pearsons_correlation_coefficient, two_tailed_p_value = pearsonr(r2eff_vs_err_ref, r2eff_vs_err) 1943 1944 # With intercept at axis. 1945 # Calculate sample correlation coefficient, measure of goodness-of-fit of linear regression 1946 x = r2eff_vs_err_ref 1947 x_m = mean(x) 1948 y = r2eff_vs_err 1949 y_m = mean(r2eff_vs_err) 1950 1951 # Solve by linear least squares. f(x) = a*x + b. 1952 n = len(y) 1953 a_int = (sum(x*y) - 1./n * sum(x) * sum(y) ) / ( sum(x**2) - 1./n * (sum(x))**2 ) 1954 b_int = 1./n * sum(y) - a_int * 1./n * sum(x) 1955 1956 r_xy_int = sum( (x - x_m)*(y - y_m) ) / sqrt( sum((x - x_m)**2) * sum((y - y_m)**2) ) 1957 1958 # Without intercept. 1959 a = sum(x*y) / sum(x**2) 1960 r_xy = sum(x*y) / sqrt(sum(x**2) * sum(y**2)) 1961 1962 # Without intercept for just R2eff. 1963 x_2 = r2eff_arr_ref 1964 y_2 = r2eff_arr 1965 a_2 = sum(x_2*y_2) / sum(x_2**2) 1966 r_xy_2 = sum(x_2*y_2) / sqrt(sum(x_2**2) * sum(y_2**2)) 1967 1968 #print(method_ref, method_cur, sampling_sparseness, glob_ini, pearsons_correlation_coefficient, r_xy**2, a, r_xy_int**2, a_int, b_int) 1969 print(method_ref, method_cur, sampling_sparseness, glob_ini, pearsons_correlation_coefficient, r_xy**2, a, r_xy_2**2, a_2) 1970 1971 # Store to result dic. 1972 res_dic[method_cur][str(glob_ini)]['pearsons_correlation_coefficient'] = pearsons_correlation_coefficient 1973 res_dic[method_cur]['pearsons_correlation_coefficient'].append(pearsons_correlation_coefficient) 1974 res_dic[method_cur][str(glob_ini)]['two_tailed_p_value'] = two_tailed_p_value 1975 res_dic[method_cur]['two_tailed_p_value'].append(two_tailed_p_value) 1976 res_dic[method_cur][str(glob_ini)]['r_xy'] = r_xy 1977 res_dic[method_cur]['r_xy'].append(r_xy) 1978 res_dic[method_cur][str(glob_ini)]['a'] = a 1979 res_dic[method_cur]['a'].append(a) 1980 1981 res_dic[method_cur][str(glob_ini)]['r_xy_2'] = r_xy_2 1982 res_dic[method_cur]['r_xy_2'].append(r_xy_2) 1983 res_dic[method_cur][str(glob_ini)]['a_2'] = a_2 1984 res_dic[method_cur]['a_2'].append(a_2) 1985 1986 res_dic[method_cur][str(glob_ini)]['r_xy_int'] = r_xy_int 1987 res_dic[method_cur]['r_xy_int'].append(r_xy_int) 1988 res_dic[method_cur][str(glob_ini)]['a_int'] = a_int 1989 res_dic[method_cur]['a_int'].append(a_int) 1990 res_dic[method_cur][str(glob_ini)]['b_int'] = b_int 1991 res_dic[method_cur]['b_int'].append(b_int) 1992 1993 res_dic[method_cur]['sampling_sparseness'] = asarray(res_dic[method_cur]['sampling_sparseness']) 1994 res_dic[method_cur]['glob_ini'] = asarray(res_dic[method_cur]['glob_ini']) 1995 res_dic[method_cur]['r2eff_norm_std'] = asarray(res_dic[method_cur]['r2eff_norm_std']) 1996 1997 res_dic[method_cur]['pearsons_correlation_coefficient'] = asarray(res_dic[method_cur]['pearsons_correlation_coefficient']) 1998 res_dic[method_cur]['two_tailed_p_value'] = asarray(res_dic[method_cur]['two_tailed_p_value']) 1999 res_dic[method_cur]['r_xy'] = asarray(res_dic[method_cur]['r_xy']) 2000 res_dic[method_cur]['a'] = asarray(res_dic[method_cur]['a']) 2001 res_dic[method_cur]['r_xy_2'] = asarray(res_dic[method_cur]['r_xy_2']) 2002 res_dic[method_cur]['a_2'] = asarray(res_dic[method_cur]['a_2']) 2003 res_dic[method_cur]['r_xy_int'] = asarray(res_dic[method_cur]['r_xy_int']) 2004 res_dic[method_cur]['a_int'] = asarray(res_dic[method_cur]['a_int']) 2005 res_dic[method_cur]['b_int'] = asarray(res_dic[method_cur]['b_int']) 2006 2007 return res_dic

2008 2009

2010 - def plot_r2eff_stat(self, r2eff_stat_dic=None, methods=[], list_glob_ini=[], show=False, write_stats=False):

2011 2012 # Define figure 2013 fig, axises = plt.subplots(nrows=2, ncols=1) 2014 fig.suptitle('Stats per NI') 2015 ax1, ax2 = axises 2016 2017 # Catch min and max values for all methods. 2018 min_a = 1.0 2019 max_a = 0.0 2020 2021 min_r_xy2 = 1.0 2022 max_r_xy2 = 0.0 2023 2024 # Prepare header for writing. 2025 selection = r2eff_stat_dic['selection'] 2026 2027 # For writing out stats. 2028 headings = [] 2029 data_dic = {} 2030 data_dic_methods = [] 2031 i_max = 0 2032 2033 for method in methods: 2034 if method not in r2eff_stat_dic: 2035 continue 2036 2037 # Use NI as x. 2038 NI = r2eff_stat_dic[method]['glob_ini'] 2039 # Use sampling_sparseness as x. 2040 SS = r2eff_stat_dic[method]['sampling_sparseness'] 2041 2042 # Add to headings. 2043 headings = headings + ['method', 'SS', 'NI', 'slope_r2eff', 'rxy2_r2eff', 'slope_r2eff_vs_err', 'rxy2_r2eff_vs_err'] 2044 2045 # Get stats. 2046 # Linear regression slope, without intercept 2047 a = r2eff_stat_dic[method]['a'] 2048 2049 if max(a) > max_a: 2050 max_a = max(a) 2051 if min(a) < min_a: 2052 min_a = min(a) 2053 2054 # sample correlation coefficient, without intercept 2055 r_xy = r2eff_stat_dic[method]['r_xy'] 2056 r_xy2 = r_xy**2 2057 2058 if max(r_xy2) > max_r_xy2: 2059 max_r_xy2 = max(r_xy2) 2060 if min(r_xy2) < min_r_xy2: 2061 min_r_xy2 = min(r_xy2) 2062 2063 # For just the R2eff values 2064 a_r2eff = r2eff_stat_dic[method]['a_2'] 2065 r_xy_r2eff = r2eff_stat_dic[method]['r_xy_2'] 2066 r_xy_r2eff2 = r_xy_r2eff**2 2067 2068 # Add to data. 2069 data_dic[method] = {} 2070 data_dic_methods.append(method) 2071 for i, NI_i in enumerate(NI): 2072 SS_i = SS[i] 2073 a_i = a[i] 2074 r_xy2_i = r_xy2[i] 2075 a_r2eff_i = a_r2eff[i] 2076 r_xy_r2eff2_i = r_xy_r2eff2[i] 2077 data_dic[method][str(i)] = ["%3.5f"%SS_i, "%i"%NI_i, "%3.5f"%a_r2eff_i, "%3.5f"%r_xy_r2eff2_i, "%3.5f"%a_i, "%3.5f"%r_xy2_i] 2078 if i > i_max: 2079 i_max = i 2080 2081 #ax1.plot(NI, a, ".-", label='%s LR'%method) 2082 #ax2.plot(NI, r_xy2, "o--", label='%s SC'%method) 2083 t = ax1.plot(SS, a_r2eff, ".--", label='%s slope R2eff'%method) 2084 color = t[0].get_color() 2085 ax1.plot(SS, a, ".-", label='%s slope'%method, color=color) 2086 2087 t = ax2.plot(SS, r_xy_r2eff2, "o--", label='%s r2 R2eff'%method) 2088 color = t[0].get_color() 2089 ax2.plot(SS, r_xy2, "o-", label='%s r2'%method, color=color) 2090 2091 # Loop over methods for writing data. 2092 data = [] 2093 2094 for i in range(0, i_max+1): 2095 data_i = [] 2096 for method in data_dic_methods: 2097 data_dic_m = data_dic[method] 2098 # Loop over all possible data points. 2099 if str(i) in data_dic_m: 2100 data_i = data_i + [method] + data_dic_m[str(i)] 2101 else: 2102 data_i = data_i + [method] + ["0", "0", "0", "0", "0", "0"] 2103 2104 data.append(data_i) 2105 2106 # Set legends. 2107 ax1.legend(loc='lower left', shadow=True, prop = fontP) 2108 #ax1.set_xlabel('NI') 2109 ax1.set_xlabel('SS') 2110 #ax1.set_ylabel(r'$\sigma ( R_{2,\mathrm{eff}} )$') 2111 ax1.set_ylabel('Linear regression slope, without intercept') 2112 #ax1.set_xticks(NI) 2113 #ax1.set_xticks(SS) 2114 ax1.set_ylim(min_a*0.95, max_a*1.05) 2115 ax1.invert_xaxis() 2116 2117 ax2.legend(loc='lower right', shadow=True, prop = fontP) 2118 ax2.set_ylabel('Sample correlation ' + r'$r_{xy}^2$') 2119 #ax2.set_xticks(NI) 2120 #ax2.set_xticks(SS) 2121 ax2.set_ylim(min_r_xy2*0.95, max_r_xy2*1.05) 2122 ax2.invert_xaxis() 2123 2124 # Determine filename. 2125 if selection == None: 2126 file_name_ini = 'r2eff_stat_all' 2127 else: 2128 file_name_ini = 'r2eff_stat_sel' 2129 2130 # Write png. 2131 png_file_name = file_name_ini + '.png' 2132 png_file_path = get_file_path(file_name=png_file_name, dir=self.results_dir) 2133 2134 # Write to file. 2135 if write_stats: 2136 # save figure 2137 plt.savefig(png_file_path, bbox_inches='tight') 2138 2139 file_name = file_name_ini + '.txt' 2140 path = self.results_dir 2141 file_obj, file_path = open_write_file(file_name=file_name, dir=path, force=True, compress_type=0, verbosity=1, return_path=True) 2142 2143 # Write data. 2144 write_data(out=file_obj, headings=headings, data=data) 2145 2146 # Close file. 2147 file_obj.close() 2148 2149 # Plot data. 2150 if show: 2151 plt.show()

2152 2153

2154 - def col_min(self, method=None, model=None, analysis=None, list_glob_ini=None, selection=None):

2155 2156 # Loop over the glob ini: 2157 res_dic = {} 2158 res_dic['method'] = method 2159 res_dic['selection'] = selection 2160 res_dic['analysis'] = analysis 2161 2162 for glob_ini in list_glob_ini: 2163 # Check previous, and get the pipe name. 2164 found, pipe_name, resfile, path = self.check_previous_result(method=method, model=model, analysis=analysis, glob_ini=glob_ini, bundle=method) 2165 2166 if pipes.cdp_name() != pipe_name: 2167 self.interpreter.pipe.switch(pipe_name) 2168 2169 # Results dictionary. 2170 res_dic[str(glob_ini)] = {} 2171 res_dic[str(glob_ini)]['params'] = {} 2172 2173 # Detect which params are in use. 2174 for cur_spin, mol_name, resi, resn, spin_id in spin_loop(selection=selection, full_info=True, return_id=True, skip_desel=True): 2175 params_list = cur_spin.params 2176 2177 # Store to dic. 2178 res_dic[str(glob_ini)]['params']['params_list'] = params_list 2179 2180 # Store individual. 2181 for param in params_list: 2182 # Store in list. 2183 res_dic[str(glob_ini)]['params'][param] = [] 2184 res_dic[str(glob_ini)]['params'][param+'_resi'] = [] 2185 2186 # Prepare to store individual per spin. 2187 res_dic[str(glob_ini)][param] = {} 2188 2189 # Break after first round. 2190 break 2191 2192 # Loop over the spins. 2193 for cur_spin, mol_name, resi, resn, spin_id in spin_loop(selection=selection, full_info=True, return_id=True, skip_desel=True): 2194 # Loop over params 2195 for param in params_list: 2196 2197 # Make spin dic. 2198 res_dic[str(glob_ini)][param][spin_id] = {} 2199 2200 # If param in PARAMS_R20, values are stored in with parameter key. 2201 param_key_list = [] 2202 if param in PARAMS_R20: 2203 # Loop over frq key. 2204 for exp_type, frq, offset, ei, mi, oi, in loop_exp_frq_offset(return_indices=True): 2205 # Get the parameter key. 2206 param_key = generate_r20_key(exp_type=exp_type, frq=frq) 2207 param_key_list.append(param_key) 2208 res_dic[str(glob_ini)][param][spin_id]['param_key_list'] = param_key_list 2209 2210 # Get the Value. 2211 param_val = deepcopy(getattr(cur_spin, param)[param_key]) 2212 # Store in list and per spin. 2213 res_dic[str(glob_ini)]['params'][param].append(param_val) 2214 res_dic[str(glob_ini)][param][spin_id][param_key] = param_val 2215 2216 res_dic[str(glob_ini)]['params'][param+'_resi'].append("%i_%3.0f"%(resi, frq/1E6)) 2217 2218 else: 2219 # Get the value. 2220 param_val = deepcopy(getattr(cur_spin, param)) 2221 # Store in list and per spin. 2222 res_dic[str(glob_ini)]['params'][param].append(param_val) 2223 res_dic[str(glob_ini)][param][spin_id][param_key] = param_val 2224 2225 res_dic[str(glob_ini)]['params'][param+'_resi'].append("%i"%(resi)) 2226 2227 # Print 2228 subtitle(file=sys.stdout, text="The minimised valus for '%s' model for pipe='%s at %s'" % (model, pipe_name, glob_ini), prespace=3) 2229 2230 # Convert to numpy array. 2231 for param in params_list: 2232 res_dic[str(glob_ini)]['params'][param] = asarray(res_dic[str(glob_ini)]['params'][param]) 2233 param_vals_str = " ".join(format(x, "2.3f") for x in res_dic[str(glob_ini)]['params'][param]) 2234 print(param, param_vals_str) 2235 2236 return res_dic

2237 2238

2239 - def plot_min_corr(self, corr_data, show=False, write_stats=False):

2240 # Define figure. 2241 # Nr of columns is number of datasets. 2242 nr_cols = len(corr_data) 2243 # Nr of rows, is the number of parameters. 2244 data_xy_0, methods_0, glob_inis_0 = corr_data[0] 2245 glob_ini_0 = glob_inis_0[0] 2246 params_list = data_xy_0[0][str(glob_ini_0)]['params']['params_list'] 2247 nr_rows = len(params_list) 2248 analysis = data_xy_0[0]['analysis'] 2249 2250 # Define figure 2251 fig, axises = plt.subplots(nrows=nr_rows, ncols=nr_cols) 2252 fig.suptitle('Correlation plot') 2253 2254 # axises is a tuple with number of elements corresponding to number of rows. 2255 # Each sub-tuple contains axis for each column. 2256 2257 # For writing out stats. 2258 data_dic = {} 2259 2260 # Loop over the rows. 2261 for i, row_axises in enumerate(axises): 2262 param = params_list[i] 2263 2264 # Loop over the columns. 2265 for j, ax in enumerate(row_axises): 2266 # Extract from lists. 2267 data, methods, glob_inis = corr_data[j] 2268 data_x, data_y = data 2269 method_x, method_y = methods 2270 glob_ini_x, glob_ini_y = glob_inis 2271 2272 x = data_x[str(glob_ini_x)]['params'][param] 2273 y = data_y[str(glob_ini_y)]['params'][param] 2274 # Relative uncertainty / fractional uncertainty / precision 2275 precision = abs(y-x) / ((x+y)/2) 2276 # Count outliers. Value differ more than than the value itself. 2277 precision_outlier = precision > 1.00 2278 #precision_outlier = precision > 0.020 2279 precision_outlier_nr = sum(precision_outlier) 2280 resis = data_x[str(glob_ini_x)]['params'][param+'_resi'] 2281 resi_outlier_arr = asarray(resis)[precision_outlier] 2282 resi_outlier_arr_str = ":"+','.join(str(e) for e in resi_outlier_arr) 2283 2284 np = len(y) 2285 2286 # Linear a, with no intercept. 2287 a = sum(x * y) / sum(x**2) 2288 min_xy = min(concatenate((x, y))) 2289 max_xy = max(concatenate((x, y))) 2290 2291 dx = (max_xy - min_xy) / np 2292 x_arange = arange(min_xy, max_xy + dx, dx) 2293 y_arange = a * x_arange 2294 2295 ax.plot(x, x, 'o', label='%s vs. %s' % (method_x, method_x)) 2296 ax.plot(x, y, '.', label='%s vs. %s' % (method_y, method_x) ) 2297 2298 #x_label = '%s'%param 2299 y_label = '%s'%param 2300 2301 #ax.set_xlabel(x_label) 2302 ax.set_ylabel(y_label) 2303 2304 ax.set_title('For %s %i vs. %s %i. np=%i' % (method_y, glob_ini_y, method_x, glob_ini_x, np), fontsize=10) 2305 ax.legend(loc='upper left', shadow=True, prop=fontP) 2306 2307 # kex has values in 1000 area. 2308 if param == 'kex': 2309 ax.ticklabel_format(style='sci', axis='x', scilimits=(0, 0)) 2310 ax.ticklabel_format(style='sci', axis='y', scilimits=(0, 0)) 2311 2312 ## If r2 or dw parameter, do a straight line: 2313 #if param in PARAMS_R20 + ['dw']: 2314 2315 # ax.plot(x_arange, x_arange, 'b--') 2316 # ax.plot(x_arange, y_arange, 'g--') 2317 2318 # Do a straight line for all. 2319 ax.plot(x_arange, x_arange, 'b--') 2320 ax.plot(x_arange, y_arange, 'g--') 2321 2322 # Store to data dic 2323 method_xy_NI = "%s_%s_%s%s_%s%s" % (analysis, param, method_x, glob_ini_x, method_y, glob_ini_y) 2324 data_dic[method_xy_NI] = [] 2325 2326 # Add to data. 2327 for k, x_k in enumerate(x): 2328 y_k = y[k] 2329 x_arange_k = x_arange[k] 2330 y_arange_k = y_arange[k] 2331 precision_k = precision[k] 2332 resi_k = resis[k] 2333 2334 data_dic[method_xy_NI].append(["%3.5f"%x_k, "%3.5f"%y_k, "%3.5f"%x_arange_k, "%3.5f"%y_arange_k, "%3.5f"%precision_k, "%s"%resi_k, "%i"%precision_outlier_nr, "%s"%resi_outlier_arr_str]) 2335 2336 plt.tight_layout() 2337 2338 # Loop over columns for writing data. 2339 # Write to file. 2340 if write_stats: 2341 # Re-order the data. 2342 headings_all = [] 2343 method_xy_NI_all = [] 2344 # Loop over the columns. 2345 for j in range(nr_cols): 2346 headings_j = [] 2347 method_xy_NI_j = [] 2348 # Loop over rows 2349 for i in range(nr_rows): 2350 # Extract from lists. 2351 data, methods, glob_inis = corr_data[j] 2352 method_x, method_y = methods 2353 glob_ini_x, glob_ini_y = glob_inis 2354 param = params_list[i] 2355 2356 # Add to headings 2357 method_x_NI = "%s_%s_%s%s" % (analysis, param, method_x, glob_ini_x) 2358 method_y_NI = "%s_%s_%s%s" % (analysis, param, method_y, glob_ini_y) 2359 method_x_NI_lin = "%s_%s_lin_%s%s" % (analysis, param, method_x, glob_ini_x) 2360 method_y_NI_lin = "%s_%s_lin_%s%s" % (analysis, param, method_y, glob_ini_y) 2361 headings_j = headings_j + [method_x_NI, method_y_NI, method_x_NI_lin, method_y_NI_lin, "abs_diff_frac", "resi", "outlier_nr", "outlier_resi"] 2362 2363 method_xy_NI = "%s_%s_%s%s_%s%s" % (analysis, param, method_x, glob_ini_x, method_y, glob_ini_y) 2364 method_xy_NI_j.append(method_xy_NI) 2365 2366 headings_all.append(headings_j) 2367 method_xy_NI_all.append(method_xy_NI_j) 2368 2369 # Loop over the columns. 2370 for j, headings_j in enumerate(headings_all): 2371 method_xy_NI_j = method_xy_NI_all[j] 2372 2373 data_w = [] 2374 method_xy_NI_r2 = method_xy_NI_j[0] 2375 data_r2 = data_dic[method_xy_NI_r2] 2376 2377 # Loop over the rows of data. 2378 for k, data_r2_k in enumerate(data_r2): 2379 data_row = data_r2_k 2380 # Loop over the columns. 2381 for method_xy_NI in method_xy_NI_j[1:]: 2382 data_param = data_dic[method_xy_NI] 2383 2384 try: 2385 data_param_row = data_param[k] 2386 except IndexError: 2387 data_param_row = len(data_param[0]) * ['0.0'] 2388 data_param_row[-1] = ':' 2389 2390 data_row = data_row + data_param_row 2391 2392 data_w.append(data_row) 2393 2394 # Define file name. 2395 data, methods, glob_inis = corr_data[j] 2396 data_x, data_y = data 2397 method_x, method_y = methods 2398 glob_ini_x, glob_ini_y = glob_inis 2399 2400 # Get the spin selection for correlation. 2401 selection = data_x['selection'] 2402 2403 file_name_ini = '%s_%s_%s_%s_%s' % (analysis, method_x, glob_ini_x, method_y, glob_ini_y) 2404 if selection == None: 2405 file_name_ini = file_name_ini + '_all' 2406 else: 2407 file_name_ini = file_name_ini + '_sel' 2408 2409 file_name = file_name_ini + '.txt' 2410 path = self.results_dir 2411 2412 # save figure 2413 # Write png. 2414 png_file_name = file_name_ini + '.png' 2415 png_file_path = get_file_path(file_name=png_file_name, dir=path) 2416 plt.savefig(png_file_path, bbox_inches='tight') 2417 2418 # Write file 2419 file_obj, file_path = open_write_file(file_name=file_name, dir=path, force=True, compress_type=0, verbosity=1, return_path=True) 2420 2421 # Write data. 2422 write_data(out=file_obj, headings=headings_j, data=data_w) 2423 2424 # Close file. 2425 file_obj.close() 2426 2427 if show: 2428 plt.show()

2429 2430

2431 - def get_min_stat_dic(self, list_r2eff_dics=None, list_glob_ini=None):

2432 2433 # Loop over the result dictionaries: 2434 res_dic = {} 2435 for i, min_dic in enumerate(list_r2eff_dics): 2436 # Let the reference dic be initial dic 2437 min_dic_ref = list_r2eff_dics[0] 2438 method_ref = min_dic_ref['method'] 2439 res_dic['method_ref'] = method_ref 2440 glob_ini_ref = list_glob_ini[0] 2441 res_dic['glob_ini_ref'] = str(glob_ini_ref) 2442 selection = min_dic_ref['selection'] 2443 res_dic['selection'] = selection 2444 params_list = min_dic_ref[str(glob_ini_ref)]['params']['params_list'] 2445 res_dic['params_list'] = params_list 2446 analysis = min_dic_ref['analysis'] 2447 res_dic['analysis'] = analysis 2448 2449 # Get the current method 2450 method_cur = min_dic['method'] 2451 res_dic[method_cur] = {} 2452 2453 # Loop over params 2454 for j, param in enumerate(params_list): 2455 res_dic[param] = {} 2456 2457 # Let the reference param array be the initial glob. 2458 param_arr_ref = min_dic_ref[str(glob_ini_ref)]['params'][param] 2459 res_dic[param]['param_arr_ref'] = param_arr_ref 2460 2461 res_dic[method_cur][param] = {} 2462 res_dic[method_cur][param]['method'] = method_cur 2463 res_dic[method_cur][param]['sampling_sparseness'] = [] 2464 res_dic[method_cur][param]['glob_ini'] = [] 2465 2466 # Other stats. 2467 res_dic[method_cur][param]['r_xy'] = [] 2468 res_dic[method_cur][param]['a'] = [] 2469 res_dic[method_cur][param]['precision_outlier_nr'] = [] 2470 res_dic[method_cur][param]['resi_outlier'] = [] 2471 2472 # Now loop over glob_ini: 2473 for glob_ini in list_glob_ini: 2474 # Get the array, if it exists. 2475 if str(glob_ini) not in min_dic: 2476 continue 2477 2478 # Get the data. 2479 param_arr = min_dic[str(glob_ini)]['params'][param] 2480 resis = min_dic[str(glob_ini)]['params'][param+'_resi'] 2481 2482 # This require that all number of points are equal. 2483 # If they are not of same length, then dont even bother to continue. 2484 if len(param_arr) != len(param_arr_ref): 2485 continue 2486 2487 # Store x 2488 sampling_sparseness = float(glob_ini) / float(glob_ini_ref) * 100. 2489 res_dic[method_cur][param]['sampling_sparseness'].append(sampling_sparseness) 2490 res_dic[method_cur][param]['glob_ini'].append(glob_ini) 2491 2492 # Store to result dic. 2493 res_dic[method_cur][param][str(glob_ini)] = {} 2494 res_dic[method_cur][param][str(glob_ini)]['sampling_sparseness'] = sampling_sparseness 2495 res_dic[method_cur][param][str(glob_ini)]['param_arr'] = param_arr 2496 2497 # With intercept at axis. 2498 # Calculate sample correlation coefficient, measure of goodness-of-fit of linear regression 2499 x = param_arr_ref 2500 x_m = mean(x) 2501 y = param_arr 2502 y_m = mean(y) 2503 2504 # Without intercept. 2505 a = sum(x*y) / sum(x**2) 2506 r_xy = sum(x*y) / sqrt(sum(x**2) * sum(y**2)) 2507 2508 # Relative uncertainty / fractional uncertainty / precision 2509 precision = abs(y-x) / ((x+y)/2) 2510 # Count outliers. Value differ more than than the value itself. 2511 precision_outlier = precision > 1.00 2512 #precision_outlier = precision > 0.02 2513 precision_outlier_nr = sum(precision_outlier) 2514 2515 resi_outlier_arr = asarray(resis)[precision_outlier] 2516 resi_outlier_arr_str = ":"+','.join(str(e) for e in resi_outlier_arr) 2517 2518 print(param, method_ref, method_cur, sampling_sparseness, glob_ini, r_xy**2, a, precision_outlier_nr, resi_outlier_arr_str) 2519 2520 # Store to result dic. 2521 res_dic[method_cur][param][str(glob_ini)]['r_xy'] = r_xy 2522 res_dic[method_cur][param]['r_xy'].append(r_xy) 2523 res_dic[method_cur][param][str(glob_ini)]['a'] = a 2524 res_dic[method_cur][param]['a'].append(a) 2525 res_dic[method_cur][param][str(glob_ini)]['precision_outlier_nr'] = precision_outlier_nr 2526 res_dic[method_cur][param]['precision_outlier_nr'].append(precision_outlier_nr) 2527 res_dic[method_cur][param]['resi_outlier'].append(resi_outlier_arr_str) 2528 2529 res_dic[method_cur][param]['sampling_sparseness'] = asarray(res_dic[method_cur][param]['sampling_sparseness']) 2530 res_dic[method_cur][param]['glob_ini'] = asarray(res_dic[method_cur][param]['glob_ini']) 2531 2532 res_dic[method_cur][param]['r_xy'] = asarray(res_dic[method_cur][param]['r_xy']) 2533 res_dic[method_cur][param]['a'] = asarray(res_dic[method_cur][param]['a']) 2534 res_dic[method_cur][param]['precision_outlier_nr'] = asarray(res_dic[method_cur][param]['precision_outlier_nr']) 2535 res_dic[method_cur][param]['resi_outlier'] = asarray(res_dic[method_cur][param]['resi_outlier']) 2536 2537 return res_dic

2538 2539

2540 - def plot_min_stat(self, min_stat_dic=None, methods=[], list_glob_ini=[], show=False, write_stats=False):

2541 2542 # Catch min and max values for all methods. 2543 min_a = 1.0 2544 max_a = 0.0 2545 2546 min_r_xy2 = 1.0 2547 max_r_xy2 = 0.0 2548 2549 # Prepare header for writing. 2550 selection = min_stat_dic['selection'] 2551 params_list = min_stat_dic['params_list'] 2552 analysis = min_stat_dic['analysis'] 2553 2554 # For writing out stats. 2555 headings = [] 2556 data_dic = {} 2557 data_dic_methods = [] 2558 2559 i_max = 0 2560 2561 for method in methods: 2562 if method not in min_stat_dic: 2563 continue 2564 2565 # Define figure 2566 fig, axises = plt.subplots(nrows=len(params_list), ncols=1) 2567 fig.suptitle('Stats per NI %s' % method) 2568 2569 # Loop over params 2570 data_dic[method] = {} 2571 data_dic_methods.append(method) 2572 2573 for j, param in enumerate(params_list): 2574 data_dic[method][param] = {} 2575 2576 # Use NI as x. 2577 NI = min_stat_dic[method][param]['glob_ini'] 2578 2579 # Use sampling_sparseness as x. 2580 SS = min_stat_dic[method][param]['sampling_sparseness'] 2581 2582 # Add to headings. 2583 headings = headings + ['method_%s'%param, 'SS', 'NI', 'slope', 'rxy2', 'outlier_nr', 'resi_outlier'] 2584 2585 # Get stats. 2586 # Linear regression slope, without intercept 2587 a = min_stat_dic[method][param]['a'] 2588 2589 if max(a) > max_a: 2590 max_a = max(a) 2591 if min(a) < min_a: 2592 min_a = min(a) 2593 2594 # sample correlation coefficient, without intercept 2595 r_xy = min_stat_dic[method][param]['r_xy'] 2596 r_xy2 = r_xy**2 2597 2598 if max(r_xy2) > max_r_xy2: 2599 max_r_xy2 = max(r_xy2) 2600 if min(r_xy2) < min_r_xy2: 2601 min_r_xy2 = min(r_xy2) 2602 2603 # Get the precision_outlier_nr, where values change more than the value itself. 2604 precision_outlier_nr = min_stat_dic[method][param]['precision_outlier_nr'] 2605 resi_outlier = min_stat_dic[method][param]['resi_outlier'] 2606 2607 # Add to data. 2608 for i, NI_i in enumerate(NI): 2609 SS_i = SS[i] 2610 a_i = a[i] 2611 r_xy2_i = r_xy2[i] 2612 precision_outlier_nr_i = precision_outlier_nr[i] 2613 resi_outlier_i = resi_outlier[i] 2614 data_dic[method][param][str(i)] = ["%3.5f"%SS_i, "%i"%NI_i, "%3.5f"%a_i, "%3.5f"%r_xy2_i, "%i"%precision_outlier_nr_i, "%s"%resi_outlier_i] 2615 if i > i_max: 2616 i_max = i 2617 2618 ax = axises[j] 2619 ax.plot(SS, a, ".-", label='%s_%s_a' % (method, param) ) 2620 ax.plot(SS, r_xy2, "o--", label='%s_%s_r_xy2' % (method, param) ) 2621 ax.legend(loc='lower left', shadow=True, prop = fontP) 2622 ax.set_xlabel('SS') 2623 ax.invert_xaxis() 2624 #ax.set_ylim(min_a*0.95, max_a*1.05) 2625 2626 2627 # Loop over methods for writing data. 2628 data = [] 2629 2630 # Loop over all lines. 2631 for i in range(0, i_max+1): 2632 data_i = [] 2633 for method in data_dic_methods: 2634 data_dic_m = data_dic[method] 2635 # Loop over all params 2636 for j, param in enumerate(params_list): 2637 # Loop over all possible data points. 2638 if str(i) in data_dic_m[param]: 2639 data_i = data_i + ["%s_%s" % (method, param)] + data_dic_m[param][str(i)] 2640 else: 2641 data_i = data_i + ["%s_%s" % (method, param)] + ["0", "0", "0", "0", "0", ":"] 2642 2643 data.append(data_i) 2644 2645 # Determine filename. 2646 if selection == None: 2647 file_name_ini = '%s_stat_all' % analysis 2648 else: 2649 file_name_ini = '%s_stat_sel' % analysis 2650 2651 # Write png. 2652 png_file_name = file_name_ini + '.png' 2653 png_file_path = get_file_path(file_name=png_file_name, dir=self.results_dir) 2654 2655 # Write to file. 2656 if write_stats: 2657 # save figure 2658 plt.savefig(png_file_path, bbox_inches='tight') 2659 2660 file_name = file_name_ini + '.txt' 2661 path = self.results_dir 2662 file_obj, file_path = open_write_file(file_name=file_name, dir=path, force=True, compress_type=0, verbosity=1, return_path=True) 2663 2664 # Write data. 2665 write_data(out=file_obj, headings=headings, data=data) 2666 2667 # Close file. 2668 file_obj.close() 2669 2670 # Plot data. 2671 if show: 2672 plt.show()

2673 2674

2675 - def write_results(self, method=None, model=None, analysis=None, list_glob_ini=None, selection=None, write_disp=True):

2676 2677 for glob_ini in list_glob_ini: 2678 # Check previous, and get the pipe name. 2679 found, pipe_name, resfile, path = self.check_previous_result(method=method, model=model, analysis=analysis, glob_ini=glob_ini, bundle=method) 2680 2681 if pipes.cdp_name() != pipe_name: 2682 self.interpreter.pipe.switch(pipe_name) 2683 2684 # Printout. 2685 section(file=sys.stdout, text="Results writing for pipe='%s"%(pipe_name), prespace=2, postspace=0) 2686 model_params = MODEL_PARAMS[model] 2687 subsection(file=sys.stdout, text="Model %s, with params='%s"%(model, model_params), prespace=0) 2688 2689 # Set path 2690 model_path = model.replace(" ", "_") 2691 analysis_path = analysis.replace(" ", "_") 2692 path = self.results_dir+sep+model_path+sep+analysis_path 2693 2694 # Dispersion curves. 2695 if write_disp: 2696 path_disp = path+sep+"disp_curves"+sep+method+sep+str(glob_ini) 2697 self.interpreter.relax_disp.plot_disp_curves(dir=path_disp, force=True) 2698 self.interpreter.relax_disp.write_disp_curves(dir=path_disp, force=True) 2699 2700 # The selected models for the final run. 2701 self.interpreter.value.write(param='model', file='model.out', dir=path, force=True) 2702 2703 models_tested = None 2704 2705 # For CPMG models. 2706 filep = str(glob_ini)+"_"+method+"_" 2707 path_par = path+sep+"r2" 2708 if has_cpmg_exp_type(): 2709 # The R20 parameter. 2710 self.write_results_test(path=path_par, model=model, models_tested=models_tested, param='r2', file_name_ini=filep+'r20') 2711 2712 # The R20A and R20B parameters. 2713 self.write_results_test(path=path_par, model=model, models_tested=models_tested, param='r2a', file_name_ini=filep+'r20a') 2714 self.write_results_test(path=path_par, model=model, models_tested=models_tested, param='r2b', file_name_ini=filep+'r20b') 2715 2716 # The pA and pB parameters. 2717 path_par = path+sep+"pop" 2718 search = method+"_"+"pA" 2719 self.write_results_test(path=path_par, model=model, models_tested=models_tested, search=search, param='pA', file_name_ini=filep+'pA') 2720 self.write_results_test(path=path_par, model=model, models_tested=models_tested, param='pB', file_name_ini=filep+'pB') 2721 2722 # The dw parameter. 2723 path_par = path+sep+"dw" 2724 search = method+"_"+"dw" 2725 self.write_results_test(path=path_par, model=model, models_tested=models_tested, search=search, param='dw', file_name_ini=filep+'dw') 2726 2727 # The k_AB, kex and tex parameters. 2728 path_par = path+sep+"rate" 2729 params = ['k_AB', 'kex', 'tex'] 2730 for param in params: 2731 search = method+"_"+param 2732 self.write_results_test(path=path_par, model=model, models_tested=models_tested, search=search, param=param, file_name_ini=filep+param) 2733 2734 # Minimisation statistics. 2735 if not (model == MODEL_R2EFF and has_fixed_time_exp_type()): 2736 path_par = path+sep+"chi2" 2737 self.interpreter.value.write(param='chi2', file=filep+'chi2.out', dir=path_par, force=True) 2738 search = method+"_"+"chi2" 2739 col_file_name="collect_%s.sh"%search 2740 self.write_convert_file(file_name=col_file_name, path=path_par, search=search) 2741 2742 self.interpreter.grace.write(y_data_type='chi2', file='chi2.agr', dir=path_par+sep+"grace", force=True)

2743 2744

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

2746 """Create a set of results, text and Grace files for the current data pipe. 2747 2748 @keyword path: The directory to place the files into. 2749 @type path: str 2750 @keyword model: The model tested. 2751 @type model: None or str 2752 @keyword model_tested: List of models tested, if the pipe is final. 2753 @type model_tested: None or list of str. 2754 @keyword param: The param to write out. 2755 @type param: None or list of str. 2756 @keyword file_name_ini: The initial part of the file name for the grace and text files. 2757 @type file_name_ini: None or str. 2758 """ 2759 2760 # If not set, use the name of the parameter. 2761 if file_name_ini == None: 2762 file_name_ini = param 2763 2764 # If the model is in the list of models which support the parameter. 2765 write_result = False 2766 if model != None: 2767 # Get the model params. 2768 model_params = MODEL_PARAMS[model] 2769 2770 if param in model_params: 2771 write_result = True 2772 2773 # If this is the final pipe, then check if the model has been tested at any time. 2774 elif model == None: 2775 # Loop through all tested models. 2776 for model_tested in models_tested: 2777 # If one of the models tested has a parameter which belong in the list of models which support the parameter, then write it out. 2778 model_params = MODEL_PARAMS[model_tested] 2779 2780 if param in model_params: 2781 write_result = True 2782 break 2783 2784 # Write results if some of the models supports the parameter. 2785 if write_result: 2786 self.interpreter.value.write(param=param, file='%s.out'%file_name_ini, dir=path, force=True) 2787 # Write convert file 2788 if search != None: 2789 col_file_name="collect_%s.sh"%search 2790 self.write_convert_file(file_name=col_file_name, path=path, search=search) 2791 2792 # Write grace 2793 self.interpreter.grace.write(x_data_type='res_num', y_data_type=param, file='%s.agr'%file_name_ini, dir=path+sep+"grace", force=True)

2794

2795 - def write_convert_file(self, file_name=None, path=None, search=None):

2796 file_obj, file_path = open_write_file(file_name=file_name, dir=path, force=True, compress_type=0, verbosity=1, return_path=True) 2797 2798 # Write file 2799 file_obj.write('#! /bin/bash' + '\n') 2800 file_obj.write('SEARCH=%s'%(search) + '\n') 2801 file_obj.write('FILES=(*_${SEARCH}.out)' + '\n') 2802 file_obj.write('readarray -t FILESSORT < <(for a in "${FILES[@]}"; do echo "$a"; done | sort -Vr)' + '\n') 2803 file_obj.write('# Skip the first two lines of header' + '\n') 2804 file_obj.write("tail -n+3 ${FILESSORT[0]} | sed 's,^# ,,' | grep -v 'None None' | awk '{print $2,$3,$5}' | column -t > collect_${SEARCH}.tmp" + '\n') 2805 file_obj.write('# Make array' + '\n') 2806 file_obj.write('ACUT=(collect_${SEARCH}.tmp)' + '\n') 2807 file_obj.write('for f in "${FILESSORT[@]}"; do' + '\n') 2808 file_obj.write(' FNAME="${f%.*}"' + '\n') 2809 file_obj.write(' NI=`echo $f | cut -d"_" -f1`' + '\n') 2810 file_obj.write(' echo "Processing $f with NI=$NI"' + '\n') 2811 file_obj.write(' tail -n+3 $f | sed "s,^# ,," | grep -v "None None" | sed "s,value,${NI}," | sed "s,error,${NI}," | awk %s{print $6,$7}%s | column -t > ${FNAME}.tmp'%("'","'") + '\n') 2812 file_obj.write(' ACUT+=(${FNAME}.tmp)' + '\n') 2813 file_obj.write('done' + '\n') 2814 file_obj.write('#paste "${ACUT[@]}" | column -t > collect_${SEARCH}.txt' + '\n') 2815 file_obj.write('paste "${ACUT[@]}" > collect_${SEARCH}.txt' + '\n') 2816 file_obj.write('rm ${ACUT[@]}' + '\n') 2817 2818 2819 # Close the batch script, then make it executable (expanding any ~ characters). 2820 file_obj.close() 2821 2822 chmod(file_path, S_IRWXU|S_IRGRP|S_IROTH)

2823 2824

2825 - def create_mc_data(self, number=500, distribution="measured", fixed_error=None, methods=None, model=None, model_from=None, analysis=None, analysis_from=None, list_glob_ini=None, force=False):

2826 """Create MC data.""" 2827 2828 # Set default 2829 if model_from == None: 2830 model_from = model 2831 if analysis_from == None: 2832 analysis_from = analysis 2833 2834 # Loop over the methods. 2835 for method in methods: 2836 # Change the self key. 2837 self.set_self(key='method', value=method) 2838 2839 # Loop over the glob ini: 2840 for glob_ini in list_glob_ini: 2841 # Check previous, and get the pipe name. 2842 found_pipe, pipe_name, resfile, path = self.check_previous_result(method=self.method, model=model, analysis=analysis, glob_ini=glob_ini, bundle=self.method) 2843 2844 # Try from analysis 2845 if not found_pipe: 2846 # Check previous, and get the pipe name. 2847 found_analysis, pipe_name, resfile, path = self.check_previous_result(method=self.method, model=model, analysis=analysis_from, glob_ini=glob_ini, bundle=self.method) 2848 2849 # Print 2850 subtitle(file=sys.stdout, text="MC data for pipe='%s'" % (pipe_name), prespace=3) 2851 2852 # Select the model. 2853 self.interpreter.relax_disp.select_model(model) 2854 2855 # Create data 2856 self.interpreter.monte_carlo.setup(number=number) 2857 self.interpreter.monte_carlo.create_data(distribution=distribution, fixed_error=fixed_error) 2858 self.interpreter.monte_carlo.initial_values() 2859 2860 # Save results, and store the current settings dic to pipe. 2861 cdp.settings = self.settings 2862 2863 # Define new pipe names. 2864 pipe_name = self.name_pipe(method=self.method, model=model, analysis=analysis, glob_ini=glob_ini) 2865 resfile = pipe_name.replace(" ", "_") 2866 model_path = model.replace(" ", "_") 2867 path = self.results_dir+sep+model_path 2868 2869 if found_pipe and not force: 2870 file_path = get_file_path(file_name=resfile, dir=path) 2871 text = "The file '%s' already exists. Set the force flag to True to overwrite." % (file_path) 2872 warn(RelaxWarning(text)) 2873 else: 2874 self.interpreter.results.write(file=resfile, dir=path, force=force)

2875 2876

2877 - def summarize_results(self, methods=None, model=None, analysis=None, list_glob_ini=None, selections=None):

2878 """ Collect all results and make statistics""" 2879 2880 ref_method, ni_method = methods 2881 ref_ni, ni_list = list_glob_ini 2882 2883 def collect_data(pipe_name=None, selections=None, ni=None): 2884 # Swith pipe 2885 self.interpreter.pipe.switch(pipe_name) 2886 2887 # Collect intensities 2888 all_intensities = [] 2889 all_intensities_err = [] 2890 2891 for selection in selections: 2892 sel_intensities = [] 2893 sel_intensities_err = [] 2894 for cur_spin, mol_name, resi, resn, spin_id in spin_loop(selection=selection, full_info=True, return_id=True, skip_desel=False): 2895 2896 # Loop over spectrum_ids. 2897 for s_id in cdp.spectrum_ids: 2898 # Check for bad data has skipped peak_intensity points 2899 if s_id in cur_spin.peak_intensity: 2900 peak_intensity_point = cur_spin.peak_intensity[s_id] 2901 peak_intensity_err_point = cur_spin.peak_intensity_err[s_id] 2902 2903 sel_intensities.append(peak_intensity_point) 2904 sel_intensities_err.append(peak_intensity_err_point) 2905 2906 # Collect 2907 all_intensities.append(sel_intensities) 2908 all_intensities_err.append(sel_intensities_err) 2909 2910 # Collect R2eff 2911 all_r2eff = [] 2912 all_r2eff_err = [] 2913 2914 for selection in selections: 2915 sel_r2eff = [] 2916 sel_r2eff_err = [] 2917 2918 # Loop over the spins. 2919 for cur_spin, mol_name, resi, resn, spin_id in spin_loop(selection=selection, full_info=True, return_id=True, skip_desel=False): 2920 # Loop over the R2eff points 2921 for exp_type, frq, offset, point, ei, mi, oi, di in loop_exp_frq_offset_point(return_indices=True): 2922 # Define the key. 2923 data_key = return_param_key_from_data(exp_type=self.exp_type, frq=frq, offset=offset, point=point) 2924 2925 # Check for bad data has skipped r2eff points 2926 if data_key in cur_spin.r2eff: 2927 r2eff_point = cur_spin.r2eff[data_key] 2928 r2eff_err_point = cur_spin.r2eff_err[data_key] 2929 2930 sel_r2eff.append(r2eff_point) 2931 sel_r2eff_err.append(r2eff_err_point) 2932 2933 # Collect 2934 all_r2eff.append(sel_r2eff) 2935 all_r2eff_err.append(sel_r2eff_err) 2936 2937 # Collect minimisations parameteers 2938 for cur_spin, mol_name, resi, resn, spin_id in spin_loop(selection=selection, full_info=True, return_id=True, skip_desel=True): 2939 params_list = cur_spin.params 2940 2941 # Break after first round. 2942 break 2943 2944 # Collect 2945 all_pars = [] 2946 all_pars_err = [] 2947 2948 # Loop over selections 2949 for selection in selections: 2950 sel_pars = [] 2951 sel_pars_err = [] 2952 2953 # Loop over params 2954 for param in params_list: 2955 # If param in PARAMS_R20, values are stored in with parameter key. 2956 sel_pars_list = [] 2957 sel_pars_list_err = [] 2958 2959 # Loop over the spins. 2960 for cur_spin, mol_name, resi, resn, spin_id in spin_loop(selection=selection, full_info=True, return_id=True, skip_desel=True): 2961 if param in PARAMS_R20: 2962 # Loop over frq key. 2963 for exp_type, frq, offset, ei, mi, oi, in loop_exp_frq_offset(return_indices=True): 2964 # Get the parameter key. 2965 param_key = generate_r20_key(exp_type=exp_type, frq=frq) 2966 2967 # Get the Value. 2968 param_val = deepcopy(getattr(cur_spin, param)[param_key]) 2969 param_err = deepcopy(getattr(cur_spin, param+"_err")[param_key]) 2970 2971 else: 2972 # Get the value. 2973 param_val = deepcopy(getattr(cur_spin, param)) 2974 param_err = deepcopy(getattr(cur_spin, param+"_err")) 2975 2976 # Collect 2977 sel_pars_list.append(param_val) 2978 sel_pars_list_err.append(param_err) 2979 2980 # Collect 2981 sel_pars.append(sel_pars_list) 2982 sel_pars_err.append(sel_pars_list_err) 2983 2984 # Collect 2985 all_pars.append(sel_pars) 2986 all_pars_err.append(sel_pars_err) 2987 2988 # Update data dictionary 2989 data_dic[pipe_name] = {} 2990 data_dic[pipe_name]['ni'] = ni 2991 data_dic[pipe_name]['all_intensities'] = all_intensities 2992 data_dic[pipe_name]['all_intensities_err'] = all_intensities_err 2993 data_dic[pipe_name]['all_r2eff'] = all_r2eff 2994 data_dic[pipe_name]['all_r2eff_err'] = all_r2eff_err 2995 data_dic['params_list'] = params_list 2996 data_dic[pipe_name]['all_pars'] = all_pars 2997 data_dic[pipe_name]['all_pars_err'] = all_pars_err

2998 2999 # Load the reference pipe 3000 found, ref_pipe_name, resfile, path = self.check_previous_result(method=ref_method, model=model, analysis=analysis, glob_ini=ref_ni, bundle=ref_method) 3001 3002 #Create data dic 3003 data_dic = {} 3004 # For ref 3005 collect_data(pipe_name=ref_pipe_name, selections=selections, ni=ref_ni) 3006 3007 # For other 3008 # Loop over ni and collect pipe names 3009 ni_pipe_names = [] 3010 for i, ni in enumerate(ni_list): 3011 # Get the ni pipe name 3012 found, ni_pipe_name, resfile, path = self.check_previous_result(method=ni_method, model=model, analysis=analysis, glob_ini=ni, bundle=ref_method) 3013 ni_pipe_names.append(ni_pipe_name) 3014 3015 # Collect data 3016 collect_data(pipe_name=ni_pipe_name, selections=selections, ni=ni) 3017 3018 # Swith pipe back to ref 3019 self.interpreter.pipe.switch(ref_pipe_name) 3020 3021 # Delete ni pipe, to free memory 3022 if ni_pipe_name != ref_pipe_name: 3023 self.interpreter.pipe.delete(pipe_name=ni_pipe_name) 3024 3025 # Now define method for stats 3026 def get_stats(x=None, y=None, x_err=None, y_err=None): 3027 # Convert to numpy 3028 x = asarray(x) 3029 y = asarray(y) 3030 3031 # Define the ratio weighted values 3032 g = x/x 3033 h = y/x 3034 3035 # Calculate the deviation 3036 d_xy = x - y 3037 d_gh = g - h 3038 3039 # Calculate the mean of the deviations 3040 mean_d_xy = mean(d_xy) 3041 mean_d_gh = mean(d_gh) 3042 3043 # Calculate the standard deviations 3044 std_d_xy = std(d_xy, ddof=1) 3045 std_d_gh = std(d_gh, ddof=1) 3046 3047 # Calculate the root mean square deviation 3048 rmsd_xy = sqrt(mean(square(d_xy))) 3049 rmsd_gh = sqrt(mean(square(d_gh))) 3050 3051 if x_err != None and y_err != None: 3052 pool_std_xy = sqrt( mean(square(x_err) + square(y_err)) ) 3053 3054 g_err = x_err / x 3055 h_err = y_err / x 3056 3057 pool_std_gh= sqrt( mean(square(g_err) + square(h_err)) ) 3058 3059 return mean_d_xy, mean_d_gh, std_d_xy, std_d_gh, pool_std_xy, pool_std_gh, rmsd_xy, rmsd_gh 3060 3061 else: 3062 return mean_d_xy, mean_d_gh, std_d_xy, std_d_gh, rmsd_xy, rmsd_gh

Source Code for Module auto_analyses.relax_disp_repeat_cpmg