1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24 """The Lipari-Szabo model-free analysis API object."""
25
26
27
28 import bmrblib
29 from copy import deepcopy
30 from math import pi
31 from minfx.grid import grid_split
32 from numpy import array, dot, float64, int32, zeros
33 from numpy.linalg import inv
34 from re import match, search
35 import string
36 from types import MethodType
37 from warnings import warn
38
39
40 from lib.arg_check import is_num_list, is_str_list
41 from lib.errors import RelaxError, RelaxFault, RelaxNoModelError, RelaxNoSequenceError, RelaxNoTensorError
42 from lib.float import isInf
43 from lib.physical_constants import h_bar, mu0, return_gyromagnetic_ratio
44 from lib.warnings import RelaxDeselectWarning, RelaxWarning
45 from multi import Processor_box
46 from pipe_control import diffusion_tensor, interatomic, mol_res_spin, pipes, relax_data, sequence
47 from pipe_control.bmrb import list_sample_conditions
48 from pipe_control.exp_info import bmrb_write_citations, bmrb_write_methods, bmrb_write_software
49 from pipe_control.interatomic import return_interatom_list
50 from pipe_control.mol_res_spin import count_spins, exists_mol_res_spin_data, find_index, get_molecule_names, return_spin, return_spin_from_index, return_spin_indices, spin_loop
51 from specific_analyses.api_base import API_base
52 from specific_analyses.api_common import API_common
53 from specific_analyses.model_free.bmrb import sf_csa_read, sf_model_free_read, to_bmrb_model
54 from specific_analyses.model_free.data import compare_objects
55 from specific_analyses.model_free.molmol import Molmol
56 from specific_analyses.model_free.model import determine_model_type
57 from specific_analyses.model_free.parameters import are_mf_params_set, assemble_param_names, assemble_param_vector, assemble_scaling_matrix, linear_constraints
58 from specific_analyses.model_free.optimisation import MF_grid_command, MF_memo, MF_minimise_command, grid_search_config, minimise_data_setup, relax_data_opt_structs, reset_min_stats
59 from specific_analyses.model_free.parameter_object import Model_free_params
60 from specific_analyses.model_free.pymol import Pymol
61 from target_functions.mf import Mf
62
63
65 """Parent class containing all the model-free specific functions."""
66
67
68 instance = None
69
89
90
91 - def back_calc_ri(self, spin_index=None, ri_id=None, ri_type=None, frq=None):
92 """Back-calculation of relaxation data from the model-free parameter values.
93
94 @keyword spin_index: The global spin index.
95 @type spin_index: int
96 @keyword ri_id: The relaxation data ID string.
97 @type ri_id: str
98 @keyword ri_type: The relaxation data type.
99 @type ri_type: str
100 @keyword frq: The field strength.
101 @type frq: float
102 @return: The back calculated relaxation data value corresponding to the index.
103 @rtype: float
104 """
105
106
107 spin, spin_id = return_spin_from_index(global_index=spin_index, return_spin_id=True)
108
109
110 if hasattr(cdp, 'diff_tensor') and (cdp.diff_tensor.type == 'spheroid' or cdp.diff_tensor.type == 'ellipsoid'):
111 interatoms = interatomic.return_interatom_list(spin_id)
112 for i in range(len(interatoms)):
113
114 if not interatoms[i].dipole_pair:
115 continue
116
117
118 if not hasattr(interatoms[i], 'vector') or interatoms[i].vector == None:
119 warn(RelaxDeselectWarning(spin_id, 'missing structural data'))
120 return
121
122
123 self.overfit_deselect(data_check=False, verbose=False)
124
125
126 value = self.minimise(min_algor='back_calc', min_options=(spin_index, ri_id, ri_type, frq))
127
128
129 return value
130
131
132 - def bmrb_read(self, file_path, version=None, sample_conditions=None):
170
171
173 """Write the model-free results to a BMRB NMR-STAR v3.1 formatted file.
174
175 @param file_path: The full file path.
176 @type file_path: str
177 @keyword version: The BMRB NMR-STAR dictionary format to output to.
178 @type version: str
179 """
180
181
182 cdp = pipes.get_pipe()
183
184
185 star = bmrblib.create_nmr_star('relax_model_free_results', file_path, version)
186
187
188 global_chi2 = None
189 if hasattr(cdp, 'chi2'):
190 global_chi2 = cdp.chi2
191
192
193 rex_frq = cdp.spectrometer_frq[cdp.ri_ids[0]]
194
195
196 mol_name_list = []
197 res_num_list = []
198 res_name_list = []
199 atom_name_list = []
200
201 csa_list = []
202 r_list = []
203 isotope_list = []
204 element_list = []
205
206 local_tm_list = []
207 s2_list = []
208 s2f_list = []
209 s2s_list = []
210 te_list = []
211 tf_list = []
212 ts_list = []
213 rex_list = []
214
215 local_tm_err_list = []
216 s2_err_list = []
217 s2f_err_list = []
218 s2s_err_list = []
219 te_err_list = []
220 tf_err_list = []
221 ts_err_list = []
222 rex_err_list = []
223
224 chi2_list = []
225 model_list = []
226
227
228 for spin, mol_name, res_num, res_name, spin_id in spin_loop(full_info=True, return_id=True):
229
230 if spin.name == 'H' or (hasattr(spin, 'element') and spin.element == 'H'):
231 warn(RelaxWarning("Skipping the proton spin '%s'." % spin_id))
232 continue
233
234
235 if res_num == None:
236 raise RelaxError("For the BMRB, the residue of spin '%s' must be numbered." % spin_id)
237 if res_name == None:
238 raise RelaxError("For the BMRB, the residue of spin '%s' must be named." % spin_id)
239 if spin.name == None:
240 raise RelaxError("For the BMRB, the spin '%s' must be named." % spin_id)
241 if not hasattr(spin, 'isotope') or spin.isotope == None:
242 raise RelaxError("For the BMRB, the spin isotope type of '%s' must be specified." % spin_id)
243 if not hasattr(spin, 'element') or spin.element == None:
244 raise RelaxError("For the BMRB, the spin element type of '%s' must be specified. Please use the spin user function for setting the element type." % spin_id)
245
246
247 mol_name_list.append(mol_name)
248 res_num_list.append(res_num)
249 res_name_list.append(res_name)
250 atom_name_list.append(spin.name)
251
252
253 if hasattr(spin, 'csa'):
254 csa_list.append(spin.csa * 1e6)
255 else:
256 csa_list.append(None)
257
258
259 interatoms = return_interatom_list(spin_id)
260 for i in range(len(interatoms)):
261
262 if not interatoms[i].dipole_pair:
263 continue
264
265
266 if hasattr(interatoms[i], 'r'):
267 r_list.append(interatoms[i].r)
268 else:
269 r_list.append(None)
270
271
272 break
273
274
275 if hasattr(spin, 'isotope'):
276 isotope_list.append(int(spin.isotope.strip(string.ascii_letters)))
277 else:
278 isotope_list.append(None)
279
280
281 if hasattr(spin, 'element'):
282 element_list.append(spin.element)
283 else:
284 element_list.append(None)
285
286
287 if hasattr(spin, 'local_tm'):
288 local_tm_list.append(spin.local_tm)
289 else:
290 local_tm_list.append(None)
291 if hasattr(spin, 'local_tm_err'):
292 local_tm_err_list.append(spin.local_tm_err)
293 else:
294 local_tm_err_list.append(None)
295
296
297 s2_list.append(spin.s2)
298 s2f_list.append(spin.s2f)
299 s2s_list.append(spin.s2s)
300 te_list.append(spin.te)
301 tf_list.append(spin.tf)
302 ts_list.append(spin.ts)
303 if spin.rex == None:
304 rex_list.append(None)
305 else:
306 rex_list.append(spin.rex / (2.0*pi*rex_frq**2))
307
308
309 params = ['s2', 's2f', 's2s', 'te', 'tf', 'ts', 'rex']
310 for param in params:
311
312 err_list = locals()[param+'_err_list']
313
314
315 if hasattr(spin, param+'_err'):
316
317 val = getattr(spin, param+'_err')
318
319
320 if param == 'rex' and val != None:
321 val = val / (2.0*pi*rex_frq**2)
322
323
324 err_list.append(val)
325
326
327 else:
328 err_list.append(None)
329
330
331
332 chi2_list.append(spin.chi2)
333
334
335 model_list.append(to_bmrb_model(spin.model))
336
337
338 entity_ids = zeros(len(mol_name_list), int32)
339 mol_names = get_molecule_names()
340 for i in range(len(mol_name_list)):
341 for j in range(len(mol_names)):
342 if mol_name_list[i] == mol_names[j]:
343 entity_ids[i] = j+1
344
345
346
347
348
349
350 bmrb_write_citations(star)
351
352
353
354
355
356
357 mol_res_spin.bmrb_write_entity(star)
358
359
360
361
362
363
364 bmrb_write_methods(star)
365
366
367 software_ids, software_labels = bmrb_write_software(star)
368
369
370
371
372
373
374 star.chem_shift_anisotropy.add(entity_ids=entity_ids, res_nums=res_num_list, res_names=res_name_list, atom_names=atom_name_list, atom_types=element_list, isotope=isotope_list, csa=csa_list)
375
376
377
378
379
380
381 relax_data.bmrb_write(star)
382
383
384
385
386
387
388 for i in range(len(software_ids)):
389 if software_labels[i] == 'relax':
390 software_id = software_ids[i]
391
392
393 star.model_free.add(global_chi2=global_chi2, software_ids=[software_id], software_labels=['relax'], entity_ids=entity_ids, res_nums=res_num_list, res_names=res_name_list, atom_names=atom_name_list, atom_types=element_list, isotope=isotope_list, bond_length=r_list, local_tc=local_tm_list, s2=s2_list, s2f=s2f_list, s2s=s2s_list, te=te_list, tf=tf_list, ts=ts_list, rex=rex_list, local_tc_err=local_tm_err_list, s2_err=s2_err_list, s2f_err=s2f_err_list, s2s_err=s2s_err_list, te_err=te_err_list, tf_err=tf_err_list, ts_err=ts_err_list, rex_err=rex_err_list, rex_frq=rex_frq, chi2=chi2_list, model_fit=model_list)
394
395
396
397
398
399
400 diffusion_tensor.bmrb_write(star)
401
402
403
404 star.write()
405
406
407 - def calculate(self, spin_id=None, verbosity=1, sim_index=None):
408 """Calculation of the model-free chi-squared value.
409
410 @keyword spin_id: The spin identification string.
411 @type spin_id: str
412 @keyword verbosity: The amount of information to print. The higher the value, the greater the verbosity.
413 @type verbosity: int
414 @keyword sim_index: The optional MC simulation index.
415 @type sim_index: int
416 """
417
418
419 if not exists_mol_res_spin_data():
420 raise RelaxNoSequenceError
421
422
423 model_type = determine_model_type()
424
425
426 if model_type != 'local_tm' and not diffusion_tensor.diff_data_exists():
427 raise RelaxNoTensorError('diffusion')
428
429
430 if model_type != 'local_tm' and cdp.diff_tensor.type != 'sphere' and not hasattr(cdp, 'structure'):
431 raise RelaxNoPdbError
432
433
434 for spin, id in spin_loop(spin_id, return_id=True):
435
436 if not spin.select:
437 continue
438
439
440 if not spin.model:
441 raise RelaxNoModelError
442
443
444 if not hasattr(spin, 'isotope'):
445 raise RelaxSpinTypeError
446
447
448 unset_param = are_mf_params_set(spin)
449 if unset_param != None:
450 raise RelaxNoValueError(unset_param)
451
452
453 if not hasattr(spin, 'csa') or spin.csa == None:
454 raise RelaxNoValueError("CSA")
455
456
457 interatoms = return_interatom_list(spin_id)
458 for interatom in interatoms:
459
460 if not interatom.dipole_pair:
461 continue
462
463
464 if model_type != 'local_tm' and cdp.diff_tensor.type != 'sphere' and not hasattr(interatom, 'vector'):
465 raise RelaxNoVectorsError
466
467
468 if model_type != 'local_tm' and cdp.diff_tensor.type != 'sphere' and hasattr(interatom, 'vector') and is_num_list(interatom.vector[0], raise_error=False):
469 raise RelaxMultiVectorError
470
471
472 if spin_id != interatom.spin_id1:
473 spin_id2 = interatom.spin_id1
474 else:
475 spin_id2 = interatom.spin_id2
476 spin2 = return_spin(spin_id2)
477
478
479 if not hasattr(spin2, 'isotope'):
480 raise RelaxSpinTypeError
481
482
483 if not hasattr(interatom, 'r') or interatom.r == None:
484 raise RelaxNoValueError("interatomic distance", spin_id=spin_id, spin_id2=spin_id2)
485
486
487 if not hasattr(spin, 'ri_data') or not hasattr(spin, 'ri_data_err'):
488 continue
489
490
491 for ri_id in cdp.ri_ids:
492
493 err = spin.ri_data_err[ri_id]
494
495
496 if err != None and err == 0.0:
497 raise RelaxError("Zero error for spin '%s' for the relaxation data ID '%s', minimisation not possible." % (id, ri_id))
498 elif err != None and err < 0.0:
499 raise RelaxError("Negative error of %s for spin '%s' for the relaxation data ID '%s', minimisation not possible." % (err, id, ri_id))
500
501
502 param_vector = assemble_param_vector(spin=spin, sim_index=sim_index)
503
504
505 data = relax_data_opt_structs(spin, sim_index=sim_index)
506
507
508 ri_data = [array(data[0])]
509 ri_data_err = [array(data[1])]
510 num_frq = [data[2]]
511 num_ri = [data[3]]
512 ri_labels = [data[4]]
513 frq = [data[5]]
514 remap_table = [data[6]]
515 noe_r1_table = [data[7]]
516 gx = [return_gyromagnetic_ratio(spin.isotope)]
517 if sim_index == None:
518 csa = [spin.csa]
519 else:
520 csa = [spin.csa_sim[sim_index]]
521
522
523 interatoms = return_interatom_list(spin_id)
524 for i in range(len(interatoms)):
525
526 if not interatoms[i].dipole_pair:
527 continue
528
529
530 if spin_id != interatoms[i].spin_id1:
531 spin_id2 = interatoms[i].spin_id1
532 else:
533 spin_id2 = interatoms[i].spin_id2
534 spin2 = return_spin(spin_id2)
535
536
537 if sim_index == None:
538 r = [interatoms[i].r]
539 else:
540 r = [interatoms[i].r_sim[sim_index]]
541
542
543 if model_type != 'local_tm' and cdp.diff_tensor.type != 'sphere':
544 xh_unit_vectors = [interatoms[i].vector]
545 else:
546 xh_unit_vectors = [None]
547
548
549 gh = [return_gyromagnetic_ratio(spin2.isotope)]
550
551
552 num_params = [len(spin.params)]
553
554
555 param_values = [assemble_param_vector(model_type='mf')]
556
557
558 if model_type == 'local_tm':
559 diff_params = [spin.local_tm]
560 diff_type = 'sphere'
561 else:
562
563 diff_type = cdp.diff_tensor.type
564
565
566 if diff_type == 'sphere':
567 diff_params = [cdp.diff_tensor.tm]
568
569
570 elif diff_type == 'spheroid':
571 diff_params = [cdp.diff_tensor.tm, cdp.diff_tensor.Da, cdp.diff_tensor.theta, cdp.diff_tensor.phi]
572
573
574 elif diff_type == 'ellipsoid':
575 diff_params = [cdp.diff_tensor.tm, cdp.diff_tensor.Da, cdp.diff_tensor.Dr, cdp.diff_tensor.alpha, cdp.diff_tensor.beta, cdp.diff_tensor.gamma]
576
577
578 mf = Mf(init_params=param_vector, model_type='mf', diff_type=diff_type, diff_params=diff_params, num_spins=1, equations=[spin.equation], param_types=[spin.params], param_values=param_values, relax_data=ri_data, errors=ri_data_err, bond_length=r, csa=csa, num_frq=num_frq, frq=frq, num_ri=num_ri, remap_table=remap_table, noe_r1_table=noe_r1_table, ri_labels=ri_labels, gx=gx, gh=gh, h_bar=h_bar, mu0=mu0, num_params=num_params, vectors=xh_unit_vectors)
579
580
581 try:
582 chi2 = mf.func(param_vector)
583 except OverflowError:
584 chi2 = 1e200
585
586
587 if model_type == 'all' or model_type == 'diff':
588 cdp.chi2 = cdp.chi2 + chi2
589 else:
590 spin.chi2 = chi2
591
592
594 """Create the Monte Carlo Ri data.
595
596 @keyword data_id: The spin identification string, as yielded by the base_data_loop() generator method.
597 @type data_id: str
598 @return: The Monte Carlo simulation data.
599 @rtype: list of floats
600 """
601
602
603 mc_data = []
604
605
606 spin = return_spin(data_id)
607 global_index = find_index(data_id)
608
609
610 if not spin.select:
611 return
612
613
614 if not hasattr(spin, 'model') or not spin.model:
615 return None
616
617
618 for ri_id in cdp.ri_ids:
619
620 value = self.back_calc_ri(spin_index=global_index, ri_id=ri_id, ri_type=cdp.ri_type[ri_id], frq=cdp.spectrometer_frq[ri_id])
621
622
623 mc_data.append(value)
624
625
626 return mc_data
627
628
630 """Initialise the spin specific data structures.
631
632 @param data_cont: The spin data container.
633 @type data_cont: SpinContainer instance
634 @keyword sim: The Monte Carlo simulation flag, which if true will initialise the simulation data structure.
635 @type sim: bool
636 """
637
638
639 for name in self._PARAMS.loop(scope='spin'):
640
641 if name in ['ri_data', 'ri_data_bc', 'ri_data_err']:
642 continue
643
644
645 list_data = [ 'params' ]
646 if name in list_data:
647 init_data = []
648
649
650 init_data = None
651 if self._PARAMS.type(name) == bool:
652 init_data = False
653 if name == 'select':
654 init_data = True
655
656
657 if not hasattr(data_cont, name):
658 setattr(data_cont, name, init_data)
659
660
661 - def deselect(self, model_info, sim_index=None):
662 """Deselect models or simulations.
663
664 @param model_info: The model index from model_loop(). This is zero for the global models or equal to the global spin index (which covers the molecule, residue, and spin indices).
665 @type model_info: int
666 @keyword sim_index: The optional Monte Carlo simulation index. If None, then models will be deselected, otherwise the given simulation will.
667 @type sim_index: None or int
668 """
669
670
671 model_type = determine_model_type()
672
673
674 if model_type == 'mf' or model_type == 'local_tm':
675
676 spin = return_spin_from_index(model_info)
677
678
679 if sim_index == None:
680 spin.select = False
681
682
683 else:
684 spin.select_sim[sim_index] = False
685
686
687 else:
688
689 if sim_index == None:
690 raise RelaxError("Cannot deselect the global model.")
691
692
693 else:
694
695 cdp.select_sim[sim_index] = False
696
697
698 - def duplicate_data(self, pipe_from=None, pipe_to=None, model_info=None, global_stats=False, verbose=True):
699 """Duplicate the data specific to a single model-free model.
700
701 @keyword pipe_from: The data pipe to copy the data from.
702 @type pipe_from: str
703 @keyword pipe_to: The data pipe to copy the data to.
704 @type pipe_to: str
705 @param model_info: The model index from model_loop(). This is zero for the global models or equal to the global spin index (which covers the molecule, residue, and spin indices).
706 @type model_info: int
707 @keyword global_stats: The global statistics flag.
708 @type global_stats: bool
709 @keyword verbose: A flag which if True will cause info about each spin to be printed out as the sequence is generated.
710 @type verbose: bool
711 """
712
713
714 if model_info == None:
715 raise RelaxError("The model_info argument cannot be None.")
716
717
718 if not pipes.has_pipe(pipe_to):
719 pipes.create(pipe_to, pipe_type='mf', switch=False)
720
721
722 dp_from = pipes.get_pipe(pipe_from)
723 dp_to = pipes.get_pipe(pipe_to)
724
725
726 for data_name in dir(dp_from):
727
728 if data_name in ['diff_tensor', 'mol', 'interatomic', 'structure', 'exp_info']:
729 continue
730
731
732 if search('^_', data_name) or data_name in list(dp_from.__class__.__dict__.keys()):
733 continue
734
735
736 data_from = getattr(dp_from, data_name)
737
738
739 if hasattr(dp_to, data_name):
740
741 data_to = getattr(dp_to, data_name)
742
743
744 if data_from != data_to:
745 raise RelaxError("The object " + repr(data_name) + " is not consistent between the pipes " + repr(pipe_from) + " and " + repr(pipe_to) + ".")
746
747
748 continue
749
750
751 setattr(dp_to, data_name, deepcopy(data_from))
752
753
754 if hasattr(dp_from, 'diff_tensor'):
755
756 if not hasattr(dp_to, 'diff_tensor'):
757 setattr(dp_to, 'diff_tensor', deepcopy(dp_from.diff_tensor))
758
759
760 else:
761
762 for data_name in dp_from.diff_tensor._mod_attr:
763
764 data_from = None
765 if hasattr(dp_from.diff_tensor, data_name):
766 data_from = getattr(dp_from.diff_tensor, data_name)
767
768
769 if data_from and not hasattr(dp_to.diff_tensor, data_name):
770 raise RelaxError("The diffusion tensor object " + repr(data_name) + " of the " + repr(pipe_from) + " data pipe is not located in the " + repr(pipe_to) + " data pipe.")
771 elif data_from:
772 data_to = getattr(dp_to.diff_tensor, data_name)
773 else:
774 continue
775
776
777 if data_from != data_to:
778 raise RelaxError("The object " + repr(data_name) + " is not consistent between the pipes " + repr(pipe_from) + " and " + repr(pipe_to) + ".")
779
780
781 if hasattr(dp_from, 'structure'):
782
783 if not hasattr(dp_to, 'structure'):
784 setattr(dp_to, 'structure', deepcopy(dp_from.structure))
785
786
787 else:
788
789 compare_objects(dp_from.structure, dp_to.structure, pipe_from, pipe_to)
790
791
792 if len(dp_from.structure.structural_data) != len(dp_from.structure.structural_data):
793 raise RelaxError("The number of structural models is not consistent between the pipes " + repr(pipe_from) + " and " + repr(pipe_to) + ".")
794
795
796 for i in range(len(dp_from.structure.structural_data)):
797
798 model_from = dp_from.structure.structural_data[i]
799 model_to = dp_to.structure.structural_data[i]
800
801
802 if model_from.num != model_to.num:
803 raise RelaxError("The structure models are not consistent between the pipes " + repr(pipe_from) + " and " + repr(pipe_to) + ".")
804
805
806 if len(model_from.mol) != len(model_to.mol):
807 raise RelaxError("The number of molecules is not consistent between the pipes " + repr(pipe_from) + " and " + repr(pipe_to) + ".")
808
809
810 for mol_index in range(len(model_from.mol)):
811
812 compare_objects(model_from.mol[mol_index], model_to.mol[mol_index], pipe_from, pipe_to)
813
814
815 if dp_from.mol.is_empty():
816 return
817
818
819 if dp_to.mol.is_empty():
820 sequence.copy(pipe_from=pipe_from, pipe_to=pipe_to, preserve_select=True, verbose=verbose)
821
822
823 if dp_to.interatomic.is_empty():
824 interatomic.copy(pipe_from=pipe_from, pipe_to=pipe_to, verbose=verbose)
825
826
827 pipes.switch(pipe_from)
828 model_type = determine_model_type()
829
830
831 spin, spin_id = return_spin_from_index(model_info, pipe=pipe_from, return_spin_id=True)
832 if model_type == 'mf' or (model_type == 'local_tm' and not global_stats):
833
834 for name in dir(spin):
835
836 if search('^__', name):
837 continue
838
839
840 obj = getattr(spin, name)
841
842
843 if isinstance(obj, MethodType):
844 continue
845
846
847 new_obj = deepcopy(getattr(spin, name))
848 setattr(dp_to.mol[spin._mol_index].res[spin._res_index].spin[spin._spin_index], name, new_obj)
849
850
851 interatoms = interatomic.return_interatom_list(spin_id)
852 for interatom in interatoms:
853
854 if not interatom.dipole_pair:
855 continue
856
857
858 if spin_id != interatom.spin_id1:
859 spin_id2 = interatom.spin_id1
860 else:
861 spin_id2 = interatom.spin_id2
862 spin2 = return_spin(spin_id2)
863
864
865 mol_index, res_index, spin_index = return_spin_indices(spin_id2)
866 dp_to.mol[mol_index].res[res_index].spin[spin_index] = deepcopy(spin2)
867
868
869 else:
870
871 dp_to.mol = deepcopy(dp_from.mol)
872
873
874 - def eliminate(self, name, value, model_info, args, sim=None):
875 """Model-free model elimination, parameter by parameter.
876
877 @param name: The parameter name.
878 @type name: str
879 @param value: The parameter value.
880 @type value: float
881 @param model_info: The model index from model_loop(). This is zero for the global models or equal to the global spin index (which covers the molecule, residue, and spin indices).
882 @type model_info: int
883 @param args: The c1 and c2 elimination constant overrides.
884 @type args: None or tuple of float
885 @keyword sim: The Monte Carlo simulation index.
886 @type sim: int
887 @return: True if the model is to be eliminated, False otherwise.
888 @rtype: bool
889 """
890
891
892 c1 = 50.0 * 1e-9
893 c2 = 1.5
894
895
896 if args != None:
897 c1, c2 = args
898
899
900 model_type = determine_model_type()
901
902
903 if model_type != 'mf' and model_type != 'local_tm':
904 raise RelaxError("Elimination of the global model is not yet supported.")
905
906
907 spin, spin_id = return_spin_from_index(model_info, return_spin_id=True)
908
909
910 if model_type == 'local_tm':
911 tm = spin.local_tm
912 else:
913 tm = cdp.diff_tensor.tm
914
915
916 if tm == None:
917 return False
918
919
920 if name == 'local_tm' and value >= c1:
921 if sim == None:
922 print("Data pipe '%s': The local tm parameter of %.5g is greater than %.5g, eliminating spin system '%s'." % (pipes.cdp_name(), value, c1, spin_id))
923 else:
924 print("Data pipe '%s': The local tm parameter of %.5g is greater than %.5g, eliminating simulation %i of spin system '%s'." % (pipes.cdp_name(), value, c1, sim, spin_id))
925 return True
926
927
928 if match('t[efs]', name) and value >= c2 * tm:
929 if sim == None:
930 print("Data pipe '%s': The %s value of %.5g is greater than %.5g, eliminating spin system '%s'." % (pipes.cdp_name(), name, value, c2*tm, spin_id))
931 else:
932 print("Data pipe '%s': The %s value of %.5g is greater than %.5g, eliminating simulation %i of spin system '%s'." % (pipes.cdp_name(), name, value, c2*tm, sim, spin_id))
933 return True
934
935
936 return False
937
938
960
961
994
995
996 - def grid_search(self, lower=None, upper=None, inc=None, constraints=True, verbosity=1, sim_index=None):
997 """The model-free grid search function.
998
999 @keyword lower: The lower bounds of the grid search which must be equal to the
1000 number of parameters in the model.
1001 @type lower: array of numbers
1002 @keyword upper: The upper bounds of the grid search which must be equal to the
1003 number of parameters in the model.
1004 @type upper: array of numbers
1005 @keyword inc: The increments for each dimension of the space for the grid search.
1006 The number of elements in the array must equal to the number of
1007 parameters in the model.
1008 @type inc: array of int
1009 @keyword constraints: If True, constraints are applied during the grid search (eliminating
1010 parts of the grid). If False, no constraints are used.
1011 @type constraints: bool
1012 @keyword verbosity: A flag specifying the amount of information to print. The higher
1013 the value, the greater the verbosity.
1014 @type verbosity: int
1015 @keyword sim_index: The index of the simulation to apply the grid search to. If None,
1016 the normal model is optimised.
1017 @type sim_index: int
1018 """
1019
1020
1021 self.minimise(min_algor='grid', lower=lower, upper=upper, inc=inc, constraints=constraints, verbosity=verbosity, sim_index=sim_index)
1022
1023
1025 """Create bounds for the OpenDX mapping function.
1026
1027 @param param: The name of the parameter to return the lower and upper bounds of.
1028 @type param: str
1029 @param spin_id: The spin identification string.
1030 @type spin_id: str
1031 @return: The upper and lower bounds of the parameter.
1032 @rtype: list of float
1033 """
1034
1035
1036 if self._PARAMS.scope(param) == 'global':
1037 return diffusion_tensor.map_bounds(param)
1038
1039
1040 spin = return_spin(spin_id)
1041
1042
1043 if search('^s2', param):
1044 return [0.0, 1.0]
1045
1046
1047 elif search('^t', param) or param == 'local_tm':
1048 return [0.0, 1e-8]
1049
1050
1051 elif param == 'rex':
1052 return [0.0, 30.0 / (2.0 * pi * cdp.spectrometer_frq[cdp.ri_ids[0]])**2]
1053
1054
1055 elif param == 'r':
1056 return [1.0 * 1e-10, 1.1 * 1e-10]
1057
1058
1059 elif param == 'csa':
1060 return [-100 * 1e-6, -300 * 1e-6]
1061
1062
1063 - def minimise(self, min_algor=None, min_options=None, func_tol=None, grad_tol=None, max_iterations=None, constraints=False, scaling=True, verbosity=0, sim_index=None, lower=None, upper=None, inc=None):
1064 """Model-free minimisation function.
1065
1066 Three categories of models exist for which the approach to minimisation is different. These
1067 are:
1068
1069 Single spin optimisations: The 'mf' and 'local_tm' model types which are the
1070 model-free parameters for single spins, optionally with a local tm parameter. These
1071 models have no global parameters.
1072
1073 Diffusion tensor optimisations: The 'diff' diffusion tensor model type. No spin
1074 specific parameters exist.
1075
1076 Optimisation of everything: The 'all' model type consisting of all model-free and all
1077 diffusion tensor parameters.
1078
1079
1080 @keyword min_algor: The minimisation algorithm to use.
1081 @type min_algor: str
1082 @keyword min_options: An array of options to be used by the minimisation algorithm.
1083 @type min_options: array of str
1084 @keyword func_tol: The function tolerance which, when reached, terminates optimisation.
1085 Setting this to None turns of the check.
1086 @type func_tol: None or float
1087 @keyword grad_tol: The gradient tolerance which, when reached, terminates optimisation.
1088 Setting this to None turns of the check.
1089 @type grad_tol: None or float
1090 @keyword max_iterations: The maximum number of iterations for the algorithm.
1091 @type max_iterations: int
1092 @keyword constraints: If True, constraints are used during optimisation.
1093 @type constraints: bool
1094 @keyword scaling: If True, diagonal scaling is enabled during optimisation to allow
1095 the problem to be better conditioned.
1096 @type scaling: bool
1097 @keyword verbosity: The amount of information to print. The higher the value, the
1098 greater the verbosity.
1099 @type verbosity: int
1100 @keyword sim_index: The index of the simulation to optimise. This should be None if
1101 normal optimisation is desired.
1102 @type sim_index: None or int
1103 @keyword lower: The lower bounds of the grid search which must be equal to the
1104 number of parameters in the model. This optional argument is only
1105 used when doing a grid search.
1106 @type lower: array of numbers
1107 @keyword upper: The upper bounds of the grid search which must be equal to the
1108 number of parameters in the model. This optional argument is only
1109 used when doing a grid search.
1110 @type upper: array of numbers
1111 @keyword inc: The increments for each dimension of the space for the grid search.
1112 The number of elements in the array must equal to the number of
1113 parameters in the model. This argument is only used when doing a
1114 grid search.
1115 @type inc: array of int
1116 """
1117
1118
1119 if not exists_mol_res_spin_data():
1120 raise RelaxNoSequenceError
1121
1122
1123 for spin in spin_loop():
1124
1125 if not spin.select:
1126 continue
1127
1128
1129 if not spin.model:
1130 raise RelaxNoModelError
1131
1132
1133 if not hasattr(spin, 'isotope'):
1134 raise RelaxSpinTypeError
1135
1136
1137 if sim_index == None and min_algor != 'back_calc':
1138 reset_min_stats()
1139
1140
1141 data_store = Data_container()
1142 opt_params = Data_container()
1143
1144
1145 data_store.h_bar = h_bar
1146 data_store.mu0 = mu0
1147 opt_params.min_algor = min_algor
1148 opt_params.min_options = min_options
1149 opt_params.func_tol = func_tol
1150 opt_params.grad_tol = grad_tol
1151 opt_params.max_iterations = max_iterations
1152
1153
1154 opt_params.verbosity = verbosity
1155
1156
1157 data_store.model_type = determine_model_type()
1158 if not data_store.model_type:
1159 return
1160
1161
1162 if min_algor == 'back_calc' and data_store.model_type != 'local_tm':
1163 data_store.model_type = 'mf'
1164
1165
1166 if data_store.model_type != 'local_tm' and not diffusion_tensor.diff_data_exists():
1167 raise RelaxNoTensorError('diffusion')
1168
1169
1170 if data_store.model_type != 'local_tm' and cdp.diff_tensor.type != 'sphere':
1171
1172 if not hasattr(cdp, 'structure'):
1173 raise RelaxNoPdbError
1174
1175
1176 for spin, spin_id in spin_loop(return_id=True):
1177
1178 if not spin.select:
1179 continue
1180
1181
1182 interatoms = return_interatom_list(spin_id)
1183
1184
1185 for i in range(len(interatoms)):
1186
1187 if not interatoms[i].dipole_pair:
1188 continue
1189
1190
1191 if not hasattr(interatoms[i], 'vector'):
1192 raise RelaxNoVectorsError
1193
1194
1195 if data_store.model_type == 'diff':
1196
1197 for spin in spin_loop():
1198 unset_param = are_mf_params_set(spin)
1199 if unset_param != None:
1200 raise RelaxNoValueError(unset_param)
1201
1202
1203 if verbosity >= 1:
1204 if data_store.model_type == 'mf':
1205 print("Only the model-free parameters for single spins will be used.")
1206 elif data_store.model_type == 'local_mf':
1207 print("Only a local tm value together with the model-free parameters for single spins will be used.")
1208 elif data_store.model_type == 'diff':
1209 print("Only diffusion tensor parameters will be used.")
1210 elif data_store.model_type == 'all':
1211 print("The diffusion tensor parameters together with the model-free parameters for all spins will be used.")
1212
1213
1214 for spin, spin_id in spin_loop(return_id=True):
1215
1216 if not spin.select:
1217 continue
1218
1219
1220 if not hasattr(spin, 'csa') or spin.csa == None:
1221 raise RelaxNoValueError("CSA")
1222
1223
1224 interatoms = return_interatom_list(spin_id)
1225
1226
1227 count = 0
1228 for i in range(len(interatoms)):
1229
1230 if not interatoms[i].dipole_pair:
1231 continue
1232
1233
1234 if not hasattr(interatoms[i], 'r') or interatoms[i].r == None:
1235 raise RelaxNoValueError("interatomic distance", spin_id=interatoms[i].spin_id1, spin_id2=interatoms[i].spin_id2)
1236
1237
1238 count += 1
1239
1240
1241 if count > 1:
1242 raise RelaxError("The spin '%s' has %s dipolar relaxation interactions defined, but only a maximum of one is currently supported." % (spin_id, count))
1243
1244
1245 if data_store.model_type == 'mf' or data_store.model_type == 'local_tm':
1246 num_instances = count_spins(skip_desel=False)
1247 num_data_sets = 1
1248 data_store.num_spins = 1
1249 elif data_store.model_type == 'diff' or data_store.model_type == 'all':
1250 num_instances = 1
1251 num_data_sets = count_spins(skip_desel=False)
1252 data_store.num_spins = count_spins()
1253
1254
1255 if min_algor == 'back_calc':
1256 num_instances = 1
1257 num_data_sets = 0
1258 data_store.num_spins = 1
1259
1260
1261 processor_box = Processor_box()
1262 processor = processor_box.processor
1263
1264
1265
1266
1267 for i in range(num_instances):
1268
1269 if data_store.model_type == 'diff' or data_store.model_type == 'all':
1270 spin_index = None
1271 spin, data_store.spin_id = None, None
1272 elif min_algor == 'back_calc':
1273 spin_index = opt_params.min_options[0]
1274 spin, data_store.spin_id = return_spin_from_index(global_index=spin_index, return_spin_id=True)
1275 else:
1276 spin_index = i
1277 spin, data_store.spin_id = return_spin_from_index(global_index=spin_index, return_spin_id=True)
1278
1279
1280 if spin and (data_store.model_type == 'mf' or data_store.model_type == 'local_tm') and not min_algor == 'back_calc':
1281
1282 if not spin.select:
1283 continue
1284
1285
1286 if not hasattr(spin, 'ri_data') or not hasattr(spin, 'ri_data_err'):
1287 continue
1288
1289
1290 if spin and (data_store.model_type == 'mf' or data_store.model_type == 'local_tm'):
1291 interatoms = return_interatom_list(data_store.spin_id)
1292 if not len(interatoms):
1293 continue
1294
1295
1296 if min_algor == 'back_calc':
1297
1298 opt_params.param_vector = assemble_param_vector(spin=spin, model_type=data_store.model_type)
1299
1300
1301 data_store.scaling_matrix = None
1302
1303 else:
1304
1305 opt_params.param_vector = assemble_param_vector(spin=spin, sim_index=sim_index)
1306
1307
1308 num_params = len(opt_params.param_vector)
1309
1310
1311 data_store.scaling_matrix = assemble_scaling_matrix(num_params, model_type=data_store.model_type, spin=spin, scaling=scaling)
1312 if len(data_store.scaling_matrix):
1313 opt_params.param_vector = dot(inv(data_store.scaling_matrix), opt_params.param_vector)
1314
1315
1316 opt_params.inc, opt_params.lower, opt_params.upper = None, None, None
1317 if match('^[Gg]rid', min_algor):
1318 opt_params.inc, opt_params.lower, opt_params.upper = grid_search_config(num_params, spin=spin, lower=lower, upper=upper, inc=inc, scaling_matrix=data_store.scaling_matrix)
1319
1320
1321 if match('^[Ss]et', min_algor):
1322 opt_params.min_options = dot(inv(data_store.scaling_matrix), opt_params.min_options)
1323
1324
1325 if constraints:
1326 opt_params.A, opt_params.b = linear_constraints(num_params, model_type=data_store.model_type, spin=spin, scaling_matrix=data_store.scaling_matrix)
1327 else:
1328 opt_params.A, opt_params.b = None, None
1329
1330
1331 minimise_data_setup(data_store, min_algor, num_data_sets, opt_params.min_options, spin=spin, sim_index=sim_index)
1332
1333
1334 if constraints and not match('^[Gg]rid', min_algor):
1335 algor = opt_params.min_options[0]
1336 else:
1337 algor = min_algor
1338
1339
1340 if min_algor == 'back_calc' or match('[Ll][Mm]$', algor) or match('[Ll]evenburg-[Mm]arquardt$', algor):
1341 mf = Mf(init_params=opt_params.param_vector, model_type=data_store.model_type, diff_type=data_store.diff_type, diff_params=data_store.diff_params, scaling_matrix=data_store.scaling_matrix, num_spins=data_store.num_spins, equations=data_store.equations, param_types=data_store.param_types, param_values=data_store.param_values, relax_data=data_store.ri_data, errors=data_store.ri_data_err, bond_length=data_store.r, csa=data_store.csa, num_frq=data_store.num_frq, frq=data_store.frq, num_ri=data_store.num_ri, remap_table=data_store.remap_table, noe_r1_table=data_store.noe_r1_table, ri_labels=data_store.ri_types, gx=data_store.gx, gh=data_store.gh, h_bar=data_store.h_bar, mu0=data_store.mu0, num_params=data_store.num_params, vectors=data_store.xh_unit_vectors)
1342
1343
1344 if match('[Ll][Mm]$', algor) or match('[Ll]evenburg-[Mm]arquardt$', algor):
1345
1346 number_ri = 0
1347 for k in range(len(ri_data_err)):
1348 number_ri = number_ri + len(ri_data_err[k])
1349
1350
1351 lm_error = zeros(number_ri, float64)
1352 index = 0
1353 for k in range(len(ri_data_err)):
1354 lm_error[index:index+len(ri_data_err[k])] = ri_data_err[k]
1355 index = index + len(ri_data_err[k])
1356
1357 opt_params.min_options = opt_params.min_options + (mf.lm_dri, lm_error)
1358
1359
1360 if min_algor == 'back_calc':
1361 return mf.calc_ri()
1362
1363
1364 if match('^[Gg]rid', min_algor) and data_store.model_type == 'diff':
1365
1366 print("Parallelised diffusion tensor grid search.")
1367
1368
1369 for subdivision in grid_split(divisions=processor.processor_size(), lower=opt_params.lower, upper=opt_params.upper, inc=opt_params.inc):
1370
1371 opt_params.subdivision = subdivision
1372
1373
1374 command = MF_grid_command()
1375
1376
1377 command.store_data(deepcopy(data_store), deepcopy(opt_params))
1378
1379
1380 memo = MF_memo(model_free=self, model_type=data_store.model_type, spin=spin, sim_index=sim_index, scaling=scaling, scaling_matrix=data_store.scaling_matrix)
1381 processor.add_to_queue(command, memo)
1382
1383
1384 processor.run_queue()
1385
1386
1387 return
1388
1389
1390 if search('^[Gg]rid', min_algor):
1391 command = MF_grid_command()
1392
1393
1394 else:
1395 command = MF_minimise_command()
1396
1397
1398 command.store_data(deepcopy(data_store), deepcopy(opt_params))
1399
1400
1401 memo = MF_memo(model_free=self, model_type=data_store.model_type, spin=spin, sim_index=sim_index, scaling=scaling, scaling_matrix=data_store.scaling_matrix)
1402 processor.add_to_queue(command, memo)
1403
1404
1405 processor.run_queue()
1406
1407
1409 """Return a description of the model.
1410
1411 @param model_info: The model index from model_loop(). This is zero for the global models or equal to the global spin index (which covers the molecule, residue, and spin indices).
1412 @type model_info: int
1413 @return: The model description.
1414 @rtype: str
1415 """
1416
1417
1418 model_type = determine_model_type()
1419
1420
1421 if model_type == 'all':
1422 return "Global model - all diffusion tensor parameters and spin specific model-free parameters."
1423 elif model_type == 'diff':
1424 return "Diffusion tensor model."
1425
1426
1427 else:
1428
1429 spin, spin_id = return_spin_from_index(model_info, return_spin_id=True)
1430
1431
1432 return "Model-free model of spin '%s'." % spin_id
1433
1434
1436 """Generator method for looping over the models (global or local).
1437
1438 If the model type is 'all' or 'diff', then this yields the single value of zero. Otherwise
1439 the global spin index is yielded.
1440
1441
1442 @return: The model index. This is zero for the global models or equal to the global spin
1443 index (which covers the molecule, residue, and spin indices).
1444 @rtype: int
1445 """
1446
1447
1448 model_type = determine_model_type()
1449
1450
1451 if model_type == 'all' or model_type == 'diff':
1452 yield 0
1453
1454
1455 else:
1456
1457 global_index = -1
1458 for spin in spin_loop():
1459
1460 global_index = global_index + 1
1461
1462
1463 yield global_index
1464
1465
1467 """Return the k, n, and chi2 model statistics.
1468
1469 k - number of parameters.
1470 n - number of data points.
1471 chi2 - the chi-squared value.
1472
1473
1474 @keyword model_info: The model index from model_loop(). This is zero for the global models or equal to the global spin index (which covers the molecule, residue, and spin indices).
1475 @type model_info: int
1476 @keyword spin_id: The spin identification string. Either this or the instance keyword argument must be supplied.
1477 @type spin_id: None or str
1478 @keyword global_stats: A parameter which determines if global or local statistics are returned. If None, then the appropriateness of global or local statistics is automatically determined.
1479 @type global_stats: None or bool
1480 @return: The optimisation statistics, in tuple format, of the number of parameters (k), the number of data points (n), and the chi-squared value (chi2).
1481 @rtype: tuple of (int, int, float)
1482 """
1483
1484
1485 if model_info == None and spin_id == None:
1486 raise RelaxError("Either the model_info or spin_id argument must be supplied.")
1487 elif model_info != None and spin_id != None:
1488 raise RelaxError("The model_info arg " + repr(model_info) + " and spin_id arg " + repr(spin_id) + " clash. Only one should be supplied.")
1489
1490
1491 model_type = determine_model_type()
1492
1493
1494 if global_stats == None:
1495 if model_type in ['mf', 'local_tm']:
1496 global_stats = False
1497 else:
1498 global_stats = True
1499
1500
1501 if not global_stats:
1502
1503 if spin_id:
1504 spin = return_spin(spin_id)
1505 else:
1506 spin = return_spin_from_index(model_info)
1507
1508
1509 if not spin.select:
1510 return None, None, None
1511
1512
1513 if not hasattr(spin, 'ri_data'):
1514 return None, None, None
1515
1516
1517 param_vector = assemble_param_vector(spin=spin)
1518 k = len(param_vector)
1519
1520
1521 n = len(spin.ri_data)
1522
1523
1524 chi2 = spin.chi2
1525
1526
1527 elif global_stats:
1528
1529 param_vector = assemble_param_vector()
1530 k = len(param_vector)
1531
1532
1533 n = 0
1534 chi2 = 0
1535 for spin in spin_loop():
1536
1537 if not spin.select:
1538 continue
1539
1540
1541 if not hasattr(spin, 'ri_data') or not len(spin.ri_data):
1542 continue
1543
1544 n = n + len(spin.ri_data)
1545
1546
1547 if model_type == 'local_tm':
1548 chi2 = chi2 + spin.chi2
1549
1550
1551 if model_type != 'local_tm':
1552 if not hasattr(cdp, 'chi2'):
1553 raise RelaxError("Global statistics are not available, most likely because the global model has not been optimised.")
1554 chi2 = cdp.chi2
1555
1556
1557 return k, n, chi2
1558
1559
1561 """Return the type of the model, either being 'local' or 'global'.
1562
1563 @return: The model type, one of 'local' or 'global'.
1564 @rtype: str
1565 """
1566
1567
1568 model_type = determine_model_type()
1569
1570
1571 if model_type in ['all', 'diff']:
1572 return 'global'
1573
1574
1575 else:
1576 return 'local'
1577
1578
1580 """Function for returning the number of instances.
1581
1582 @return: The number of instances used for optimisation. Either the number of spins if
1583 the local optimisations are setup ('mf' and 'local_tm'), or 1 for the global
1584 models.
1585 @rtype: int
1586 """
1587
1588
1589 if not exists_mol_res_spin_data():
1590 return 0
1591
1592
1593 model_type = determine_model_type()
1594
1595
1596 if model_type == 'mf' or model_type == 'local_tm':
1597 return count_spins()
1598
1599
1600 elif model_type == 'diff' or model_type == 'all':
1601 return 1
1602
1603
1604 else:
1605 raise RelaxFault
1606
1607
1609 """Deselect spins which have insufficient data to support minimisation.
1610
1611 @keyword data_check: A flag to signal if the presence of base data is to be checked for.
1612 @type data_check: bool
1613 @keyword verbose: A flag which if True will allow printouts.
1614 @type verbose: bool
1615 """
1616
1617
1618 if verbose:
1619 print("\nOver-fit spin deselection:")
1620
1621
1622 if not exists_mol_res_spin_data():
1623 raise RelaxNoSequenceError
1624
1625
1626 need_vect = False
1627 if hasattr(cdp, 'diff_tensor') and (cdp.diff_tensor.type == 'spheroid' or cdp.diff_tensor.type == 'ellipsoid'):
1628 need_vect = True
1629
1630
1631 deselect_flag = False
1632 spin_count = 0
1633 for spin, spin_id in spin_loop(return_id=True):
1634
1635 if not spin.select:
1636 continue
1637
1638
1639 interatoms = interatomic.return_interatom_list(spin_id)
1640
1641
1642 dipole_relax = False
1643 for i in range(len(interatoms)):
1644
1645 if not interatoms[i].dipole_pair:
1646 continue
1647
1648
1649 if spin_id != interatoms[i].spin_id1:
1650 spin_id2 = interatoms[i].spin_id1
1651 else:
1652 spin_id2 = interatoms[i].spin_id2
1653 spin2 = return_spin(spin_id2)
1654
1655
1656 dipole_relax = True
1657
1658
1659 if not dipole_relax or not hasattr(spin, 'csa') or spin.csa == None:
1660 warn(RelaxDeselectWarning(spin_id, 'an absence of relaxation mechanisms'))
1661 spin.select = False
1662 deselect_flag = True
1663 continue
1664
1665
1666 if data_check:
1667
1668 data_points = 0
1669 inf_data = False
1670 if hasattr(cdp, 'ri_ids') and hasattr(spin, 'ri_data'):
1671 for id in cdp.ri_ids:
1672 if id in spin.ri_data and spin.ri_data[id] != None:
1673 data_points += 1
1674
1675
1676 if isInf(spin.ri_data[id]):
1677 inf_data = True
1678
1679
1680 if inf_data:
1681 warn(RelaxDeselectWarning(spin_id, 'infinite relaxation data'))
1682 spin.select = False
1683 deselect_flag = True
1684 continue
1685
1686
1687 if not hasattr(spin, 'ri_data'):
1688 warn(RelaxDeselectWarning(spin_id, 'missing relaxation data'))
1689 spin.select = False
1690 deselect_flag = True
1691 continue
1692
1693
1694 elif data_points < 3:
1695 warn(RelaxDeselectWarning(spin_id, 'insufficient relaxation data, 3 or more data points are required'))
1696 spin.select = False
1697 deselect_flag = True
1698 continue
1699
1700
1701 elif hasattr(spin, 'params') and spin.params and len(spin.params) > data_points:
1702 warn(RelaxDeselectWarning(spin_id, 'over-fitting - more parameters than data points'))
1703 spin.select = False
1704 deselect_flag = True
1705 continue
1706
1707
1708 for i in range(len(interatoms)):
1709
1710 if not interatoms[i].dipole_pair:
1711 continue
1712
1713
1714 if need_vect:
1715 if not hasattr(interatoms[i], 'vector') or interatoms[i].vector == None:
1716 warn(RelaxDeselectWarning(spin_id, 'missing structural data'))
1717 spin.select = False
1718 deselect_flag = True
1719 continue
1720
1721
1722 spin_count += 1
1723
1724
1725 if spin_count == 0:
1726 warn(RelaxWarning("No spins are selected therefore the optimisation or calculation cannot proceed."))
1727
1728
1729 if verbose and not deselect_flag:
1730 print("No spins have been deselected.")
1731
1732
1733 - def set_error(self, model_info, index, error):
1734 """Set the parameter errors.
1735
1736 @param model_info: The model index from model_loop(). This is zero for the global models or equal to the global spin index (which covers the molecule, residue, and spin indices).
1737 @type model_info: int
1738 @param index: The index of the parameter to set the errors for.
1739 @type index: int
1740 @param error: The error value.
1741 @type error: float
1742 """
1743
1744
1745 inc = 0
1746
1747
1748 model_type = determine_model_type()
1749
1750
1751 param_names = self.data_names(set='params', scope='spin')
1752
1753
1754
1755
1756
1757 if model_type == 'diff' or model_type == 'all':
1758
1759 if cdp.diff_tensor.type == 'sphere':
1760
1761 if index == 0:
1762 cdp.diff_tensor.set(param='tm', value=error, category='err')
1763
1764
1765 inc = inc + 1
1766
1767
1768 elif cdp.diff_tensor.type == 'spheroid':
1769
1770 if index == 0:
1771 cdp.diff_tensor.set(param='tm', value=error, category='err')
1772 elif index == 1:
1773 cdp.diff_tensor.set(param='Da', value=error, category='err')
1774 elif index == 2:
1775 cdp.diff_tensor.set(param='theta', value=error, category='err')
1776 elif index == 3:
1777 cdp.diff_tensor.set(param='phi', value=error, category='err')
1778
1779
1780 inc = inc + 4
1781
1782
1783 elif cdp.diff_tensor.type == 'ellipsoid':
1784
1785 if index == 0:
1786 cdp.diff_tensor.set(param='tm', value=error, category='err')
1787 elif index == 1:
1788 cdp.diff_tensor.set(param='Da', value=error, category='err')
1789 elif index == 2:
1790 cdp.diff_tensor.set(param='Dr', value=error, category='err')
1791 elif index == 3:
1792 cdp.diff_tensor.set(param='alpha', value=error, category='err')
1793 elif index == 4:
1794 cdp.diff_tensor.set(param='beta', value=error, category='err')
1795 elif index == 5:
1796 cdp.diff_tensor.set(param='gamma', value=error, category='err')
1797
1798
1799 inc = inc + 6
1800
1801
1802
1803
1804
1805 if model_type == 'all':
1806
1807 for spin in spin_loop():
1808
1809 if not spin.select:
1810 continue
1811
1812
1813 for param in param_names:
1814
1815 if index == inc:
1816 setattr(spin, param + "_err", error)
1817
1818
1819 inc = inc + 1
1820
1821
1822
1823
1824
1825 if model_type == 'mf' or model_type == 'local_tm':
1826
1827 spin = return_spin_from_index(model_info)
1828
1829
1830 if not spin.select:
1831 return
1832
1833
1834 for param in param_names:
1835
1836 if index == inc:
1837 setattr(spin, param + "_err", error)
1838
1839
1840 inc = inc + 1
1841
1842
1843 - def set_param_values(self, param=None, value=None, index=None, spin_id=None, error=False, force=True):
1844 """Set the model-free parameter values.
1845
1846 @keyword param: The parameter name list.
1847 @type param: list of str
1848 @keyword value: The parameter value list.
1849 @type value: list
1850 @keyword index: The index for parameters which are of the list-type. This is unused.
1851 @type index: None or int
1852 @keyword spin_id: The spin identification string, only used for spin specific parameters.
1853 @type spin_id: None or str
1854 @keyword error: A flag which if True will allow the parameter errors to be set instead of the values.
1855 @type error: bool
1856 @keyword force: A flag which if True will cause current values to be overwritten. If False, a RelaxError will raised if the parameter value is already set.
1857 @type force: bool
1858 """
1859
1860
1861 is_str_list(param, 'parameter name')
1862
1863
1864 diff_params = []
1865 diff_vals = []
1866 mf_params = []
1867 mf_vals = []
1868 for i in range(len(param)):
1869
1870 if self._PARAMS.scope(param[i]) == 'global':
1871 if error:
1872 diff_params.append(param[i] + '_err')
1873 else:
1874 diff_params.append(param[i])
1875 diff_vals.append(value[i])
1876
1877
1878 else:
1879 mf_params.append(param[i])
1880 mf_vals.append(value[i])
1881
1882
1883 if diff_params:
1884 diffusion_tensor.set(value=diff_vals, param=diff_params)
1885
1886
1887 for i in range(len(mf_params)):
1888
1889 if mf_params[i] not in self.data_names(set='params', scope='spin') and mf_params[i] not in self.data_names(set='generic', scope='spin'):
1890 raise RelaxError("The parameter '%s' is unknown. It should be one of %s or %s" % (mf_params[i], self.data_names(set='params', scope='spin'), self.data_names(set='generic', scope='spin')))
1891
1892
1893 if error:
1894 mf_params[i] += '_err'
1895
1896
1897 for spin in spin_loop(spin_id):
1898 setattr(spin, mf_params[i], mf_vals[i])
1899
1900
1902 """Set all simulation selection flags.
1903
1904 @param model_info: The model index from model_loop(). This is zero for the global models or equal to the global spin index (which covers the molecule, residue, and spin indices).
1905 @type model_info: int
1906 @param select_sim: The selection flags.
1907 @type select_sim: bool
1908 """
1909
1910
1911 model_type = determine_model_type()
1912
1913
1914 if model_type == 'all' or model_type == 'diff':
1915 cdp.select_sim = select_sim
1916
1917
1918 else:
1919
1920 spin = return_spin_from_index(model_info)
1921
1922
1923 if not spin.select:
1924 return
1925
1926
1927 spin.select_sim = deepcopy(select_sim)
1928
1929
1931 """Function to update the other model-free parameters.
1932
1933 @param param: The name of the parameter which has been changed.
1934 @type param: str
1935 @param spin: The SpinContainer object.
1936 @type spin: SpinContainer
1937 """
1938
1939
1940 if param == 's2f':
1941
1942 if hasattr(spin, 's2s') and spin.s2s != None:
1943 spin.s2 = spin.s2f * spin.s2s
1944
1945
1946
1947 if param == 's2s':
1948
1949 if hasattr(spin, 's2f') and spin.s2f != None:
1950 spin.s2 = spin.s2f * spin.s2s
1951
1952
1954 """Initialise the Monte Carlo parameter values."""
1955
1956
1957 model_type = determine_model_type()
1958
1959
1960 param_names = self.data_names(set='params', scope='spin')
1961
1962
1963 min_names = self.data_names(set='min', scope='spin')
1964
1965
1966 if model_type == 'diff' or model_type == 'all':
1967
1968 if cdp.diff_tensor.type == 'sphere':
1969 diff_params = ['tm']
1970
1971
1972 elif cdp.diff_tensor.type == 'spheroid':
1973 diff_params = ['tm', 'Da', 'theta', 'phi']
1974
1975
1976 elif cdp.diff_tensor.type == 'ellipsoid':
1977 diff_params = ['tm', 'Da', 'Dr', 'alpha', 'beta', 'gamma']
1978
1979
1980
1981
1982
1983
1984 if model_type == 'diff' or model_type == 'all':
1985
1986 for object_name in diff_params:
1987
1988 sim_object_name = object_name + '_sim'
1989
1990
1991 if hasattr(cdp.diff_tensor, sim_object_name):
1992 raise RelaxError("Monte Carlo parameter values have already been set.")
1993
1994
1995 for object_name in min_names:
1996
1997 sim_object_name = object_name + '_sim'
1998
1999
2000 if hasattr(cdp, sim_object_name):
2001 raise RelaxError("Monte Carlo parameter values have already been set.")
2002
2003
2004 if model_type != 'diff':
2005 for spin in spin_loop():
2006
2007 if not spin.select:
2008 continue
2009
2010
2011 for object_name in param_names:
2012
2013 sim_object_name = object_name + '_sim'
2014
2015
2016 if hasattr(spin, sim_object_name):
2017 raise RelaxError("Monte Carlo parameter values have already been set.")
2018
2019
2020
2021
2022
2023
2024 for object_name in min_names:
2025
2026 if not hasattr(cdp, object_name):
2027 continue
2028
2029
2030 sim_object_name = object_name + '_sim'
2031
2032
2033 setattr(cdp, sim_object_name, [])
2034
2035
2036 sim_object = getattr(cdp, sim_object_name)
2037
2038
2039 for j in range(cdp.sim_number):
2040
2041 object = getattr(cdp, object_name)
2042
2043
2044 sim_object.append(deepcopy(object))
2045
2046
2047 if model_type == 'diff' or model_type == 'all':
2048
2049 cdp.diff_tensor.set_sim_num(cdp.sim_number)
2050
2051
2052 for object_name in diff_params:
2053 for j in range(cdp.sim_number):
2054 cdp.diff_tensor.set(param=object_name, value=deepcopy(getattr(cdp.diff_tensor, object_name)), category='sim', sim_index=j)
2055
2056
2057 if model_type != 'diff':
2058 for spin in spin_loop():
2059
2060 if not spin.select:
2061 continue
2062
2063
2064 for object_name in param_names:
2065
2066 sim_object_name = object_name + '_sim'
2067
2068
2069 setattr(spin, sim_object_name, [])
2070
2071
2072 sim_object = getattr(spin, sim_object_name)
2073
2074
2075 for j in range(cdp.sim_number):
2076
2077 sim_object.append(deepcopy(getattr(spin, object_name)))
2078
2079
2080 for object_name in min_names:
2081
2082 sim_object_name = object_name + '_sim'
2083
2084
2085 setattr(spin, sim_object_name, [])
2086
2087
2088 sim_object = getattr(spin, sim_object_name)
2089
2090
2091 for j in range(cdp.sim_number):
2092
2093 sim_object.append(deepcopy(getattr(spin, object_name)))
2094
2095
2097 """Return the simulation chi-squared values.
2098
2099 @param model_info: The model index from model_loop(). This is zero for the global models or equal to the global spin index (which covers the molecule, residue, and spin indices).
2100 @type model_info: int
2101 @keyword index: The optional simulation index.
2102 @type index: int
2103 @return: The list of simulation chi-squared values. If the index is supplied, only a single value will be returned.
2104 @rtype: list of float or float
2105 """
2106
2107
2108 model_type = determine_model_type()
2109
2110
2111 if model_type == 'all' or model_type == 'diff':
2112 return cdp.chi2_sim
2113
2114
2115 else:
2116
2117 spin = return_spin_from_index(model_info)
2118
2119
2120 return spin.chi2_sim
2121
2122
2124 """Return the array of simulation parameter values.
2125
2126 @param model_info: The model index from model_loop(). This is zero for the global models or equal to the global spin index (which covers the molecule, residue, and spin indices).
2127 @type model_info: int
2128 @param index: The index of the parameter to return the array of values for.
2129 @type index: int
2130 @return: The array of simulation parameter values.
2131 @rtype: list of float
2132 """
2133
2134
2135 inc = 0
2136
2137
2138 model_type = determine_model_type()
2139
2140
2141 param_names = self.data_names(set='params', scope='spin')
2142
2143
2144
2145
2146
2147 if model_type == 'diff' or model_type == 'all':
2148
2149 if cdp.diff_tensor.type == 'sphere':
2150
2151 if index == 0:
2152 return cdp.diff_tensor.tm_sim
2153
2154
2155 inc = inc + 1
2156
2157
2158 elif cdp.diff_tensor.type == 'spheroid':
2159
2160 if index == 0:
2161 return cdp.diff_tensor.tm_sim
2162 elif index == 1:
2163 return cdp.diff_tensor.Da_sim
2164 elif index == 2:
2165 return cdp.diff_tensor.theta_sim
2166 elif index == 3:
2167 return cdp.diff_tensor.phi_sim
2168
2169
2170 inc = inc + 4
2171
2172
2173 elif cdp.diff_tensor.type == 'ellipsoid':
2174
2175 if index == 0:
2176 return cdp.diff_tensor.tm_sim
2177 elif index == 1:
2178 return cdp.diff_tensor.Da_sim
2179 elif index == 2:
2180 return cdp.diff_tensor.Dr_sim
2181 elif index == 3:
2182 return cdp.diff_tensor.alpha_sim
2183 elif index == 4:
2184 return cdp.diff_tensor.beta_sim
2185 elif index == 5:
2186 return cdp.diff_tensor.gamma_sim
2187
2188
2189 inc = inc + 6
2190
2191
2192
2193
2194
2195 if model_type == 'all':
2196
2197 for spin in spin_loop():
2198
2199 if not spin.select:
2200 continue
2201
2202
2203 for param in param_names:
2204
2205 if index == inc:
2206 return getattr(spin, param + "_sim")
2207
2208
2209 inc = inc + 1
2210
2211
2212
2213
2214
2215 if model_type == 'mf' or model_type == 'local_tm':
2216
2217 spin = return_spin_from_index(model_info)
2218
2219
2220 if not spin.select:
2221 return
2222
2223
2224 for param in param_names:
2225
2226 if index == inc:
2227 return getattr(spin, param + "_sim")
2228
2229
2230 inc = inc + 1
2231
2232
2234 """Return the array of selected simulation flags for the spin.
2235
2236 @param model_info: The model index from model_loop(). This is zero for the global models or equal to the global spin index (which covers the molecule, residue, and spin indices).
2237 @type model_info: int
2238 @return: The array of selected simulation flags.
2239 @rtype: list of int
2240 """
2241
2242
2243 model_type = determine_model_type()
2244
2245
2246 if model_type == 'all' or model_type == 'diff':
2247 return cdp.select_sim
2248
2249
2250 else:
2251
2252 spin = return_spin_from_index(model_info)
2253
2254
2255 if not spin.select:
2256 return
2257
2258
2259 return spin.select_sim
2260
2261
2263 """Skip certain data.
2264
2265 @param model_info: The model index from model_loop(). This is zero for the global models or equal to the global spin index (which covers the molecule, residue, and spin indices).
2266 @type model_info: int
2267 @return: True if the data should be skipped, False otherwise.
2268 @rtype: bool
2269 """
2270
2271
2272 model_type = determine_model_type()
2273
2274
2275 if (model_type == 'mf' or model_type == 'local_tm') and not return_spin_from_index(model_info).select:
2276 return True
2277
2278
2279 return False
2280
2281
2282
2284 """Empty class to be used for data storage."""
2285