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, rdc_vect=None, T_flags=None, j_couplings=None, rdc_pseudo_flags=None, pcs_pseudo_flags=None, temp=None, frq=None, dip_const=None, absolute_rdc=None, atomic_pos=None, paramag_centre=None, scaling_matrix=None, centre_fixed=True)
(Constructor)

Set up the class instance for optimisation.

The N-state models

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:

model, the type of N-state model. This can be '2-domain', 'population', or 'fixed'.
N, the number of states (or structures).
init_params, the initial parameter values.
scaling_matrix, the matrix used for parameter scaling during optimisation.

2-domain N-state model

If the model type is set to '2-domain', then the following data structures should be supplied:

full_tensors, the alignment tensors in matrix form.
red_data, the alignment tensors in 5D form in a rank-1 array.
red_errors, the alignment tensor errors in 5D form in a rank-1 array. This data is not obligatory.
full_in_ref_frame, an array of flags specifying if the tensor in the reference frame is the full or reduced tensor.

The population N-state model

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).

PCS base data

If pseudocontact shift data is to be used for optimisation, then the following should be supplied:

pcs, the pseudocontact shifts.
pcs_errors, the optional pseudocontact shift error values.
temp, the temperatures of the experiments.
frq, the frequencies of the experiments.

PCS and PRE base data

If either pseudocontact shift or PRE data is to be used for optimisation, then the following should be supplied:

atomic_pos, the positions of all atoms.
paramag_centre, the paramagnetic centre position.

RDC base data

If residual dipolar coupling data is to be used for optimisation, then the following should be supplied:

rdcs, the residual dipolar couplings.
rdc_errors, the optional residual dipolar coupling errors.
rdc_vect, the interatomic vectors.
T_flags, the array of flags specifying if the data for the given alignment is of the T = J+D type.
j_couplings, the J couplings for the case when RDC data is of the type T = J+D.
dip_const, the dipolar constants.
absolute_rdc, the flags specifying whether each RDC is signless.

Parameters:

model (str) - The N-state model type. This can be one of '2-domain', 'population' or 'fixed'.
N (int) - The number of states.
init_params (numpy float64 array) - The initial parameter values. Optimisation must start at some point!
probs (list of float) - The probabilities for each state c. The length of this list should be equal to N.
full_tensors (list of rank-2, 3D numpy arrays) - 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]
red_data (numpy float64 array) - An array of the {Sxx, Syy, Sxy, Sxz, Syz} values for all reduced tensors. The format is the same as for full_tensors.
red_errors (numpy float64 array) - An array of the {Sxx, Syy, Sxy, Sxz, Syz} errors for all reduced tensors. The array format is the same as for full_tensors.
full_in_ref_frame (numpy rank-1 array) - An array of flags specifying if the tensor in the reference frame is the full or reduced tensor.
fixed_tensors (list of bool) - An array of flags specifying if the tensor is fixed or will be optimised.
pcs (numpy rank-2 array) - The PCS lists. The first index must correspond to the different alignment media i and the second index to the spin systems j.
pcs_errors (numpy rank-2 array) - The PCS error lists. The dimensions of this argument are the same as for 'pcs'.
pcs_weights (numpy rank-2 array) - The PCS weight lists. The dimensions of this argument are the same as for 'pcs'.
rdcs (numpy rank-2 array) - The RDC lists. The first index must correspond to the different alignment media i and the second index to the spin pairs k.
rdc_errors (numpy rank-2 array) - The RDC error lists. The dimensions of this argument are the same as for 'rdcs'.
rdc_weights (numpy rank-2 array) - The RDC weight lists. The dimensions of this argument are the same as for 'rdcs'.
rdc_vect (list of numpy rank-3 array) - The unit vector lists for the RDC. The first index must correspond to the spin pair, the second index to each structure (its size being equal to the number of states) and the third to each pseudo-atom present.
T_flags (numpy rank-2 int32 array) - The array of flags specifying if the data for the given alignment is of the T = J+D type.
j_couplings (numpy rank-1 array) - The J couplings list, used when the RDC data is of the type T = J+D. The number of elements must match the second dimension of the rdcs argument.
rdc_pseudo_flags (numpy rank-1 int32 array) - The array of flags specifying if one of the atoms of the interatomic pair for the RDC are pseudo-atoms. The indices correspond to the interatomic pairs.
pcs_pseudo_flags (numpy rank-1 int32 array) - The array of flags specifying if one of the atoms of the interatomic pair for the PCS are pseudo-atoms. The indices correspond to the interatomic pairs.
temp (numpy rank-1 array) - The temperature of each experiment, used for the PCS.
frq (numpy rank-1 array) - The frequency of each alignment, used for the PCS.
dip_const (list of lists of floats) - The dipolar constants for each spin pair. The first index corresponds to the spin pairs k and the second to each pseudo-atom present.
absolute_rdc (numpy rank-2 int32 array) - The signless or absolute value flags for the RDCs.
atomic_pos (numpy rank-3 array) - The atomic positions of all spins for the PCS and PRE data. The first index is the spin systems j and the second is the structure or state c.
paramag_centre (numpy rank-1, 3D array or rank-2, Nx3 array) - The paramagnetic centre position (or positions).
scaling_matrix (numpy rank-2 array) - The square and diagonal scaling matrix.
centre_fixed (bool) - A flag which if False will cause the paramagnetic centre to be optimised.

func_standard(self, params)

source code

The target function for optimisation of the standard N-state model.

Description

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 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).

Indices

For this calculation, five indices are looped over and used in the various data structures. These include:

i, the index over alignments,
j, the index over spin systems,
c, the index over the N-states (or over the structures),
n, the index over the first dimension of the alignment tensor n = {x, y, z},
m, the index over the second dimension of the alignment tensor m = {x, y, z}.

Equations

To calculate the function value, a chain of equations are used. This includes the chi-squared equation and the RDC and PCS equations.

The chi-squared equation

The equations are:

                ___
                \    (Dij - Dij(theta)) ** 2
chi^2(theta)  =  >   ----------------------- ,
                /__       sigma_ij ** 2
                 ij

                ___
                \    (delta_ij - delta_ij(theta)) ** 2
chi^2(theta)  =  >   --------------------------------- ,
                /__             sigma_ij ** 2
                 ij

where:

theta is the parameter vector,
Dij are the measured RDCs for alignment i, spin j,
Dij(theta) are the back calculated RDCs for alignment i, spin j,
delta_ij are the measured PCSs for alignment i, spin j,
delta_ij(theta) are the back calculated PCSs for alignment i, spin j,
sigma_ij are the RDC or PCS errors.

Both chi-squared values sum.

The RDC equation

The RDC equation is:

                  _N_
                  \              T
Dij(theta)  =  dj  >   pc . mu_jc . Ai . mu_jc,
                  /__
                  c=1

where:

dj is the dipolar constant for spin j,
N is the total number of states or structures,
pc is the weight or probability associated with state c,
mu_jc is the unit vector corresponding to spin j and state c,
Ai is the alignment tensor.

In the fixed and equal probability case, the equation is:

                  _N_
               dj \         T
Dij(theta)  =  --  >   mu_jc . Ai . mu_jc,
               N  /__
                  c=1

The dipolar constant is henceforth defined as:

   dj = 3 / (2pi) d',

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:

          mu0 gI.gS.h_bar
   d' = - --- ----------- ,
          4pi    r**3

where:

mu0 is the permeability of free space,
gI and gS are the gyromagnetic ratios of the I and S spins,
h_bar is Dirac's constant which is equal to Planck's constant divided by 2pi,
r is the distance between the two spins.

The PCS equation

The PCS equation is:

                      _N_
                      \                    T
   delta_ij(theta)  =  >  pc . dijc . mu_jc . Ai . mu_jc,
                      /__
                      c=1

where:

djci is the PCS constant for spin j, state c and experiment or alignment i,
N is the total number of states or structures,
pc is the weight or probability associated with state c,
mu_jc is the unit vector corresponding to spin j and state c,
Ai is the alignment tensor.

In the fixed and equal probability case, the equation is:

                        _N_
                      1 \               T
   delta_ij(theta)  = -  >  dijc . mu_jc . Ai . mu_jc,
                      N /__
                        c=1

The PCS constant is defined as:

          mu0 15kT   1
   dijc = --- ----- ---- ,
          4pi Bo**2 r**3

where:

mu0 is the permeability of free space,
k is Boltzmann's constant,
T is the absolute temperature (different for each experiment),
Bo is the magnetic field strength (different for each experiment),
r is the distance between the paramagnetic centre (electron spin) and the nuclear spin (different for each spin and state).

Stored data structures

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 PCSs and the alignment tensors.

Dij(theta)

The back calculated RDCs. This is a rank-2 tensor with indices {i, j}.

delta_ij(theta)

The back calculated PCS. This is a rank-2 tensor with indices {i, j}.

Ai

The alignment tensors. This is a rank-3 tensor with indices {i, n, m}.

Parameters:

params (numpy rank-1 array) - The vector of parameter values.

Returns: float

The chi-squared or SSE value.

dfunc_standard(self, params)

source code

The gradient function for optimisation of the standard N-state model.

Description

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 or PCS 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).

Indices

For this calculation, six indices are looped over and used in the various data structures. These include:

k, the index over all parameters,
i, the index over alignments,
j, the index over spin systems,
c, the index over the N-states (or over the structures),
m, the index over the first dimension of the alignment tensor m = {x, y, z}.
n, the index over the second dimension of the alignment tensor n = {x, y, z},

Equations

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.

The chi-squared gradient

The equation is:

                     ___
dchi^2(theta)        \   / Dij - Dij(theta)     dDij(theta) \ 
-------------  =  -2  >  | ----------------  .  ----------- |
   dthetak           /__ \   sigma_ij**2         dthetak    /
                     ij

where:

theta is the parameter vector,
Dij are the measured RDCs or PCSs,
Dij(theta) are the back calculated RDCs or PCSs,
sigma_ij are the RDC or PCS errors,
dDij(theta)/dthetak is the RDC or PCS gradient for parameter k.

The RDC gradient

This gradient is different for the various parameter types.

pc partial derivative

The population parameter partial derivative is:

    dDij(theta)               T
    -----------  =  dj . mu_jc . Ai . mu_jc,
        dpc

where:

dj is the dipolar constant for spin j,
mu_jc is the unit vector corresponding to spin j and state c,
Ai is the alignment tensor.

Amn partial derivative

The alignment tensor element partial derivative is:

                   _N_
dDij(theta)        \              T   dAi
-----------  =  dj  >   pc . mu_jc . ---- . mu_jc,
   dAmn            /__               dAmn
                   c=1

where:

dj is the dipolar constant for spin j,
pc is the weight or probability associated with state c,
mu_jc is the unit vector corresponding to spin j and state c,
dAi/dAmn is the partial derivative of the alignment tensor with respect to element Amn.

In the case of fixed and equal populations, the equation is:

                   _N_
dDij(theta)     dj \         T   dAi
-----------  =  --  >   mu_jc . ---- . mu_jc,
   dAmn         N  /__          dAmn
                   c=1

The PCS gradient

This gradient is also different for the various parameter types.

pc partial derivative

The population parameter partial derivative is:

    ddeltaij(theta)                 T
    ---------------  =  dijc . mu_jc . Ai . mu_jc,
         dpc

where:

djc is the pseudocontact shift constant for spin j and state c,
mu_jc is the unit vector corresponding to spin j and state c,
Ai is the alignment tensor.

Amn partial derivative

The alignment tensor element partial derivative is:

                       _N_
   ddelta_ij(theta)    \                   T   dAi
   ----------------  =  >  pc . djc . mu_jc . ---- . mu_jc,
         dAmn          /__                    dAmn
                       c=1

where:

djc is the pseudocontact shift constant for spin j and state c,
pc is the weight or probability associated with state c,
mu_jc is the unit vector corresponding to spin j and state c,
dAi/dAmn is the partial derivative of the alignment tensor with respect to element Amn.

In the case of fixed and equal populations, the equation is:

                         _N_
   ddelta_ij(theta)    1 \              T   dAi
   ----------------  = -  >  djc . mu_jc . ---- . mu_jc,
         dAmn          N /__               dAmn
                         c=1

xi partial derivative

The paramagnetic position partial derivative is:

                       _N_
   ddelta_ij(theta)    \        / ddjc                       dr_jcT                          dr_jc \ 
   ----------------  =  >  pc . | ----.r_jcT.Ai.r_jc  +  djc.------.Ai.r_jc  +  djc.r_jcT.Ai.----- | ,
         dxi           /__      \ dxi                         dxi                             dxi  /
                       c=1

where xi are the paramagnetic position coordinates {x0, x1, x2} and the last two terms in the sum are equal due to the symmetry of the alignment tensor, and:

   ddjc    mu0 15kT                 5 (si - xi)
   ----  = --- ----- ---------------------------------------------  ,
   dxi     4pi Bo**2 ((sx-x0)**2 + (sy-x1)**2 + (sz-x2)**2)**(7/2)

and:

   dr      | 1 |   dr      | 0 |   dr      | 0 |
   --  = - | 0 | , --  = - | 1 | , --  = - | 0 | .
   dx      | 0 |   dy      | 0 |   dy      | 1 |

The pseudocontact shift constant is defined here as:

         mu0 15kT    1
   djc = --- ----- ------ ,
         4pi Bo**2 rjc**5

The alignment tensor gradient

The five unique elements of the tensor {Axx, Ayy, Axy, Axz, Ayz} give five different partial derivatives. These are:

    dAi   | 1  0  0 |
   ---- = | 0  0  0 |,
   dAxx   | 0  0 -1 |

    dAi   | 0  0  0 |
   ---- = | 0  1  0 |,
   dAyy   | 0  0 -1 |

    dAi   | 0  1  0 |
   ---- = | 1  0  0 |,
   dAxy   | 0  0  0 |

    dAi   | 0  0  1 |
   ---- = | 0  0  0 |,
   dAxz   | 1  0  0 |

    dAi   | 0  0  0 |
   ---- = | 0  0  1 |.
   dAyz   | 0  1  0 |

As these are invariant, they can be pre-calculated.

Stored data structures

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.

dDij(theta)/dthetak

The back calculated RDC gradient. This is a rank-3 tensor with indices {k, i, j}.

ddeltaij(theta)/dthetak

The back calculated PCS gradient. This is a rank-3 tensor with indices {k, i, j}.

dAi/dAmn

The alignment tensor gradients. This is a rank-3 tensor with indices {5, n, m}.

Parameters:

params (numpy rank-1 array) - The vector of parameter values. This is unused as it is assumed that func() was called first.

Returns: numpy rank-1 array

The chi-squared or SSE gradient.

d2func_standard(self, params)

source code

The Hessian function for optimisation of the standard N-state model.

Description

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).

Indices

For this calculation, six indices are looped over and used in the various data structures. These include:

k, the index over all parameters,
i, the index over alignments,
j, the index over spin systems,
c, the index over the N-states (or over the structures),
m, the index over the first dimension of the alignment tensor m = {x, y, z}.
n, the index over the second dimension of the alignment tensor n = {x, y, z},

Equations

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.

The chi-squared Hessian

The equation is:

                      ___
d2chi^2(theta)        \       1      / dDij(theta)   dDij(theta)                         d2Dij(theta)   \ 
---------------  =  2  >  ---------- | ----------- . -----------  -  (Dij-Dij(theta)) . --------------- |.
dthetaj.dthetak       /__ sigma_i**2 \  dthetaj       dthetak                           dthetaj.dthetak /
                      ij

where:

theta is the parameter vector,
Dij are the measured RDCs or PCSs,
Dij(theta) are the back calculated RDCs or PCSs,
sigma_ij are the RDC or PCS errors,
dDij(theta)/dthetak is the RDC or PCS gradient for parameter k.
d2Dij(theta)/dthetaj.dthetak is the RDC or PCS Hessian for parameters j and k.

The RDC Hessian

pc-pd second partial derivatives

The probability parameter second partial derivative is:

d2Dij(theta)
------------  =  0.
  dpc.dpd

pc-Anm second partial derivatives

The probability parameter-tensor element second partial derivative is:

d2Dij(theta)               T   dAi
------------  =  dj . mu_jc . ---- . mu_jc.
  dpc.dAmn                    dAmn

Amn-Aop second partial derivatives

The alignment tensor element second partial derivative is:

d2Dij(theta)
------------  =  0.
 dAmn.dAop

The PCS Hessian

pc-pd second partial derivatives

The probability parameter second partial derivative is:

d2delta_ij(theta)
-----------------  =  0.
     dpc.dpd

pc-Anm second partial derivatives

The probability parameter-tensor element second partial derivative is:

d2delta_ij(theta)                T   dAi
-----------------  =  djc . mu_jc . ---- . mu_jc.
    dpc.dAmn                        dAmn

Amn-Aop second partial derivatives

The alignment tensor element second partial derivative is:

   d2delta_ij(theta)
   -----------------  =  0
       dAmn.dAop

The alignment tensor Hessian

The five unique elements of the tensor {Axx, Ayy, Axy, Axz, Ayz} all have the same second partial derivative of:

     d2Ai      | 0  0  0 |
   --------- = | 0  0  0 |.
   dAmn.dAop   | 0  0  0 |

Parameters:

params (numpy rank-1 array) - The vector of parameter values. This is unused as it is assumed that func() was called first.

Returns: numpy rank-2 array

The chi-squared or SSE Hessian.

Class N_state_opt

The N-state models

2-domain N-state model

The population N-state model

PCS base data

PCS and PRE base data

RDC base data

func_2domain(self, params)

func_standard(self, params)

Description

Indices

Equations

The chi-squared equation

The RDC equation

The PCS equation

Stored data structures

Dij(theta)

delta_ij(theta)

Ai

dfunc_standard(self, params)

Description

Indices

Equations

The chi-squared gradient

The RDC gradient

pc partial derivative

Amn partial derivative

The PCS gradient

pc partial derivative

Amn partial derivative

xi partial derivative

The alignment tensor gradient

Stored data structures

dDij(theta)/dthetak

ddeltaij(theta)/dthetak

dAi/dAmn

d2func_standard(self, params)

Description

Indices

Equations

The chi-squared Hessian

The RDC Hessian

pc-pd second partial derivatives

pc-Anm second partial derivatives

Amn-Aop second partial derivatives

The PCS Hessian

pc-pd second partial derivatives

pc-Anm second partial derivatives

Amn-Aop second partial derivatives

The alignment tensor Hessian