Author: bugman
Date: Thu Nov 21 22:50:28 2013
New Revision: 21593
URL: http://svn.gna.org/viewcvs/relax?rev=21593&view=rev
Log:
Removed a latent bug in the 'MMQ 2-site' dispersion model.
This was not being seen but might have caused problems in the future.
Modified:
branches/relax_disp/target_functions/relax_disp.py
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=21593&r1=21592&r2=21593&view=diff
==============================================================================
--- branches/relax_disp/target_functions/relax_disp.py (original)
+++ branches/relax_disp/target_functions/relax_disp.py Thu Nov 21 22:50:28
2013
@@ -936,7 +936,6 @@
self.M0[1] = pB
# Initialise.
- aliased_dwH = 0.0
chi2_sum = 0.0
# Loop over the experiment types.
@@ -947,78 +946,6 @@
for frq_index in range(self.num_frq):
# The R20 index.
r20_index = frq_index + exp_index*self.num_frq +
spin_index*self.num_frq*self.num_exp
-
- # Convert dw from ppm to rad/s.
- dw_frq = dw[spin_index] *
self.frqs[exp_index][spin_index][frq_index]
- dwH_frq = dwH[spin_index] *
self.frqs_H[exp_index][spin_index][frq_index]
-
- # Alias the dw frequency combinations.
- if self.exp_types[exp_index] == EXP_TYPE_CPMG_SQ:
- aliased_dw = dw_frq
- elif self.exp_types[exp_index] ==
EXP_TYPE_CPMG_PROTON_SQ:
- aliased_dw = dwH_frq
- elif self.exp_types[exp_index] == EXP_TYPE_CPMG_DQ:
- aliased_dw = dwH_frq + dw_frq
- elif self.exp_types[exp_index] == EXP_TYPE_CPMG_ZQ:
- aliased_dw = dwH_frq - dw_frq
- elif self.exp_types[exp_index] == EXP_TYPE_CPMG_MQ:
- aliased_dw = dw_frq
- aliased_dwH = dwH_frq
- elif self.exp_types[exp_index] ==
EXP_TYPE_CPMG_PROTON_MQ:
- aliased_dw = dwH_frq
- aliased_dwH = dw_frq
-
- # Back calculate the R2eff values for each experiment
type.
- self.r2eff_mmq[exp_index](M0=self.M0, m1=self.m1,
m2=self.m2, R20A=R20[r20_index], R20B=R20[r20_index], pA=pA, pB=pB,
dw=aliased_dw, dwH=aliased_dwH, k_AB=k_AB, k_BA=k_BA,
inv_tcpmg=self.inv_relax_times[exp_index][frq_index],
tcp=self.tau_cpmg[exp_index][frq_index],
back_calc=self.back_calc[exp_index][spin_index][frq_index],
num_points=self.num_disp_points[exp_index][frq_index],
power=self.power[exp_index][frq_index], n=self.n[exp_index][frq_index])
-
- # 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[exp_index][frq_index]):
- if
self.missing[exp_index][spin_index][frq_index][point_index]:
-
self.back_calc[exp_index][spin_index][frq_index][point_index] =
self.values[exp_index][spin_index][frq_index][point_index]
-
- # Calculate and return the chi-squared value.
- chi2_sum +=
chi2(self.values[exp_index][spin_index][frq_index],
self.back_calc[exp_index][spin_index][frq_index],
self.errors[exp_index][spin_index][frq_index])
-
- # Return the total chi-squared value.
- return chi2_sum
-
-
- def func_mq_CR72(self, params):
- """Target function for the CR72 model extended for MQ CPMG data.
-
- @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]]
- dwH = params[self.end_index[1]:self.end_index[2]]
- pA = params[self.end_index[2]]
- kex = params[self.end_index[2]+1]
-
- # Once off parameter conversions.
- pB = 1.0 - pA
- k_BA = pA * kex
- k_AB = pB * kex
-
- # Initialise.
- chi2_sum = 0.0
-
- # Loop over the experiment types.
- for exp_index in range(self.num_exp):
- # Loop over the spins.
- for spin_index in range(self.num_spins):
- # Loop over the spectrometer frequencies.
- for frq_index in range(self.num_frq):
- # The R20 index.
- r20_index = frq_index + spin_index*self.num_frq
# Convert dw from ppm to rad/s.
dw_frq = dw[spin_index] *
self.frqs[exp_index][spin_index][frq_index]
@@ -1041,6 +968,79 @@
aliased_dw = dwH_frq
aliased_dwH = dw_frq
+ # Back calculate the R2eff values for each experiment
type.
+ self.r2eff_mmq[exp_index](M0=self.M0, m1=self.m1,
m2=self.m2, R20A=R20[r20_index], R20B=R20[r20_index], pA=pA, pB=pB,
dw=aliased_dw, dwH=aliased_dwH, k_AB=k_AB, k_BA=k_BA,
inv_tcpmg=self.inv_relax_times[exp_index][frq_index],
tcp=self.tau_cpmg[exp_index][frq_index],
back_calc=self.back_calc[exp_index][spin_index][frq_index],
num_points=self.num_disp_points[exp_index][frq_index],
power=self.power[exp_index][frq_index], n=self.n[exp_index][frq_index])
+
+ # 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[exp_index][frq_index]):
+ if
self.missing[exp_index][spin_index][frq_index][point_index]:
+
self.back_calc[exp_index][spin_index][frq_index][point_index] =
self.values[exp_index][spin_index][frq_index][point_index]
+
+ # Calculate and return the chi-squared value.
+ chi2_sum +=
chi2(self.values[exp_index][spin_index][frq_index],
self.back_calc[exp_index][spin_index][frq_index],
self.errors[exp_index][spin_index][frq_index])
+
+ # Return the total chi-squared value.
+ return chi2_sum
+
+
+ def func_mq_CR72(self, params):
+ """Target function for the CR72 model extended for MQ CPMG data.
+
+ @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]]
+ dwH = params[self.end_index[1]:self.end_index[2]]
+ pA = params[self.end_index[2]]
+ kex = params[self.end_index[2]+1]
+
+ # Once off parameter conversions.
+ pB = 1.0 - pA
+ k_BA = pA * kex
+ k_AB = pB * kex
+
+ # Initialise.
+ chi2_sum = 0.0
+
+ # Loop over the experiment types.
+ for exp_index in range(self.num_exp):
+ # Loop over the spins.
+ for spin_index in range(self.num_spins):
+ # Loop over the spectrometer frequencies.
+ for frq_index in range(self.num_frq):
+ # The R20 index.
+ r20_index = frq_index + spin_index*self.num_frq
+
+ # Convert dw from ppm to rad/s.
+ dw_frq = dw[spin_index] *
self.frqs[exp_index][spin_index][frq_index]
+ dwH_frq = dwH[spin_index] *
self.frqs_H[exp_index][spin_index][frq_index]
+
+ # Alias the dw frequency combinations.
+ aliased_dwH = 0.0
+ if self.exp_types[exp_index] == EXP_TYPE_CPMG_SQ:
+ aliased_dw = dw_frq
+ elif self.exp_types[exp_index] ==
EXP_TYPE_CPMG_PROTON_SQ:
+ aliased_dw = dwH_frq
+ elif self.exp_types[exp_index] == EXP_TYPE_CPMG_DQ:
+ aliased_dw = dwH_frq + dw_frq
+ elif self.exp_types[exp_index] == EXP_TYPE_CPMG_ZQ:
+ aliased_dw = dwH_frq - dw_frq
+ elif self.exp_types[exp_index] == EXP_TYPE_CPMG_MQ:
+ aliased_dw = dw_frq
+ aliased_dwH = dwH_frq
+ elif self.exp_types[exp_index] ==
EXP_TYPE_CPMG_PROTON_MQ:
+ aliased_dw = dwH_frq
+ aliased_dwH = dw_frq
+
# Back calculate the R2eff values.
r2eff_mq_cr72(r20=R20[r20_index], pA=pA, pB=pB,
dw=aliased_dw, dwH=aliased_dwH, kex=kex, k_AB=k_AB, k_BA=k_BA,
cpmg_frqs=self.cpmg_frqs[exp_index][frq_index],
inv_tcpmg=self.inv_relax_times[exp_index][frq_index],
tcp=self.tau_cpmg[exp_index][frq_index],
back_calc=self.back_calc[exp_index][spin_index][frq_index],
num_points=self.num_disp_points[exp_index][frq_index],
power=self.power[exp_index][frq_index])
r21593 - /branches/relax_disp/target_functions/relax_disp.py