Author: bugman
Date: Wed Sep 24 17:01:32 2014
New Revision: 26015
URL: http://svn.gna.org/viewcvs/relax?rev=26015&view=rev
Log:
Implemented the SciPy quadratic integration target function for the double
rotor frame order model.
This simply follows from what all the other quadratic integration target
functions and
lib.frame_order module functions do.
Modified:
branches/frame_order_cleanup/lib/frame_order/double_rotor.py
branches/frame_order_cleanup/target_functions/frame_order.py
Modified: branches/frame_order_cleanup/lib/frame_order/double_rotor.py
URL:
http://svn.gna.org/viewcvs/relax/branches/frame_order_cleanup/lib/frame_order/double_rotor.py?rev=26015&r1=26014&r2=26015&view=diff
==============================================================================
--- branches/frame_order_cleanup/lib/frame_order/double_rotor.py
(original)
+++ branches/frame_order_cleanup/lib/frame_order/double_rotor.py Wed
Sep 24 17:01:32 2014
@@ -23,8 +23,12 @@
"""Module for the double rotor frame order model."""
# Python module imports.
-from math import pi
+from math import cos, pi, sin
from numpy import add, divide, dot, eye, float64, multiply, sinc, swapaxes,
tensordot
+try:
+ from scipy.integrate import dblquad
+except ImportError:
+ pass
# relax module imports.
from lib.compat import norm
@@ -181,6 +185,57 @@
divide(pcs_theta, float(num), pcs_theta)
+def pcs_numeric_quad_int_double_rotor(sigma_max=None, sigma_max_2=None,
c=None, r_pivot_atom=None, r_ln_pivot=None, r_inter_pivot=None, A=None,
R_eigen=None, RT_eigen=None, Ri_prime=None, Ri2_prime=None):
+ """Determine the averaged PCS value via SciPy quadratic numerical
integration.
+
+ @keyword sigma_max: The maximum opening angle for the first
rotor.
+ @type sigma_max: float
+ @keyword sigma_max_2: The maximum opening angle for the second
rotor.
+ @type sigma_max_2: float
+ @keyword c: The PCS constant (without the interatomic
distance and in Angstrom units).
+ @type c: numpy rank-1 array
+ @keyword r_pivot_atom: The pivot point to atom vector.
+ @type r_pivot_atom: numpy rank-1, 3D array
+ @keyword r_ln_pivot: The lanthanide position to pivot point
vector.
+ @type r_ln_pivot: numpy rank-1, 3D array
+ @keyword r_inter_pivot: The vector between the two pivots.
+ @type r_inter_pivot: numpy rank-1, 3D array
+ @keyword A: The full alignment tensor of the non-moving
domain.
+ @type A: numpy rank-2, 3D array
+ @keyword R_eigen: The eigenframe rotation matrix.
+ @type R_eigen: numpy rank-2, 3D array
+ @keyword RT_eigen: The transpose of the eigenframe rotation
matrix (for faster calculations).
+ @type RT_eigen: numpy rank-2, 3D array
+ @keyword Ri_prime: The array of pre-calculated rotation
matrices for the in-frame double rotor motion for the 1st mode of motion,
used to calculate the PCS for each state i in the numerical integration.
+ @type Ri_prime: numpy rank-2, 3D array
+ @keyword Ri2_prime: The array of pre-calculated rotation
matrices for the in-frame double rotor motion for the 2nd mode of motion,
used to calculate the PCS for each state i in the numerical integration.
+ @type Ri2_prime: numpy rank-2, 3D array
+ """
+
+ # Preset the 1st rotation matrix elements for state i.
+ Ri_prime[0, 1] = 0.0
+ Ri_prime[1, 0] = 0.0
+ Ri_prime[1, 1] = 1.0
+ Ri_prime[1, 2] = 0.0
+ Ri_prime[2, 1] = 0.0
+
+ # Preset the 2nd rotation matrix elements for state i.
+ Ri2_prime[0, 0] = 1.0
+ Ri2_prime[0, 1] = 0.0
+ Ri2_prime[0, 2] = 0.0
+ Ri2_prime[1, 0] = 0.0
+ Ri2_prime[2, 0] = 0.0
+
+ # Perform numerical integration.
+ result = dblquad(pcs_pivot_motion_double_rotor_quad_int, -sigma_max,
sigma_max, lambda sigma2: -sigma_max_2, lambda sigma2: sigma_max_2,
args=(r_pivot_atom, r_ln_pivot, r_inter_pivot, A, R_eigen, RT_eigen,
Ri_prime, Ri2_prime))
+
+ # The surface area normalisation factor.
+ SA = 4.0 * sigma_max * sigma_max_2
+
+ # Return the value.
+ return c * result[0] / SA
+
+
def pcs_pivot_motion_double_rotor_qr_int(full_in_ref_frame=None,
r_pivot_atom=None, r_pivot_atom_rev=None, r_ln_pivot=None,
r_inter_pivot=None, A=None, Ri=None, Ri2=None, pcs_theta=None,
pcs_theta_err=None, missing_pcs=None):
"""Calculate the PCS value after a pivoted motion for the double rotor
model.
@@ -249,3 +304,75 @@
# The PCS.
pcs_theta[i, j] += proj * length_i
+
+
+def pcs_pivot_motion_double_rotor_quad_int(sigma_i, sigma2_i, r_pivot_atom,
r_ln_pivot, r_inter_pivot, A, R_eigen, RT_eigen, Ri_prime, Ri2_prime):
+ """Calculate the PCS value after a pivoted motion for the double rotor
model.
+
+ @param sigma_i: The 1st torsion angle for state i.
+ @type sigma_i: float
+ @param sigma2_i: The 1st torsion angle for state i.
+ @type sigma2_i: float
+ @param r_pivot_atom: The pivot point to atom vector.
+ @type r_pivot_atom: numpy rank-2, 3D array
+ @param r_ln_pivot: The lanthanide position to pivot point
vector.
+ @type r_ln_pivot: numpy rank-2, 3D array
+ @param r_inter_pivot: The vector between the two pivots.
+ @type r_inter_pivot: numpy rank-1, 3D array
+ @param A: The full alignment tensor of the non-moving
domain.
+ @type A: numpy rank-2, 3D array
+ @param R_eigen: The eigenframe rotation matrix.
+ @type R_eigen: numpy rank-2, 3D array
+ @param RT_eigen: The transpose of the eigenframe rotation
matrix (for faster calculations).
+ @type RT_eigen: numpy rank-2, 3D array
+ @param Ri_prime: The empty rotation matrix for state i.
+ @type Ri_prime: numpy rank-2, 3D array
+ @param Ri2_prime: The 2nd empty rotation matrix for state i.
+ @type Ri2_prime: numpy rank-2, 3D array
+ @return: The PCS value for the changed position.
+ @rtype: float
+ """
+
+ # The 1st rotation matrix.
+ c_sigma = cos(sigma_i)
+ s_sigma = sin(sigma_i)
+ Ri_prime[0, 0] = c_sigma
+ Ri_prime[0, 2] = s_sigma
+ Ri_prime[2, 0] = -s_sigma
+ Ri_prime[2, 2] = c_sigma
+
+ # The 2nd rotation matrix.
+ c_sigma = cos(sigma2_i)
+ s_sigma = sin(sigma2_i)
+ Ri2_prime[1, 1] = c_sigma
+ Ri2_prime[1, 2] = -s_sigma
+ Ri2_prime[2, 1] = s_sigma
+ Ri2_prime[2, 2] = c_sigma
+
+ # The rotations.
+ Ri = dot(R_eigen, dot(Ri_prime, RT_eigen))
+ Ri2 = dot(R_eigen, dot(Ri2_prime, RT_eigen))
+
+ # Rotate the first pivot to atomic position vectors.
+ rot_vect = dot(r_pivot_atom, Ri)
+
+ # Add the inter-pivot vector to obtain the 2nd pivot to atomic position
vectors.
+ add(r_inter_pivot, rot_vect, rot_vect)
+
+ # Rotate the 2nd pivot to atomic position vectors.
+ rot_vect = dot(rot_vect, Ri2)
+
+ # Add the lanthanide to pivot vector.
+ add(rot_vect, r_ln_pivot, rot_vect)
+
+ # The vector length.
+ length = norm(rot_vect)
+
+ # The projection.
+ proj = dot(rot_vect, dot(A, rot_vect))
+
+ # The PCS.
+ pcs = proj / length**5
+
+ # Return the PCS value (without the PCS constant).
+ return pcs
Modified: branches/frame_order_cleanup/target_functions/frame_order.py
URL:
http://svn.gna.org/viewcvs/relax/branches/frame_order_cleanup/target_functions/frame_order.py?rev=26015&r1=26014&r2=26015&view=diff
==============================================================================
--- branches/frame_order_cleanup/target_functions/frame_order.py
(original)
+++ branches/frame_order_cleanup/target_functions/frame_order.py Wed
Sep 24 17:01:32 2014
@@ -36,7 +36,7 @@
from lib.errors import RelaxError
from lib.float import isNaN
from lib.frame_order.conversions import create_rotor_axis_alpha
-from lib.frame_order.double_rotor import compile_2nd_matrix_double_rotor,
pcs_numeric_qr_int_double_rotor
+from lib.frame_order.double_rotor import compile_2nd_matrix_double_rotor,
pcs_numeric_qr_int_double_rotor, pcs_numeric_quad_int_double_rotor
from lib.frame_order.free_rotor import compile_2nd_matrix_free_rotor
from lib.frame_order.iso_cone import compile_2nd_matrix_iso_cone,
pcs_numeric_quad_int_iso_cone, pcs_numeric_qr_int_iso_cone
from lib.frame_order.iso_cone_free_rotor import
compile_2nd_matrix_iso_cone_free_rotor
@@ -365,6 +365,7 @@
self.R_eigen_2 = zeros((3, 3), float64)
self.R_ave = zeros((3, 3), float64)
self.Ri_prime = zeros((3, 3), float64)
+ self.Ri2_prime = zeros((3, 3), float64)
self.tensor_3D = zeros((3, 3), float64)
# The cone axis storage and molecular frame z-axis.
@@ -457,6 +458,98 @@
# Calculate and sum the single alignment chi-squared value (for
the PCS).
for align_index in range(self.num_align):
+ chi2_sum = chi2_sum + chi2(self.pcs[align_index],
self.pcs_theta[align_index], self.pcs_error[align_index])
+
+ # Return the chi-squared value.
+ return chi2_sum
+
+
+ def func_double_rotor_quad_int(self, params):
+ """SciPy quadratic integration target function for the double rotor
model.
+
+ This function optimises the model parameters using the RDC and PCS
base data. Quasi-random, Sobol' sequence based, numerical integration is
used for the PCS.
+
+
+ @param params: The vector of parameter values. These can include
{pivot_x, pivot_y, pivot_z, pivot_disp, ave_pos_x, ave_pos_y, ave_pos_z,
ave_pos_alpha, ave_pos_beta, ave_pos_gamma, eigen_alpha, eigen_beta,
eigen_gamma, cone_sigma_max, cone_sigma_max_2}.
+ @type params: list of float
+ @return: The chi-squared or SSE value.
+ @rtype: float
+ """
+
+ # Scaling.
+ if self.scaling_flag:
+ params = dot(params, self.scaling_matrix)
+
+ # Unpack the parameters.
+ if self.pivot_opt:
+ pivot2 = outer(self.spin_ones_struct, params[:3])
+ param_disp = params[3]
+ self._translation_vector = params[4:7]
+ ave_pos_alpha, ave_pos_beta, ave_pos_gamma, eigen_alpha,
eigen_beta, eigen_gamma, sigma_max, sigma_max_2 = params[7:]
+ else:
+ pivot2 = self.pivot
+ param_disp = params[0]
+ self._translation_vector = params[1:4]
+ ave_pos_alpha, ave_pos_beta, ave_pos_gamma, eigen_alpha,
eigen_beta, eigen_gamma, sigma_max, sigma_max_2 = params[4:]
+
+ # Reconstruct the full eigenframe of the motion.
+ euler_to_R_zyz(eigen_alpha, eigen_beta, eigen_gamma, self.R_eigen)
+
+ # The Kronecker product of the eigenframe.
+ Rx2_eigen = kron_prod(self.R_eigen, self.R_eigen)
+
+ # Generate the 2nd degree Frame Order super matrix.
+ frame_order_2nd =
compile_2nd_matrix_double_rotor(self.frame_order_2nd, Rx2_eigen, sigma_max,
sigma_max_2)
+
+ # The average frame rotation matrix (and reduce and rotate the
tensors).
+ self.reduce_and_rot(ave_pos_alpha, ave_pos_beta, ave_pos_gamma,
frame_order_2nd)
+
+ # Pre-transpose matrices for faster calculations.
+ RT_eigen = transpose(self.R_eigen)
+ RT_ave = transpose(self.R_ave)
+
+ # Pre-calculate all the necessary vectors.
+ if self.pcs_flag:
+ # The 1st pivot point (sum of the 2nd pivot and the displacement
along the eigenframe z-axis).
+ pivot = pivot2 + param_disp * self.R_eigen[:, 2]
+
+ # Calculate the vectors.
+ self.calc_vectors(pivot=pivot, pivot2=pivot2, R_ave=self.R_ave,
RT_ave=RT_ave)
+
+ # Initial chi-squared (or SSE) value.
+ chi2_sum = 0.0
+
+ # RDCs.
+ if self.rdc_flag:
+ # Loop over each alignment.
+ for align_index in range(self.num_align):
+ # Loop over the RDCs.
+ for j in range(self.num_interatom):
+ # The back calculated RDC.
+ if not self.missing_rdc[align_index, j]:
+ self.rdc_theta[align_index, j] =
rdc_tensor(self.dip_const[j], self.rdc_vect[j], self.A_3D_bc[align_index])
+
+ # Calculate and sum the single alignment chi-squared value
(for the RDC).
+ chi2_sum = chi2_sum + chi2(self.rdc[align_index],
self.rdc_theta[align_index], self.rdc_error[align_index])
+
+ # PCS via numerical integration.
+ if self.pcs_flag:
+ # Loop over each alignment.
+ for align_index in range(self.num_align):
+ # Loop over the PCSs.
+ for j in range(self.num_spins):
+ # The back calculated PCS.
+ if not self.missing_pcs[align_index, j]:
+ # Forwards and reverse rotations.
+ if self.full_in_ref_frame[align_index]:
+ r_pivot_atom = self.r_pivot_atom[j]
+ else:
+ r_pivot_atom = self.r_pivot_atom_rev[j]
+
+ # The numerical integration.
+ self.pcs_theta[align_index, j] =
pcs_numeric_quad_int_double_rotor(sigma_max=sigma_max,
sigma_max_2=sigma_max_2, c=self.pcs_const[align_index, j],
r_pivot_atom=r_pivot_atom, r_ln_pivot=self.r_ln_pivot[0],
r_inter_pivot=self.r_inter_pivot[j], A=self.A_3D[align_index],
R_eigen=self.R_eigen, RT_eigen=RT_eigen, Ri_prime=self.Ri_prime,
Ri2_prime=self.Ri2_prime)
+
+ # Calculate and sum the single alignment chi-squared value
(for the PCS).
chi2_sum = chi2_sum + chi2(self.pcs[align_index],
self.pcs_theta[align_index], self.pcs_error[align_index])
# Return the chi-squared value.
r26015 - in /branches/frame_order_cleanup: lib/frame_order/double_rotor.py target_functions/frame_order.py