1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23 """The N-state model or structural ensemble analysis optimisation functions."""
24
25
26 from numpy import array, dot, float64, ones, zeros
27 from numpy.linalg import inv
28
29
30 from lib.errors import RelaxError, RelaxNoModelError
31 from pipe_control import align_tensor
32 from pipe_control.align_tensor import opt_uses_align_data, opt_uses_tensor
33 from pipe_control.interatomic import interatomic_loop
34 from pipe_control.mol_res_spin import return_spin, spin_loop
35 from pipe_control.pcs import return_pcs_data
36 from pipe_control.rdc import check_rdcs, return_rdc_data
37 from specific_analyses.n_state_model.data import base_data_types, tensor_loop
38 from specific_analyses.n_state_model.parameters import assemble_param_vector, assemble_scaling_matrix, update_model
39 from target_functions.n_state_model import N_state_opt
40
41
43 """Extract and unpack the back calculated data.
44
45 @param model: The instantiated class containing the target function.
46 @type model: class instance
47 """
48
49
50 if not hasattr(cdp, 'align_tensors'):
51 return
52
53
54 align_index = 0
55 for i in range(len(cdp.align_ids)):
56
57 if not opt_uses_tensor(cdp.align_tensors[i]):
58 continue
59
60
61 align_id = cdp.align_ids[i]
62
63
64 rdc_flag = False
65 if hasattr(cdp, 'rdc_ids') and align_id in cdp.rdc_ids:
66 rdc_flag = True
67 pcs_flag = False
68 if hasattr(cdp, 'pcs_ids') and align_id in cdp.pcs_ids:
69 pcs_flag = True
70
71
72 pcs_index = 0
73 for spin in spin_loop():
74
75 if not spin.select:
76 continue
77
78
79 if pcs_flag and hasattr(spin, 'pcs'):
80
81 if not hasattr(spin, 'pcs_bc'):
82 spin.pcs_bc = {}
83
84
85 spin.pcs_bc[align_id] = model.deltaij_theta[align_index, pcs_index] * 1e6
86
87
88 pcs_index = pcs_index + 1
89
90
91 rdc_index = 0
92 for interatom in interatomic_loop():
93
94 spin1 = return_spin(interatom.spin_id1)
95 spin2 = return_spin(interatom.spin_id2)
96
97
98 if not check_rdcs(interatom):
99 continue
100
101
102 if rdc_flag and hasattr(interatom, 'rdc'):
103
104 if not hasattr(interatom, 'rdc_bc'):
105 interatom.rdc_bc = {}
106
107
108 interatom.rdc_bc[align_id] = model.rdc_theta[align_index, rdc_index]
109
110
111 rdc_index = rdc_index + 1
112
113
114 align_index += 1
115
116
118 """Set up the atomic position data structures for optimisation using PCSs and PREs as base data sets.
119
120 @keyword sim_index: The index of the simulation to optimise. This should be None if normal optimisation is desired.
121 @type sim_index: None or int
122 @return: The atomic positions (the first index is the spins, the second is the structures, and the third is the atomic coordinates) and the paramagnetic centre.
123 @rtype: numpy rank-3 array, numpy rank-1 array.
124 """
125
126
127 atomic_pos = []
128
129
130 for spin in spin_loop():
131
132 if not spin.select:
133 continue
134
135
136 if not hasattr(spin, 'pcs') and not hasattr(spin, 'pre'):
137 continue
138
139
140 if type(spin.pos[0]) in [float, float64]:
141 atomic_pos.append([spin.pos])
142 else:
143 atomic_pos.append(spin.pos)
144
145
146 atomic_pos = array(atomic_pos, float64)
147
148
149 if not hasattr(cdp, 'paramagnetic_centre'):
150 paramag_centre = zeros(3, float64)
151 elif sim_index != None and not cdp.paramag_centre_fixed:
152 if not hasattr(cdp, 'paramagnetic_centre_sim') or cdp.paramagnetic_centre_sim[sim_index] == None:
153 paramag_centre = zeros(3, float64)
154 else:
155 paramag_centre = array(cdp.paramagnetic_centre_sim[sim_index])
156 else:
157 paramag_centre = array(cdp.paramagnetic_centre)
158
159
160 return atomic_pos, paramag_centre
161
162
164 """Set up the data structures for optimisation using alignment tensors as base data sets.
165
166 @keyword sim_index: The index of the simulation to optimise. This should be None if
167 normal optimisation is desired.
168 @type sim_index: None or int
169 @return: The assembled data structures for using alignment tensors as the base
170 data for optimisation. These include:
171 - full_tensors, the data of the full alignment tensors.
172 - red_tensor_elem, the tensors as concatenated rank-1 5D arrays.
173 - red_tensor_err, the tensor errors as concatenated rank-1 5D
174 arrays.
175 - full_in_ref_frame, flags specifying if the tensor in the reference
176 frame is the full or reduced tensor.
177 @rtype: tuple of (list, numpy rank-1 array, numpy rank-1 array, numpy rank-1
178 array)
179 """
180
181
182 n = len(cdp.align_tensors.reduction)
183 full_tensors = zeros(n*5, float64)
184 red_tensors = zeros(n*5, float64)
185 red_err = ones(n*5, float64) * 1e-5
186 full_in_ref_frame = zeros(n, float64)
187
188
189 for i, tensor in tensor_loop(red=False):
190
191 full_tensors[5*i + 0] = tensor.Axx
192 full_tensors[5*i + 1] = tensor.Ayy
193 full_tensors[5*i + 2] = tensor.Axy
194 full_tensors[5*i + 3] = tensor.Axz
195 full_tensors[5*i + 4] = tensor.Ayz
196
197
198 if cdp.ref_domain == tensor.domain:
199 full_in_ref_frame[i] = 1
200
201
202 for i, tensor in tensor_loop(red=True):
203
204 if sim_index != None:
205 red_tensors[5*i + 0] = tensor.Axx_sim[sim_index]
206 red_tensors[5*i + 1] = tensor.Ayy_sim[sim_index]
207 red_tensors[5*i + 2] = tensor.Axy_sim[sim_index]
208 red_tensors[5*i + 3] = tensor.Axz_sim[sim_index]
209 red_tensors[5*i + 4] = tensor.Ayz_sim[sim_index]
210
211
212 else:
213 red_tensors[5*i + 0] = tensor.Axx
214 red_tensors[5*i + 1] = tensor.Ayy
215 red_tensors[5*i + 2] = tensor.Axy
216 red_tensors[5*i + 3] = tensor.Axz
217 red_tensors[5*i + 4] = tensor.Ayz
218
219
220 if hasattr(tensor, 'Axx_err'):
221 red_err[5*i + 0] = tensor.Axx_err
222 red_err[5*i + 1] = tensor.Ayy_err
223 red_err[5*i + 2] = tensor.Axy_err
224 red_err[5*i + 3] = tensor.Axz_err
225 red_err[5*i + 4] = tensor.Ayz_err
226
227
228 return full_tensors, red_tensors, red_err, full_in_ref_frame
229
230
232 """Set up the data structures for the fixed alignment tensors.
233
234 @return: The assembled data structures for the fixed alignment tensors.
235 @rtype: numpy rank-1 array.
236 """
237
238
239 n = align_tensor.num_tensors(skip_fixed=False) - align_tensor.num_tensors(skip_fixed=True)
240 tensors = zeros(n*5, float64)
241
242
243 if n == 0:
244 return None
245
246
247 index = 0
248 for i in range(len(cdp.align_tensors)):
249
250 if not opt_uses_align_data(cdp.align_tensors[i].name):
251 continue
252
253
254 tensors[5*index + 0] = cdp.align_tensors[i].Axx
255 tensors[5*index + 1] = cdp.align_tensors[i].Ayy
256 tensors[5*index + 2] = cdp.align_tensors[i].Axy
257 tensors[5*index + 3] = cdp.align_tensors[i].Axz
258 tensors[5*index + 4] = cdp.align_tensors[i].Ayz
259
260
261 index += 1
262
263
264 return tensors
265
266
268 """Initialise the target function for optimisation or direct calculation.
269
270 @param sim_index: The index of the simulation to optimise. This should be None if normal optimisation is desired.
271 @type sim_index: None or int
272 @param scaling: If True, diagonal scaling is enabled during optimisation to allow the problem to be better conditioned.
273 @type scaling: bool
274 """
275
276
277 if not hasattr(cdp, 'model'):
278 raise RelaxNoModelError('N-state')
279
280
281 if cdp.model == '2-domain':
282
283 if not hasattr(cdp, 'N'):
284 raise RelaxError("The number of states has not been set.")
285
286
287 if not hasattr(cdp, 'ref_domain'):
288 raise RelaxError("The reference domain has not been set.")
289
290
291 update_model()
292
293
294 param_vector = assemble_param_vector(sim_index=sim_index)
295
296
297 data_types = base_data_types()
298
299
300 probs = None
301 if hasattr(cdp, 'probs') and len(cdp.probs) and cdp.probs[0] != None:
302 probs = cdp.probs
303
304
305 scaling_matrix = None
306 if len(param_vector):
307 scaling_matrix = assemble_scaling_matrix(data_types=data_types, scaling=scaling)
308 param_vector = dot(inv(scaling_matrix), param_vector)
309
310
311 full_tensors, red_tensor_elem, red_tensor_err, full_in_ref_frame = None, None, None, None
312 if 'tensor' in data_types:
313 full_tensors, red_tensor_elem, red_tensor_err, full_in_ref_frame = minimise_setup_tensors(sim_index=sim_index)
314
315
316 pcs, pcs_err, pcs_weight, temp, frq, pcs_pseudo_flags = None, None, None, None, None, None
317 if 'pcs' in data_types:
318 pcs, pcs_err, pcs_weight, temp, frq, pcs_pseudo_flags = return_pcs_data(sim_index=sim_index)
319
320
321 rdcs, rdc_err, rdc_weight, rdc_vector, rdc_dj, absolute_rdc, T_flags, j_couplings, rdc_pseudo_flags = None, None, None, None, None, None, None, None, None
322 if 'rdc' in data_types:
323 rdcs, rdc_err, rdc_weight, rdc_vector, rdc_dj, absolute_rdc, T_flags, j_couplings, rdc_pseudo_flags = return_rdc_data(sim_index=sim_index)
324
325
326 fixed_tensors = None
327 if 'rdc' in data_types or 'pcs' in data_types:
328 full_tensors = minimise_setup_fixed_tensors()
329
330
331 fixed_tensors = []
332 for i in range(len(cdp.align_tensors)):
333
334 if not opt_uses_align_data(cdp.align_tensors[i].name):
335 continue
336
337 if cdp.align_tensors[i].fixed:
338 fixed_tensors.append(True)
339 else:
340 fixed_tensors.append(False)
341
342
343 atomic_pos, paramag_centre, centre_fixed = None, None, True
344 if 'pcs' in data_types or 'pre' in data_types:
345 atomic_pos, paramag_centre = minimise_setup_atomic_pos(sim_index=sim_index)
346
347
348 if hasattr(cdp, 'paramag_centre_fixed'):
349 centre_fixed = cdp.paramag_centre_fixed
350
351
352 model = N_state_opt(model=cdp.model, N=cdp.N, init_params=param_vector, probs=probs, full_tensors=full_tensors, red_data=red_tensor_elem, red_errors=red_tensor_err, full_in_ref_frame=full_in_ref_frame, fixed_tensors=fixed_tensors, pcs=pcs, rdcs=rdcs, pcs_errors=pcs_err, rdc_errors=rdc_err, T_flags=T_flags, j_couplings=j_couplings, rdc_pseudo_flags=rdc_pseudo_flags, pcs_pseudo_flags=pcs_pseudo_flags, pcs_weights=pcs_weight, rdc_weights=rdc_weight, rdc_vect=rdc_vector, temp=temp, frq=frq, dip_const=rdc_dj, absolute_rdc=absolute_rdc, atomic_pos=atomic_pos, paramag_centre=paramag_centre, scaling_matrix=scaling_matrix, centre_fixed=centre_fixed)
353
354
355 return model, param_vector, data_types, scaling_matrix
356