1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23 """Analysis specific code for the Frame Order theory of domain dynamics."""
24
25
26 from copy import deepcopy
27 from math import cos, pi
28 from minfx.generic import generic_minimise
29 from minfx.grid import grid_point_array
30 from numpy import arccos, array, dot, eye, float64, ones, transpose, zeros
31 from re import search
32 from warnings import warn
33
34
35 from float import isNaN, isInf
36 from generic_fns import align_tensor, pipes
37 from generic_fns.angles import wrap_angles
38 from generic_fns.structure.cones import Iso_cone, Pseudo_elliptic
39 from generic_fns.structure.geometric import create_cone_pdb, generate_vector_dist, generate_vector_residues
40 from generic_fns.structure.internal import Internal
41 from maths_fns import frame_order, order_parameters
42 from maths_fns.coord_transform import spherical_to_cartesian
43 from maths_fns.rotation_matrix import euler_to_R_zyz, two_vect_to_R
44 from relax_errors import RelaxError, RelaxInfError, RelaxModelError, RelaxNaNError, RelaxNoModelError
45 from relax_io import open_write_file
46 from relax_warnings import RelaxWarning, RelaxDeselectWarning
47 from specific_fns.api_base import API_base
48 from specific_fns.api_common import API_common
49
50
52 """Class containing the specific methods of the Frame Order theories."""
53
55 """Initialise the class by placing API_common methods into the API."""
56
57
58 super(Frame_order, self).__init__()
59
60
61 self.eliminate = self._eliminate_false
62 self.overfit_deselect = self._overfit_deselect_dummy
63 self.return_conversion_factor = self._return_no_conversion_factor
64 self.set_param_values = self._set_param_values_global
65
66
67 self.PARAMS.add('ave_pos_alpha', scope='global', units='rad', desc='The average position alpha Euler angle', py_type=float, set='params', err=True, sim=True)
68 self.PARAMS.add('ave_pos_beta', scope='global', units='rad', desc='The average position beta Euler angle', py_type=float, set='params', err=True, sim=True)
69 self.PARAMS.add('ave_pos_gamma', scope='global', units='rad', desc='The average position gamma Euler angle', py_type=float, set='params', err=True, sim=True)
70 self.PARAMS.add('eigen_alpha', scope='global', units='rad', desc='The Eigenframe alpha Euler angle', py_type=float, set='params', err=True, sim=True)
71 self.PARAMS.add('eigen_beta', scope='global', units='rad', desc='The Eigenframe beta Euler angle', py_type=float, set='params', err=True, sim=True)
72 self.PARAMS.add('eigen_gamma', scope='global', units='rad', desc='The Eigenframe gamma Euler angle', py_type=float, set='params', err=True, sim=True)
73 self.PARAMS.add('axis_theta', scope='global', units='rad', desc='The cone axis polar angle (for the isotropic cone model)', py_type=float, set='params', err=True, sim=True)
74 self.PARAMS.add('axis_phi', scope='global', units='rad', desc='The cone axis azimuthal angle (for the isotropic cone model)', py_type=float, set='params', err=True, sim=True)
75 self.PARAMS.add('cone_theta_x', scope='global', units='rad', desc='The pseudo-ellipse cone opening half-angle for the x-axis', py_type=float, set='params', err=True, sim=True)
76 self.PARAMS.add('cone_theta_y', scope='global', units='rad', desc='The pseudo-ellipse cone opening half-angle for the y-axis', py_type=float, set='params', err=True, sim=True)
77 self.PARAMS.add('cone_theta', scope='global', units='rad', desc='The isotropic cone opening half-angle', py_type=float, set='params', err=True, sim=True)
78 self.PARAMS.add('cone_s1', scope='global', units='', desc='The isotropic cone order parameter', py_type=float, set='params', err=True, sim=True)
79 self.PARAMS.add('cone_sigma_max', scope='global', units='rad', desc='The torsion angle', py_type=float, set='params', err=True, sim=True)
80 self.PARAMS.add('params', scope='global', desc='The model parameters', py_type=list)
81
82
83 self.PARAMS.add_min_data(min_stats_global=True)
84
85
87 """Assemble and return the limit vectors.
88
89 @return: The lower and upper limit vectors.
90 @rtype: numpy rank-1 array, numpy rank-1 array
91 """
92
93
94 lower = zeros(len(cdp.params), float64)
95 upper = 2.0*pi * ones(len(cdp.params), float64)
96
97
98 return lower, upper
99
100
102 """Assemble and return the parameter vector.
103
104 @return: The parameter vector.
105 @rtype: numpy rank-1 array
106 @keyword sim_index: The Monte Carlo simulation index.
107 @type sim_index: int
108 """
109
110
111 if sim_index == None:
112
113 if cdp.model in ['free rotor', 'iso cone, torsionless', 'iso cone, free rotor']:
114 param_vect = [cdp.ave_pos_beta, cdp.ave_pos_gamma]
115 else:
116 param_vect = [cdp.ave_pos_alpha, cdp.ave_pos_beta, cdp.ave_pos_gamma]
117
118
119 if cdp.model in ['iso cone', 'pseudo-ellipse', 'pseudo-ellipse, torsionless', 'pseudo-ellipse, free rotor']:
120 param_vect.append(cdp.eigen_alpha)
121 param_vect.append(cdp.eigen_beta)
122 param_vect.append(cdp.eigen_gamma)
123
124
125 elif cdp.model in ['free rotor', 'iso cone, torsionless', 'iso cone, free rotor', 'rotor']:
126 param_vect.append(cdp.axis_theta)
127 param_vect.append(cdp.axis_phi)
128
129
130 if cdp.model in ['pseudo-ellipse', 'pseudo-ellipse, torsionless', 'pseudo-ellipse, free rotor']:
131 param_vect.append(cdp.cone_theta_x)
132 param_vect.append(cdp.cone_theta_y)
133
134
135 elif cdp.model in ['iso cone', 'iso cone, torsionless']:
136 param_vect.append(cdp.cone_theta)
137 elif cdp.model in ['iso cone, free rotor']:
138 param_vect.append(cdp.cone_s1)
139
140
141 if cdp.model in ['rotor', 'line', 'iso cone', 'pseudo-ellipse']:
142 param_vect.append(cdp.cone_sigma_max)
143
144
145 else:
146
147 if cdp.model in ['free rotor', 'iso cone, torsionless', 'iso cone, free rotor']:
148 param_vect = [cdp.ave_pos_beta_sim[sim_index], cdp.ave_pos_gamma_sim[sim_index]]
149 else:
150 param_vect = [cdp.ave_pos_alpha_sim[sim_index], cdp.ave_pos_beta_sim[sim_index], cdp.ave_pos_gamma_sim[sim_index]]
151
152
153 if cdp.model in ['iso cone', 'pseudo-ellipse', 'pseudo-ellipse, torsionless', 'pseudo-ellipse, free rotor']:
154 param_vect.append(cdp.eigen_alpha_sim[sim_index])
155 param_vect.append(cdp.eigen_beta_sim[sim_index])
156 param_vect.append(cdp.eigen_gamma_sim[sim_index])
157
158
159 elif cdp.model in ['free rotor', 'iso cone, torsionless', 'iso cone, free rotor', 'rotor']:
160 param_vect.append(cdp.axis_theta_sim[sim_index])
161 param_vect.append(cdp.axis_phi_sim[sim_index])
162
163
164 if cdp.model in ['pseudo-ellipse', 'pseudo-ellipse, torsionless', 'pseudo-ellipse, free rotor']:
165 param_vect.append(cdp.cone_theta_x_sim[sim_index])
166 param_vect.append(cdp.cone_theta_y_sim[sim_index])
167
168
169 elif cdp.model in ['iso cone', 'iso cone, torsionless']:
170 param_vect.append(cdp.cone_theta_sim[sim_index])
171 elif cdp.model in ['iso cone, free rotor']:
172 param_vect.append(cdp.cone_s1_sim[sim_index])
173
174
175 if cdp.model in ['rotor', 'line', 'iso cone', 'pseudo-ellipse']:
176 param_vect.append(cdp.cone_sigma_max_sim[sim_index])
177
178
179 return array(param_vect, float64)
180
181
183 """Back-calculation of the reduced alignment tensor.
184
185 @return: The reduced alignment tensors.
186 @rtype: numpy array
187 """
188
189
190 param_vector = self._assemble_param_vector()
191
192
193 full_tensors, red_tensors, red_tensor_err, full_in_ref_frame = self._minimise_setup_tensors()
194
195
196 target = frame_order.Frame_order(model=cdp.model, full_tensors=full_tensors, red_tensors=red_tensors, red_errors=red_tensor_err, full_in_ref_frame=full_in_ref_frame)
197
198
199 target.func(param_vector)
200
201
202 self._store_bc_tensors(target)
203
204
205 return target.red_tensors_bc
206
207
208 - def _cone_pdb(self, size=30.0, file=None, dir=None, inc=40, force=False):
209 """Create a PDB file containing a geometric object representing the Frame Order cone models.
210
211 @param size: The size of the geometric object in Angstroms.
212 @type size: float
213 @param inc: The number of increments for the filling of the cone objects.
214 @type inc: int
215 @param file: The name of the PDB file to create.
216 @type file: str
217 @param dir: The name of the directory to place the PDB file into.
218 @type dir: str
219 @param force: Flag which if set to True will cause any pre-existing file to be
220 overwritten.
221 @type force: bool
222 """
223
224
225 pipes.test()
226
227
228 if cdp.model == 'rigid':
229 raise RelaxError("The 'rigid' frame order model has no cone representation.")
230
231
232 if not hasattr(cdp, 'pivot'):
233 raise RelaxError("The pivot point for the domain motion has not been set.")
234
235
236 neg_cone = True
237
238
239 sim = False
240 num_sim = 0
241 if hasattr(cdp, 'sim_number'):
242 sim = True
243 num_sim = cdp.sim_number
244
245
246 inv_mat = -eye(3)
247
248
249 structure = Internal()
250
251
252 model = structure.add_model(model=1)
253 if neg_cone:
254 model_neg = structure.add_model(model=2)
255
256
257 structure.add_molecule(name=cdp.model)
258
259
260 mol = model.mol[0]
261 if neg_cone:
262 mol_neg = model_neg.mol[0]
263
264
265
266
267
268
269 structure.add_atom(mol_name=cdp.model, pdb_record='HETATM', atom_num=1, atom_name='R', res_name='PIV', res_num=1, pos=cdp.pivot, element='C')
270
271
272
273
274
275
276 print("\nGenerating the central axis.")
277
278
279 if cdp.model in ['free rotor', 'iso cone, torsionless', 'iso cone, free rotor', 'rotor']:
280 theta_name = 'axis_theta'
281 phi_name = 'axis_phi'
282 else:
283 theta_name = 'eigen_beta'
284 phi_name = 'eigen_gamma'
285
286
287 axis = zeros(3, float64)
288 spherical_to_cartesian([1.0, getattr(cdp, theta_name), getattr(cdp, phi_name)], axis)
289 print("Central axis: %s." % axis)
290
291
292 axis_pos = axis
293 axis_neg = dot(inv_mat, axis)
294
295
296 axis_sim_pos = None
297 axis_sim_neg = None
298 if sim:
299
300 axis_sim = zeros((cdp.sim_number, 3), float64)
301
302
303 for i in range(cdp.sim_number):
304 spherical_to_cartesian([1.0, getattr(cdp, theta_name+'_sim')[i], getattr(cdp, phi_name+'_sim')[i]], axis_sim[i])
305
306
307 axis_sim_pos = axis_sim
308 axis_sim_neg = transpose(dot(inv_mat, transpose(axis_sim_pos)))
309
310
311 print("\nGenerating the axis vectors.")
312 res_num = generate_vector_residues(mol=mol, vector=axis_pos, atom_name='z-ax', res_name_vect='ZAX', sim_vectors=axis_sim_pos, res_num=2, origin=cdp.pivot, scale=size)
313
314
315 if neg_cone:
316 res_num = generate_vector_residues(mol=mol_neg, vector=axis_neg, atom_name='z-ax', res_name_vect='ZAX', sim_vectors=axis_sim_neg, res_num=2, origin=cdp.pivot, scale=size)
317
318
319
320
321
322
323 if cdp.model not in ['free rotor', 'iso cone, torsionless', 'iso cone, free rotor', 'rotor']:
324
325 print("\nGenerating the x and y axes.")
326
327
328 axes = zeros((3, 3), float64)
329 euler_to_R_zyz(cdp.eigen_alpha, cdp.eigen_beta, cdp.eigen_gamma, axes)
330 print("Axis system:\n%s" % axes)
331
332
333 axes_pos = axes
334 axes_neg = dot(inv_mat, axes)
335
336
337 axes_sim_pos = None
338 axes_sim_neg = None
339 if sim:
340
341 axes_sim = zeros((3, cdp.sim_number, 3), float64)
342
343
344 for i in range(cdp.sim_number):
345 euler_to_R_zyz(cdp.eigen_alpha_sim[i], cdp.eigen_beta_sim[i], cdp.eigen_gamma_sim[i], axes_sim[:, i])
346
347
348 axes_sim_pos = axes_sim
349 axes_sim_neg = dot(inv_mat, axes_sim_pos)
350
351
352 print("\nGenerating the axis vectors.")
353 label = ['x', 'y']
354 for i in range(2):
355
356 if sim:
357 axis_sim_pos = axes_sim_pos[:, i]
358 axis_sim_neg = axes_sim_neg[:, i]
359
360
361 res_num = generate_vector_residues(mol=mol, vector=axes_pos[:, i], atom_name='%s-ax'%label[i], res_name_vect='%sAX'%label[i].upper(), sim_vectors=axis_sim_pos, res_num=res_num+1, origin=cdp.pivot, scale=size)
362 if neg_cone:
363 res_num = generate_vector_residues(mol=mol_neg, vector=axes_neg[:, i], atom_name='%s-ax'%label[i], res_name_vect='%sAX'%label[i].upper(), sim_vectors=axis_sim_neg, res_num=res_num, origin=cdp.pivot, scale=size)
364
365
366
367
368
369
370 if cdp.model not in ['rotor', 'free rotor']:
371
372 if cdp.model not in ['iso cone, torsionless', 'iso cone, free rotor']:
373 R = axes
374 else:
375 R = zeros((3, 3), float64)
376 two_vect_to_R(array([0, 0, 1], float64), axis, R)
377
378
379 R_pos = R
380 R_neg = dot(inv_mat, R)
381
382
383 if cdp.model in ['pseudo-ellipse', 'pseudo-ellipse, torsionless', 'pseudo-ellipse, free rotor']:
384 cone = Pseudo_elliptic(cdp.cone_theta_x, cdp.cone_theta_y)
385
386
387 else:
388 cone = Iso_cone(cdp.cone_theta)
389
390
391 create_cone_pdb(mol=mol, cone=cone, start_res=mol.res_num[-1]+1, apex=cdp.pivot, R=R_pos, inc=inc, distribution='regular')
392
393
394 if neg_cone:
395 create_cone_pdb(mol=mol_neg, cone=cone, start_res=mol_neg.res_num[-1]+1, apex=cdp.pivot, R=R_neg, inc=inc, distribution='regular')
396
397
398
399
400
401
402 print("\nGenerating the PDB file.")
403
404
405 pdb_file = open_write_file(file, dir, force=force)
406 structure.write_pdb(pdb_file)
407 pdb_file.close()
408
409
410 - def _domain_to_pdb(self, domain=None, pdb=None):
411 """Match domains to PDB files.
412
413 @keyword domain: The domain to associate the PDB file to.
414 @type domain: str
415 @keyword pdb: The PDB file to associate the domain to.
416 @type pdb: str
417 """
418
419
420 exists = False
421 for i in range(len(cdp.align_tensors)):
422 if hasattr(cdp.align_tensors[i], 'domain') and domain == cdp.align_tensors[i].domain:
423 exists = True
424 if not exists:
425 raise RelaxError("The domain '%s' cannot be found" % domain)
426
427
428 if not hasattr(cdp, 'domain_to_pdb'):
429 cdp.domain_to_pdb = []
430
431
432 if search('.pdb$', pdb):
433 pdb = pdb[:-4]
434
435
436 cdp.domain_to_pdb.append([domain, pdb])
437
438
439 - def _grid_row(self, incs, lower, upper, dist_type=None, end_point=True):
440 """Set up a row of the grid search for a given parameter.
441
442 @param incs: The number of increments.
443 @type incs: int
444 @param lower: The lower bounds.
445 @type lower: float
446 @param upper: The upper bounds.
447 @type upper: float
448 @keyword dist_type: The spacing or distribution type between grid nodes. If None, then a linear grid row is returned. If 'acos', then an inverse cos distribution of points is returned (e.g. for uniform sampling in angular space).
449 @type dist_type: None or str
450 @keyword end_point: A flag which if False will cause the end point to be removed.
451 @type end_point: bool
452 @return: The row of the grid.
453 @rtype: list of float
454 """
455
456
457 row = []
458
459
460 if dist_type == None:
461
462 for i in range(incs):
463
464 row.append(lower + i * (upper - lower) / (incs - 1.0))
465
466
467 elif dist_type == 'acos':
468
469 v = zeros(incs, float64)
470 val = (cos(lower) - cos(upper)) / (incs - 1.0)
471 for i in range(incs):
472 v[-i-1] = cos(upper) + float(i) * val
473
474
475 row = arccos(v)
476
477
478 if not end_point:
479 row = row[:-1]
480
481
482 return list(row)
483
484
486 """Set up the data structures for optimisation using alignment tensors as base data sets.
487
488 @keyword sim_index: The simulation index. This should be None if normal optimisation is
489 desired.
490 @type sim_index: None or int
491 @return: The assembled data structures for using alignment tensors as the base
492 data for optimisation. These include:
493 - full_tensors, the full tensors as concatenated arrays.
494 - red_tensors, the reduced tensors as concatenated arrays.
495 - red_err, the reduced tensor errors as concatenated arrays.
496 @rtype: tuple of 3 numpy nx5D, rank-1 arrays
497 """
498
499
500 if not hasattr(cdp, 'ref_domain'):
501 raise RelaxError("The reference domain has not been set up.")
502 if not hasattr(cdp.align_tensors, 'reduction'):
503 raise RelaxError("The tensor reductions have not been specified.")
504 for i, tensor in self._tensor_loop():
505 if not hasattr(tensor, 'domain'):
506 raise RelaxError("The domain that the '%s' tensor is attached to has not been set" % tensor.name)
507
508
509 n = len(cdp.align_tensors.reduction)
510 full_tensors = zeros(n*5, float64)
511 red_tensors = zeros(n*5, float64)
512 red_err = ones(n*5, float64) * 1e-5
513 full_in_ref_frame = zeros(n, float64)
514
515
516 for i, tensor in self._tensor_loop(red=False):
517
518 full_tensors[5*i + 0] = tensor.Axx
519 full_tensors[5*i + 1] = tensor.Ayy
520 full_tensors[5*i + 2] = tensor.Axy
521 full_tensors[5*i + 3] = tensor.Axz
522 full_tensors[5*i + 4] = tensor.Ayz
523
524
525 if cdp.ref_domain == tensor.domain:
526 full_in_ref_frame[i] = 1
527
528
529 for i, tensor in self._tensor_loop(red=True):
530
531 if sim_index != None:
532 red_tensors[5*i + 0] = tensor.Axx_sim[sim_index]
533 red_tensors[5*i + 1] = tensor.Ayy_sim[sim_index]
534 red_tensors[5*i + 2] = tensor.Axy_sim[sim_index]
535 red_tensors[5*i + 3] = tensor.Axz_sim[sim_index]
536 red_tensors[5*i + 4] = tensor.Ayz_sim[sim_index]
537
538
539 else:
540 red_tensors[5*i + 0] = tensor.Axx
541 red_tensors[5*i + 1] = tensor.Ayy
542 red_tensors[5*i + 2] = tensor.Axy
543 red_tensors[5*i + 3] = tensor.Axz
544 red_tensors[5*i + 4] = tensor.Ayz
545
546
547 if hasattr(tensor, 'Axx_err'):
548 red_err[5*i + 0] = tensor.Axx_err
549 red_err[5*i + 1] = tensor.Ayy_err
550 red_err[5*i + 2] = tensor.Axy_err
551 red_err[5*i + 3] = tensor.Axz_err
552 red_err[5*i + 4] = tensor.Ayz_err
553
554
555 return full_tensors, red_tensors, red_err, full_in_ref_frame
556
557
558 - def _pivot(self, pivot=None):
559 """Set the pivot point for the 2 body motion.
560
561 @param pivot: The pivot point of the two bodies (domains, etc.) in the structural
562 coordinate system.
563 @type pivot: list of num
564 """
565
566
567 pipes.test()
568
569
570 cdp.pivot = pivot
571
572
573 for i in range(3):
574 cdp.pivot[i] = float(cdp.pivot[i])
575
576
577 - def _ref_domain(self, ref=None):
578 """Set the reference domain for the frame order, multi-domain models.
579
580 @param ref: The reference domain.
581 @type ref: str
582 """
583
584
585 pipes.test()
586
587
588 exists = False
589 for tensor_cont in cdp.align_tensors:
590 if hasattr(tensor_cont, 'domain') and tensor_cont.domain == ref:
591 exists = True
592 if not exists:
593 raise RelaxError("The reference domain cannot be found within any of the loaded tensors.")
594
595
596 cdp.ref_domain = ref
597
598
599 if hasattr(cdp, 'model'):
600 self._update_model()
601
602
604 """Select the Frame Order model.
605
606 @param model: The Frame Order model. This can be one of 'pseudo-ellipse', 'pseudo-ellipse, torsionless', 'pseudo-ellipse, free rotor', 'iso cone', 'iso cone, torsionless', 'iso cone, free rotor', 'line', 'line, torsionless', 'line, free rotor', 'rotor', 'rigid', 'free rotor'.
607 @type model: str
608 """
609
610
611 pipes.test()
612
613
614 if not model in ['pseudo-ellipse', 'pseudo-ellipse, torsionless', 'pseudo-ellipse, free rotor', 'iso cone', 'iso cone, torsionless', 'iso cone, free rotor', 'line', 'line, torsionless', 'line, free rotor', 'rotor', 'rigid', 'free rotor']:
615 raise RelaxError("The model name " + repr(model) + " is invalid.")
616
617
618 cdp.model = model
619
620
621 cdp.params = []
622
623
624 self._update_model()
625
626
628 """Store the back-calculated reduced alignment tensors.
629
630 @param target_fn: The frame-order target function class.
631 @type target_fn: class instance
632 """
633
634
635 for i, tensor in self._tensor_loop(red=True):
636
637 name = tensor.name + ' bc'
638
639
640 align_tensor.init(tensor=name, params=(target_fn.red_tensors_bc[5*i + 0], target_fn.red_tensors_bc[5*i + 1], target_fn.red_tensors_bc[5*i + 2], target_fn.red_tensors_bc[5*i + 3], target_fn.red_tensors_bc[5*i + 4]), param_types=2)
641
642
644 """Generator method for looping over the full or reduced tensors.
645
646 @keyword red: A flag which if True causes the reduced tensors to be returned, and if False
647 the full tensors are returned.
648 @type red: bool
649 @return: The tensor index and the tensor.
650 @rtype: (int, AlignTensorData instance)
651 """
652
653
654 n = len(cdp.align_tensors.reduction)
655
656
657 data = cdp.align_tensors
658 list = data.reduction
659
660
661 if red:
662 index = 1
663 else:
664 index = 0
665
666
667 for i in range(n):
668 yield i, data[list[i][index]]
669
670
672 """Update the model parameters as necessary."""
673
674
675 if not hasattr(cdp, 'params'):
676 cdp.params = []
677
678
679 if not len(cdp.params):
680
681 if cdp.model not in ['free rotor', 'iso cone, torsionless', 'iso cone, free rotor']:
682 cdp.params.append('ave_pos_alpha')
683 cdp.params.append('ave_pos_beta')
684 cdp.params.append('ave_pos_gamma')
685
686
687 if cdp.model in ['iso cone', 'pseudo-ellipse', 'pseudo-ellipse, torsionless', 'pseudo-ellipse, free rotor']:
688 cdp.params.append('eigen_alpha')
689 cdp.params.append('eigen_beta')
690 cdp.params.append('eigen_gamma')
691
692
693 elif cdp.model in ['free rotor', 'iso cone, torsionless', 'iso cone, free rotor', 'rotor']:
694 cdp.params.append('axis_theta')
695 cdp.params.append('axis_phi')
696
697
698 if cdp.model in ['pseudo-ellipse', 'pseudo-ellipse, torsionless', 'pseudo-ellipse, free rotor']:
699 cdp.params.append('cone_theta_x')
700 cdp.params.append('cone_theta_y')
701
702
703 elif cdp.model in ['iso cone', 'iso cone, torsionless']:
704 cdp.params.append('cone_theta')
705 elif cdp.model in ['iso cone, free rotor']:
706 cdp.params.append('cone_s1')
707
708
709 if cdp.model in ['rotor', 'line', 'iso cone', 'pseudo-ellipse']:
710 cdp.params.append('cone_sigma_max')
711
712
713 for param in cdp.params:
714 if not hasattr(cdp, param):
715 setattr(cdp, param, 0.0)
716
717
719 """Unpack and store the Frame Order optimisation results.
720
721 @param results: The results tuple returned by the minfx generic_minimise() function.
722 @type results: tuple
723 @param sim_index: The index of the simulation to optimise. This should be None for normal
724 optimisation.
725 @type sim_index: None or int
726 """
727
728
729 if len(results) == 4:
730 param_vector, func, iter_count, warning = results
731 f_count = iter_count
732 g_count = 0.0
733 h_count = 0.0
734 else:
735 param_vector, func, iter_count, f_count, g_count, h_count, warning = results
736
737
738 if isInf(func):
739 raise RelaxInfError('chi-squared')
740
741
742 if isNaN(func):
743 raise RelaxNaNError('chi-squared')
744
745
746 ave_pos_alpha, ave_pos_beta, ave_pos_gamma = None, None, None
747 eigen_alpha, eigen_beta, eigen_gamma = None, None, None
748 axis_theta, axis_phi = None, None
749 cone_theta_x, cone_theta_y = None, None
750 cone_theta = None
751 cone_s1 = None
752 cone_sigma_max = None
753 if cdp.model == 'pseudo-ellipse':
754 ave_pos_alpha, ave_pos_beta, ave_pos_gamma, eigen_alpha, eigen_beta, eigen_gamma, cone_theta_x, cone_theta_y, cone_sigma_max = param_vector
755 elif cdp.model in ['pseudo-ellipse, torsionless', 'pseudo-ellipse, free rotor']:
756 ave_pos_alpha, ave_pos_beta, ave_pos_gamma, eigen_alpha, eigen_beta, eigen_gamma, cone_theta_x, cone_theta_y = param_vector
757 elif cdp.model == 'iso cone':
758 ave_pos_alpha, ave_pos_beta, ave_pos_gamma, eigen_alpha, eigen_beta, eigen_gamma, cone_theta, cone_sigma_max = param_vector
759 elif cdp.model in ['iso cone, torsionless']:
760 ave_pos_beta, ave_pos_gamma, axis_theta, axis_phi, cone_theta = param_vector
761 ave_pos_alpha = None
762 elif cdp.model in ['iso cone, free rotor']:
763 ave_pos_beta, ave_pos_gamma, axis_theta, axis_phi, cone_s1 = param_vector
764 ave_pos_alpha = None
765 elif cdp.model == 'line':
766 ave_pos_alpha, ave_pos_beta, ave_pos_gamma, eigen_alpha, eigen_beta, eigen_gamma, cone_theta_x, cone_sigma_max = param_vector
767 elif cdp.model in ['line, torsionless', 'line, free rotor']:
768 ave_pos_alpha, ave_pos_beta, ave_pos_gamma, eigen_alpha, eigen_beta, eigen_gamma, cone_theta_x, cone_sigma_max = param_vector
769 elif cdp.model in ['rotor']:
770 ave_pos_alpha, ave_pos_beta, ave_pos_gamma, axis_theta, axis_phi, cone_sigma_max = param_vector
771 elif cdp.model in ['free rotor']:
772 ave_pos_beta, ave_pos_gamma, axis_theta, axis_phi = param_vector
773 ave_pos_alpha = None
774 elif cdp.model == 'rigid':
775 ave_pos_alpha, ave_pos_beta, ave_pos_gamma = param_vector
776
777
778 if sim_index != None:
779
780 if ave_pos_alpha != None:
781 cdp.ave_pos_alpha_sim[sim_index] = wrap_angles(ave_pos_alpha, 0.0, 2.0*pi)
782 if ave_pos_beta != None:
783 cdp.ave_pos_beta_sim[sim_index] = wrap_angles(ave_pos_beta, 0.0, 2.0*pi)
784 if ave_pos_gamma != None:
785 cdp.ave_pos_gamma_sim[sim_index] = wrap_angles(ave_pos_gamma, 0.0, 2.0*pi)
786
787
788 if eigen_alpha != None:
789 cdp.eigen_alpha_sim[sim_index] = wrap_angles(eigen_alpha, 0.0, 2.0*pi)
790 if eigen_beta != None:
791 cdp.eigen_beta_sim[sim_index] = wrap_angles(eigen_beta, 0.0, 2.0*pi)
792 if eigen_gamma != None:
793 cdp.eigen_gamma_sim[sim_index] = wrap_angles(eigen_gamma, 0.0, 2.0*pi)
794 if axis_theta != None:
795 cdp.axis_theta_sim[sim_index] = axis_theta
796 if axis_phi != None:
797 cdp.axis_phi_sim[sim_index] = axis_phi
798
799
800 if cone_theta != None:
801 cdp.cone_theta_sim[sim_index] = cone_theta
802 if cone_s1 != None:
803 cdp.cone_s1_sim[sim_index] = cone_s1
804 cdp.cone_theta_sim[sim_index] = order_parameters.iso_cone_S_to_theta(cone_s1)
805 if cone_theta_x != None:
806 cdp.cone_theta_x_sim[sim_index] = cone_theta_x
807 if cone_theta_y != None:
808 cdp.cone_theta_y_sim[sim_index] = cone_theta_y
809 if cone_sigma_max != None:
810 cdp.cone_sigma_max_sim[sim_index] = abs(cone_sigma_max)
811
812
813 cdp.chi2_sim[sim_index] = func
814 cdp.iter_sim[sim_index] = iter_count
815 cdp.f_count_sim[sim_index] = f_count
816 cdp.g_count_sim[sim_index] = g_count
817 cdp.h_count_sim[sim_index] = h_count
818 cdp.warning_sim[sim_index] = warning
819
820
821 else:
822
823 if ave_pos_alpha != None:
824 cdp.ave_pos_alpha = wrap_angles(ave_pos_alpha, 0.0, 2.0*pi)
825 if ave_pos_beta != None:
826 cdp.ave_pos_beta = wrap_angles(ave_pos_beta, 0.0, 2.0*pi)
827 if ave_pos_gamma != None:
828 cdp.ave_pos_gamma = wrap_angles(ave_pos_gamma, 0.0, 2.0*pi)
829
830
831 if eigen_alpha != None:
832 cdp.eigen_alpha = wrap_angles(eigen_alpha, 0.0, 2.0*pi)
833 if eigen_beta != None:
834 cdp.eigen_beta = wrap_angles(eigen_beta, 0.0, 2.0*pi)
835 if eigen_gamma != None:
836 cdp.eigen_gamma = wrap_angles(eigen_gamma, 0.0, 2.0*pi)
837 if axis_theta != None:
838 cdp.axis_theta = axis_theta
839 if axis_phi != None:
840 cdp.axis_phi = axis_phi
841
842
843 if cone_theta != None:
844 cdp.cone_theta = cone_theta
845 if cone_s1 != None:
846 cdp.cone_s1 = cone_s1
847 cdp.cone_theta = order_parameters.iso_cone_S_to_theta(cone_s1)
848 if cone_theta_x != None:
849 cdp.cone_theta_x = cone_theta_x
850 if cone_theta_y != None:
851 cdp.cone_theta_y = cone_theta_y
852 if cone_sigma_max != None:
853 cdp.cone_sigma_max = abs(cone_sigma_max)
854
855
856 cdp.chi2 = func
857 cdp.iter = iter_count
858 cdp.f_count = f_count
859 cdp.g_count = g_count
860 cdp.h_count = h_count
861 cdp.warning = warning
862
863
865 """Generator method for looping nothing.
866
867 The loop essentially consists of a single element.
868
869 @return: Nothing.
870 @rtype: None
871 """
872
873 yield None
874
875
876 - def calculate(self, spin_id=None, verbosity=1, sim_index=None):
877 """Calculate the chi-squared value for the current parameter values.
878
879 @keyword spin_id: The spin identification string (unused).
880 @type spin_id: None
881 @keyword verbosity: The amount of information to print. The higher the value, the greater the verbosity.
882 @type verbosity: int
883 @keyword sim_index: The optional MC simulation index (unused).
884 @type sim_index: None or int
885 """
886
887
888 param_vector = self._assemble_param_vector()
889
890
891 full_tensors, red_tensors, red_tensor_err, full_in_ref_frame = self._minimise_setup_tensors()
892
893
894 target = frame_order.Frame_order(model=cdp.model, full_tensors=full_tensors, red_tensors=red_tensors, red_errors=red_tensor_err, full_in_ref_frame=full_in_ref_frame)
895
896
897 chi2 = target.func(param_vector)
898
899
900 cdp.chi2 = chi2
901
902
903 self._store_bc_tensors(target)
904
905
907 """Create the Monte Carlo data by back calculating the reduced tensor data.
908
909 @keyword data_id: Unused.
910 @type data_id: None
911 @return: The Monte Carlo simulation data.
912 @rtype: list of floats
913 """
914
915
916 red_tensors_bc = self._back_calc()
917
918
919 return red_tensors_bc
920
921
923 """Return a vector of parameter names.
924
925 @keyword model_info: The model index from model_info().
926 @type model_info: int
927 @return: The vector of parameter names.
928 @rtype: list of str
929 """
930
931
932 self._update_model()
933
934
935 return cdp.params
936
937
939 """Return a vector of parameter values.
940
941 @keyword model_info: The model index from model_info(). This is zero for the global models or equal to the global spin index (which covers the molecule, residue, and spin indices).
942 @type model_info: int
943 @keyword sim_index: The Monte Carlo simulation index.
944 @type sim_index: int
945 @return: The vector of parameter values.
946 @rtype: list of str
947 """
948
949
950 return self._assemble_param_vector(sim_index=sim_index)
951
952
953 - def grid_search(self, lower=None, upper=None, inc=None, constraints=False, verbosity=0, sim_index=None):
954 """Perform a grid search.
955
956 @keyword lower: The lower bounds of the grid search which must be equal to the
957 number of parameters in the model.
958 @type lower: list of float
959 @keyword upper: The upper bounds of the grid search which must be equal to the
960 number of parameters in the model.
961 @type upper: list of float
962 @keyword inc: The increments for each dimension of the space for the grid search.
963 The number of elements in the array must equal to the number of
964 parameters in the model.
965 @type inc: int or list of int
966 @keyword constraints: If True, constraints are applied during the grid search (eliminating
967 parts of the grid). If False, no constraints are used.
968 @type constraints: bool
969 @keyword verbosity: A flag specifying the amount of information to print. The higher
970 the value, the greater the verbosity.
971 @type verbosity: int
972 @keyword sim_index: The Monte Carlo simulation index.
973 @type sim_index: None or int
974 """
975
976
977 if not hasattr(cdp, 'model'):
978 raise RelaxNoModelError('Frame Order')
979
980
981 n = len(cdp.params)
982
983
984 if isinstance(inc, int):
985 incs = [inc]*n
986 else:
987 incs = inc
988
989
990 if len(incs) != n:
991 raise RelaxError("The size of the increment list %s does not match the number of parameters in %s." % (incs, cdp.params))
992
993
994 lower_list = lower
995 upper_list = upper
996 grid = []
997 """This structure is a list of lists. The first dimension corresponds to the model
998 parameter. The second dimension are the grid node positions."""
999
1000
1001 for i in range(n):
1002
1003 if incs[i] == None:
1004 grid.append(None)
1005 continue
1006
1007
1008 dist_type = None
1009 end_point = True
1010
1011
1012 if cdp.params[i] in ['ave_pos_alpha', 'ave_pos_gamma', 'eigen_alpha', 'eigen_gamma', 'axis_phi']:
1013 lower = 0.0
1014 upper = 2*pi * (1.0 - 1.0/incs[i])
1015
1016
1017 if cdp.params[i] in ['ave_pos_beta', 'eigen_beta', 'axis_theta']:
1018
1019 if not isinstance(inc, list):
1020 incs[i] = int(incs[i] / 2) + 1
1021
1022
1023 dist_type = 'acos'
1024 end_point = False
1025
1026
1027 lower = 0.0
1028 upper = pi
1029
1030
1031 if cdp.params[i] == 'cone_s1':
1032 lower = -0.125
1033 upper = 1.0
1034
1035
1036 if cdp.params[i] in ['cone_theta', 'cone_theta_x', 'cone_theta_y', 'cone_sigma_max']:
1037 lower = pi * (1.0/incs[i])
1038 upper = pi * (1.0 - 1.0/incs[i])
1039
1040
1041 if lower_list:
1042 lower = lower_list[i]
1043 if upper_list:
1044 upper = upper_list[i]
1045
1046
1047 row = self._grid_row(incs[i], lower, upper, dist_type=dist_type, end_point=end_point)
1048 grid.append(row)
1049
1050
1051 if not end_point:
1052 incs[i] -= 1
1053
1054
1055 total_pts = 1
1056 for i in range(n):
1057
1058 if grid[i] == None:
1059 continue
1060
1061 total_pts = total_pts * len(grid[i])
1062
1063
1064 max_pts = 50e6
1065 if total_pts > max_pts:
1066 raise RelaxError("The total number of grid points '%s' exceeds the maximum of '%s'." % (total_pts, int(max_pts)))
1067
1068
1069 pts = zeros((total_pts, n), float64)
1070 indices = zeros(n, int)
1071 for i in range(total_pts):
1072
1073 for j in range(n):
1074
1075 if grid[j] == None:
1076
1077 pts[i, j] = getattr(cdp, cdp.params[j])
1078
1079
1080 else:
1081 pts[i, j] = grid[j][indices[j]]
1082
1083
1084 for j in range(n):
1085 if incs[j] != None and indices[j] < incs[j]-1:
1086 indices[j] += 1
1087 break
1088 else:
1089 indices[j] = 0
1090
1091
1092 self.minimise(min_algor='grid', min_options=pts, constraints=constraints, verbosity=verbosity, sim_index=sim_index)
1093
1094
1096 """State that the parameter is not spin specific.
1097
1098 @param name: The name of the parameter.
1099 @type name: str
1100 @return: False.
1101 @rtype: bool
1102 """
1103
1104
1105 return False
1106
1107
1109 """Create bounds for the OpenDX mapping function.
1110
1111 @param param: The name of the parameter to return the lower and upper bounds of.
1112 @type param: str
1113 @param spin_id: The spin identification string (unused).
1114 @type spin_id: None
1115 @return: The upper and lower bounds of the parameter.
1116 @rtype: list of float
1117 """
1118
1119
1120 if param in ['ave_pos_alpha', 'ave_pos_beta', 'ave_pos_gamma']:
1121 return [0.0, 2*pi]
1122
1123
1124 if param == 'axis_theta':
1125 return [0.0, pi]
1126
1127
1128 if param == 'axis_phi':
1129 return [0.0, 2*pi]
1130
1131
1132 if param == 'cone_theta':
1133 return [0.0, pi]
1134
1135
1136 - 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):
1137 """Minimisation function.
1138
1139 @param min_algor: The minimisation algorithm to use.
1140 @type min_algor: str
1141 @param min_options: An array of options to be used by the minimisation algorithm.
1142 @type min_options: array of str
1143 @param func_tol: The function tolerance which, when reached, terminates optimisation.
1144 Setting this to None turns of the check.
1145 @type func_tol: None or float
1146 @param grad_tol: The gradient tolerance which, when reached, terminates optimisation.
1147 Setting this to None turns of the check.
1148 @type grad_tol: None or float
1149 @param max_iterations: The maximum number of iterations for the algorithm.
1150 @type max_iterations: int
1151 @param constraints: If True, constraints are used during optimisation.
1152 @type constraints: bool
1153 @param scaling: If True, diagonal scaling is enabled during optimisation to allow
1154 the problem to be better conditioned.
1155 @type scaling: bool
1156 @param verbosity: A flag specifying the amount of information to print. The higher
1157 the value, the greater the verbosity.
1158 @type verbosity: int
1159 @param sim_index: The index of the simulation to optimise. This should be None if
1160 normal optimisation is desired.
1161 @type sim_index: None or int
1162 @keyword lower: The 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.
1163 @type lower: array of numbers
1164 @keyword upper: The 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.
1165 @type upper: array of numbers
1166 @keyword inc: The 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.
1167 @type inc: array of int
1168 """
1169
1170
1171 if constraints:
1172
1173 constraints = False
1174
1175
1176 if not search('^[Gg]rid', min_algor):
1177 min_algor = min_options[0]
1178 min_options = min_options[1:]
1179
1180
1181 warn(RelaxWarning("Constraints are as of yet not implemented - turning this option off."))
1182
1183
1184 lower, upper = self._assemble_limit_arrays()
1185
1186
1187 param_vector = self._assemble_param_vector()
1188
1189
1190 full_tensors, red_tensors, red_tensor_err, full_in_ref_frame = self._minimise_setup_tensors(sim_index)
1191
1192
1193 target = frame_order.Frame_order(model=cdp.model, full_tensors=full_tensors, red_tensors=red_tensors, red_errors=red_tensor_err, full_in_ref_frame=full_in_ref_frame)
1194
1195
1196 if search('^[Gg]rid', min_algor):
1197 results = grid_point_array(func=target.func, args=(), points=min_options, verbosity=verbosity)
1198
1199
1200 else:
1201 results = generic_minimise(func=target.func, args=(), x0=param_vector, min_algor=min_algor, min_options=min_options, func_tol=func_tol, grad_tol=grad_tol, maxiter=max_iterations, l=lower, u=upper, full_output=True, print_flag=verbosity)
1202
1203
1204 self._unpack_opt_results(results, sim_index)
1205
1206
1207 self._store_bc_tensors(target)
1208
1209
1210
1212 """Dummy generator method.
1213
1214 In this case only a single model per spin system is assumed. Hence the yielded data is the
1215 spin container object.
1216
1217
1218 @return: Information about the model which for this analysis is the spin container.
1219 @rtype: SpinContainer instance
1220 """
1221
1222
1223 yield None
1224
1225
1227 """Return the k, n, and chi2 model statistics.
1228
1229 k - number of parameters.
1230 n - number of data points.
1231 chi2 - the chi-squared value.
1232
1233
1234 @keyword model_info: Unused.
1235 @type model_info: None
1236 @keyword spin_id: The spin identification string (unused).
1237 @type spin_id: None
1238 @keyword global_stats: Unused.
1239 @type global_stats: None
1240 @return: The optimisation statistics, in tuple format, of the number of
1241 parameters (k), the number of data points (n), and the chi-squared
1242 value (chi2).
1243 @rtype: tuple of (int, int, float)
1244 """
1245
1246
1247 param_vector = self._assemble_param_vector()
1248 k = len(param_vector)
1249
1250
1251 n = len(cdp.align_tensors.reduction)
1252
1253
1254 if not hasattr(cdp, 'chi2'):
1255 raise RelaxError("Statistics are not available, most likely because the model has not been optimised.")
1256 chi2 = cdp.chi2
1257
1258
1259 return k, n, chi2
1260
1261
1263 """Return the alignment tensor error structure.
1264
1265 @param data_id: The information yielded by the base_data_loop() generator method.
1266 @type data_id: None
1267 @return: The array of tensor error values.
1268 @rtype: list of float
1269 """
1270
1271
1272 full_tensors, red_tensors, red_tensor_err, full_in_ref_frame = self._minimise_setup_tensors()
1273
1274
1275 return red_tensor_err
1276
1277
1279 """Return a string representing the parameters units.
1280
1281 @param param: The name of the parameter to return the units string for.
1282 @type param: str
1283 @return: The parameter units string.
1284 @rtype: str
1285 """
1286
1287
1288 if param in ['ave_pos_alpha', 'ave_pos_beta', 'ave_pos_gamma']:
1289 return 'rad'
1290
1291
1292 if param in ['eigen_alpha', 'eigen_beta', 'eigen_gamma', 'axis_theta', 'axis_phi']:
1293 return 'rad'
1294
1295
1296 if param in ['cone_theta_x', 'cone_theta_y', 'cone_sigma_max', 'cone_theta']:
1297 return 'rad'
1298
1299
1300 - def set_error(self, model_info, index, error):
1301 """Set the parameter errors.
1302
1303 @param model_info: The model information originating from model_loop() (unused).
1304 @type model_info: None
1305 @param index: The index of the parameter to set the errors for.
1306 @type index: int
1307 @param error: The error value.
1308 @type error: float
1309 """
1310
1311
1312 inc = 0
1313
1314
1315 for param in self.data_names(set='params'):
1316
1317 if index == inc:
1318 setattr(cdp, param + "_err", error)
1319
1320
1321 inc = inc + 1
1322
1323
1324 if cdp.model == 'iso cone, free rotor' and inc == index:
1325 setattr(cdp, 'cone_theta_err', error)
1326
1327
1328
1330 """Set the simulation selection flag for the spin.
1331
1332 @param model_info: The model information originating from model_loop() (unused).
1333 @type model_info: None
1334 @param select_sim: The selection flag for the simulations.
1335 @type select_sim: bool
1336 """
1337
1338
1339 cdp.select_sim = deepcopy(select_sim)
1340
1341
1343 """Initialise the Monte Carlo parameter values."""
1344
1345
1346 param_names = self.data_names(set='params')
1347
1348
1349 if cdp.model == 'iso cone, free rotor':
1350 param_names.append('cone_theta')
1351
1352
1353 min_names = self.data_names(set='min')
1354
1355
1356
1357
1358
1359
1360 for object_name in param_names:
1361
1362 sim_object_name = object_name + '_sim'
1363
1364
1365 if hasattr(cdp, sim_object_name):
1366 raise RelaxError("Monte Carlo parameter values have already been set.")
1367
1368
1369
1370
1371
1372
1373 for object_name in param_names:
1374
1375 if not hasattr(cdp, object_name):
1376 continue
1377
1378
1379 sim_object_name = object_name + '_sim'
1380
1381
1382 setattr(cdp, sim_object_name, [])
1383
1384
1385 sim_object = getattr(cdp, sim_object_name)
1386
1387
1388 for j in range(cdp.sim_number):
1389
1390 sim_object.append(deepcopy(getattr(cdp, object_name)))
1391
1392
1393 for object_name in min_names:
1394
1395 sim_object_name = object_name + '_sim'
1396
1397
1398 setattr(cdp, sim_object_name, [])
1399
1400
1401 sim_object = getattr(cdp, sim_object_name)
1402
1403
1404 for j in range(cdp.sim_number):
1405
1406 sim_object.append(deepcopy(getattr(cdp, object_name)))
1407
1408
1410 """Pack the Monte Carlo simulation data.
1411
1412 @param data_id: The spin identification string, as yielded by the base_data_loop() generator method.
1413 @type data_id: None
1414 @param sim_data: The Monte Carlo simulation data.
1415 @type sim_data: list of float
1416 """
1417
1418
1419 sim_data = transpose(sim_data)
1420
1421
1422 for i, tensor in self._tensor_loop(red=True):
1423
1424 tensor.set_sim_num(cdp.sim_number)
1425
1426
1427 for j in range(cdp.sim_number):
1428
1429 tensor.set(param='Axx', value=sim_data[5*i + 0][j], category='sim', sim_index=j)
1430 tensor.set(param='Ayy', value=sim_data[5*i + 1][j], category='sim', sim_index=j)
1431 tensor.set(param='Axy', value=sim_data[5*i + 2][j], category='sim', sim_index=j)
1432 tensor.set(param='Axz', value=sim_data[5*i + 3][j], category='sim', sim_index=j)
1433 tensor.set(param='Ayz', value=sim_data[5*i + 4][j], category='sim', sim_index=j)
1434
1435
1437 """Return the array of simulation parameter values.
1438
1439 @param model_info: The model information originating from model_loop() (unused).
1440 @type model_info: None
1441 @param index: The index of the parameter to return the array of values for.
1442 @type index: int
1443 """
1444
1445
1446 inc = 0
1447
1448
1449 param_names = self.data_names(set='params')
1450
1451
1452 for param in param_names:
1453
1454 if not hasattr(cdp, param):
1455 continue
1456
1457
1458 if index == inc:
1459 return getattr(cdp, param + "_sim")
1460
1461
1462 inc = inc + 1
1463
1464
1465 if cdp.model == 'iso cone, free rotor' and inc == index:
1466 return getattr(cdp, 'cone_theta_sim')
1467
1468
1470 """Return the array of selected simulation flags for the spin.
1471
1472 @param model_info: The model information originating from model_loop() (unused).
1473 @type model_info: None
1474 @return: The array of selected simulation flags.
1475 @rtype: list of int
1476 """
1477
1478
1479 return cdp.select_sim
1480