Perfect! This will make it easy to have two 'B14' models. The http://wiki.nmr-relax.com/Tutorial_for_adding_relaxation_dispersion_models_to_relax#Debugging link in the commit message is a little off though.
Cheers, Edward On 6 May 2014 17:25, <[email protected]> wrote: > 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 > [email protected] > > 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 [email protected] 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
