1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23 """Module for handling the frame order model parameters."""
24
25
26 from math import pi
27 from numpy import array, float64, identity, zeros
28
29
30 from specific_analyses.frame_order.data import pivot_fixed, translation_fixed
31
32
34 """Assemble and return the parameter vector.
35
36 @return: The parameter vector.
37 @rtype: numpy rank-1 array
38 @keyword sim_index: The Monte Carlo simulation index.
39 @type sim_index: int
40 """
41
42
43 param_vect = []
44
45
46 ext = ''
47 if sim_index != None:
48 ext = '_sim'
49
50
51 for param_name in cdp.params:
52
53 obj = getattr(cdp, param_name+ext)
54
55
56 if sim_index == None:
57 param_vect.append(obj)
58 else:
59 param_vect.append(obj[sim_index])
60
61
62 return array(param_vect, float64)
63
64
66 """Create and return the scaling matrix.
67
68 @keyword scaling: If False, then the identity matrix will be returned.
69 @type scaling: bool
70 @return: The square and diagonal scaling matrix.
71 @rtype: numpy rank-2 array
72 """
73
74
75 scaling_matrix = identity(param_num(), float64)
76
77
78 if not scaling:
79 return scaling_matrix
80
81
82 if not pivot_fixed():
83 for i in range(3):
84 scaling_matrix[i, i] = 1e2
85
86
87 return scaling_matrix
88
89
91 """Create the linear constraint matrices A and b.
92
93 Standard notation
94 =================
95
96 The parameter constraints for the motional amplitude parameters are::
97
98 0 <= theta <= pi,
99 0 <= theta_x <= theta_y <= pi,
100 -0.125 <= S <= 1,
101 0 <= sigma_max <= pi,
102
103 The pivot point parameter, domain position parameters, and eigenframe parameters are unconstrained.
104
105
106 Matrix notation
107 ===============
108
109 In the notation A.x >= b, where A is an matrix of coefficients, x is an array of parameter values, and b is a vector of scalars, these inequality constraints are::
110
111 | 1 0 0 0 0 | | 0 |
112 | | | |
113 |-1 0 0 0 0 | | -pi |
114 | | | |
115 | 0 1 0 0 0 | | 0 |
116 | | | |
117 | 0 -1 0 0 0 | | theta | | -pi |
118 | | | | | |
119 | 0 -1 1 0 0 | | theta_x | | 0 |
120 | | | | | |
121 | 0 0 1 0 0 | . | theta_y | >= | 0 |
122 | | | | | |
123 | 0 0 -1 0 0 | | S | | -pi |
124 | | | | | |
125 | 0 0 0 1 0 | | sigma_max | | -0.125 |
126 | | | |
127 | 0 0 0 -1 0 | | -1 |
128 | | | |
129 | 0 0 0 0 1 | | 0 |
130 | | | |
131 | 0 0 0 0 -1 | | -pi |
132
133
134 @keyword scaling_matrix: The diagonal, square scaling matrix.
135 @type scaling_matrix: numpy rank-2 square matrix
136 @return: The matrices A and b.
137 @rtype: numpy rank-2 NxM array, numpy rank-1 N array
138 """
139
140
141 A = []
142 b = []
143 n = param_num()
144 zero_array = zeros(n, float64)
145 i = 0
146 j = 0
147
148
149 for i in range(n):
150
151 if cdp.params[i] in ['cone_theta', 'cone_theta_x', 'cone_theta_y', 'cone_sigma_max']:
152
153 A.append(zero_array * 0.0)
154 A.append(zero_array * 0.0)
155 A[j][i] = 1.0
156 A[j+1][i] = -1.0
157 b.append(0.0)
158 b.append(-pi / scaling_matrix[i, i])
159 j = j + 2
160
161
162 if cdp.params[i] == 'cone_theta_y':
163 for m in range(n):
164 if cdp.params[m] == 'cone_theta_x':
165 A.append(zero_array * 0.0)
166 A[j][i] = 1.0
167 A[j][m] = -1.0
168 b.append(0.0)
169 j = j + 1
170
171
172
173 if cdp.params[i] == 'cone_s1':
174
175 A.append(zero_array * 0.0)
176 A.append(zero_array * 0.0)
177 A[j][i] = 1.0
178 A[j+1][i] = -1.0
179 b.append(-0.125 / scaling_matrix[i, i])
180 b.append(-1 / scaling_matrix[i, i])
181 j = j + 2
182
183
184 A = array(A, float64)
185 b = array(b, float64)
186
187
188 return A, b
189
190
192 """Determine the number of parameters in the model.
193
194 @return: The number of model parameters.
195 @rtype: int
196 """
197
198
199 update_model()
200
201
202 return len(cdp.params)
203
204
206 """Update the model parameters as necessary."""
207
208
209 cdp.params = []
210
211
212 if not pivot_fixed():
213 cdp.params.append('pivot_x')
214 cdp.params.append('pivot_y')
215 cdp.params.append('pivot_z')
216
217
218 if cdp.model == 'double rotor':
219 cdp.params.append('pivot_x_2')
220 cdp.params.append('pivot_y_2')
221 cdp.params.append('pivot_z_2')
222
223
224 if not translation_fixed():
225 cdp.params.append('ave_pos_x')
226 cdp.params.append('ave_pos_y')
227 cdp.params.append('ave_pos_z')
228
229
230 if cdp.model not in ['free rotor', 'iso cone, free rotor']:
231 cdp.params.append('ave_pos_alpha')
232 cdp.params.append('ave_pos_beta')
233 cdp.params.append('ave_pos_gamma')
234
235
236 if cdp.model in ['pseudo-ellipse', 'pseudo-ellipse, torsionless', 'pseudo-ellipse, free rotor']:
237 cdp.params.append('eigen_alpha')
238 cdp.params.append('eigen_beta')
239 cdp.params.append('eigen_gamma')
240
241
242 if cdp.model in ['iso cone', 'free rotor', 'iso cone, torsionless', 'iso cone, free rotor', 'double rotor']:
243 cdp.params.append('axis_theta')
244 cdp.params.append('axis_phi')
245
246
247 if cdp.model in ['double rotor']:
248 cdp.params.append('axis_theta_2')
249 cdp.params.append('axis_phi_2')
250
251
252 if cdp.model in ['rotor']:
253 cdp.params.append('axis_alpha')
254
255
256 if cdp.model in ['pseudo-ellipse', 'pseudo-ellipse, torsionless', 'pseudo-ellipse, free rotor']:
257 cdp.params.append('cone_theta_x')
258 cdp.params.append('cone_theta_y')
259
260
261 if cdp.model in ['iso cone', 'iso cone, torsionless']:
262 cdp.params.append('cone_theta')
263 if cdp.model in ['iso cone, free rotor']:
264 cdp.params.append('cone_s1')
265
266
267 if cdp.model in ['double rotor', 'rotor', 'line', 'iso cone', 'pseudo-ellipse']:
268 cdp.params.append('cone_sigma_max')
269
270
271 if cdp.model in ['double rotor']:
272 cdp.params.append('cone_sigma_max_2')
273
274
275 for param in cdp.params:
276 if not param in ['pivot_x', 'pivot_y', 'pivot_z', 'pivot_x_2', 'pivot_y_2', 'pivot_z_2'] and not hasattr(cdp, param):
277 setattr(cdp, param, 0.0)
278