Hi Troels,
Thanks to the lightning quick infrastructure you are putting into
place, we can also simplify the target_functions.relax_disp to
lib.dispersion API. You will notice a lot of code like in this TP02
model:
+ # Once off parameter conversions.
+ pB = 1.0 - pA
This was originally in lib.dispersion (well at least for some of the
models), but I shifted it into the func_*() target functions to speed
the code up, as then the calculation would happen only once rather
than once for each iteration of that massive looping beast you have
killed.
So now that the lib.dispersion functions are only called once per
target function call, all of these 'Once off parameter conversions'
can be shifted back into the lib.dispersion functions. Then the
number of arguments for these functions will drop significantly, as
the {pB, k_AB, k_BA} parameters will no longer need to be passed in.
This will be far more significant for the 3-site models where there
are many, many parameter conversions. This will have the added
benefit of simplifying the use of the lib.dispersion modules outside
of the dispersion target functions, for example with the unit testing.
Cheers,
Edward
On 13 June 2014 07:21, <tlinnet@xxxxxxxxxxxxx> wrote:
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.
_______________________________________________
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
_______________________________________________
relax (http://www.nmr-relax.com)
This is the relax-devel mailing list
relax-devel@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-devel
Re: r23901 - /branches/disp_spin_speed/target_functions/relax_disp.py