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