Author: tlinnet
Date: Sun Jun 15 15:15:11 2014
New Revision: 23957
URL: http://svn.gna.org/viewcvs/relax?rev=23957&view=rev
Log:
Moved the bulk operation of model CPMG 2site 3d into the lib file.
This is to keep the API clean.
Task #7807 (https://gna.org/task/index.php?7807): Speed-up of dispersion
models for Clustered analysis.
Modified:
branches/disp_spin_speed/lib/dispersion/ns_cpmg_2site_3d.py
Modified: branches/disp_spin_speed/lib/dispersion/ns_cpmg_2site_3d.py
URL:
http://svn.gna.org/viewcvs/relax/branches/disp_spin_speed/lib/dispersion/ns_cpmg_2site_3d.py?rev=23957&r1=23956&r2=23957&view=diff
==============================================================================
--- branches/disp_spin_speed/lib/dispersion/ns_cpmg_2site_3d.py (original)
+++ branches/disp_spin_speed/lib/dispersion/ns_cpmg_2site_3d.py Sun Jun 15
15:15:11 2014
@@ -54,9 +54,9 @@
"""
# Python module imports.
-from numpy import dot, fabs, isfinite, log, min, ones, ndarray
-from numpy.ma import fix_invalid, masked_less_equal, masked_where
-import numpy as np
+from numpy import dot, fabs, isfinite, log, min, sum
+from numpy.ma import fix_invalid, masked_where
+
# relax module imports.
from lib.dispersion.ns_matrices import rcpmg_3d
@@ -64,7 +64,7 @@
from lib.linear_algebra.matrix_exponential import matrix_exponential
-def r2eff_ns_cpmg_2site_3D(r180x=None, M0=None, r10a=0.0, r10b=0.0,
r20a=None, r20b=None, pA=None, dw=None, kex=None, inv_tcpmg=None, tcp=None,
back_calc=None, num_points=None, power=None):
+def r2eff_ns_cpmg_2site_3D(r180x=None, M0=None, r10a=0.0, r10b=0.0,
r20a=None, r20b=None, pA=None, dw=None, dw_orig=None, kex=None,
inv_tcpmg=None, tcp=None, back_calc=None, num_points=None, power=None):
"""The 2-site numerical solution to the Bloch-McConnell equation.
This function calculates and stores the R2eff values.
@@ -79,43 +79,41 @@
@keyword r10b: The R1 value for state B.
@type r10b: float
@keyword r20a: The R2 value for state A in the absence of
exchange.
- @type r20a: float
+ @type r20a: numpy float array of rank [NE][NS][[NM][NO][ND]
@keyword r20b: The R2 value for state B in the absence of
exchange.
- @type r20b: float
+ @type r20b: numpy float array of rank [NE][NS][[NM][NO][ND]
@keyword pA: The population of state A.
@type pA: float
@keyword dw: The chemical exchange difference between
states A and B in rad/s.
- @type dw: float
+ @type dw: numpy float array of rank [NE][NS][[NM][NO][ND]
+ @keyword dw_orig: The chemical exchange difference between
states A and B in ppm. This is only for faster checking of zero value,
which result in no exchange.
+ @type dw_orig: numpy float array of rank-1
@keyword kex: The kex parameter value (the exchange rate in
rad/s).
@type kex: float
@keyword inv_tcpmg: The inverse of the total duration of the CPMG
element (in inverse seconds).
- @type inv_tcpmg: float
+ @type inv_tcpmg: numpy float array of rank [NE][NS][[NM][NO][ND]
@keyword tcp: The tau_CPMG times (1 / 4.nu1).
- @type tcp: numpy rank-1 float array
+ @type tcp: numpy float array of rank [NE][NS][[NM][NO][ND]
@keyword back_calc: The array for holding the back calculated
R2eff values. Each element corresponds to one of the CPMG nu1 frequencies.
- @type back_calc: numpy rank-1 float array
+ @type back_calc: numpy float array of rank [NE][NS][[NM][NO][ND]
@keyword num_points: The number of points on the dispersion curve,
equal to the length of the tcp and back_calc arguments.
- @type num_points: int
+ @type num_points: numpy int array of rank [NE][NS][[NM][NO][ND]
@keyword power: The matrix exponential power array.
- @type power: numpy int16, rank-1 array
+ @type power: numpy int array of rank [NE][NS][[NM][NO][ND]
"""
-
- # This is temporary hack between rank 1 and multi rank.
- dw_a = ones([num_points]) * dw
- r20a_a = ones([num_points]) * r20a
# Flag to tell if values should be replaced if math function is
violated.
t_dw_zero = False
# Catch parameter values that will result in no exchange, returning
flat R2eff = R20 lines (when kex = 0.0, k_AB = 0.0).
if pA == 1.0 or kex == 0.0:
- back_calc[:] = r20a_a
+ back_calc[:] = r20a
return
# Test if dw is zero. Wait for replacement, since this is spin
specific.
- if min(fabs(dw_a)) == 0.0:
+ if min(fabs(dw_orig)) == 0.0:
t_dw_zero = True
- mask_dw_zero = masked_where(dw_a == 0.0, dw_a)
+ mask_dw_zero = masked_where(dw == 0.0, dw)
# Once off parameter conversions.
pB = 1.0 - pA
@@ -126,40 +124,58 @@
M0[1] = pA
M0[4] = pB
- # The matrix R that contains all the contributions to the evolution,
i.e. relaxation, exchange and chemical shift evolution.
- R = rcpmg_3d(R1A=r10a, R1B=r10b, R2A=r20a, R2B=r20b, pA=pA, pB=pB,
dw=dw, k_AB=k_AB, k_BA=k_BA)
+ # Extract the total numbers of experiments, number of spins, number of
magnetic field strength, number of offsets, maximum number of dispersion
point.
+ NE, NS, NM, NO, ND = back_calc.shape
- # Loop over the time points, back calculating the R2eff values.
- for i in range(num_points):
- # Initial magnetisation.
- Mint = M0.reshape(7, 1)
+ # Loop over the spins
+ for si in range(NS):
+ # Loop over the spectrometer frequencies.
+ for mi in range(NM):
- # This matrix is a propagator that will evolve the magnetization
with the matrix R for a delay tcp.
- Rexpo = matrix_exponential(R*tcp[i])
+ # Extract the values from the higher dimensional arrays.
+ R2A_si_mi=r20a[0][si][mi][0][0]
+ R2B_si_mi=r20b[0][si][mi][0][0]
+ dw_si_mi = dw[0][si][mi][0][0]
+ num_points_si_mi = int(num_points[0][si][mi][0][0])
- # Temp matrix.
- t_mat =
Rexpo.dot(r180x).dot(Rexpo).dot(Rexpo).dot(r180x).dot(Rexpo).dot(Rexpo).dot(r180x).dot(Rexpo).dot(Rexpo).dot(r180x).dot(Rexpo)
+ # The matrix R that contains all the contributions to the
evolution, i.e. relaxation, exchange and chemical shift evolution.
+ R = rcpmg_3d(R1A=r10a, R1B=r10b, R2A=R2A_si_mi, R2B=R2B_si_mi,
pA=pA, pB=pB, dw=dw_si_mi, k_AB=k_AB, k_BA=k_BA)
- # Loop over the CPMG elements, propagating the magnetisation.
- for j in range(power[i]/2):
- Mint = t_mat.dot(Mint)
+ # Loop over the time points, back calculating the R2eff values.
+ for di in range(num_points_si_mi):
+ # Extract the values from the higher dimensional arrays.
+ tcp_si_mi_di = tcp[0][si][mi][0][di]
+ inv_tcpmg_si_mi_di = inv_tcpmg[0][si][mi][0][di]
+ power_si_mi_di = int(power[0][si][mi][0][di])
+ r20a_si_mi_di = r20a[0][si][mi][0][di]
- # The next lines calculate the R2eff using a two-point
approximation, i.e. assuming that the decay is mono-exponential.
- Mx = Mint[1][0] / pA
- if Mx <= 0.0 or isNaN(Mx):
- back_calc[i] = r20a
- else:
- back_calc[i]= -inv_tcpmg * log(Mx)
+ # Initial magnetisation.
+ Mint = M0
+
+ # This matrix is a propagator that will evolve the
magnetization with the matrix R for a delay tcp.
+ Rexpo = matrix_exponential(R*tcp_si_mi_di)
+
+ # Temp matrix.
+ t_mat =
Rexpo.dot(r180x).dot(Rexpo).dot(Rexpo).dot(r180x).dot(Rexpo).dot(Rexpo).dot(r180x).dot(Rexpo).dot(Rexpo).dot(r180x).dot(Rexpo)
+
+ # Loop over the CPMG elements, propagating the
magnetisation.
+ for j in range(power_si_mi_di/2):
+ Mint = t_mat.dot(Mint)
+
+ # The next lines calculate the R2eff using a two-point
approximation, i.e. assuming that the decay is mono-exponential.
+ Mx = Mint[1] / pA
+ #print back_calc[0][si][mi][0]
+ #print lkhj
+
+ if Mx <= 0.0 or isNaN(Mx):
+ back_calc[0][si][mi][0][di] = r20a_si_mi_di
+ else:
+ back_calc[0][si][mi][0][di] = - inv_tcpmg_si_mi_di *
log(Mx)
# Replace data in array.
# If dw is zero.
if t_dw_zero:
- back_calc[mask_dw_zero.mask] = r20a_a[mask_dw_zero.mask]
-
- # If Mx is less than 0.0
- if min(Mx) <= 0.0:
- mask_min_mx = masked_less_equal(Mx, 0.0)
- back_calc[mask_min_mx.mask] = r20a_a[mask_min_mx.mask]
+ back_calc[mask_dw_zero.mask] = r20a[mask_dw_zero.mask]
# Catch errors, taking a sum over array is the fastest way to check for
# +/- inf (infinity) and nan (not a number).
_______________________________________________
relax (http://www.nmr-relax.com)
This is the relax-commits mailing list
relax-commits@xxxxxxx
To unsubscribe from this list, get a password
reminder, or change your subscription options,
visit the list information page at
https://mail.gna.org/listinfo/relax-commits
Re: r23957 - /branches/disp_spin_speed/lib/dispersion/ns_cpmg_2site_3d.py