Author: tlinnet
Date: Tue May 6 17:25:23 2014
New Revision: 23013
URL: http://svn.gna.org/viewcvs/relax?rev=23013&view=rev
Log:
Split the func_B14 into full, with a calc function.
This is to prepare for the splitting up of B14, into a full:R2a!=R2b, and
"normal" which is r2a=r2b.
sr #3154: (https://gna.org/support/?3154) Implementation of Baldwin (2014)
B14 model - 2-site exact solution model for all time scales.
This follows the tutorial for adding relaxation dispersion models at:
http://wiki.nmr-relax.com/Tutorial_for_adding_relaxation_dispersion_models_to_relax#Debugging
Modified:
trunk/target_functions/relax_disp.py
Modified: trunk/target_functions/relax_disp.py
URL:
http://svn.gna.org/viewcvs/relax/trunk/target_functions/relax_disp.py?rev=23013&r1=23012&r2=23013&view=diff
==============================================================================
--- trunk/target_functions/relax_disp.py (original)
+++ trunk/target_functions/relax_disp.py Tue May 6 17:25:23 2014
@@ -352,7 +352,7 @@
if model == MODEL_TSMFK01:
self.func = self.func_TSMFK01
if model == MODEL_B14:
- self.func = self.func_B14
+ self.func = self.func_B14_full
if model == MODEL_NS_CPMG_2SITE_3D_FULL:
self.func = self.func_ns_cpmg_2site_3D_full
if model == MODEL_NS_CPMG_2SITE_3D:
@@ -391,6 +391,76 @@
self.func = self.func_ns_mmq_3site_linear
+ def calc_B14_chi2(self, R20A=None, R20B=None, dw=None, pA=None,
kex=None):
+ """Calculate the chi-squared value of the Baldwin (2014) 2-site
exact solution model for all time scales.
+
+
+ @keyword R20A: The R2 value for state A in the absence of
exchange.
+ @type R20A: list of float
+ @keyword R20B: The R2 value for state B in the absence of
exchange.
+ @type R20B: list of float
+ @keyword dw: The chemical shift differences in ppm for each
spin.
+ @type dw: list of float
+ @keyword pA: The population of state A.
+ @type pA: float
+ @keyword kex: The rate of exchange.
+ @type kex: float
+ @return: The chi-squared value.
+ @rtype: float
+ """
+
+ # Once off parameter conversions.
+ pB = 1.0 - pA
+ k_BA = pA * kex
+ k_AB = pB * kex
+ g_fact = 1/sqrt(2)
+
+ # Initialise.
+ chi2_sum = 0.0
+
+ # Loop over the experiment types.
+ for ei in range(self.num_exp):
+ # Loop over the spins.
+ for si in range(self.num_spins):
+ # Loop over the spectrometer frequencies.
+ for mi in range(self.num_frq):
+ # Convert dw from ppm to rad/s.
+ dw_frq = dw[si] * self.frqs[ei][si][mi]
+
+ # Alias the dw frequency combinations.
+ if self.exp_types[ei] == EXP_TYPE_CPMG_SQ:
+ aliased_dw = dw_frq
+ elif self.exp_types[ei] == EXP_TYPE_CPMG_PROTON_SQ:
+ aliased_dw = dw_frq
+
+ # The R20 index.
+ r20_index = mi + si*self.num_frq
+
+ # Repetitive calculations (to speed up calculations).
+ r20a = R20A[r20_index]
+ r20b = R20B[r20_index]
+ deltaR2 = r20a - r20b
+
+ # The Carver and Richards (1972) alpha_minus short
notation.
+ alpha_m = deltaR2 + k_AB - k_BA
+ zeta = 2 * aliased_dw * alpha_m
+ Psi = alpha_m**2 + 4 * k_BA * k_AB - aliased_dw**2
+
+ # Back calculate the R2eff values.
+ r2eff_B14(r20a=r20a, r20b=r20b, deltaR2=deltaR2,
alpha_m=alpha_m, pA=pA, pB=pB, dw=dw_frq, zeta=zeta, Psi=Psi,
g_fact=g_fact, kex=kex, k_AB=k_AB, k_BA=k_BA, ncyc=self.power[ei][mi],
inv_tcpmg=self.inv_relax_times[ei][mi], tcp=self.tau_cpmg[ei][mi],
back_calc=self.back_calc[ei][si][mi][0],
num_points=self.num_disp_points[ei][si][mi][0])
+
+ # 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[ei][si][mi][0]):
+ if self.missing[ei][si][mi][0][di]:
+ self.back_calc[ei][si][mi][0][di] =
self.values[ei][si][mi][0][di]
+
+ # Calculate and return the chi-squared value.
+ chi2_sum += chi2(self.values[ei][si][mi][0],
self.back_calc[ei][si][mi][0], self.errors[ei][si][mi][0])
+
+ # Return the total chi-squared value.
+ return chi2_sum
+
+
def calc_CR72_chi2(self, R20A=None, R20B=None, dw=None, pA=None,
kex=None):
"""Calculate the chi-squared value of the 'NS CPMG 2-site star'
models.
@@ -763,7 +833,7 @@
raise RelaxError("The '%s' CPMG model is not compatible
with the '%s' experiment type." % (self.model, self.exp_types[0]))
- def func_B14(self, params):
+ def func_B14_full(self, params):
"""Target function for the Baldwin (2014) 2-site exact solution
model for all time scales.
This assumes that pA > pB, and hence this must be implemented as a
constraint.
@@ -786,56 +856,8 @@
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
- g_fact = 1/sqrt(2)
-
- # Initialise.
- chi2_sum = 0.0
-
- # Loop over the experiment types.
- for ei in range(self.num_exp):
- # Loop over the spins.
- for si in range(self.num_spins):
- # Loop over the spectrometer frequencies.
- for mi in range(self.num_frq):
- # Convert dw from ppm to rad/s.
- dw_frq = dw[si] * self.frqs[ei][si][mi]
-
- # Alias the dw frequency combinations.
- if self.exp_types[ei] == EXP_TYPE_CPMG_SQ:
- aliased_dw = dw_frq
- elif self.exp_types[ei] == EXP_TYPE_CPMG_PROTON_SQ:
- aliased_dw = dw_frq
-
- # The R20 index.
- r20_index = mi + si*self.num_frq
-
- # Repetitive calculations (to speed up calculations).
- r20a = R20A[r20_index]
- r20b = R20B[r20_index]
- deltaR2 = r20a - r20b
-
- # The Carver and Richards (1972) alpha_minus short
notation.
- alpha_m = deltaR2 + k_AB - k_BA
- zeta = 2 * aliased_dw * alpha_m
- Psi = alpha_m**2 + 4 * k_BA * k_AB - aliased_dw**2
-
- # Back calculate the R2eff values.
- r2eff_B14(r20a=r20a, r20b=r20b, deltaR2=deltaR2,
alpha_m=alpha_m, pA=pA, pB=pB, dw=dw_frq, zeta=zeta, Psi=Psi,
g_fact=g_fact, kex=kex, k_AB=k_AB, k_BA=k_BA, ncyc=self.power[ei][mi],
inv_tcpmg=self.inv_relax_times[ei][mi], tcp=self.tau_cpmg[ei][mi],
back_calc=self.back_calc[ei][si][mi][0],
num_points=self.num_disp_points[ei][si][mi][0])
-
- # 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[ei][si][mi][0]):
- if self.missing[ei][si][mi][0][di]:
- self.back_calc[ei][si][mi][0][di] =
self.values[ei][si][mi][0][di]
-
- # Calculate and return the chi-squared value.
- chi2_sum += chi2(self.values[ei][si][mi][0],
self.back_calc[ei][si][mi][0], self.errors[ei][si][mi][0])
-
- # Return the total chi-squared value.
- return chi2_sum
+ # Calculate and return the chi-squared value.
+ return self.calc_B14_chi2(R20A=R20A, R20B=R20B, dw=dw, pA=pA,
kex=kex)
def func_CR72(self, params):
_______________________________________________
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: r23013 - /trunk/target_functions/relax_disp.py