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 parameter handling."""
24
25
26 from numpy import array, float64, zeros
27 from re import search
28 from warnings import warn
29
30
31 from lib.errors import RelaxNoModelError
32 from lib.warnings import RelaxWarning
33 from pipe_control import align_tensor
34 from pipe_control.pipes import check_pipe
35 from pipe_control.align_tensor import opt_uses_align_data, opt_uses_tensor
36 from specific_analyses.n_state_model.data import base_data_types
37
38
40 """Assemble all the parameters of the model into a single array.
41
42 @param sim_index: The index of the simulation to optimise. This should be None if normal optimisation is desired.
43 @type sim_index: None or int
44 @return: The parameter vector used for optimisation.
45 @rtype: numpy array
46 """
47
48
49 if not hasattr(cdp, 'model') or not isinstance(cdp.model, str):
50 raise RelaxNoModelError
51
52
53 data_types = base_data_types()
54
55
56 param_vector = []
57
58
59 if opt_uses_align_data():
60 for i in range(len(cdp.align_tensors)):
61
62 if not opt_uses_tensor(cdp.align_tensors[i]):
63 continue
64
65
66 if not hasattr(cdp.align_tensors[i], 'A_5D'):
67 param_vector += [None, None, None, None, None]
68
69
70 else:
71 param_vector += list(cdp.align_tensors[i].A_5D)
72
73
74 if sim_index != None:
75
76 if cdp.model in ['2-domain', 'population']:
77 probs = cdp.probs_sim[sim_index]
78
79
80 if cdp.model == '2-domain':
81 alpha = cdp.alpha_sim[sim_index]
82 beta = cdp.beta_sim[sim_index]
83 gamma = cdp.gamma_sim[sim_index]
84
85
86 else:
87
88 if cdp.model in ['2-domain', 'population']:
89 probs = cdp.probs
90
91
92 if cdp.model == '2-domain':
93 alpha = cdp.alpha
94 beta = cdp.beta
95 gamma = cdp.gamma
96
97
98 if cdp.model in ['2-domain', 'population']:
99 param_vector = param_vector + probs[0:-1]
100
101
102 if cdp.model == '2-domain':
103 for i in range(cdp.N):
104 param_vector.append(alpha[i])
105 param_vector.append(beta[i])
106 param_vector.append(gamma[i])
107
108
109 if hasattr(cdp, 'paramag_centre_fixed') and not cdp.paramag_centre_fixed:
110 if not hasattr(cdp, 'paramagnetic_centre'):
111 for i in range(3):
112 param_vector.append(None)
113 elif sim_index != None:
114 if cdp.paramagnetic_centre_sim[sim_index] is None:
115 for i in range(3):
116 param_vector.append(None)
117 else:
118 for i in range(3):
119 param_vector.append(cdp.paramagnetic_centre_sim[sim_index][i])
120 else:
121 for i in range(3):
122 param_vector.append(cdp.paramagnetic_centre[i])
123
124
125 return array(param_vector, float64)
126
127
129 """Disassemble the parameter vector and place the values into the relevant variables.
130
131 For the 2-domain N-state model, the parameters are stored in the probability and Euler angle data structures. For the population N-state model, only the probabilities are stored. If RDCs are present and alignment tensors are optimised, then these are stored as well.
132
133 @keyword data_types: The base data types used in the optimisation. This list can contain the elements 'rdc', 'pcs' or 'tensor'.
134 @type data_types: list of str
135 @keyword param_vector: The parameter vector returned from optimisation.
136 @type param_vector: numpy array
137 @keyword sim_index: The index of the simulation to optimise. This should be None if normal optimisation is desired.
138 @type sim_index: None or int
139 """
140
141
142 if not hasattr(cdp, 'model') or not isinstance(cdp.model, str):
143 raise RelaxNoModelError
144
145
146 if ('rdc' in data_types or 'pcs' in data_types) and not align_tensor.all_tensors_fixed():
147
148 tensor_num = 0
149 for i in range(len(cdp.align_tensors)):
150
151 if not opt_uses_tensor(cdp.align_tensors[i]):
152 continue
153
154
155 if sim_index == None:
156 cdp.align_tensors[i].set(param='Axx', value=param_vector[5*tensor_num])
157 cdp.align_tensors[i].set(param='Ayy', value=param_vector[5*tensor_num+1])
158 cdp.align_tensors[i].set(param='Axy', value=param_vector[5*tensor_num+2])
159 cdp.align_tensors[i].set(param='Axz', value=param_vector[5*tensor_num+3])
160 cdp.align_tensors[i].set(param='Ayz', value=param_vector[5*tensor_num+4])
161
162
163 else:
164 cdp.align_tensors[i].set(param='Axx', value=param_vector[5*tensor_num], category='sim', sim_index=sim_index)
165 cdp.align_tensors[i].set(param='Ayy', value=param_vector[5*tensor_num+1], category='sim', sim_index=sim_index)
166 cdp.align_tensors[i].set(param='Axy', value=param_vector[5*tensor_num+2], category='sim', sim_index=sim_index)
167 cdp.align_tensors[i].set(param='Axz', value=param_vector[5*tensor_num+3], category='sim', sim_index=sim_index)
168 cdp.align_tensors[i].set(param='Ayz', value=param_vector[5*tensor_num+4], category='sim', sim_index=sim_index)
169
170
171 tensor_num += 1
172
173
174 param_vector = param_vector[5*tensor_num:]
175
176
177 if sim_index != None:
178
179 if cdp.model in ['2-domain', 'population']:
180 probs = cdp.probs_sim[sim_index]
181
182
183 if cdp.model == '2-domain':
184 alpha = cdp.alpha_sim[sim_index]
185 beta = cdp.beta_sim[sim_index]
186 gamma = cdp.gamma_sim[sim_index]
187
188
189 else:
190
191 if cdp.model in ['2-domain', 'population']:
192 probs = cdp.probs
193
194
195 if cdp.model == '2-domain':
196 alpha = cdp.alpha
197 beta = cdp.beta
198 gamma = cdp.gamma
199
200
201 if cdp.model in ['2-domain', 'population']:
202 for i in range(cdp.N-1):
203 probs[i] = param_vector[i]
204
205
206 probs[-1] = 1 - sum(probs[0:-1])
207
208
209 if cdp.model == '2-domain':
210 for i in range(cdp.N):
211 alpha[i] = param_vector[cdp.N-1 + 3*i]
212 beta[i] = param_vector[cdp.N-1 + 3*i + 1]
213 gamma[i] = param_vector[cdp.N-1 + 3*i + 2]
214
215
216 if hasattr(cdp, 'paramag_centre_fixed') and not cdp.paramag_centre_fixed:
217
218 if not hasattr(cdp, 'paramagnetic_centre'):
219 cdp.paramagnetic_centre = zeros(3, float64)
220
221
222 if sim_index == None:
223 cdp.paramagnetic_centre[0] = param_vector[-3]
224 cdp.paramagnetic_centre[1] = param_vector[-2]
225 cdp.paramagnetic_centre[2] = param_vector[-1]
226
227
228 else:
229 if cdp.paramagnetic_centre_sim[sim_index] is None:
230 cdp.paramagnetic_centre_sim[sim_index] = [None, None, None]
231 cdp.paramagnetic_centre_sim[sim_index][0] = param_vector[-3]
232 cdp.paramagnetic_centre_sim[sim_index][1] = param_vector[-2]
233 cdp.paramagnetic_centre_sim[sim_index][2] = param_vector[-1]
234
235
237 """Remove all structures or states which have no probability."""
238
239
240 check_pipe()
241
242
243 if not hasattr(cdp, 'model'):
244 raise RelaxNoModelError('N-state')
245
246
247 if not hasattr(cdp, 'probs'):
248 raise RelaxError("The N-state model populations do not exist.")
249
250
251 if not hasattr(cdp, 'select_models'):
252 cdp.state_select = [True] * cdp.N
253
254
255 for i in range(len(cdp.N)):
256
257 if cdp.probs[i] < 1e-5:
258 cdp.state_select[i] = False
259
260
262 """Function for setting up the linear constraint matrices A and b.
263
264 Standard notation
265 =================
266
267 The N-state model constraints are::
268
269 0 <= pc <= 1,
270
271 where p is the probability and c corresponds to state c.
272
273
274 Matrix notation
275 ===============
276
277 In the notation A.x >= b, where A is an matrix of coefficients, x is an array of parameter
278 values, and b is a vector of scalars, these inequality constraints are::
279
280 | 1 0 0 | | 0 |
281 | | | |
282 |-1 0 0 | | -1 |
283 | | | |
284 | 0 1 0 | | 0 |
285 | | | p0 | | |
286 | 0 -1 0 | | | | -1 |
287 | | . | p1 | >= | |
288 | 0 0 1 | | | | 0 |
289 | | | p2 | | |
290 | 0 0 -1 | | -1 |
291 | | | |
292 |-1 -1 -1 | | -1 |
293 | | | |
294 | 1 1 1 | | 0 |
295
296 This example is for a 4-state model, the last probability pn is not included as this
297 parameter does not exist (because the sum of pc is equal to 1). The Euler angle parameters
298 have been excluded here but will be included in the returned A and b objects. These
299 parameters simply add columns of zero to the A matrix and have no effect on b. The last two
300 rows correspond to the inequality::
301
302 0 <= pN <= 1.
303
304 As::
305 N-1
306 \
307 pN = 1 - > pc,
308 /__
309 c=1
310
311 then::
312
313 -p1 - p2 - ... - p(N-1) >= -1,
314
315 p1 + p2 + ... + p(N-1) >= 0.
316
317
318 @keyword data_types: The base data types used in the optimisation. This list can
319 contain the elements 'rdc', 'pcs' or 'tensor'.
320 @type data_types: list of str
321 @keyword scaling_matrix: The diagonal scaling matrix.
322 @type scaling_matrix: numpy rank-2 square matrix
323 @return: The matrices A and b.
324 @rtype: tuple of len 2 of a numpy rank-2, size NxM matrix and numpy
325 rank-1, size N array
326 """
327
328
329 pop_start = 0
330 if ('rdc' in data_types or 'pcs' in data_types) and not align_tensor.all_tensors_fixed():
331
332 for i in range(len(cdp.align_tensors)):
333
334 if not opt_uses_tensor(cdp.align_tensors[i]):
335 continue
336
337
338 pop_start += 5
339
340
341 A = []
342 b = []
343 zero_array = zeros(param_num(), float64)
344 i = pop_start
345 j = 0
346
347
348 if cdp.model in ['2-domain', 'population']:
349
350 for k in range(cdp.N - 1):
351
352 A.append(zero_array * 0.0)
353 A.append(zero_array * 0.0)
354 A[j][i] = 1.0
355 A[j+1][i] = -1.0
356 b.append(0.0)
357 b.append(-1.0 / scaling_matrix[i, i])
358 j = j + 2
359
360
361 i = i + 1
362
363
364 A.append(zero_array * 0.0)
365 A.append(zero_array * 0.0)
366 for i in range(pop_start, pop_start+cdp.N-1):
367 A[-2][i] = -1.0
368 A[-1][i] = 1.0
369 b.append(-1.0 / scaling_matrix[i, i])
370 b.append(0.0)
371
372
373 A = array(A, float64)
374 b = array(b, float64)
375
376
377 return A, b
378
379
381 """Set the number of states in the N-state model.
382
383 @param N: The number of states.
384 @type N: int
385 """
386
387
388 check_pipe()
389
390
391 cdp.N = N
392
393
394 if hasattr(cdp, 'model'):
395 update_model()
396
397
399 """Return the N-state model index for the given parameter.
400
401 @param param: The N-state model parameter.
402 @type param: str
403 @return: The N-state model index.
404 @rtype: str
405 """
406
407
408 if search('^p[0-9]*$', param):
409 return int(param[1:])
410
411
412 if search('^alpha', param):
413 return int(param[5:])
414
415
416 if search('^beta', param):
417 return int(param[4:])
418
419
420 if search('^gamma', param):
421 return int(param[5:])
422
423
424 return None
425
426
428 """Determine the number of parameters in the model.
429
430 @return: The number of model parameters.
431 @rtype: int
432 """
433
434
435 data_types = base_data_types()
436
437
438 num = 0
439
440
441 if ('rdc' in data_types or 'pcs' in data_types) and not align_tensor.all_tensors_fixed():
442
443 for i in range(len(cdp.align_tensors)):
444
445 if not opt_uses_tensor(cdp.align_tensors[i]):
446 continue
447
448
449 num += 5
450
451
452 if cdp.model in ['2-domain', 'population']:
453 num = num + (cdp.N - 1)
454
455
456 if cdp.model == '2-domain':
457 num = num + 3*cdp.N
458
459
460 if hasattr(cdp, 'paramag_centre_fixed') and not cdp.paramag_centre_fixed:
461 num = num + 3
462
463
464 return num
465
466
467 -def ref_domain(ref=None):
468 """Set the reference domain for the '2-domain' N-state model.
469
470 @param ref: The reference domain.
471 @type ref: str
472 """
473
474
475 check_pipe()
476
477
478 if not hasattr(cdp, 'model'):
479 raise RelaxNoModelError('N-state')
480
481
482 if cdp.model != '2-domain':
483 raise RelaxError("Setting the reference domain is only possible for the '2-domain' N-state model.")
484
485
486 exists = False
487 for tensor_cont in cdp.align_tensors:
488 if tensor_cont.domain == ref:
489 exists = True
490 if not exists:
491 raise RelaxError("The reference domain cannot be found within any of the loaded tensors.")
492
493
494 cdp.ref_domain = ref
495
496
497 update_model()
498
499
501 """Select the N-state model type.
502
503 @param model: The N-state model type. Can be one of '2-domain', 'population', or 'fixed'.
504 @type model: str
505 """
506
507
508 check_pipe()
509
510
511 if not model in ['2-domain', 'population', 'fixed']:
512 raise RelaxError("The model name " + repr(model) + " is invalid.")
513
514
515 if hasattr(cdp, 'model'):
516 warn(RelaxWarning("The N-state model has already been set up. Switching from model '%s' to '%s'." % (cdp.model, model)))
517
518
519 cdp.model = model
520
521
522 cdp.params = []
523
524
525 update_model()
526
527
529 """Update the model parameters as necessary."""
530
531
532 if not hasattr(cdp, 'params'):
533 cdp.params = []
534
535
536 if not hasattr(cdp, 'N'):
537
538 if hasattr(cdp, 'structure'):
539 cdp.N = cdp.structure.num_models()
540
541
542 else:
543 return
544
545
546 if not cdp.params:
547
548 if cdp.model in ['2-domain', 'population']:
549 for i in range(cdp.N-1):
550 cdp.params.append('p' + repr(i))
551
552
553 if cdp.model == '2-domain':
554 for i in range(cdp.N):
555 cdp.params.append('alpha' + repr(i))
556 cdp.params.append('beta' + repr(i))
557 cdp.params.append('gamma' + repr(i))
558
559
560 if cdp.model in ['2-domain', 'population']:
561 if not hasattr(cdp, 'probs'):
562 cdp.probs = [None] * cdp.N
563 if cdp.model == '2-domain':
564 if not hasattr(cdp, 'alpha'):
565 cdp.alpha = [None] * cdp.N
566 if not hasattr(cdp, 'beta'):
567 cdp.beta = [None] * cdp.N
568 if not hasattr(cdp, 'gamma'):
569 cdp.gamma = [None] * cdp.N
570
571
572 data_types = base_data_types()
573
574
575 if hasattr(cdp, 'align_ids'):
576 for id in cdp.align_ids:
577
578 if not hasattr(cdp, 'align_tensors'):
579 align_tensor.init(tensor=id, align_id=id)
580
581
582 exists = False
583 for tensor in cdp.align_tensors:
584 if id == tensor.align_id:
585 exists = True
586
587
588 if not exists:
589 align_tensor.init(tensor=id, align_id=id)
590