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