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