1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24 from math import sqrt
25 from numpy import array, dot, float64, ones, rank, transpose, zeros
26
27
28 from alignment_tensor import dAi_dAxx, dAi_dAyy, dAi_dAxy, dAi_dAxz, dAi_dAyz, to_tensor
29 from chi2 import chi2, dchi2_element, d2chi2_element
30 from float import isNaN
31 from paramag_centre import vectors_single_centre, vectors_centre_per_state
32 from pcs import ave_pcs_tensor, ave_pcs_tensor_ddeltaij_dAmn, pcs_tensor
33 from physical_constants import pcs_constant
34 from rdc import ave_rdc_tensor, ave_rdc_tensor_dDij_dAmn, rdc_tensor
35 from relax_errors import RelaxError, RelaxImplementError
36 from rotation_matrix import euler_to_R_zyz
37
38
40 """Class containing the target function of the optimisation of the N-state model."""
41
42 - def __init__(self, model=None, N=None, init_params=None, probs=None, full_tensors=None, red_data=None, red_errors=None, full_in_ref_frame=None, fixed_tensors=None, pcs=None, pcs_errors=None, pcs_weights=None, rdcs=None, rdc_errors=None, rdc_weights=None, xh_vect=None, temp=None, frq=None, dip_const=None, atomic_pos=None, paramag_centre=None, scaling_matrix=None, centre_fixed=True):
43 """Set up the class instance for optimisation.
44
45 The N-state models
46 ==================
47
48 All constant data required for the N-state model are initialised here. Depending on the base data used for optimisation, different data structures can be supplied. However a number of structures must be provided for the N-state model. These are:
49
50 - model, the type of N-state model. This can be '2-domain', 'population', or 'fixed'.
51 - N, the number of states (or structures).
52 - init_params, the initial parameter values.
53 - scaling_matrix, the matrix used for parameter scaling during optimisation.
54
55
56 2-domain N-state model
57 ----------------------
58
59 If the model type is set to '2-domain', then the following data structures should be supplied:
60
61 - full_tensors, the alignment tensors in matrix form.
62 - red_data, the alignment tensors in 5D form in a rank-1 array.
63 - red_errors, the alignment tensor errors in 5D form in a rank-1 array. This data is not obligatory.
64 - full_in_ref_frame, an array of flags specifying if the tensor in the reference frame is the full or reduced tensor.
65
66
67 The population N-state model
68 ============================
69
70 In this model, populations are optimised for each state. Additionally the alignment tensors for anisotropic data can also be optimised if they have not been supplied (through the full_tensors arg).
71
72
73 PCS base data
74 -------------
75
76 If pseudocontact shift data is to be used for optimisation, then the following should be supplied:
77
78 - pcs, the pseudocontact shifts.
79 - pcs_errors, the optional pseudocontact shift error values.
80 - temp, the temperatures of the experiments.
81 - frq, the frequencies of the experiments.
82
83
84 PCS and PRE base data
85 ---------------------
86
87 If either pseudocontact shift or PRE data is to be used for optimisation, then the following should be supplied:
88
89 - atomic_pos, the positions of all atoms.
90 - paramag_centre, the paramagnetic centre position.
91
92
93 RDC base data
94 -------------
95
96 If residual dipolar coupling data is to be used for optimisation, then the following should be supplied:
97
98 - rdcs, the residual dipolar couplings.
99 - rdc_errors, the optional residual dipolar coupling errors.
100 - xh_vect, the heteronucleus to proton unit vectors.
101 - dip_const, the dipolar constants.
102
103
104 @keyword model: The N-state model type. This can be one of '2-domain', 'population' or 'fixed'.
105 @type model: str
106 @keyword N: The number of states.
107 @type N: int
108 @keyword init_params: The initial parameter values. Optimisation must start at some point!
109 @type init_params: numpy float64 array
110 @keyword probs: The probabilities for each state c. The length of this list should be equal to N.
111 @type probs: list of float
112 @keyword full_tensors: An array of the {Sxx, Syy, Sxy, Sxz, Syz} values for all full tensors. The format is [Sxx1, Syy1, Sxy1, Sxz1, Syz1, Sxx2, Syy2, Sxy2, Sxz2, Syz2, ..., Sxxn, Syyn, Sxyn, Sxzn, Syzn]
113 @type full_tensors: list of rank-2, 3D numpy arrays
114 @keyword red_data: An array of the {Sxx, Syy, Sxy, Sxz, Syz} values for all reduced tensors. The format is the same as for full_tensors.
115 @type red_data: numpy float64 array
116 @keyword red_errors: An array of the {Sxx, Syy, Sxy, Sxz, Syz} errors for all reduced tensors. The array format is the same as for full_tensors.
117 @type red_errors: numpy float64 array
118 @keyword full_in_ref_frame: An array of flags specifying if the tensor in the reference frame is the full or reduced tensor.
119 @type full_in_ref_frame: numpy rank-1 array
120 @keyword fixed_tensors: An array of flags specifying if the tensor is fixed or will be optimised.
121 @type fixed_tensors: list of bool
122 @keyword pcs: The PCS lists. The first index must correspond to the different alignment media i and the second index to the spin systems j.
123 @type pcs: numpy rank-2 array
124 @keyword pcs_errors: The PCS error lists. The dimensions of this argument are the same as for 'pcs'.
125 @type pcs_errors: numpy rank-2 array
126 @keyword pcs_weights: The PCS weight lists. The dimensions of this argument are the same as for 'pcs'.
127 @type pcs_weights: numpy rank-2 array
128 @keyword rdcs: The RDC lists. The first index must correspond to the different alignment media i and the second index to the spin systems j.
129 @type rdcs: numpy rank-2 array
130 @keyword rdc_errors: The RDC error lists. The dimensions of this argument are the same as for 'rdcs'.
131 @type rdc_errors: numpy rank-2 array
132 @keyword rdc_weights: The RDC weight lists. The dimensions of this argument are the same as for 'rdcs'.
133 @type rdc_weights: numpy rank-2 array
134 @keyword xh_vect: The unit XH vector lists. The first index must correspond to the spin systems and the second index to each structure (its size being equal to the number of states).
135 @type xh_vect: numpy rank-2 array
136 @keyword temp: The temperature of each experiment, used for the PCS.
137 @type temp: numpy rank-1 array
138 @keyword frq: The frequency of each alignment, used for the PCS.
139 @type frq: numpy rank-1 array
140 @keyword dip_const: The dipolar constants for each XH vector. The indices correspond to the spin systems j.
141 @type dip_const: numpy rank-1 array
142 @keyword atomic_pos: The atomic positions of all spins. The first index is the spin systems j and the second is the structure or state c.
143 @type atomic_pos: numpy rank-3 array
144 @keyword paramag_centre: The paramagnetic centre position (or positions).
145 @type paramag_centre: numpy rank-1, 3D array or rank-2, Nx3 array
146 @keyword scaling_matrix: The square and diagonal scaling matrix.
147 @type scaling_matrix: numpy rank-2 array
148 @keyword centre_fixed: A flag which if False will cause the paramagnetic centre to be optimised.
149 @type centre_fixed: bool
150 """
151
152
153 self.N = N
154 self.params = 1.0 * init_params
155 self.fixed_tensors = fixed_tensors
156 self.deltaij = pcs
157 self.Dij = rdcs
158 self.dip_vect = xh_vect
159 self.dip_const = dip_const
160 self.temp = temp
161 self.frq = frq
162 self.atomic_pos = atomic_pos
163 self.paramag_centre = paramag_centre
164 self.centre_fixed = centre_fixed
165 self.total_num_params = len(init_params)
166
167
168 self.chi2 = 0.0
169 self.dchi2 = zeros((self.total_num_params), float64)
170 self.d2chi2 = zeros((self.total_num_params, self.total_num_params), float64)
171
172
173 self.scaling_matrix = scaling_matrix
174 if self.scaling_matrix != None:
175 self.scaling_flag = True
176 else:
177 self.scaling_flag = False
178
179
180 if model == '2-domain':
181
182 if red_data == None or not len(red_data):
183 raise RelaxError("The red_data argument " + repr(red_data) + " must be supplied.")
184 if red_errors == None or not len(red_errors):
185 raise RelaxError("The red_errors argument " + repr(red_errors) + " must be supplied.")
186 if full_in_ref_frame == None or not len(full_in_ref_frame):
187 raise RelaxError("The full_in_ref_frame argument " + repr(full_in_ref_frame) + " must be supplied.")
188
189
190 self.full_tensors = array(full_tensors, float64)
191 self.num_tensors = len(self.full_tensors) / 5
192 self.red_data = red_data
193 self.red_errors = red_errors
194 self.full_in_ref_frame = full_in_ref_frame
195
196
197 self.A = zeros((self.num_tensors, 3, 3), float64)
198 for i in range(self.num_tensors):
199 to_tensor(self.A[i], self.full_tensors[5*i:5*i+5])
200
201
202
203
204
205
206 self.R = zeros((self.N, 3, 3), float64)
207 self.RT = zeros((self.N, 3, 3), float64)
208 self.red_bc = zeros((self.num_tensors, 3, 3), float64)
209 self.red_bc_vector = zeros(self.num_tensors*5, float64)
210
211
212 self.func = self.func_2domain
213 self.dfunc = None
214 self.d2func = None
215
216
217 elif model in ['population', 'fixed']:
218
219 self.rdc_flag = True
220 self.pcs_flag = True
221 if rdcs == None or len(rdcs) == 0:
222 self.rdc_flag = False
223 if pcs == None or len(pcs) == 0:
224 self.pcs_flag = False
225
226
227 if self.rdc_flag and (xh_vect == None or not len(xh_vect)):
228 raise RelaxError("The xh_vect argument " + repr(xh_vect) + " must be supplied.")
229 if self.pcs_flag and (atomic_pos == None or not len(atomic_pos)):
230 raise RelaxError("The atomic_pos argument " + repr(atomic_pos) + " must be supplied.")
231
232
233 self.num_spins = 0
234 if self.rdc_flag:
235 self.num_spins = len(rdcs[0])
236 elif self.pcs_flag:
237 self.num_spins = len(pcs[0])
238
239
240 self.num_align = 0
241 if self.rdc_flag:
242 self.num_align = len(rdcs)
243 elif self.pcs_flag:
244 self.num_align = len(pcs)
245
246
247 self.A = zeros((self.num_align, 3, 3), float64)
248 self.dA = zeros((5, 3, 3), float64)
249
250
251 if full_tensors != None:
252
253 self.full_tensors = array(full_tensors, float64)
254
255
256 self.num_align_params = 0
257 index = 0
258 for i in range(self.num_align):
259
260 if fixed_tensors[i]:
261 to_tensor(self.A[i], self.full_tensors[5*index:5*index+5])
262 index += 1
263
264
265 if not fixed_tensors[i]:
266 self.num_align_params += 5
267
268
269 else:
270
271 dAi_dAxx(self.dA[0])
272 dAi_dAyy(self.dA[1])
273 dAi_dAxy(self.dA[2])
274 dAi_dAxz(self.dA[3])
275 dAi_dAyz(self.dA[4])
276
277
278 if self.pcs_flag:
279 err = False
280 for i in xrange(len(pcs_errors)):
281 for j in xrange(len(pcs_errors[i])):
282 if not isNaN(pcs_errors[i, j]):
283 err = True
284 if err:
285 self.pcs_sigma_ij = pcs_errors
286 else:
287
288 self.pcs_sigma_ij = 0.03 * 1e-6 * ones((self.num_align, self.num_spins), float64)
289
290
291 if self.rdc_flag:
292 err = False
293 for i in xrange(len(rdc_errors)):
294 for j in xrange(len(rdc_errors[i])):
295 if not isNaN(rdc_errors[i, j]):
296 err = True
297 if err:
298 self.rdc_sigma_ij = rdc_errors
299 else:
300
301 self.rdc_sigma_ij = ones((self.num_align, self.num_spins), float64)
302
303
304 if self.rdc_flag:
305 self.missing_Dij = zeros((self.num_align, self.num_spins), float64)
306
307
308 if self.pcs_flag:
309 self.missing_deltaij = zeros((self.num_align, self.num_spins), float64)
310
311
312 if self.rdc_flag or self.pcs_flag:
313 for i in xrange(self.num_align):
314 for j in xrange(self.num_spins):
315 if self.rdc_flag:
316 if isNaN(self.Dij[i, j]):
317
318 self.missing_Dij[i, j] = 1
319
320
321 self.Dij[i, j] = 0.0
322
323
324 self.rdc_sigma_ij[i, j] = 1.0
325
326
327 rdc_weights[i, j] = 1.0
328
329 if self.pcs_flag:
330 if isNaN(self.deltaij[i, j]):
331
332 self.missing_deltaij[i, j] = 1
333
334
335 self.deltaij[i, j] = 0.0
336
337
338 self.pcs_sigma_ij[i, j] = 1.0
339
340
341 pcs_weights[i, j] = 1.0
342
343
344 if self.rdc_flag:
345 self.rdc_sigma_ij[i, j] = self.rdc_sigma_ij[i, j] / sqrt(rdc_weights[i, j])
346
347
348 if self.pcs_flag:
349 self.pcs_sigma_ij[i, j] = self.pcs_sigma_ij[i, j] / sqrt(pcs_weights[i, j])
350
351
352
353 if self.pcs_flag:
354
355 self.paramag_unit_vect = zeros(atomic_pos.shape, float64)
356 self.paramag_dist = zeros((self.num_spins, self.N), float64)
357 self.pcs_const = zeros((self.num_align, self.num_spins, self.N), float64)
358 if self.paramag_centre == None:
359 self.paramag_centre = zeros(3, float64)
360
361
362 self.paramag_info()
363
364
365 self.deltaij_theta = zeros((self.num_align, self.num_spins), float64)
366 self.ddeltaij_theta = zeros((self.total_num_params, self.num_align, self.num_spins), float64)
367 self.d2deltaij_theta = zeros((self.total_num_params, self.total_num_params, self.num_align, self.num_spins), float64)
368
369
370 self.Dij_theta = zeros((self.num_align, self.num_spins), float64)
371 self.dDij_theta = zeros((self.total_num_params, self.num_align, self.num_spins), float64)
372 self.d2Dij_theta = zeros((self.total_num_params, self.total_num_params, self.num_align, self.num_spins), float64)
373
374
375 if not self.centre_fixed:
376 self.func = self.func_population
377 self.dfunc = None
378 self.d2func = None
379
380
381 else:
382 self.func = self.func_population
383 self.dfunc = self.dfunc_population
384 self.d2func = self.d2func_population
385
386
387 if model == 'fixed':
388
389 self.func = self.func_tensor_opt
390 self.dfunc = self.dfunc_tensor_opt
391 self.d2func = self.d2func_tensor_opt
392
393
394 self.zero_hessian = zeros(self.num_spins, float64)
395
396
397 if probs:
398 self.probs = probs
399
400
401 else:
402 self.probs = ones(self.N, float64) / self.N
403
404
405 - def func_2domain(self, params):
406 """The target function for optimisation of the 2-domain N-state model.
407
408 This function should be passed to the optimisation algorithm. It accepts, as an array, a vector of parameter values and, using these, returns the single chi-squared value corresponding to that coordinate in the parameter space. If no tensor errors are supplied, then the SSE (the sum of squares error) value is returned instead. The chi-squared is simply the SSE normalised to unit variance (the SSE divided by the error squared).
409
410 @param params: The vector of parameter values.
411 @type params: list of float
412 @return: The chi-squared or SSE value.
413 @rtype: float
414 """
415
416
417 if self.scaling_flag:
418 params = dot(params, self.scaling_matrix)
419
420
421 self.red_bc = self.red_bc * 0.0
422
423
424 for c in xrange(self.N):
425
426 euler_to_R_zyz(params[self.N-1+3*c], params[self.N-1+3*c+1], params[self.N-1+3*c+2], self.R[c])
427
428
429 self.RT[c] = transpose(self.R[c])
430
431
432 if c < self.N-1:
433 pc = params[c]
434
435
436 else:
437 pc = 1.0
438 for c2 in xrange(self.N-1):
439 pc = pc - params[c2]
440
441
442 for i in xrange(self.num_tensors):
443
444 if self.full_in_ref_frame[i]:
445 self.red_bc[i] = self.red_bc[i] + pc * dot(self.RT[c], dot(self.A[i], self.R[c]))
446
447
448 else:
449 self.red_bc[i] = self.red_bc[i] + pc * dot(self.R[c], dot(self.A[i], self.RT[c]))
450
451
452 for i in xrange(self.num_tensors):
453 self.red_bc_vector[5*i] = self.red_bc[i, 0, 0]
454 self.red_bc_vector[5*i+1] = self.red_bc[i, 1, 1]
455 self.red_bc_vector[5*i+2] = self.red_bc[i, 0, 1]
456 self.red_bc_vector[5*i+3] = self.red_bc[i, 0, 2]
457 self.red_bc_vector[5*i+4] = self.red_bc[i, 1, 2]
458
459
460 return chi2(self.red_data, self.red_bc_vector, self.red_errors)
461
462
464 """The target function for optimisation of the flexible population N-state model.
465
466 Description
467 ===========
468
469 This function should be passed to the optimisation algorithm. It accepts, as an array, a vector of parameter values and, using these, returns the single chi-squared value corresponding to that coordinate in the parameter space. If no RDC errors are supplied, then the SSE (the sum of squares error) value is returned instead. The chi-squared is simply the SSE normalised to unit variance (the SSE divided by the error squared).
470
471
472 Indices
473 =======
474
475 For this calculation, five indices are looped over and used in the various data structures. These include:
476 - i, the index over alignments,
477 - j, the index over spin systems,
478 - c, the index over the N-states (or over the structures),
479 - n, the index over the first dimension of the alignment tensor n = {x, y, z},
480 - m, the index over the second dimension of the alignment tensor m = {x, y, z}.
481
482
483 Equations
484 =========
485
486 To calculate the function value, a chain of equations are used. This includes the chi-squared equation and the RDC equation.
487
488
489 The chi-squared equation
490 ------------------------
491
492 The equations are::
493
494 ___
495 \ (Dij - Dij(theta)) ** 2
496 chi^2(theta) = > ----------------------- ,
497 /__ sigma_ij ** 2
498 ij
499
500 ___
501 \ (delta_ij - delta_ij(theta)) ** 2
502 chi^2(theta) = > --------------------------------- ,
503 /__ sigma_ij ** 2
504 ij
505
506 where:
507 - theta is the parameter vector,
508 - Dij are the measured RDCs for alignment i, spin j,
509 - Dij(theta) are the back calculated RDCs for alignment i, spin j,
510 - delta_ij are the measured PCSs for alignment i, spin j,
511 - delta_ij(theta) are the back calculated PCSs for alignment i, spin j,
512 - sigma_ij are the RDC or PCS errors.
513
514 Both chi-squared values sum.
515
516
517 The RDC equation
518 ----------------
519
520 The RDC equation is::
521
522 _N_
523 \ T
524 Dij(theta) = dj > pc . mu_jc . Ai . mu_jc,
525 /__
526 c=1
527
528 where:
529 - dj is the dipolar constant for spin j,
530 - N is the total number of states or structures,
531 - pc is the weight or probability associated with state c,
532 - mu_jc is the unit vector corresponding to spin j and state c,
533 - Ai is the alignment tensor.
534
535 The dipolar constant is henceforth defined as::
536
537 dj = 3 / (2pi) d',
538
539 where the factor of 2pi is to convert from units of rad.s^-1 to Hertz, the factor of 3 is associated with the alignment tensor and the pure dipolar constant in SI units is::
540
541 mu0 gI.gS.h_bar
542 d' = - --- ----------- ,
543 4pi r**3
544
545 where:
546 - mu0 is the permeability of free space,
547 - gI and gS are the gyromagnetic ratios of the I and S spins,
548 - h_bar is Dirac's constant which is equal to Planck's constant divided by 2pi,
549 - r is the distance between the two spins.
550
551
552 The PCS equation
553 ----------------
554
555 The PCS equation is::
556
557 _N_
558 \ T
559 delta_ij(theta) = > pc . dijc . mu_jc . Ai . mu_jc,
560 /__
561 c=1
562
563 where:
564 - djci is the PCS constant for spin j, state c and experiment or alignment i,
565 - N is the total number of states or structures,
566 - pc is the weight or probability associated with state c,
567 - mu_jc is the unit vector corresponding to spin j and state c,
568 - Ai is the alignment tensor.
569
570 The PCS constant is defined as::
571
572 mu0 15kT 1
573 dijc = --- ----- ---- ,
574 4pi Bo**2 r**3
575
576 where:
577 - mu0 is the permeability of free space,
578 - k is Boltzmann's constant,
579 - T is the absolute temperature (different for each experiment),
580 - Bo is the magnetic field strength (different for each experiment),
581 - r is the distance between the paramagnetic centre (electron spin) and the nuclear spin (different for each spin and state).
582
583
584 Stored data structures
585 ======================
586
587 There are a number of data structures calculated by this function and stored for subsequent use in the gradient and Hessian functions. This include the back calculated RDCs and the alignment tensors.
588
589 Dij(theta)
590 ----------
591
592 The back calculated RDCs. This is a rank-2 tensor with indices {i, j}.
593
594 delta_ij(theta)
595 ---------------
596
597 The back calculated PCS. This is a rank-2 tensor with indices {i, j}.
598
599 Ai
600 --
601
602 The alignment tensors. This is a rank-3 tensor with indices {i, n, m}.
603
604
605 @param params: The vector of parameter values.
606 @type params: numpy rank-1 array
607 @return: The chi-squared or SSE value.
608 @rtype: float
609 """
610
611
612 if self.scaling_flag:
613 params = dot(params, self.scaling_matrix)
614
615
616 chi2_sum = 0.0
617
618
619 if self.N > 1:
620 self.probs = params[-(self.N-1):]
621
622
623 if not self.centre_fixed:
624
625 self.paramag_centre = params[-3:]
626
627
628 self.paramag_info()
629
630
631 index = 0
632 for i in xrange(self.num_align):
633
634 if not self.fixed_tensors[i]:
635 to_tensor(self.A[i], params[5*index:5*index + 5])
636 index += 1
637
638
639 for j in xrange(self.num_spins):
640
641 if self.rdc_flag:
642
643 if not self.missing_Dij[i, j]:
644 self.Dij_theta[i, j] = ave_rdc_tensor(self.dip_const[j], self.dip_vect[j], self.N, self.A[i], weights=self.probs)
645
646
647 if self.pcs_flag:
648
649 if not self.missing_deltaij[i, j]:
650 self.deltaij_theta[i, j] = ave_pcs_tensor(self.pcs_const[i, j], self.paramag_unit_vect[j], self.N, self.A[i], weights=self.probs)
651
652
653 if self.rdc_flag:
654 chi2_sum = chi2_sum + chi2(self.Dij[i], self.Dij_theta[i], self.rdc_sigma_ij[i])
655
656
657 if self.pcs_flag:
658 chi2_sum = chi2_sum + chi2(self.deltaij[i], self.deltaij_theta[i], self.pcs_sigma_ij[i])
659
660
661 return chi2_sum
662
663
665 """The target function for optimisation of the alignment tensor from RDC and/or PCS data.
666
667 Description
668 ===========
669
670 This function should be passed to the optimisation algorithm. It accepts, as an array, a vector of parameter values and, using these, returns the single chi-squared value corresponding to that coordinate in the parameter space. If no RDC or PCS errors are supplied, then the SSE (the sum of squares error) value is returned instead. The chi-squared is simply the SSE normalised to unit variance (the SSE divided by the error squared).
671
672
673 Indices
674 =======
675
676 For this calculation, five indices are looped over and used in the various data structures. These include:
677 - i, the index over alignments,
678 - j, the index over spin systems,
679 - c, the index over the N-states (or over the structures),
680 - n, the index over the first dimension of the alignment tensor n = {x, y, z},
681 - m, the index over the second dimension of the alignment tensor m = {x, y, z}.
682
683
684 Equations
685 =========
686
687 To calculate the function value, a chain of equations are used. This includes the chi-squared equation and the RDC and PCS equations.
688
689
690 The chi-squared equation
691 ------------------------
692
693 The equations are::
694
695 ___
696 \ (Dij - Dij(theta)) ** 2
697 chi^2(theta) = > ----------------------- ,
698 /__ sigma_ij ** 2
699 ij
700
701 ___
702 \ (delta_ij - delta_ij(theta)) ** 2
703 chi^2(theta) = > --------------------------------- ,
704 /__ sigma_ij ** 2
705 ij
706
707 where:
708 - theta is the parameter vector,
709 - Dij are the measured RDCs for alignment i, spin j,
710 - Dij(theta) are the back calculated RDCs for alignment i, spin j,
711 - delta_ij are the measured PCSs for alignment i, spin j,
712 - delta_ij(theta) are the back calculated PCSs for alignment i, spin j,
713 - sigma_ij are the RDC or PCS errors.
714
715 Both chi-squared values sum.
716
717
718 The RDC equation
719 ----------------
720
721 The RDC equation is::
722
723 _N_
724 dj \ T
725 Dij(theta) = -- > mu_jc . Ai . mu_jc,
726 N /__
727 c=1
728
729 where:
730 - dj is the dipolar constant for spin j,
731 - N is the total number of states or structures,
732 - mu_jc is the unit vector corresponding to spin j and state c,
733 - Ai is the alignment tensor.
734
735 The dipolar constant is henceforth defined as::
736
737 dj = 3 / (2pi) d',
738
739 where the factor of 2pi is to convert from units of rad.s^-1 to Hertz, the factor of 3 is associated with the alignment tensor and the pure dipolar constant in SI units is::
740
741 mu0 gI.gS.h_bar
742 d' = - --- ----------- ,
743 4pi r**3
744
745 where:
746 - mu0 is the permeability of free space,
747 - gI and gS are the gyromagnetic ratios of the I and S spins,
748 - h_bar is Dirac's constant which is equal to Planck's constant divided by 2pi,
749 - r is the distance between the two spins.
750
751
752 The PCS equation
753 ----------------
754
755 The PCS equation is::
756
757 _N_
758 1 \ T
759 delta_ij(theta) = - > dijc . mu_jc . Ai . mu_jc,
760 N /__
761 c=1
762
763 where:
764 - djci is the PCS constant for spin j, state c and experiment or alignment i,
765 - N is the total number of states or structures,
766 - mu_jc is the unit vector corresponding to spin j and state c,
767 - Ai is the alignment tensor.
768
769 The PCS constant is defined as::
770
771 mu0 15kT 1
772 dijc = --- ----- ---- ,
773 4pi Bo**2 r**3
774
775 where:
776 - mu0 is the permeability of free space,
777 - k is Boltzmann's constant,
778 - T is the absolute temperature (different for each experiment),
779 - Bo is the magnetic field strength (different for each experiment),
780 - r is the distance between the paramagnetic centre (electron spin) and the nuclear spin (different for each spin and state).
781
782
783 Stored data structures
784 ======================
785
786 There are a number of data structures calculated by this function and stored for subsequent use in the gradient and Hessian functions. This include the back calculated RDCs and the alignment tensors.
787
788 Dij(theta)
789 ----------
790
791 The back calculated RDCs. This is a rank-2 tensor with indices {i, j}.
792
793 delta_ij(theta)
794 ---------------
795
796 The back calculated PCS. This is a rank-2 tensor with indices {i, j}.
797
798 Ai
799 --
800
801 The alignment tensors. This is a rank-3 tensor with indices {i, n, m}.
802
803
804 @param params: The vector of parameter values.
805 @type params: numpy rank-1 array
806 @return: The chi-squared or SSE value.
807 @rtype: float
808 """
809
810
811 if self.scaling_flag:
812 params = dot(params, self.scaling_matrix)
813
814
815 chi2_sum = 0.0
816
817
818 if not self.centre_fixed:
819
820 self.paramag_centre = params[-3:]
821
822
823 self.paramag_info()
824
825
826 index = 0
827 for i in xrange(self.num_align):
828
829 if not self.fixed_tensors[i]:
830 to_tensor(self.A[i], params[5*index:5*index + 5])
831 index += 1
832
833
834 for j in xrange(self.num_spins):
835
836 if self.rdc_flag:
837
838 if not self.missing_Dij[i, j]:
839 self.Dij_theta[i, j] = ave_rdc_tensor(self.dip_const[j], self.dip_vect[j], self.N, self.A[i], weights=self.probs)
840
841
842 if self.pcs_flag:
843
844 if not self.missing_deltaij[i, j]:
845 self.deltaij_theta[i, j] = ave_pcs_tensor(self.pcs_const[i, j], self.paramag_unit_vect[j], self.N, self.A[i], weights=self.probs)
846
847
848 if self.rdc_flag:
849 chi2_sum = chi2_sum + chi2(self.Dij[i], self.Dij_theta[i], self.rdc_sigma_ij[i])
850
851
852 if self.pcs_flag:
853 chi2_sum = chi2_sum + chi2(self.deltaij[i], self.deltaij_theta[i], self.pcs_sigma_ij[i])
854
855
856 return chi2_sum
857
858
860 """The gradient function for optimisation of the flexible population N-state model.
861
862 Description
863 ===========
864
865 This function should be passed to the optimisation algorithm. It accepts, as an array, a vector of parameter values and, using these, returns the chi-squared gradient corresponding to that coordinate in the parameter space. If no RDC errors are supplied, then the SSE (the sum of squares error) gradient is returned instead. The chi-squared gradient is simply the SSE gradient normalised to unit variance (the SSE divided by the error squared).
866
867
868 Indices
869 =======
870
871 For this calculation, six indices are looped over and used in the various data structures. These include:
872 - k, the index over all parameters,
873 - i, the index over alignments,
874 - j, the index over spin systems,
875 - c, the index over the N-states (or over the structures),
876 - m, the index over the first dimension of the alignment tensor m = {x, y, z}.
877 - n, the index over the second dimension of the alignment tensor n = {x, y, z},
878
879
880 Equations
881 =========
882
883 To calculate the chi-squared gradient, a chain of equations are used. This includes the chi-squared gradient, the RDC gradient and the alignment tensor gradient.
884
885
886 The chi-squared gradient
887 ------------------------
888
889 The equation is::
890 ___
891 dchi^2(theta) \ / Dij - Dij(theta) dDij(theta) \
892 ------------- = -2 > | ---------------- . ----------- |
893 dthetak /__ \ sigma_ij**2 dthetak /
894 ij
895
896 where:
897 - theta is the parameter vector,
898 - Dij are the measured RDCs or PCSs,
899 - Dij(theta) are the back calculated RDCs or PCSs,
900 - sigma_ij are the RDC or PCS errors,
901 - dDij(theta)/dthetak is the RDC or PCS gradient for parameter k.
902
903
904 The RDC gradient
905 ----------------
906
907 This gradient is different for the various parameter types.
908
909 pc partial derivative
910 ~~~~~~~~~~~~~~~~~~~~~
911
912 The population parameter partial derivative is::
913
914 dDij(theta) T
915 ----------- = dj . mu_jc . Ai . mu_jc,
916 dpc
917
918 where:
919 - dj is the dipolar constant for spin j,
920 - mu_jc is the unit vector corresponding to spin j and state c,
921 - Ai is the alignment tensor.
922
923 Amn partial derivative
924 ~~~~~~~~~~~~~~~~~~~~~~
925
926 The alignment tensor element partial derivative is::
927
928 _N_
929 dDij(theta) \ T dAi
930 ----------- = dj > pc . mu_jc . ---- . mu_jc,
931 dAmn /__ dAmn
932 c=1
933
934 where:
935 - dj is the dipolar constant for spin j,
936 - pc is the weight or probability associated with state c,
937 - mu_jc is the unit vector corresponding to spin j and state c,
938 - dAi/dAmn is the partial derivative of the alignment tensor with respect to element Amn.
939
940
941 The PCS gradient
942 ----------------
943
944 This gradient is also different for the various parameter types.
945
946 pc partial derivative
947 ~~~~~~~~~~~~~~~~~~~~~
948
949 The population parameter partial derivative is::
950
951 ddeltaij(theta) T
952 --------------- = dijc . mu_jc . Ai . mu_jc,
953 dpc
954
955 where:
956 - djc is the pseudocontact shift constant for spin j and state c,
957 - mu_jc is the unit vector corresponding to spin j and state c,
958 - Ai is the alignment tensor.
959
960 Amn partial derivative
961 ~~~~~~~~~~~~~~~~~~~~~~
962
963 The alignment tensor element partial derivative is::
964
965 _N_
966 ddelta_ij(theta) \ T dAi
967 ---------------- = > pc . djc . mu_jc . ---- . mu_jc,
968 dAmn /__ dAmn
969 c=1
970
971 where:
972 - djc is the pseudocontact shift constant for spin j and state c,
973 - pc is the weight or probability associated with state c,
974 - mu_jc is the unit vector corresponding to spin j and state c,
975 - dAi/dAmn is the partial derivative of the alignment tensor with respect to element Amn.
976
977
978 The alignment tensor gradient
979 -----------------------------
980
981 The five unique elements of the tensor {Axx, Ayy, Axy, Axz, Ayz} give five different partial derivatives. These are::
982
983 dAi | 1 0 0 |
984 ---- = | 0 0 0 |,
985 dAxx | 0 0 -1 |
986
987 dAi | 0 0 0 |
988 ---- = | 0 1 0 |,
989 dAyy | 0 0 -1 |
990
991 dAi | 0 1 0 |
992 ---- = | 1 0 0 |,
993 dAxy | 0 0 0 |
994
995 dAi | 0 0 1 |
996 ---- = | 0 0 0 |,
997 dAxz | 1 0 0 |
998
999 dAi | 0 0 0 |
1000 ---- = | 0 0 1 |.
1001 dAyz | 0 1 0 |
1002
1003 As these are invariant, they can be pre-calculated.
1004
1005
1006 Stored data structures
1007 ======================
1008
1009 There are a number of data structures calculated by this function and stored for subsequent use in the Hessian function. This include the back calculated RDC and PCS gradients and the alignment tensor gradients.
1010
1011 dDij(theta)/dthetak
1012 -------------------
1013
1014 The back calculated RDC gradient. This is a rank-3 tensor with indices {k, i, j}.
1015
1016 ddeltaij(theta)/dthetak
1017 -----------------------
1018
1019 The back calculated PCS gradient. This is a rank-3 tensor with indices {k, i, j}.
1020
1021 dAi/dAmn
1022 --------
1023
1024 The alignment tensor gradients. This is a rank-3 tensor with indices {5, n, m}.
1025
1026
1027 @param params: The vector of parameter values. This is unused as it is assumed that
1028 func_population() was called first.
1029 @type params: numpy rank-1 array
1030 @return: The chi-squared or SSE gradient.
1031 @rtype: numpy rank-1 array
1032 """
1033
1034
1035 if self.scaling_flag:
1036 params = dot(params, self.scaling_matrix)
1037
1038
1039 self.dchi2 = self.dchi2 * 0.0
1040
1041
1042 index = 0
1043 for i in xrange(self.num_align):
1044
1045 if self.fixed_tensors[i]:
1046 continue
1047
1048
1049 for j in xrange(self.num_spins):
1050
1051 if self.fixed_tensors[i] and self.rdc_flag and not self.missing_Dij[i, j]:
1052 self.dDij_theta[i*5, i, j] = ave_rdc_tensor_dDij_dAmn(self.dip_const[j], self.dip_vect[j], self.N, self.dA[0], weights=self.probs)
1053 self.dDij_theta[i*5+1, i, j] = ave_rdc_tensor_dDij_dAmn(self.dip_const[j], self.dip_vect[j], self.N, self.dA[1], weights=self.probs)
1054 self.dDij_theta[i*5+2, i, j] = ave_rdc_tensor_dDij_dAmn(self.dip_const[j], self.dip_vect[j], self.N, self.dA[2], weights=self.probs)
1055 self.dDij_theta[i*5+3, i, j] = ave_rdc_tensor_dDij_dAmn(self.dip_const[j], self.dip_vect[j], self.N, self.dA[3], weights=self.probs)
1056 self.dDij_theta[i*5+4, i, j] = ave_rdc_tensor_dDij_dAmn(self.dip_const[j], self.dip_vect[j], self.N, self.dA[4], weights=self.probs)
1057
1058
1059 if self.fixed_tensors[i] and self.pcs_flag and not self.missing_deltaij[i, j]:
1060 self.ddeltaij_theta[i*5, i, j] = ave_pcs_tensor_ddeltaij_dAmn(self.pcs_const[i, j], self.paramag_unit_vect[j], self.N, self.dA[0], weights=self.probs)
1061 self.ddeltaij_theta[i*5+1, i, j] = ave_pcs_tensor_ddeltaij_dAmn(self.pcs_const[i, j], self.paramag_unit_vect[j], self.N, self.dA[1], weights=self.probs)
1062 self.ddeltaij_theta[i*5+2, i, j] = ave_pcs_tensor_ddeltaij_dAmn(self.pcs_const[i, j], self.paramag_unit_vect[j], self.N, self.dA[2], weights=self.probs)
1063 self.ddeltaij_theta[i*5+3, i, j] = ave_pcs_tensor_ddeltaij_dAmn(self.pcs_const[i, j], self.paramag_unit_vect[j], self.N, self.dA[3], weights=self.probs)
1064 self.ddeltaij_theta[i*5+4, i, j] = ave_pcs_tensor_ddeltaij_dAmn(self.pcs_const[i, j], self.paramag_unit_vect[j], self.N, self.dA[4], weights=self.probs)
1065
1066
1067 for c in xrange(self.N - 1):
1068
1069 param_index = self.num_align_params + c
1070
1071
1072 for j in xrange(self.num_spins):
1073
1074 if self.rdc_flag and not self.missing_Dij[i, j]:
1075 self.dDij_theta[param_index, i, j] = rdc_tensor(self.dip_const[j], self.dip_vect[j, c], self.A[i])
1076
1077
1078 if self.pcs_flag and not self.missing_deltaij[i, j]:
1079 self.ddeltaij_theta[param_index, i, j] = pcs_tensor(self.pcs_const[i, j, c], self.paramag_unit_vect[j, c], self.A[i])
1080
1081
1082 for k in xrange(self.total_num_params):
1083
1084 if self.rdc_flag:
1085 self.dchi2[k] = self.dchi2[k] + dchi2_element(self.Dij[i], self.Dij_theta[i], self.dDij_theta[k, i], self.rdc_sigma_ij[i])
1086
1087
1088 if self.pcs_flag:
1089 self.dchi2[k] = self.dchi2[k] + dchi2_element(self.deltaij[i], self.deltaij_theta[i], self.ddeltaij_theta[k, i], self.pcs_sigma_ij[i])
1090
1091
1092 index += 1
1093
1094
1095 if self.scaling_flag:
1096 self.dchi2 = dot(self.dchi2, self.scaling_matrix)
1097
1098
1099 return self.dchi2
1100
1101
1103 """The gradient function for optimisation of the alignment tensor from RDC and/or PCS data.
1104
1105 Description
1106 ===========
1107
1108 This function should be passed to the optimisation algorithm. It accepts, as an array, a vector of parameter values and, using these, returns the chi-squared gradient corresponding to that coordinate in the parameter space. If no RDC errors are supplied, then the SSE (the sum of squares error) gradient is returned instead. The chi-squared gradient is simply the SSE gradient normalised to unit variance (the SSE divided by the error squared).
1109
1110 Indices
1111 =======
1112
1113 For this calculation, six indices are looped over and used in the various data structures. These include:
1114 - k, the index over all parameters,
1115 - i, the index over alignments,
1116 - j, the index over spin systems,
1117 - c, the index over the N-states (or over the structures),
1118 - m, the index over the first dimension of the alignment tensor m = {x, y, z}.
1119 - n, the index over the second dimension of the alignment tensor n = {x, y, z},
1120
1121
1122 Equations
1123 =========
1124
1125 To calculate the chi-squared gradient, a chain of equations are used. This includes the chi-squared gradient, the RDC gradient and the alignment tensor gradient.
1126
1127
1128 The chi-squared gradient
1129 ------------------------
1130
1131 The equation is::
1132 ___
1133 dchi^2(theta) \ / Dij - Dij(theta) dDij(theta) \
1134 ------------- = -2 > | ---------------- . ----------- |
1135 dthetak /__ \ sigma_ij**2 dthetak /
1136 ij
1137
1138 where:
1139 - theta is the parameter vector,
1140 - Dij are the measured RDCs or PCSs,
1141 - Dij(theta) are the back calculated RDCs or PCSs,
1142 - sigma_ij are the RDC or PCS errors,
1143 - dDij(theta)/dthetak is the RDC or PCS gradient for parameter k.
1144
1145
1146 The RDC gradient
1147 ----------------
1148
1149 The only parameters are the tensor components.
1150
1151 Amn partial derivative
1152 ~~~~~~~~~~~~~~~~~~~~~~
1153
1154 The alignment tensor element partial derivative is::
1155
1156 _N_
1157 dDij(theta) dj \ T dAi
1158 ----------- = -- > mu_jc . ---- . mu_jc,
1159 dAmn N /__ dAmn
1160 c=1
1161
1162 where:
1163 - dj is the dipolar constant for spin j,
1164 - N is the total number of states or structures,
1165 - mu_jc is the unit vector corresponding to spin j and state c,
1166 - dAi/dAmn is the partial derivative of the alignment tensor with respect to element Amn.
1167
1168
1169 The PCS gradient
1170 ----------------
1171
1172 Amn partial derivative
1173 ~~~~~~~~~~~~~~~~~~~~~~
1174
1175 The alignment tensor element partial derivative is::
1176
1177 _N_
1178 ddelta_ij(theta) 1 \ T dAi
1179 ---------------- = - > djc . mu_jc . ---- . mu_jc,
1180 dAmn N /__ dAmn
1181 c=1
1182
1183 where:
1184 - djc is the pseudocontact shift constant for spin j and state c,
1185 - N is the total number of states or structures,
1186 - mu_jc is the unit vector corresponding to spin j and state c,
1187 - dAi/dAmn is the partial derivative of the alignment tensor with respect to element Amn.
1188
1189
1190 The alignment tensor gradient
1191 -----------------------------
1192
1193 The five unique elements of the tensor {Axx, Ayy, Axy, Axz, Ayz} give five different partial derivatives. These are::
1194
1195 dAi | 1 0 0 |
1196 ---- = | 0 0 0 |,
1197 dAxx | 0 0 -1 |
1198
1199 dAi | 0 0 0 |
1200 ---- = | 0 1 0 |,
1201 dAyy | 0 0 -1 |
1202
1203 dAi | 0 1 0 |
1204 ---- = | 1 0 0 |,
1205 dAxy | 0 0 0 |
1206
1207 dAi | 0 0 1 |
1208 ---- = | 0 0 0 |,
1209 dAxz | 1 0 0 |
1210
1211 dAi | 0 0 0 |
1212 ---- = | 0 0 1 |.
1213 dAyz | 0 1 0 |
1214
1215 As these are invariant, they can be pre-calculated.
1216
1217
1218 Stored data structures
1219 ======================
1220
1221 There are a number of data structures calculated by this function and stored for subsequent use in the Hessian function. This include the back calculated RDC and PCS gradients and the alignment tensor gradients.
1222
1223 dDij(theta)/dthetak
1224 -------------------
1225
1226 The back calculated RDC gradient. This is a rank-3 tensor with indices {k, i, j}.
1227
1228 ddeltaij(theta)/dthetak
1229 -----------------------
1230
1231 The back calculated PCS gradient. This is a rank-3 tensor with indices {k, i, j}.
1232
1233 dAi/dAmn
1234 --------
1235
1236 The alignment tensor gradients. This is a rank-3 tensor with indices {5, n, m}.
1237
1238
1239 @param params: The vector of parameter values. This is unused as it is assumed that func_population() was called first.
1240 @type params: numpy rank-1 array
1241 @return: The chi-squared or SSE gradient.
1242 @rtype: numpy rank-1 array
1243 """
1244
1245
1246 if self.scaling_flag:
1247 params = dot(params, self.scaling_matrix)
1248
1249
1250 self.dchi2 = self.dchi2 * 0.0
1251
1252
1253 index = 0
1254 for i in xrange(self.num_align):
1255
1256 if self.fixed_tensors[i]:
1257 continue
1258
1259
1260 for j in xrange(self.num_spins):
1261
1262 if self.rdc_flag and not self.missing_Dij[i, j]:
1263 self.dDij_theta[index*5, index, j] = ave_rdc_tensor_dDij_dAmn(self.dip_const[j], self.dip_vect[j], self.N, self.dA[0], weights=self.probs)
1264 self.dDij_theta[index*5+1, index, j] = ave_rdc_tensor_dDij_dAmn(self.dip_const[j], self.dip_vect[j], self.N, self.dA[1], weights=self.probs)
1265 self.dDij_theta[index*5+2, index, j] = ave_rdc_tensor_dDij_dAmn(self.dip_const[j], self.dip_vect[j], self.N, self.dA[2], weights=self.probs)
1266 self.dDij_theta[index*5+3, index, j] = ave_rdc_tensor_dDij_dAmn(self.dip_const[j], self.dip_vect[j], self.N, self.dA[3], weights=self.probs)
1267 self.dDij_theta[index*5+4, index, j] = ave_rdc_tensor_dDij_dAmn(self.dip_const[j], self.dip_vect[j], self.N, self.dA[4], weights=self.probs)
1268
1269
1270 if self.pcs_flag and not self.missing_deltaij[i, j]:
1271 self.ddeltaij_theta[index*5, index, j] = ave_pcs_tensor_ddeltaij_dAmn(self.pcs_const[index, j], self.paramag_unit_vect[j], self.N, self.dA[0], weights=self.probs)
1272 self.ddeltaij_theta[index*5+1, index, j] = ave_pcs_tensor_ddeltaij_dAmn(self.pcs_const[index, j], self.paramag_unit_vect[j], self.N, self.dA[1], weights=self.probs)
1273 self.ddeltaij_theta[index*5+2, index, j] = ave_pcs_tensor_ddeltaij_dAmn(self.pcs_const[index, j], self.paramag_unit_vect[j], self.N, self.dA[2], weights=self.probs)
1274 self.ddeltaij_theta[index*5+3, index, j] = ave_pcs_tensor_ddeltaij_dAmn(self.pcs_const[index, j], self.paramag_unit_vect[j], self.N, self.dA[3], weights=self.probs)
1275 self.ddeltaij_theta[index*5+4, index, j] = ave_pcs_tensor_ddeltaij_dAmn(self.pcs_const[index, j], self.paramag_unit_vect[j], self.N, self.dA[4], weights=self.probs)
1276
1277
1278 for k in xrange(self.total_num_params):
1279
1280 if self.rdc_flag:
1281 self.dchi2[k] = self.dchi2[k] + dchi2_element(self.Dij[index], self.Dij_theta[index], self.dDij_theta[k, index], self.rdc_sigma_ij[index])
1282
1283
1284 if self.pcs_flag:
1285 self.dchi2[k] = self.dchi2[k] + dchi2_element(self.deltaij[index], self.deltaij_theta[index], self.ddeltaij_theta[k, index], self.pcs_sigma_ij[index])
1286
1287
1288 index += 1
1289
1290
1291 if self.scaling_flag:
1292 self.dchi2 = dot(self.dchi2, self.scaling_matrix)
1293
1294
1295 return self.dchi2
1296
1297
1299 """The Hessian function for optimisation of the flexible population N-state model.
1300
1301 Description
1302 ===========
1303
1304 This function should be passed to the optimisation algorithm. It accepts, as an array, a vector of parameter values and, using these, returns the chi-squared Hessian corresponding to that coordinate in the parameter space. If no RDC/PCS errors are supplied, then the SSE (the sum of squares error) Hessian is returned instead. The chi-squared Hessian is simply the SSE Hessian normalised to unit variance (the SSE divided by the error squared).
1305
1306
1307 Indices
1308 =======
1309
1310 For this calculation, six indices are looped over and used in the various data structures. These include:
1311 - k, the index over all parameters,
1312 - i, the index over alignments,
1313 - j, the index over spin systems,
1314 - c, the index over the N-states (or over the structures),
1315 - m, the index over the first dimension of the alignment tensor m = {x, y, z}.
1316 - n, the index over the second dimension of the alignment tensor n = {x, y, z},
1317
1318
1319 Equations
1320 =========
1321
1322 To calculate the chi-squared gradient, a chain of equations are used. This includes the chi-squared gradient, the RDC gradient and the alignment tensor gradient.
1323
1324
1325 The chi-squared Hessian
1326 -----------------------
1327
1328 The equation is::
1329 ___
1330 d2chi^2(theta) \ 1 / dDij(theta) dDij(theta) d2Dij(theta) \
1331 --------------- = 2 > ---------- | ----------- . ----------- - (Dij-Dij(theta)) . --------------- |.
1332 dthetaj.dthetak /__ sigma_i**2 \ dthetaj dthetak dthetaj.dthetak /
1333 ij
1334
1335 where:
1336 - theta is the parameter vector,
1337 - Dij are the measured RDCs or PCSs,
1338 - Dij(theta) are the back calculated RDCs or PCSs,
1339 - sigma_ij are the RDC or PCS errors,
1340 - dDij(theta)/dthetak is the RDC or PCS gradient for parameter k.
1341 - d2Dij(theta)/dthetaj.dthetak is the RDC or PCS Hessian for parameters j and k.
1342
1343
1344 The RDC Hessian
1345 ---------------
1346
1347 pc-pd second partial derivatives
1348 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1349
1350 The probability parameter second partial derivative is::
1351
1352 d2Dij(theta)
1353 ------------ = 0.
1354 dpc.dpd
1355
1356
1357 pc-Anm second partial derivatives
1358 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1359
1360 The probability parameter-tensor element second partial derivative is::
1361
1362 d2Dij(theta) T dAi
1363 ------------ = dj . mu_jc . ---- . mu_jc.
1364 dpc.dAmn dAmn
1365
1366
1367 Amn-Aop second partial derivatives
1368 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1369
1370 The alignment tensor element second partial derivative is::
1371
1372 d2Dij(theta)
1373 ------------ = 0.
1374 dAmn.dAop
1375
1376
1377 The PCS Hessian
1378 ---------------
1379
1380 pc-pd second partial derivatives
1381 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1382
1383 The probability parameter second partial derivative is::
1384
1385 d2delta_ij(theta)
1386 ----------------- = 0.
1387 dpc.dpd
1388
1389
1390 pc-Anm second partial derivatives
1391 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1392
1393 The probability parameter-tensor element second partial derivative is::
1394
1395 d2delta_ij(theta) T dAi
1396 ----------------- = djc . mu_jc . ---- . mu_jc.
1397 dpc.dAmn dAmn
1398
1399
1400 Amn-Aop second partial derivatives
1401 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1402
1403 The alignment tensor element second partial derivative is::
1404
1405 d2delta_ij(theta)
1406 ----------------- = 0
1407 dAmn.dAop
1408
1409
1410 The alignment tensor Hessian
1411 ----------------------------
1412
1413 The five unique elements of the tensor {Axx, Ayy, Axy, Axz, Ayz} all have the same second partial derivative of::
1414
1415 d2Ai | 0 0 0 |
1416 --------- = | 0 0 0 |.
1417 dAmn.dAop | 0 0 0 |
1418
1419
1420 @param params: The vector of parameter values. This is unused as it is assumed that func_population() was called first.
1421 @type params: numpy rank-1 array
1422 @return: The chi-squared or SSE Hessian.
1423 @rtype: numpy rank-2 array
1424 """
1425
1426
1427 if self.scaling_flag:
1428 params = dot(params, self.scaling_matrix)
1429
1430
1431 self.d2chi2 = self.d2chi2 * 0.0
1432
1433
1434 index = 0
1435 for i in xrange(self.num_align):
1436
1437 if self.fixed_tensors[i]:
1438 continue
1439
1440
1441 for c in xrange(self.N - 1):
1442
1443 pc_index = self.num_align_params + c
1444
1445
1446 for j in xrange(self.num_spins):
1447
1448 if self.fixed_tensors[i] and self.rdc_flag and not self.missing_Dij[i, j]:
1449 self.d2Dij_theta[pc_index, i*5+0, i, j] = self.d2Dij_theta[i*5+0, pc_index, i, j] = rdc_tensor(self.dip_const[j], self.dip_vect[j, c], self.dA[0])
1450 self.d2Dij_theta[pc_index, i*5+1, i, j] = self.d2Dij_theta[i*5+1, pc_index, i, j] = rdc_tensor(self.dip_const[j], self.dip_vect[j, c], self.dA[1])
1451 self.d2Dij_theta[pc_index, i*5+2, i, j] = self.d2Dij_theta[i*5+2, pc_index, i, j] = rdc_tensor(self.dip_const[j], self.dip_vect[j, c], self.dA[2])
1452 self.d2Dij_theta[pc_index, i*5+3, i, j] = self.d2Dij_theta[i*5+3, pc_index, i, j] = rdc_tensor(self.dip_const[j], self.dip_vect[j, c], self.dA[3])
1453 self.d2Dij_theta[pc_index, i*5+4, i, j] = self.d2Dij_theta[i*5+4, pc_index, i, j] = rdc_tensor(self.dip_const[j], self.dip_vect[j, c], self.dA[4])
1454
1455
1456 if self.fixed_tensors[i] and self.pcs_flag and not self.missing_deltaij[i, j]:
1457 self.d2deltaij_theta[pc_index, i*5+0, i, j] = self.d2deltaij_theta[i*5+0, pc_index, i, j] = pcs_tensor(self.pcs_const[i, j, c], self.paramag_unit_vect[j, c], self.dA[0])
1458 self.d2deltaij_theta[pc_index, i*5+1, i, j] = self.d2deltaij_theta[i*5+1, pc_index, i, j] = pcs_tensor(self.pcs_const[i, j, c], self.paramag_unit_vect[j, c], self.dA[1])
1459 self.d2deltaij_theta[pc_index, i*5+2, i, j] = self.d2deltaij_theta[i*5+2, pc_index, i, j] = pcs_tensor(self.pcs_const[i, j, c], self.paramag_unit_vect[j, c], self.dA[2])
1460 self.d2deltaij_theta[pc_index, i*5+3, i, j] = self.d2deltaij_theta[i*5+3, pc_index, i, j] = pcs_tensor(self.pcs_const[i, j, c], self.paramag_unit_vect[j, c], self.dA[3])
1461 self.d2deltaij_theta[pc_index, i*5+4, i, j] = self.d2deltaij_theta[i*5+4, pc_index, i, j] = pcs_tensor(self.pcs_const[i, j, c], self.paramag_unit_vect[j, c], self.dA[4])
1462
1463
1464 index += 1
1465
1466
1467 for i in xrange(self.num_align):
1468
1469 for j in xrange(self.total_num_params):
1470 for k in xrange(self.total_num_params):
1471
1472 if self.fixed_tensors[i] and self.rdc_flag:
1473 self.d2chi2[j, k] = self.d2chi2[j, k] + d2chi2_element(self.Dij[i], self.Dij_theta[i], self.dDij_theta[j, i], self.dDij_theta[k, i], self.d2Dij_theta[j, k, i], self.rdc_sigma_ij[i])
1474
1475
1476 if self.fixed_tensors[i] and self.pcs_flag:
1477 self.d2chi2[j, k] = self.d2chi2[j, k] + d2chi2_element(self.deltaij[i], self.deltaij_theta[i], self.ddeltaij_theta[j, i], self.ddeltaij_theta[k, i], self.d2deltaij_theta[j, k, i], self.pcs_sigma_ij[i])
1478
1479
1480 if self.scaling_flag:
1481 self.d2chi2 = dot(self.d2chi2, self.scaling_matrix)
1482
1483
1484 return self.d2chi2
1485
1486
1488 """The Hessian function for optimisation of the alignment tensor from RDC and/or PCS data.
1489
1490 Description
1491 ===========
1492
1493 This function should be passed to the optimisation algorithm. It accepts, as an array, a vector of parameter values and, using these, returns the chi-squared Hessian corresponding to that coordinate in the parameter space. If no RDC/PCS errors are supplied, then the SSE (the sum of squares error) Hessian is returned instead. The chi-squared Hessian is simply the SSE Hessian normalised to unit variance (the SSE divided by the error squared).
1494
1495 Indices
1496 =======
1497
1498 For this calculation, six indices are looped over and used in the various data structures. These include:
1499 - k, the index over all parameters,
1500 - i, the index over alignments,
1501 - j, the index over spin systems,
1502 - c, the index over the N-states (or over the structures),
1503 - m, the index over the first dimension of the alignment tensor m = {x, y, z}.
1504 - n, the index over the second dimension of the alignment tensor n = {x, y, z},
1505
1506
1507 Equations
1508 =========
1509
1510 To calculate the chi-squared gradient, a chain of equations are used. This includes the chi-squared gradient, the RDC gradient and the alignment tensor gradient.
1511
1512
1513 The chi-squared Hessian
1514 -----------------------
1515
1516 The equation is::
1517 ___
1518 d2chi^2(theta) \ 1 / dDij(theta) dDij(theta) d2Dij(theta) \
1519 --------------- = 2 > ---------- | ----------- . ----------- - (Dij-Dij(theta)) . --------------- |.
1520 dthetaj.dthetak /__ sigma_i**2 \ dthetaj dthetak dthetaj.dthetak /
1521 ij
1522
1523 where:
1524 - theta is the parameter vector,
1525 - Dij are the measured RDCs or PCSs,
1526 - Dij(theta) are the back calculated RDCs or PCSs,
1527 - sigma_ij are the RDC or PCS errors,
1528 - dDij(theta)/dthetak is the RDC or PCS gradient for parameter k.
1529 - d2Dij(theta)/dthetaj.dthetak is the RDC or PCS Hessian for parameters j and k.
1530
1531
1532 The RDC Hessian
1533 ---------------
1534
1535 The only parameters are the tensor components.
1536
1537
1538 Amn-Aop second partial derivatives
1539 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1540
1541 The alignment tensor element second partial derivative is::
1542
1543 d2Dij(theta)
1544 ------------ = 0.
1545 dAmn.dAop
1546
1547
1548 The PCS Hessian
1549 ---------------
1550
1551 Amn-Aop second partial derivatives
1552 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1553
1554 The alignment tensor element second partial derivative is::
1555
1556 d2delta_ij(theta)
1557 ----------------- = 0
1558 dAmn.dAop
1559
1560
1561 The alignment tensor Hessian
1562 ----------------------------
1563
1564 The five unique elements of the tensor {Axx, Ayy, Axy, Axz, Ayz} all have the same second partial derivative of::
1565
1566 d2Ai | 0 0 0 |
1567 --------- = | 0 0 0 |.
1568 dAmn.dAop | 0 0 0 |
1569
1570
1571 @param params: The vector of parameter values. This is unused as it is assumed that func_population() was called first.
1572 @type params: numpy rank-1 array
1573 @return: The chi-squared or SSE Hessian.
1574 @rtype: numpy rank-2 array
1575 """
1576
1577
1578 if self.scaling_flag:
1579 params = dot(params, self.scaling_matrix)
1580
1581
1582 self.d2chi2 = self.d2chi2 * 0.0
1583
1584
1585 index = 0
1586 for i in xrange(self.num_align):
1587
1588 if self.fixed_tensors[i]:
1589 continue
1590
1591
1592 for j in xrange(self.total_num_params):
1593 for k in xrange(self.total_num_params):
1594
1595 if self.rdc_flag:
1596 self.d2chi2[j, k] = self.d2chi2[j, k] + d2chi2_element(self.Dij[i], self.Dij_theta[i], self.dDij_theta[j, i], self.dDij_theta[k, i], self.zero_hessian, self.rdc_sigma_ij[i])
1597
1598
1599 if self.pcs_flag:
1600 self.d2chi2[j, k] = self.d2chi2[j, k] + d2chi2_element(self.deltaij[i], self.deltaij_theta[i], self.ddeltaij_theta[j, i], self.ddeltaij_theta[k, i], self.zero_hessian, self.pcs_sigma_ij[i])
1601
1602
1603 index += 1
1604
1605
1606 if self.scaling_flag:
1607 self.d2chi2 = dot(self.d2chi2, self.scaling_matrix)
1608
1609
1610 return self.d2chi2
1611
1612
1614 """Calculate the paramagnetic centre to spin vectors, distances and constants."""
1615
1616
1617 if rank(self.paramag_centre) == 1:
1618 vectors_single_centre(self.atomic_pos, self.paramag_centre, self.paramag_unit_vect, self.paramag_dist)
1619 else:
1620 vectors_centre_per_state(self.atomic_pos, self.paramag_centre, self.paramag_unit_vect, self.paramag_dist)
1621
1622
1623 for i in range(self.num_align):
1624 for j in range(self.num_spins):
1625 for c in range(self.N):
1626 self.pcs_const[i, j, c] = pcs_constant(self.temp[i], self.frq[i], self.paramag_dist[j, c])
1627