Author: tlinnet
Date: Fri Jun 13 07:21:02 2014
New Revision: 23901
URL: http://svn.gna.org/viewcvs/relax?rev=23901&view=rev
Log:
Replaced the loop structure in target function of TAP03 with numpy arrays.
This makes the model faster.
Task #7807 (https://gna.org/task/index.php?7807): Speed-up of dispersion
models for Clustered analysis.
Modified:
branches/disp_spin_speed/target_functions/relax_disp.py
Modified: branches/disp_spin_speed/target_functions/relax_disp.py
URL:
http://svn.gna.org/viewcvs/relax/branches/disp_spin_speed/target_functions/relax_disp.py?rev=23901&r1=23900&r2=23901&view=diff
==============================================================================
--- branches/disp_spin_speed/target_functions/relax_disp.py (original)
+++ branches/disp_spin_speed/target_functions/relax_disp.py Fri Jun 13
07:21:02 2014
@@ -395,7 +395,7 @@
self.func = self.func_ns_mmq_3site_linear
# Setup special numpy array structures, for higher dimensional
computation.
- if model in [MODEL_B14, MODEL_B14_FULL, MODEL_CR72, MODEL_CR72_FULL,
MODEL_DPL94, MODEL_TSMFK01]:
+ if model in [MODEL_B14, MODEL_B14_FULL, MODEL_CR72, MODEL_CR72_FULL,
MODEL_DPL94, MODEL_TAP03, MODEL_TSMFK01]:
# Get the shape of back_calc structure.
# If using just one field, or having the same number of
dispersion points, the shape would extend to that number.
# Shape has to be: [ei][si][mi][oi].
@@ -435,10 +435,12 @@
self.power_a = ones(self.numpy_array_shape, int16)
# For R1rho data.
- if model in [MODEL_DPL94]:
+ if model in [MODEL_DPL94, MODEL_TAP03]:
self.tilt_angles_a = deepcopy(zeros_a)
self.spin_lock_omega1_squared_a = deepcopy(zeros_a)
+ self.spin_lock_omega1_a = deepcopy(zeros_a)
self.phi_ex_struct = deepcopy(zeros_a)
+ self.offset_a = deepcopy(zeros_a)
self.frqs_a = deepcopy(zeros_a)
self.disp_struct = deepcopy(zeros_a)
@@ -447,6 +449,7 @@
# Create special numpy structures.
# Structure of dw. The full and the outer dimensions structures.
self.dw_struct = deepcopy(zeros_a)
+ self.no_nd_struct = ones([self.NO, self.ND], float64)
self.nm_no_nd_struct = ones([self.NM, self.NO, self.ND], float64)
# Structure of r20a and r20b. The full and outer dimensions
structures.
@@ -459,10 +462,11 @@
# Expand relax times.
self.inv_relax_times_a = 1.0 / multiply.outer(
tile(self.relax_times[:,None],(1, 1, self.NS)).reshape(self.NE, self.NS,
self.NM), self.no_nd_struct )
- if model in [MODEL_DPL94]:
+ if model in [MODEL_DPL94, MODEL_TAP03]:
self.r1_a = multiply.outer( self.r1.reshape(self.NE,
self.NS, self.NM), self.no_nd_struct )
-
- # Extract the frequencies to numpy array.
+ self.chemical_shifts_a = multiply.outer(
self.chemical_shifts, self.no_nd_struct )
+
+ # Extract the frequencies to numpy array.
self.frqs_a = multiply.outer(
asarray(self.frqs).reshape(self.NE, self.NS, self.NM), self.no_nd_struct )
# Loop over the experiment types.
@@ -476,7 +480,7 @@
# Extract number of dispersion points.
num_disp_points =
self.num_disp_points[ei][si][mi][oi]
- if model not in [MODEL_DPL94]:
+ if model not in [MODEL_DPL94, MODEL_TAP03]:
# Extract cpmg_frqs and num_disp_points from
lists.
self.cpmg_frqs_a[ei][si][mi][oi][:num_disp_points] =
self.cpmg_frqs[ei][mi][oi]
self.num_disp_points_a[ei][si][mi][oi][:num_disp_points] =
self.num_disp_points[ei][si][mi][oi]
@@ -498,12 +502,14 @@
self.power_a[ei][si][mi][oi][di] =
int(round(self.cpmg_frqs[ei][mi][0][di] * self.relax_times[ei][mi]))
self.tau_cpmg_a[ei][si][mi][oi][di] =
0.25 / self.cpmg_frqs[ei][mi][0][di]
# For R1rho data.
- if model in [MODEL_DPL94]:
+ if model in [MODEL_DPL94, MODEL_TAP03]:
self.disp_struct[ei][si][mi][oi][di] =
1.0
# Extract the frequencies to numpy array.
self.tilt_angles_a[ei][si][mi][oi][di] =
self.tilt_angles[ei][si][mi][oi][di]
self.spin_lock_omega1_squared_a[ei][si][mi][oi][di] =
self.spin_lock_omega1_squared[ei][mi][oi][di]
+
self.spin_lock_omega1_a[ei][si][mi][oi][di] =
self.spin_lock_omega1[ei][mi][oi][di]
+ self.offset_a[ei][si][mi][oi] =
self.offset[ei][si][mi][oi]
if spin_lock_nu1 != None and
len(spin_lock_nu1[ei][mi][oi]):
self.num_disp_points_a[ei][si][mi][oi][di] = num_disp_points
@@ -1908,6 +1914,49 @@
# Once off parameter conversions.
pB = 1.0 - pA
+ # Convert dw from ppm to rad/s. Use the out argument, to pass
directly to structure.
+ multiply( multiply.outer( dw.reshape(self.NE, self.NS),
self.nm_no_nd_struct ), self.frqs_struct, out=self.dw_struct )
+
+ # Reshape R20 to per experiment, spin and frequency.
+ self.r20_struct[:] = multiply.outer( R20.reshape(self.NE, self.NS,
self.NM), self.no_nd_struct )
+
+ # Back calculate the R1rho values.
+ r1rho_TAP03(r1rho_prime=self.r20_struct,
omega=self.chemical_shifts_a, offset=self.offset_a, pA=pA, pB=pB,
dw=self.dw_struct, kex=kex, R1=self.r1_a,
spin_lock_fields=self.spin_lock_omega1_a,
spin_lock_fields2=self.spin_lock_omega1_squared_a,
back_calc=self.back_calc_a, num_points=self.num_disp_points_a)
+
+ # Clean the data for all values, which is left over at the end of
arrays.
+ self.back_calc_a = self.back_calc_a*self.disp_struct
+
+ ## 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.
+ if self.has_missing:
+ # Replace with values.
+ self.back_calc_a[self.mask_replace_blank.mask] =
self.values_a[self.mask_replace_blank.mask]
+
+ # Return the total chi-squared value.
+ return chi2_rankN(self.values_a, self.back_calc_a, self.errors_a)
+
+
+ def func_TP02(self, params):
+ """Target function for the Trott and Palmer (2002) R1rho
off-resonance 2-site model.
+
+ @param params: The vector of parameter values.
+ @type params: numpy rank-1 float array
+ @return: The chi-squared value.
+ @rtype: float
+ """
+
+ # Scaling.
+ if self.scaling_flag:
+ params = dot(params, self.scaling_matrix)
+
+ # Unpack the parameter values.
+ R20 = params[:self.end_index[0]]
+ dw = params[self.end_index[0]:self.end_index[1]]
+ pA = params[self.end_index[1]]
+ kex = params[self.end_index[1]+1]
+
+ # Once off parameter conversions.
+ pB = 1.0 - pA
+
# Initialise.
chi2_sum = 0.0
@@ -1924,7 +1973,7 @@
# Loop over the offsets.
for oi in range(self.num_offsets[0][si][mi]):
# Back calculate the R1rho values.
- r1rho_TAP03(r1rho_prime=R20[r20_index],
omega=self.chemical_shifts[0][si][mi], offset=self.offset[0][si][mi][oi],
pA=pA, pB=pB, dw=dw_frq, kex=kex, R1=self.r1[si, mi],
spin_lock_fields=self.spin_lock_omega1[0][mi][oi],
spin_lock_fields2=self.spin_lock_omega1_squared[0][mi][oi],
back_calc=self.back_calc[0][si][mi][oi],
num_points=self.num_disp_points[0][si][mi][oi])
+ r1rho_TP02(r1rho_prime=R20[r20_index],
omega=self.chemical_shifts[0][si][mi], offset=self.offset[0][si][mi][oi],
pA=pA, pB=pB, dw=dw_frq, kex=kex, R1=self.r1[si, mi],
spin_lock_fields=self.spin_lock_omega1[0][mi][oi],
spin_lock_fields2=self.spin_lock_omega1_squared[0][mi][oi],
back_calc=self.back_calc[0][si][mi][oi],
num_points=self.num_disp_points[0][si][mi][oi])
# 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 di in range(self.num_disp_points[0][si][mi][oi]):
@@ -1938,58 +1987,6 @@
return chi2_sum
- def func_TP02(self, params):
- """Target function for the Trott and Palmer (2002) R1rho
off-resonance 2-site model.
-
- @param params: The vector of parameter values.
- @type params: numpy rank-1 float array
- @return: The chi-squared value.
- @rtype: float
- """
-
- # Scaling.
- if self.scaling_flag:
- params = dot(params, self.scaling_matrix)
-
- # Unpack the parameter values.
- R20 = params[:self.end_index[0]]
- dw = params[self.end_index[0]:self.end_index[1]]
- pA = params[self.end_index[1]]
- kex = params[self.end_index[1]+1]
-
- # Once off parameter conversions.
- pB = 1.0 - pA
-
- # Initialise.
- chi2_sum = 0.0
-
- # Loop over the spins.
- for si in range(self.num_spins):
- # Loop over the spectrometer frequencies.
- for mi in range(self.num_frq):
- # The R20 index.
- r20_index = mi + si*self.num_frq
-
- # Convert dw from ppm to rad/s.
- dw_frq = dw[si] * self.frqs[0][si][mi]
-
- # Loop over the offsets.
- for oi in range(self.num_offsets[0][si][mi]):
- # Back calculate the R1rho values.
- r1rho_TP02(r1rho_prime=R20[r20_index],
omega=self.chemical_shifts[0][si][mi], offset=self.offset[0][si][mi][oi],
pA=pA, pB=pB, dw=dw_frq, kex=kex, R1=self.r1[si, mi],
spin_lock_fields=self.spin_lock_omega1[0][mi][oi],
spin_lock_fields2=self.spin_lock_omega1_squared[0][mi][oi],
back_calc=self.back_calc[0][si][mi][oi],
num_points=self.num_disp_points[0][si][mi][oi])
-
- # 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 di in range(self.num_disp_points[0][si][mi][oi]):
- if self.missing[0][si][mi][oi][di]:
- self.back_calc[0][si][mi][oi][di] =
self.values[0][si][mi][oi][di]
-
- # Calculate and return the chi-squared value.
- chi2_sum += chi2(self.values[0][si][mi][oi],
self.back_calc[0][si][mi][oi], self.errors[0][si][mi][oi])
-
- # Return the total chi-squared value.
- return chi2_sum
-
-
def func_TSMFK01(self, params):
"""Target function for the the Tollinger et al. (2001) 2-site
very-slow exchange model, range of microsecond to second time scale.
r23901 - /branches/disp_spin_speed/target_functions/relax_disp.py