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):