Author: tlinnet
Date: Mon Jun 16 22:11:38 2014
New Revision: 24003
URL: http://svn.gna.org/viewcvs/relax?rev=24003&view=rev
Log:
Modified lib function for NS MMQ 2site, to have looping over spins and
frequencies inside lib function.
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_mmq_2site.py
Modified: branches/disp_spin_speed/lib/dispersion/ns_mmq_2site.py
URL:
http://svn.gna.org/viewcvs/relax/branches/disp_spin_speed/lib/dispersion/ns_mmq_2site.py?rev=24003&r1=24002&r2=24003&view=diff
==============================================================================
--- branches/disp_spin_speed/lib/dispersion/ns_mmq_2site.py (original)
+++ branches/disp_spin_speed/lib/dispersion/ns_mmq_2site.py Mon Jun 16
22:11:38 2014
@@ -132,91 +132,107 @@
@type power: numpy int16, rank-1 array
"""
- # Populate the m1 and m2 matrices (only once per function call for
speed).
- populate_matrix(matrix=m1, R20A=R20A, R20B=R20B, dw=-dw-dwH, k_AB=k_AB,
k_BA=k_BA) # D+ matrix component.
- populate_matrix(matrix=m2, R20A=R20A, R20B=R20B, dw=dw-dwH, k_AB=k_AB,
k_BA=k_BA) # Z- matrix component.
-
- # Loop over the time points, back calculating the R2eff values.
- for i in range(num_points):
- # The M1 and M2 matrices.
- M1 = matrix_exponential(m1*tcp[i]) # Equivalent to D+.
- M2 = matrix_exponential(m2*tcp[i]) # Equivalent to Z-.
-
- # The complex conjugates M1* and M2*
- M1_star = conj(M1) # Equivalent to D+*.
- M2_star = conj(M2) # Equivalent to Z-*.
-
- # Repetitive dot products (minimised for speed).
- M1_M2 = dot(M1, M2)
- M2_M1 = dot(M2, M1)
- M1_M2_M2_M1 = dot(M1_M2, M2_M1)
- M2_M1_M1_M2 = dot(M2_M1, M1_M2)
- M1_M2_star = dot(M1_star, M2_star)
- M2_M1_star = dot(M2_star, M1_star)
- M1_M2_M2_M1_star = dot(M1_M2_star, M2_M1_star)
- M2_M1_M1_M2_star = dot(M2_M1_star, M1_M2_star)
-
- # Special case of 1 CPMG block - the power is zero.
- if power[i] == 1:
- # M1.M2.
- A = M1_M2
-
- # M1*.M2*.
- B = M1_M2_star
-
- # M2.M1.
- C = M2_M1
-
- # M2*.M1*.
- D = M2_M1_star
-
- # Matrices for even number of CPMG blocks.
- elif power[i] % 2 == 0:
- # The power factor (only calculate once).
- fact = int(floor(power[i] / 2))
-
- # (M1.M2.M2.M1)^(n/2).
- A = square_matrix_power(M1_M2_M2_M1, fact)
-
- # (M2*.M1*.M1*.M2*)^(n/2).
- B = square_matrix_power(M2_M1_M1_M2_star, fact)
-
- # (M2.M1.M1.M2)^(n/2).
- C = square_matrix_power(M2_M1_M1_M2, fact)
-
- # (M1*.M2*.M2*.M1*)^(n/2).
- D = square_matrix_power(M1_M2_M2_M1_star, fact)
-
- # Matrices for odd number of CPMG blocks.
- else:
- # The power factor (only calculate once).
- fact = int(floor((power[i] - 1) / 2))
-
- # (M1.M2.M2.M1)^((n-1)/2).M1.M2.
- A = square_matrix_power(M1_M2_M2_M1, fact)
- A = dot(A, M1_M2)
-
- # (M1*.M2*.M2*.M1*)^((n-1)/2).M1*.M2*.
- B = square_matrix_power(M1_M2_M2_M1_star, fact)
- B = dot(B, M1_M2_star)
-
- # (M2.M1.M1.M2)^((n-1)/2).M2.M1.
- C = square_matrix_power(M2_M1_M1_M2, fact)
- C = dot(C, M2_M1)
-
- # (M2*.M1*.M1*.M2*)^((n-1)/2).M2*.M1*.
- D = square_matrix_power(M2_M1_M1_M2_star, fact)
- D = dot(D, M2_M1_star)
-
- # The next lines calculate the R2eff using a two-point
approximation, i.e. assuming that the decay is mono-exponential.
- A_B = dot(A, B)
- C_D = dot(C, D)
- Mx = dot(dot(F_vector, (A_B + C_D)), M0)
- Mx = Mx.real / 2.0
- if Mx <= 0.0 or isNaN(Mx):
- back_calc[i] = 1e99
- else:
- back_calc[i]= -inv_tcpmg[i] * log(Mx / pA)
+ # Extract shape of experiment.
+ NS, NM, NO = num_points.shape
+
+ # Loop over spins.
+ for si in range(NS):
+ # Loop over the spectrometer frequencies.
+ for mi in range(NM):
+ # Loop over offsets:
+ for oi in range(NO):
+
+ r20a_si_mi_oi = R20A[si][mi][oi][0]
+ r20b_si_mi_oi = R20B[si][mi][oi][0]
+ dw_si_mi_oi = dw[si][mi][oi][0]
+ dwH_si_mi_oi = dwH[si][mi][oi][0]
+ num_points_si_mi_oi = num_points[si][mi][oi]
+
+ # Populate the m1 and m2 matrices (only once per function
call for speed).
+ populate_matrix(matrix=m1, R20A=r20a_si_mi_oi,
R20B=r20b_si_mi_oi, dw=-dw_si_mi_oi - dwH_si_mi_oi, k_AB=k_AB, k_BA=k_BA)
# D+ matrix component.
+ populate_matrix(matrix=m2, R20A=r20a_si_mi_oi,
R20B=r20b_si_mi_oi, dw=dw_si_mi_oi - dwH_si_mi_oi, k_AB=k_AB, k_BA=k_BA) #
Z- matrix component.
+
+ # Loop over the time points, back calculating the R2eff
values.
+ for i in range(num_points_si_mi_oi):
+ # The M1 and M2 matrices.
+ M1 = matrix_exponential(m1*tcp[si][mi][oi][i]) #
Equivalent to D+.
+ M2 = matrix_exponential(m2*tcp[si][mi][oi][i]) #
Equivalent to Z-.
+
+ # The complex conjugates M1* and M2*
+ M1_star = conj(M1) # Equivalent to D+*.
+ M2_star = conj(M2) # Equivalent to Z-*.
+
+ # Repetitive dot products (minimised for speed).
+ M1_M2 = dot(M1, M2)
+ M2_M1 = dot(M2, M1)
+ M1_M2_M2_M1 = dot(M1_M2, M2_M1)
+ M2_M1_M1_M2 = dot(M2_M1, M1_M2)
+ M1_M2_star = dot(M1_star, M2_star)
+ M2_M1_star = dot(M2_star, M1_star)
+ M1_M2_M2_M1_star = dot(M1_M2_star, M2_M1_star)
+ M2_M1_M1_M2_star = dot(M2_M1_star, M1_M2_star)
+
+ # Special case of 1 CPMG block - the power is zero.
+ if power[si][mi][oi][i] == 1:
+ # M1.M2.
+ A = M1_M2
+
+ # M1*.M2*.
+ B = M1_M2_star
+
+ # M2.M1.
+ C = M2_M1
+
+ # M2*.M1*.
+ D = M2_M1_star
+
+ # Matrices for even number of CPMG blocks.
+ elif power[si][mi][oi][i] % 2 == 0:
+ # The power factor (only calculate once).
+ fact = int(floor(power[si][mi][oi][i] / 2))
+
+ # (M1.M2.M2.M1)^(n/2).
+ A = square_matrix_power(M1_M2_M2_M1, fact)
+
+ # (M2*.M1*.M1*.M2*)^(n/2).
+ B = square_matrix_power(M2_M1_M1_M2_star, fact)
+
+ # (M2.M1.M1.M2)^(n/2).
+ C = square_matrix_power(M2_M1_M1_M2, fact)
+
+ # (M1*.M2*.M2*.M1*)^(n/2).
+ D = square_matrix_power(M1_M2_M2_M1_star, fact)
+
+ # Matrices for odd number of CPMG blocks.
+ else:
+ # The power factor (only calculate once).
+ fact = int(floor((power[si][mi][oi][i] - 1) / 2))
+
+ # (M1.M2.M2.M1)^((n-1)/2).M1.M2.
+ A = square_matrix_power(M1_M2_M2_M1, fact)
+ A = dot(A, M1_M2)
+
+ # (M1*.M2*.M2*.M1*)^((n-1)/2).M1*.M2*.
+ B = square_matrix_power(M1_M2_M2_M1_star, fact)
+ B = dot(B, M1_M2_star)
+
+ # (M2.M1.M1.M2)^((n-1)/2).M2.M1.
+ C = square_matrix_power(M2_M1_M1_M2, fact)
+ C = dot(C, M2_M1)
+
+ # (M2*.M1*.M1*.M2*)^((n-1)/2).M2*.M1*.
+ D = square_matrix_power(M2_M1_M1_M2_star, fact)
+ D = dot(D, M2_M1_star)
+
+ # The next lines calculate the R2eff using a two-point
approximation, i.e. assuming that the decay is mono-exponential.
+ A_B = dot(A, B)
+ C_D = dot(C, D)
+ Mx = dot(dot(F_vector, (A_B + C_D)), M0)
+ Mx = Mx.real / 2.0
+ if Mx <= 0.0 or isNaN(Mx):
+ back_calc[si][mi][oi][i] = 1e99
+ else:
+ back_calc[si][mi][oi][i]= -inv_tcpmg[si][mi][oi][i]
* log(Mx / pA)
def r2eff_ns_mmq_2site_sq_dq_zq(M0=None, F_vector=array([1, 0], float64),
m1=None, m2=None, R20A=None, R20B=None, pA=None, pB=None, dw=None, dwH=None,
k_AB=None, k_BA=None, inv_tcpmg=None, tcp=None, back_calc=None,
num_points=None, power=None):
r24003 - /branches/disp_spin_speed/lib/dispersion/ns_mmq_2site.py