Author: bugman
Date: Tue Sep 10 14:07:02 2013
New Revision: 20961
URL: http://svn.gna.org/viewcvs/relax?rev=20961&view=rev
Log:
Speed ups for the optimisation of all of the R1rho dispersion models.
The spin-lock field strength data structure is now converted from Hz to
rad.s^-1 in the dispersion
target function initialisation. Previously the conversion was happening
multiple times per target
function call. This has a noticeable effect on the test suite timings.
Modified:
branches/relax_disp/lib/dispersion/dpl94.py
branches/relax_disp/lib/dispersion/m61.py
branches/relax_disp/lib/dispersion/m61b.py
branches/relax_disp/lib/dispersion/ns_r1rho_2site.py
branches/relax_disp/lib/dispersion/tp02.py
branches/relax_disp/target_functions/relax_disp.py
Modified: branches/relax_disp/lib/dispersion/dpl94.py
URL:
http://svn.gna.org/viewcvs/relax/branches/relax_disp/lib/dispersion/dpl94.py?rev=20961&r1=20960&r2=20961&view=diff
==============================================================================
--- branches/relax_disp/lib/dispersion/dpl94.py (original)
+++ branches/relax_disp/lib/dispersion/dpl94.py Tue Sep 10 14:07:02 2013
@@ -60,7 +60,7 @@
@type theta: numpy rank-1 float array
@keyword R1: The R1 relaxation rate.
@type R1: float
- @keyword spin_lock_fields: The CPMG nu1 frequencies.
+ @keyword spin_lock_fields: The R1rho spin-lock field strengths (in
rad.s^-1).
@type spin_lock_fields: numpy rank-1 float array
@keyword back_calc: The array for holding the back calculated
R1rho values. Each element corresponds to one of the spin-lock fields.
@type back_calc: numpy rank-1 float array
@@ -86,7 +86,7 @@
continue
# Denominator.
- denom = kex2 + (2.0*pi*spin_lock_fields[i])**2
+ denom = kex2 + spin_lock_fields[i]**2
# Avoid divide by zero.
if denom == 0.0:
Modified: branches/relax_disp/lib/dispersion/m61.py
URL:
http://svn.gna.org/viewcvs/relax/branches/relax_disp/lib/dispersion/m61.py?rev=20961&r1=20960&r2=20961&view=diff
==============================================================================
--- branches/relax_disp/lib/dispersion/m61.py (original)
+++ branches/relax_disp/lib/dispersion/m61.py Tue Sep 10 14:07:02 2013
@@ -56,7 +56,7 @@
@type phi_ex: float
@keyword kex: The kex parameter value (the exchange rate
in rad/s).
@type kex: float
- @keyword spin_lock_fields: The CPMG nu1 frequencies.
+ @keyword spin_lock_fields: The R1rho spin-lock field strengths (in
rad.s^-1).
@type spin_lock_fields: numpy rank-1 float array
@keyword back_calc: The array for holding the back calculated
R1rho values. Each element corresponds to one of the spin-lock fields.
@type back_calc: numpy rank-1 float array
@@ -78,7 +78,7 @@
continue
# Denominator.
- denom = kex2 + (2.0*pi*spin_lock_fields[i])**2
+ denom = kex2 + spin_lock_fields[i]**2
# Avoid divide by zero.
if denom == 0.0:
Modified: branches/relax_disp/lib/dispersion/m61b.py
URL:
http://svn.gna.org/viewcvs/relax/branches/relax_disp/lib/dispersion/m61b.py?rev=20961&r1=20960&r2=20961&view=diff
==============================================================================
--- branches/relax_disp/lib/dispersion/m61b.py (original)
+++ branches/relax_disp/lib/dispersion/m61b.py Tue Sep 10 14:07:02 2013
@@ -54,7 +54,7 @@
@type dw: float
@keyword kex: The kex parameter value (the exchange rate
in rad/s).
@type kex: float
- @keyword spin_lock_fields: The spin-lock field strengths (Hz).
+ @keyword spin_lock_fields: The R1rho spin-lock field strengths (in
rad.s^-1).
@type spin_lock_fields: numpy rank-1 float array
@keyword back_calc: The array for holding the back calculated
R1rho values. Each element corresponds to one of the spin-lock fields.
@type back_calc: numpy rank-1 float array
@@ -80,7 +80,7 @@
continue
# Denominator.
- denom = kex2_pA2dw2 + (2.0*pi*spin_lock_fields[i])**2
+ denom = kex2_pA2dw2 + spin_lock_fields[i]**2
# Avoid divide by zero.
if denom == 0.0:
Modified: branches/relax_disp/lib/dispersion/ns_r1rho_2site.py
URL:
http://svn.gna.org/viewcvs/relax/branches/relax_disp/lib/dispersion/ns_r1rho_2site.py?rev=20961&r1=20960&r2=20961&view=diff
==============================================================================
--- branches/relax_disp/lib/dispersion/ns_r1rho_2site.py (original)
+++ branches/relax_disp/lib/dispersion/ns_r1rho_2site.py Tue Sep 10 14:07:02
2013
@@ -67,7 +67,7 @@
@type k_AB: float
@keyword k_BA: The rate of exchange from site B to A
(rad/s).
@type k_BA: float
- @keyword spin_lock_fields: The R1rho spin-lock field strengths in Hz.
+ @keyword spin_lock_fields: The R1rho spin-lock field strengths (in
rad.s^-1).
@type spin_lock_fields: numpy rank-1 float array
@keyword relax_time: The total relaxation time period for each
spin-lock field strength (in seconds).
@type relax_time: float
@@ -86,17 +86,16 @@
# Loop over the time points, back calculating the R2eff values.
for i in range(num_points):
Wsl = offset[i] # Larmor frequency of spin lock
[s^-1].
- w1 = spin_lock_fields[i] * 2.0 * pi # Spin-lock field strength.
dA = Wa - Wsl # Offset of spin-lock from A.
dB = Wb - Wsl # Offset of spin-lock from B.
# The matrix R that contains all the contributions to the evolution,
i.e. relaxation, exchange and chemical shift evolution.
- R = rr1rho_3d(R1=r1, Rinf=r1rho_prime, pA=pA, pB=pB, wA=dA, wB=dB,
w1=w1, k_AB=k_AB, k_BA=k_BA)
+ R = rr1rho_3d(R1=r1, Rinf=r1rho_prime, pA=pA, pB=pB, wA=dA, wB=dB,
w1=spin_lock_fields[i], k_AB=k_AB, k_BA=k_BA)
# The following lines rotate the magnetization previous to spin-lock
into the weff frame.
# A new M0 is obtained: M0 = [sin(thetaA)*pA sin(thetaB)*pB 0 0
cos(thetaA)*pA cos(thetaB)*pB].
- thetaA = atan(w1/dA)
- thetaB = atan(w1/dB)
+ thetaA = atan(spin_lock_fields[i]/dA)
+ thetaB = atan(spin_lock_fields[i]/dB)
M0[0] = sin(thetaA) * pA
M0[1] = sin(thetaB) * pB
M0[4] = cos(thetaA) * pA
Modified: branches/relax_disp/lib/dispersion/tp02.py
URL:
http://svn.gna.org/viewcvs/relax/branches/relax_disp/lib/dispersion/tp02.py?rev=20961&r1=20960&r2=20961&view=diff
==============================================================================
--- branches/relax_disp/lib/dispersion/tp02.py (original)
+++ branches/relax_disp/lib/dispersion/tp02.py Tue Sep 10 14:07:02 2013
@@ -58,7 +58,7 @@
@type kex: float
@keyword R1: The R1 relaxation rate.
@type R1: float
- @keyword spin_lock_fields: The CPMG nu1 frequencies.
+ @keyword spin_lock_fields: The R1rho spin-lock field strengths (in
rad.s^-1).
@type spin_lock_fields: numpy rank-1 float array
@keyword back_calc: The array for holding the back calculated
R1rho values. Each element corresponds to one of the spin-lock fields.
@type back_calc: numpy rank-1 float array
@@ -77,18 +77,17 @@
# Loop over the dispersion points, back calculating the R1rho values.
for i in range(num_points):
# We assume that A resonates at 0 [s^-1], without loss of generality.
- Wsl = offset[i] # Larmor frequency of spin lock
[s^-1].
- w1 = spin_lock_fields[i] * 2.0 * pi # Spin-lock field strength.
- W = pA*Wa + pB*Wb # Pop-averaged Larmor frequency
[s^-1].
- da = Wa - Wsl # Offset of spin-lock from A.
- db = Wb - Wsl # Offset of spin-lock from B.
- d = W - Wsl # Offset of spin-lock from
pop-average.
- waeff2 = w1**2 + da**2 # Effective field at A.
- wbeff2 = w1**2 + db**2 # Effective field at B.
- weff2 = w1**2 + d**2 # Effective field at pop-average.
+ Wsl = offset[i] # Larmor frequency of
spin lock [s^-1].
+ W = pA*Wa + pB*Wb # Pop-averaged Larmor
frequency [s^-1].
+ da = Wa - Wsl # Offset of spin-lock
from A.
+ db = Wb - Wsl # Offset of spin-lock
from B.
+ d = W - Wsl # Offset of spin-lock
from pop-average.
+ waeff2 = spin_lock_fields[i]**2 + da**2 # Effective field at A.
+ wbeff2 = spin_lock_fields[i]**2 + db**2 # Effective field at B.
+ weff2 = spin_lock_fields[i]**2 + d**2 # Effective field at
pop-average.
# The rotating frame flip angle.
- theta = atan(w1 / d)
+ theta = atan(spin_lock_fields[i] / d)
# Repetitive calculations (to speed up calculations).
sin_theta2 = sin(theta)**2
Modified: branches/relax_disp/target_functions/relax_disp.py
URL:
http://svn.gna.org/viewcvs/relax/branches/relax_disp/target_functions/relax_disp.py?rev=20961&r1=20960&r2=20961&view=diff
==============================================================================
--- branches/relax_disp/target_functions/relax_disp.py (original)
+++ branches/relax_disp/target_functions/relax_disp.py Tue Sep 10 14:07:02
2013
@@ -24,6 +24,7 @@
"""Target functions for relaxation dispersion."""
# Python module imports.
+from math import pi
from numpy import complex64, dot, float64, int16, zeros
# relax module imports.
@@ -209,6 +210,10 @@
self.tau_cpmg[i] = 0.25 / self.cpmg_frqs[i]
self.power[i] = int(round(self.cpmg_frqs[i] *
self.relax_time))
+ # Convert the spin-lock data to rad.s^-1.
+ if spin_lock_nu1 != None:
+ self.spin_lock_omega1 = 2.0 * pi * self.spin_lock_nu1
+
# The inverted relaxation delay.
if model in [MODEL_NS_CPMG_2SITE_3D, MODEL_NS_CPMG_2SITE_3D_FULL,
MODEL_NS_CPMG_2SITE_EXPANDED, MODEL_NS_CPMG_2SITE_STAR,
MODEL_NS_CPMG_2SITE_STAR_FULL, MODEL_NS_R1RHO_2SITE]:
self.inv_relax_time = 1.0 / relax_time
@@ -514,7 +519,7 @@
phi_ex_scaled = phi_ex[spin_index] * self.frqs[spin_index,
frq_index]**2
# Back calculate the R2eff values.
- r1rho_DPL94(r1rho_prime=R20[r20_index],
phi_ex=phi_ex_scaled, kex=kex, theta=self.tilt_angles[spin_index, frq_index],
R1=self.r1[spin_index, frq_index], spin_lock_fields=self.spin_lock_nu1,
back_calc=self.back_calc[spin_index, frq_index],
num_points=self.num_disp_points)
+ r1rho_DPL94(r1rho_prime=R20[r20_index],
phi_ex=phi_ex_scaled, kex=kex, theta=self.tilt_angles[spin_index, frq_index],
R1=self.r1[spin_index, frq_index], spin_lock_fields=self.spin_lock_omega1,
back_calc=self.back_calc[spin_index, frq_index],
num_points=self.num_disp_points)
# For all missing data points, set the back-calculated value
to the measured values so that it has no effect on the chi-squared value.
for point_index in range(self.num_disp_points):
@@ -710,7 +715,7 @@
phi_ex_scaled = phi_ex[spin_index] * self.frqs[spin_index,
frq_index]**2
# Back calculate the R2eff values.
- r1rho_M61(r1rho_prime=R20[r20_index], phi_ex=phi_ex_scaled,
kex=kex, spin_lock_fields=self.spin_lock_nu1,
back_calc=self.back_calc[spin_index, frq_index],
num_points=self.num_disp_points)
+ r1rho_M61(r1rho_prime=R20[r20_index], phi_ex=phi_ex_scaled,
kex=kex, spin_lock_fields=self.spin_lock_omega1,
back_calc=self.back_calc[spin_index, frq_index],
num_points=self.num_disp_points)
# For all missing data points, set the back-calculated value
to the measured values so that it has no effect on the chi-squared value.
for point_index in range(self.num_disp_points):
@@ -757,7 +762,7 @@
dw_frq = dw[spin_index] * self.frqs[spin_index, frq_index]
# Back calculate the R1rho values.
- r1rho_M61b(r1rho_prime=R20[r20_index], pA=pA, dw=dw_frq,
kex=kex, spin_lock_fields=self.spin_lock_nu1,
back_calc=self.back_calc[spin_index, frq_index],
num_points=self.num_disp_points)
+ r1rho_M61b(r1rho_prime=R20[r20_index], pA=pA, dw=dw_frq,
kex=kex, spin_lock_fields=self.spin_lock_omega1,
back_calc=self.back_calc[spin_index, frq_index],
num_points=self.num_disp_points)
# For all missing data points, set the back-calculated value
to the measured values so that it has no effect on the chi-squared value.
for point_index in range(self.num_disp_points):
@@ -1004,7 +1009,7 @@
dw_frq = dw[spin_index] * self.frqs[spin_index, frq_index]
# Back calculate the R2eff values.
- ns_r1rho_2site(M0=self.M0,
r1rho_prime=r1rho_prime[r20_index], omega=self.chemical_shifts[spin_index,
frq_index], offset=self.spin_lock_offsets[spin_index, frq_index],
r1=self.r1[spin_index, frq_index], pA=pA, pB=pB, dw=dw_frq, k_AB=k_AB,
k_BA=k_BA, spin_lock_fields=self.spin_lock_nu1, relax_time=self.relax_time,
inv_relax_time=self.inv_relax_time, back_calc=self.back_calc[spin_index,
frq_index], num_points=self.num_disp_points)
+ ns_r1rho_2site(M0=self.M0,
r1rho_prime=r1rho_prime[r20_index], omega=self.chemical_shifts[spin_index,
frq_index], offset=self.spin_lock_offsets[spin_index, frq_index],
r1=self.r1[spin_index, frq_index], pA=pA, pB=pB, dw=dw_frq, k_AB=k_AB,
k_BA=k_BA, spin_lock_fields=self.spin_lock_omega1,
relax_time=self.relax_time, inv_relax_time=self.inv_relax_time,
back_calc=self.back_calc[spin_index, frq_index],
num_points=self.num_disp_points)
# For all missing data points, set the back-calculated value
to the measured values so that it has no effect on the chi-squared value.
for point_index in range(self.num_disp_points):
@@ -1054,7 +1059,7 @@
dw_frq = dw[spin_index] * self.frqs[spin_index, frq_index]
# Back calculate the R1rho values.
- r1rho_TP02(r1rho_prime=R20[r20_index],
omega=self.chemical_shifts[spin_index, frq_index],
offset=self.spin_lock_offsets[spin_index, frq_index], pA=pA, pB=pB,
dw=dw_frq, kex=kex, R1=self.r1[spin_index, frq_index],
spin_lock_fields=self.spin_lock_nu1, back_calc=self.back_calc[spin_index,
frq_index], num_points=self.num_disp_points)
+ r1rho_TP02(r1rho_prime=R20[r20_index],
omega=self.chemical_shifts[spin_index, frq_index],
offset=self.spin_lock_offsets[spin_index, frq_index], pA=pA, pB=pB,
dw=dw_frq, kex=kex, R1=self.r1[spin_index, frq_index],
spin_lock_fields=self.spin_lock_omega1, back_calc=self.back_calc[spin_index,
frq_index], num_points=self.num_disp_points)
# For all missing data points, set the back-calculated value
to the measured values so that it has no effect on the chi-squared value.
for point_index in range(self.num_disp_points):
r20961 - in /branches/relax_disp: lib/dispersion/ target_functions/