Author: bugman
Date: Fri Feb 3 11:34:37 2012
New Revision: 15294
URL: http://svn.gna.org/viewcvs/relax?rev=15294&view=rev
Log:
Converted all of the PCS numerical integration methods to use the
quasi-random Sobol' sequence data.
This is to convert from Monte Carlo to quasi-random numerical integration.
Modified:
branches/frame_order_testing/maths_fns/frame_order.py
branches/frame_order_testing/maths_fns/frame_order_matrix_ops.py
Modified: branches/frame_order_testing/maths_fns/frame_order.py
URL:
http://svn.gna.org/viewcvs/relax/branches/frame_order_testing/maths_fns/frame_order.py?rev=15294&r1=15293&r2=15294&view=diff
==============================================================================
--- branches/frame_order_testing/maths_fns/frame_order.py (original)
+++ branches/frame_order_testing/maths_fns/frame_order.py Fri Feb 3 11:34:37
2012
@@ -36,7 +36,7 @@
from maths_fns.alignment_tensor import to_5D, to_tensor
from maths_fns.chi2 import chi2
from maths_fns.coord_transform import spherical_to_cartesian
-from maths_fns.frame_order_matrix_ops import compile_2nd_matrix_free_rotor,
compile_2nd_matrix_iso_cone, compile_2nd_matrix_iso_cone_free_rotor,
compile_2nd_matrix_iso_cone_torsionless, compile_2nd_matrix_pseudo_ellipse,
compile_2nd_matrix_pseudo_ellipse_free_rotor,
compile_2nd_matrix_pseudo_ellipse_torsionless, compile_2nd_matrix_rotor,
reduce_alignment_tensor, pcs_numeric_int_iso_cone,
pcs_numeric_int_iso_cone_mcint, pcs_numeric_int_iso_cone_torsionless,
pcs_numeric_int_iso_cone_torsionless_mcint, pcs_numeric_int_pseudo_ellipse,
pcs_numeric_int_pseudo_ellipse_qrint,
pcs_numeric_int_pseudo_ellipse_torsionless,
pcs_numeric_int_pseudo_ellipse_torsionless_mcint, pcs_numeric_int_rotor,
pcs_numeric_int_rotor_qrint
+from maths_fns.frame_order_matrix_ops import compile_2nd_matrix_free_rotor,
compile_2nd_matrix_iso_cone, compile_2nd_matrix_iso_cone_free_rotor,
compile_2nd_matrix_iso_cone_torsionless, compile_2nd_matrix_pseudo_ellipse,
compile_2nd_matrix_pseudo_ellipse_free_rotor,
compile_2nd_matrix_pseudo_ellipse_torsionless, compile_2nd_matrix_rotor,
reduce_alignment_tensor, pcs_numeric_int_iso_cone,
pcs_numeric_int_iso_cone_qrint, pcs_numeric_int_iso_cone_torsionless,
pcs_numeric_int_iso_cone_torsionless_qrint, pcs_numeric_int_pseudo_ellipse,
pcs_numeric_int_pseudo_ellipse_qrint,
pcs_numeric_int_pseudo_ellipse_torsionless,
pcs_numeric_int_pseudo_ellipse_torsionless_qrint, pcs_numeric_int_rotor,
pcs_numeric_int_rotor_qrint
from maths_fns.kronecker_product import kron_prod
from maths_fns import order_parameters
from maths_fns.rotation_matrix import euler_to_R_zyz
Modified: branches/frame_order_testing/maths_fns/frame_order_matrix_ops.py
URL:
http://svn.gna.org/viewcvs/relax/branches/frame_order_testing/maths_fns/frame_order_matrix_ops.py?rev=15294&r1=15293&r2=15294&view=diff
==============================================================================
--- branches/frame_order_testing/maths_fns/frame_order_matrix_ops.py
(original)
+++ branches/frame_order_testing/maths_fns/frame_order_matrix_ops.py Fri Feb
3 11:34:37 2012
@@ -1256,11 +1256,11 @@
return c * result[0] / SA
-def pcs_numeric_int_iso_cone_mcint(N=1000, theta_max=None, sigma_max=None,
c=None, full_in_ref_frame=None, r_pivot_atom=None, r_pivot_atom_rev=None,
r_ln_pivot=None, A=None, R_eigen=None, RT_eigen=None, Ri_prime=None,
pcs_theta=None, pcs_theta_err=None, missing_pcs=None, error_flag=False):
+def pcs_numeric_int_iso_cone_qrint(points=None, theta_max=None,
sigma_max=None, c=None, full_in_ref_frame=None, r_pivot_atom=None,
r_pivot_atom_rev=None, r_ln_pivot=None, A=None, R_eigen=None, RT_eigen=None,
Ri_prime=None, pcs_theta=None, pcs_theta_err=None, missing_pcs=None,
error_flag=False):
"""Determine the averaged PCS value via numerical integration.
- @keyword N: The number of Monte Carlo samples to use.
- @type N: int
+ @keyword points: The Sobol points in the torsion-tilt angle
space.
+ @type points: numpy rank-2, 3D array
@keyword theta_max: The half cone angle.
@type theta_max: float
@keyword sigma_max: The maximum torsion angle.
@@ -1300,27 +1300,32 @@
pcs_theta_err[i, j] = 0.0
# Loop over the samples.
- for i in range(N):
- # Sigma and phi randomisation.
- sigma_i = uniform(-sigma_max, sigma_max)
- phi_i = uniform(-pi, pi)
-
- # Theta randomisation.
- v = uniform(cos(theta_max), 1.0)
- theta_i = acos(v)
+ num = 0
+ for i in range(len(points)):
+ # Unpack the point.
+ phi, theta, sigma = points[i]
+
+ # Outside of the distribution, so skip the point.
+ if theta > theta_max:
+ continue
+ if sigma > sigma_max or sigma < -sigma_max:
+ continue
# Calculate the PCSs for this state.
- pcs_pivot_motion_full_mcint(theta_i=theta_i, phi_i=phi_i,
sigma_i=sigma_i, full_in_ref_frame=full_in_ref_frame,
r_pivot_atom=r_pivot_atom, r_pivot_atom_rev=r_pivot_atom_rev,
r_ln_pivot=r_ln_pivot, A=A, R_eigen=R_eigen, RT_eigen=RT_eigen,
Ri_prime=Ri_prime, pcs_theta=pcs_theta, pcs_theta_err=pcs_theta_err,
missing_pcs=missing_pcs)
+ pcs_pivot_motion_full_qrint(theta_i=theta, phi_i=phi, sigma_i=sigma,
full_in_ref_frame=full_in_ref_frame, r_pivot_atom=r_pivot_atom,
r_pivot_atom_rev=r_pivot_atom_rev, r_ln_pivot=r_ln_pivot, A=A,
R_eigen=R_eigen, RT_eigen=RT_eigen, Ri_prime=Ri_prime, pcs_theta=pcs_theta,
pcs_theta_err=pcs_theta_err, missing_pcs=missing_pcs)
+
+ # Increment the number of points.
+ num += 1
# Calculate the PCS and error.
for i in range(len(pcs_theta)):
for j in range(len(pcs_theta[i])):
# The average PCS.
- pcs_theta[i, j] = c[i] * pcs_theta[i, j] / float(N)
+ pcs_theta[i, j] = c[i] * pcs_theta[i, j] / float(num)
# The error.
if error_flag:
- pcs_theta_err[i, j] = abs(pcs_theta_err[i, j] / float(N) -
pcs_theta[i, j]**2) / float(N)
+ pcs_theta_err[i, j] = abs(pcs_theta_err[i, j] / float(num)
- pcs_theta[i, j]**2) / float(num)
pcs_theta_err[i, j] = c[i] * sqrt(pcs_theta_err[i, j])
print "%8.3f +/- %-8.3f" % (pcs_theta[i, j]*1e6,
pcs_theta_err[i, j]*1e6)
@@ -1358,11 +1363,11 @@
return c * result[0] / SA
-def pcs_numeric_int_iso_cone_torsionless_mcint(N=1000, theta_max=None,
c=None, full_in_ref_frame=None, r_pivot_atom=None, r_pivot_atom_rev=None,
r_ln_pivot=None, A=None, R_eigen=None, RT_eigen=None, Ri_prime=None,
pcs_theta=None, pcs_theta_err=None, missing_pcs=None, error_flag=False):
+def pcs_numeric_int_iso_cone_torsionless_qrint(points=None, theta_max=None,
c=None, full_in_ref_frame=None, r_pivot_atom=None, r_pivot_atom_rev=None,
r_ln_pivot=None, A=None, R_eigen=None, RT_eigen=None, Ri_prime=None,
pcs_theta=None, pcs_theta_err=None, missing_pcs=None, error_flag=False):
"""Determine the averaged PCS value via numerical integration.
- @keyword N: The number of Monte Carlo samples to use.
- @type N: int
+ @keyword points: The Sobol points in the torsion-tilt angle
space.
+ @type points: numpy rank-2, 3D array
@keyword theta_max: The half cone angle.
@type theta_max: float
@keyword c: The PCS constant (without the interatomic
distance and in Angstrom units).
@@ -1400,26 +1405,30 @@
pcs_theta_err[i, j] = 0.0
# Loop over the samples.
- for i in range(N):
- # Phi randomisation.
- phi_i = uniform(-pi, pi)
-
- # Theta randomisation.
- v = uniform(cos(theta_max), 1.0)
- theta_i = acos(v)
+ num = 0
+ for i in range(len(points)):
+ # Unpack the point.
+ phi, theta = points[i]
+
+ # Outside of the distribution, so skip the point.
+ if theta > theta_max:
+ continue
# Calculate the PCSs for this state.
- pcs_pivot_motion_torsionless_mcint(theta_i=theta_i, phi_i=phi_i,
full_in_ref_frame=full_in_ref_frame, r_pivot_atom=r_pivot_atom,
r_pivot_atom_rev=r_pivot_atom_rev, r_ln_pivot=r_ln_pivot, A=A,
R_eigen=R_eigen, RT_eigen=RT_eigen, Ri_prime=Ri_prime, pcs_theta=pcs_theta,
pcs_theta_err=pcs_theta_err, missing_pcs=missing_pcs)
+ pcs_pivot_motion_torsionless_qrint(theta_i=theta, phi_i=phi,
full_in_ref_frame=full_in_ref_frame, r_pivot_atom=r_pivot_atom,
r_pivot_atom_rev=r_pivot_atom_rev, r_ln_pivot=r_ln_pivot, A=A,
R_eigen=R_eigen, RT_eigen=RT_eigen, Ri_prime=Ri_prime, pcs_theta=pcs_theta,
pcs_theta_err=pcs_theta_err, missing_pcs=missing_pcs)
+
+ # Increment the number of points.
+ num += 1
# Calculate the PCS and error.
for i in range(len(pcs_theta)):
for j in range(len(pcs_theta[i])):
# The average PCS.
- pcs_theta[i, j] = c[i] * pcs_theta[i, j] / float(N)
+ pcs_theta[i, j] = c[i] * pcs_theta[i, j] / float(num)
# The error.
if error_flag:
- pcs_theta_err[i, j] = abs(pcs_theta_err[i, j] / float(N) -
pcs_theta[i, j]**2) / float(N)
+ pcs_theta_err[i, j] = abs(pcs_theta_err[i, j] / float(num)
- pcs_theta[i, j]**2) / float(num)
pcs_theta_err[i, j] = c[i] * sqrt(pcs_theta_err[i, j])
print "%8.3f +/- %-8.3f" % (pcs_theta[i, j]*1e6,
pcs_theta_err[i, j]*1e6)
@@ -1527,7 +1536,7 @@
continue
# Calculate the PCSs for this state.
- pcs_pivot_motion_full_mcint(theta_i=theta, phi_i=phi, sigma_i=sigma,
full_in_ref_frame=full_in_ref_frame, r_pivot_atom=r_pivot_atom,
r_pivot_atom_rev=r_pivot_atom_rev, r_ln_pivot=r_ln_pivot, A=A,
R_eigen=R_eigen, RT_eigen=RT_eigen, Ri_prime=Ri_prime, pcs_theta=pcs_theta,
pcs_theta_err=pcs_theta_err, missing_pcs=missing_pcs)
+ pcs_pivot_motion_full_qrint(theta_i=theta, phi_i=phi, sigma_i=sigma,
full_in_ref_frame=full_in_ref_frame, r_pivot_atom=r_pivot_atom,
r_pivot_atom_rev=r_pivot_atom_rev, r_ln_pivot=r_ln_pivot, A=A,
R_eigen=R_eigen, RT_eigen=RT_eigen, Ri_prime=Ri_prime, pcs_theta=pcs_theta,
pcs_theta_err=pcs_theta_err, missing_pcs=missing_pcs)
# Increment the number of points.
num += 1
@@ -1585,11 +1594,11 @@
return c * result[0] / SA
-def pcs_numeric_int_pseudo_ellipse_torsionless_mcint(N=1000, theta_x=None,
theta_y=None, c=None, full_in_ref_frame=None, r_pivot_atom=None,
r_pivot_atom_rev=None, r_ln_pivot=None, A=None, R_eigen=None, RT_eigen=None,
Ri_prime=None, pcs_theta=None, pcs_theta_err=None, missing_pcs=None,
error_flag=False):
+def pcs_numeric_int_pseudo_ellipse_torsionless_qrint(points=None,
theta_x=None, theta_y=None, c=None, full_in_ref_frame=None,
r_pivot_atom=None, r_pivot_atom_rev=None, r_ln_pivot=None, A=None,
R_eigen=None, RT_eigen=None, Ri_prime=None, pcs_theta=None,
pcs_theta_err=None, missing_pcs=None, error_flag=False):
"""Determine the averaged PCS value via numerical integration.
- @keyword N: The number of Monte Carlo samples to use.
- @type N: int
+ @keyword points: The Sobol points in the torsion-tilt angle
space.
+ @type points: numpy rank-2, 3D array
@keyword theta_x: The x-axis half cone angle.
@type theta_x: float
@keyword theta_y: The y-axis half cone angle.
@@ -1629,29 +1638,33 @@
pcs_theta_err[i, j] = 0.0
# Loop over the samples.
- for i in range(N):
- # Phi randomisation.
- phi_i = uniform(-pi, pi)
+ num = 0
+ for i in range(len(points)):
+ # Unpack the point.
+ phi, theta = points[i]
# Calculate theta_max.
theta_max = tmax_pseudo_ellipse(phi_i, theta_x, theta_y)
- # Theta randomisation.
- v = uniform(cos(theta_max), 1.0)
- theta_i = acos(v)
+ # Outside of the distribution, so skip the point.
+ if theta > theta_max:
+ continue
# Calculate the PCSs for this state.
- pcs_pivot_motion_torsionless_mcint(theta_i=theta_i, phi_i=phi_i,
full_in_ref_frame=full_in_ref_frame, r_pivot_atom=r_pivot_atom,
r_pivot_atom_rev=r_pivot_atom_rev, r_ln_pivot=r_ln_pivot, A=A,
R_eigen=R_eigen, RT_eigen=RT_eigen, Ri_prime=Ri_prime, pcs_theta=pcs_theta,
pcs_theta_err=pcs_theta_err, missing_pcs=missing_pcs)
+ pcs_pivot_motion_torsionless_qrint(theta_i=theta, phi_i=phi,
full_in_ref_frame=full_in_ref_frame, r_pivot_atom=r_pivot_atom,
r_pivot_atom_rev=r_pivot_atom_rev, r_ln_pivot=r_ln_pivot, A=A,
R_eigen=R_eigen, RT_eigen=RT_eigen, Ri_prime=Ri_prime, pcs_theta=pcs_theta,
pcs_theta_err=pcs_theta_err, missing_pcs=missing_pcs)
+
+ # Increment the number of points.
+ num += 1
# Calculate the PCS and error.
for i in range(len(pcs_theta)):
for j in range(len(pcs_theta[i])):
# The average PCS.
- pcs_theta[i, j] = c[i] * pcs_theta[i, j] / float(N)
+ pcs_theta[i, j] = c[i] * pcs_theta[i, j] / float(num)
# The error.
if error_flag:
- pcs_theta_err[i, j] = abs(pcs_theta_err[i, j] / float(N) -
pcs_theta[i, j]**2) / float(N)
+ pcs_theta_err[i, j] = abs(pcs_theta_err[i, j] / float(num)
- pcs_theta[i, j]**2) / float(num)
pcs_theta_err[i, j] = c[i] * sqrt(pcs_theta_err[i, j])
print "%8.3f +/- %-8.3f" % (pcs_theta[i, j]*1e6,
pcs_theta_err[i, j]*1e6)
@@ -1829,7 +1842,7 @@
return pcs
-def pcs_pivot_motion_full_mcint(theta_i=None, phi_i=None, sigma_i=None,
full_in_ref_frame=None, r_pivot_atom=None, r_pivot_atom_rev=None,
r_ln_pivot=None, A=None, R_eigen=None, RT_eigen=None, Ri_prime=None,
pcs_theta=None, pcs_theta_err=None, missing_pcs=None, error_flag=False):
+def pcs_pivot_motion_full_qrint(theta_i=None, phi_i=None, sigma_i=None,
full_in_ref_frame=None, r_pivot_atom=None, r_pivot_atom_rev=None,
r_ln_pivot=None, A=None, R_eigen=None, RT_eigen=None, Ri_prime=None,
pcs_theta=None, pcs_theta_err=None, missing_pcs=None, error_flag=False):
"""Calculate the PCS value after a pivoted motion for the isotropic cone
model.
@keyword theta_i: The half cone opening angle (polar angle).
@@ -2104,7 +2117,7 @@
return pcs
-def pcs_pivot_motion_torsionless_mcint(theta_i=None, phi_i=None,
full_in_ref_frame=None, r_pivot_atom=None, r_pivot_atom_rev=None,
r_ln_pivot=None, A=None, R_eigen=None, RT_eigen=None, Ri_prime=None,
pcs_theta=None, pcs_theta_err=None, missing_pcs=None, error_flag=False):
+def pcs_pivot_motion_torsionless_qrint(theta_i=None, phi_i=None,
full_in_ref_frame=None, r_pivot_atom=None, r_pivot_atom_rev=None,
r_ln_pivot=None, A=None, R_eigen=None, RT_eigen=None, Ri_prime=None,
pcs_theta=None, pcs_theta_err=None, missing_pcs=None, error_flag=False):
"""Calculate the PCS value after a pivoted motion for the isotropic cone
model.
@keyword theta_i: The half cone opening angle (polar angle).
r15294 - in /branches/frame_order_testing/maths_fns: frame_order.py frame_order_matrix_ops.py