Re: numerical cpmg fitRe: numerical cpmg fitHi Edward, I understand your point of fitting pA and kex instead of keg and kge (or kAB/kBA). Please find attached a program with the convention of kex and pB. The states A and B are called G and E (ground/excited), so the parameters are generally associated with E and G instead of A and B. |
# -*- coding: ISO-8859-1 -*-
from math import*
import shutil
import sys
import pylab
from pylab import *
from scipy import *
from scipy.linalg import expm
from scipy.optimize import fmin
from scipy.optimize import leastsq
from random import gauss
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.patches as patches
import matplotlib.path as path
import matplotlib.pylab as pyl
import numpy as np
import scipy.integrate as itg
import random
import scipy.stats as stats
import math
import xlrd
import os
import matplotlib.gridspec as gridspec
import cmath
import math
print ('------IN------')
wb = xlrd.open_workbook('myworkbook.xls')
sh = wb.sheet_by_index(0)
PathOutputFile=str(sh.cell(25,1).value)
if not os.path.exists(PathOutputFile):
os.mkdir(PathOutputFile)
nb_fields=int(sh.cell(0,1).value)
in_nlin={}
in_freq={}
my_fields=np.zeros(nb_fields)
for kk in range(nb_fields):
my_fields[kk]=float(sh.cell(2,1+kk).value)
in_nlin[kk]=str(sh.cell(3,1+kk).value)
in_freq[kk]=str(sh.cell(4,1+kk).value)
nbCALLfunc=int(sh.cell(6,1).value)
nbMC=int(sh.cell(7,1).value)
print ('------CALCULATION START------')
'''residues listing'''
paramCPMG0={} #paramCPMG0 contains the initial parameters
which will be fitted for all residues
my_residues=[] #my_residues contains the list of the names
of the residues
for res in range(int(sh.cell(9,1).value)):
my_residues.append(str(sh.cell(10,1+res).value))
paramCPMG0[res]=[]
for param in range(3+nb_fields):
paramCPMG0[res].append(float(sh.cell(12+param,1+res).value))
Tc=float(sh.cell(1,1).value) #Tc is the constant-time of the CPMG
experiment
f={} #f contains the name of the nlin.tab for each
field which contain the data (HEIGHT, DHEIGHT, A_Z0...)
for kk in range(nb_fields):
f[kk]='f'+str(kk)
for kk in range(nb_fields):
f[kk] = open(in_nlin[kk],'r')
'''end of the header number line'''
lignes = f[0].readlines()
l_vars=lignes[0].split() #header of the nlin.tab file, read once from
the nlin.tab file of the first field
for kk in range(nb_fields):
f[kk].close()
'''colunm nulber of the height and the intensities ratio'''
for j in range(len(l_vars)): #header reading, colunm data recovery
if l_vars[j] == 'HEIGHT':
num_col_HEIGHT=j
if l_vars[j] == 'Z_A0':
ncol_Z_A0=j
'''nlin.tab file reading'''
my_intensities={} #contains the arrays of the intensities of
all residues
my_noises={} #contains the arrays of the calculated noises
from experiment noise of all residues
my_R2eff={} #contains the arrays of the calculated R2eff
from my_intensities of all residues
for kk in range(nb_fields):
f[kk] = open(in_nlin[kk],'r')
for kk in range(nb_fields):
lignes = f[kk].readlines()[1:] ## reading of the kk field nlin.tab
ii=0
for line in lignes:
intensities=[]
noises=[]
R2eff=[]
liste = line.split()
for k in range(len(liste)-ncol_Z_A0+1):
intensity=float(liste[k-1+ncol_Z_A0])*float(liste[num_col_HEIGHT-1]) #there
is a lag
noise=float(liste[num_col_HEIGHT])
#num col of DHEIGHT (because of the lag)
intensities.append(intensity)
noises.append(noise)
R=-1./Tc*log(abs(float(liste[k-1+ncol_Z_A0])))
R2eff.append(R)
my_intensities[kk, ii]=intensities
my_noises[kk, ii]=noises
my_R2eff[kk, ii]=R2eff
ii +=1
freq={} #contains the arrays of the numbers of pulses
of each CPMG experiment
'''reading of the frequencies'''
for jj in range(nb_fields):
f_freqCPMG=open(in_freq[jj], 'r')
lignes = f_freqCPMG.readlines()
freq[jj]=np.zeros(len(lignes))
kk=0
for line in lignes:
liste = line.split()
freq[jj][kk]=float(liste[1]) #the number of pulses must be
on the second column of the freq.txt file
kk+=1
f_freqCPMG.close()
'''Error bars calculation'''
my_errors={}
for jj in range(nb_fields):
for k in range (len(my_residues)):
error=[]
for kk in range(len(freq[jj])):
error.append(abs(my_noises[jj,
k][kk]/(Tc*float(my_intensities[jj, k][kk]))))
my_errors[jj,k]=error
for kk in range(nb_fields):
f[kk].close()
print('------files reading: ok------')
'''definition de la matrice d echange'''
def RCPMG_3D(R1E,R1G, R2E, R2G,df, kGE, kEG):
fG=0; fE=df; IGeq=kEG/(kEG+kGE); IEeq=kGE/(kEG+kGE)
temp=np.matrix([
[ 0, 0, 0, 0, 0, 0,
0],
[ 0, -R2G-kGE, -fG, 0, kEG, 0,
0],
[ 0, fG, -R2G-kGE, 0, 0, kEG,
0],
[2.0*R1G*IGeq, 0, 0, -R1G-kGE, 0, 0,
kEG],
[ 0, kGE, 0, 0, -R2E-kEG, -fE,
0],
[ 0, 0, kGE, 0, fE, -R2E-kEG,
0],
[2.0*R1E*IEeq, 0, 0, kGE, 0, 0,
-R1E-kEG]])
return temp
def RCPMG_2D(R2E, R2G,df, kGE, kEG):
#fG=0; fE=fG+df;
fG=0 ; fE=fG+df
temp=np.matrix([
[-R2G-kGE, -fG, kEG, 0],
[ fG, -R2G-kGE, 0, kEG],
[ kGE, 0, -R2E-kEG, -fE],
[ 0, kGE, fE, -R2E-kEG]])
return temp
'''rotational matrix '''
def R180_3Dx(a):
theta=180.0
R=np.matrix([[1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0], [0.0, 1.0, 0.0, 0.0, 0.0,
0.0, 0.0], [0.0, 0.0,cos(2*math.pi*(theta+a)/360),
-sin(2*math.pi*(theta+a)/360), 0.0, 0.0, 0.0], [0.0, 0.0,
sin(2*math.pi*(theta+a)/360),cos(2*math.pi*(theta+a)/360), 0.0, 0.0, 0.0],
[0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0], [0.0, 0.0, 0.0, 0.0,
0.0,cos(2*math.pi*(theta+a)/360), -sin(2*math.pi*(theta+a)/360)], [0.0, 0.0,
0.0, 0.0, 0.0, sin(2*math.pi*(theta+a)/360),cos(2*math.pi*(theta+a)/360)]])
return R
def R180x_2D(a):
R=np.matrix([
[1, 0, 0.0, 0.0],
[0,-1, 0.0, 0.0],
[0.0, 0.0,1, 0],
[0.0, 0.0, 0,-1]])
return R
def mpower(x,y):
'''Compute x raised to the power y when x is a square matrix and y
is a scalar.'''
s = shape(x)
if len(s) != 2 or s[0] != s[1]:
raise ValueError('matrix must be square')
if y == 0:
return eye(s[0])
[e,v] = eig(x)
d = diag(e)
return dot(dot(v,d**y),inv(v))
numbercallofe=0
def incremente():
global numbercallofe
numbercallofe +=1
'''CPMG process for R calculation'''
fit_function={}
for kk in range(5):
fit_function[kk]='func'+str(kk)
'''function which simulates a R2eff profile (exchange model) according to the
number of pulses during the constant time:'''
def MakeFunction(a):
if a==1:
print('optimization with the function 3D_slow')
def func1(r1e,r1g,r2e,r2g,df,kex,pg, nb_pulse): #R2E, R2G,df, kGE,
kEG,nb_pulse,tcpmg
kge=(1-pg)*kex; keg=pg*kex
Rcalc_array=np.zeros(len(nb_pulse))
Rtemp=RCPMG_3D(r1e,r1g,r2e,r2g,df,kge,keg)
M0=np.matrix([[0.5],[pg],[0.0],[0.0],[1-pg],[0.0],[0.0]])
for k in range(len(nb_pulse)):
tcp=Tc/(4.0*nb_pulse[k])
nb_total = 2*int(nb_pulse[k])
Mint=M0
Rexpo=expm(Rtemp*tcp)
for kk in range(0, nb_total):
Mint=Rexpo*Mint
Mint=R180_3Dx(0)*Mint
Mint=Rexpo*Mint
Rcalc_array[k]=-(1.0/Tc)*log(fabs((Mint[1])/pg));
return Rcalc_array
if a==2:
# This function is actually not required. The above "slow" version is
enough.
# The difference is just that one evaluates the magnetizations on both
states, while
# this one only considers the major state. Taking only the major state
is always correct, while
# taking both states is only correct in the case of fast exchange. In
other words, the "fast" version is not needed, because in cases where this is
correct the other one is correct, too.
print('optimization with the function 3D_fast')
def func1(r1e,r1g,r2e,r2g,df,kex,pg, nb_pulse): #R2E, R2G,df, kGE,
kEG,nb_pulse,tcpmg
kge=(1-pg)*kex; keg=pg*kex
Rcalc_array=np.zeros(len(nb_pulse))
Rtemp=RCPMG_3D(r1e,r1g,r2e,r2g,df,kge,keg)
M0=np.matrix([[0.5],[pg],[0.0],[0.0],[1-pg],[0.0],[0.0]])
for k in range(len(nb_pulse)):
tcp=Tc/(4.0*nb_pulse[k])
nb_total = 2*int(nb_pulse[k])
Mint=M0
Rexpo=expm(Rtemp*tcp)
for kk in range(0, nb_total):
Mint=Rexpo*Mint
Mint=R180_3Dx(0)*Mint
Mint=Rexpo*Mint
Rcalc_array[k]=-(1.0/Tc)*log(fabs((Mint[1]+Mint[4])));
return Rcalc_array
if a==3:
print('optimization with the function 2D_slow')
def func1(r1e,r1g,r2e,r2g,df,kex,pg, nb_pulse): #ATTENTION nb_pulse (not
nu_cpmg)
kge=(1-pg)*kex; keg=pg*kex
Rcalc_array=np.zeros(len(nb_pulse))
Rtemp=RCPMG_2D(r2e, r2g,df, kge, keg)
M0=np.matrix([[pg],[0.0],[1-pg],[0.0]])
for k in range(len(nb_pulse)):
tcp=Tc/(4.0*nb_pulse[k])
nb_total = 2*int(nb_pulse[k])
Mint=M0
Rexpo=expm(Rtemp*tcp)
for kk in range(0, nb_total):
Mint=Rexpo*Mint
Mint=R180x_2D(0)*Mint
Mint=Rexpo*Mint
Rcalc_array[k]=-(1.0/Tc)*log(fabs(Mint[0]/pg));
return Rcalc_array
if a==4:
# This function is not required. the above "slow" version is enough.
print('optimization with the function 2D_fast')
def func1(r1e,r1g,r2e,r2g,df,kex,pg, nb_pulse): #ATTENTION nb_pulse (not
nu_cpmg)
kge=(1-pg)*kex; keg=pg*kex
Rcalc_array=np.zeros(len(nb_pulse))
Rtemp=RCPMG_2D(r2e, r2g,df, kge, keg)
M0=np.matrix([[pg],[0.0],[1-pg],[0.0]])
for k in range(len(nb_pulse)):
tcp=Tc/(4.0*nb_pulse[k])
nb_total = 2*int(nb_pulse[k])
Mint=M0
Rexpo=expm(Rtemp*tcp)
for kk in range(0, nb_total):
Mint=Rexpo*Mint
Mint=R180x_2D(0)*Mint
Mint=Rexpo*Mint
Rcalc_array[k]=-(1.0/Tc)*log(fabs((Mint[0]+Mint[2])));
return Rcalc_array
if a==5:
#This function is EXACT, just as the explicit Bloch-McConnell numerical
treatments above.
# It comes from a Maple derivation based on the Bloch-McConnell
equations. It is much faster
# than the numerical Bloch-McConnell solution.
# It was derived by Nikolai Skrynnikov, and is provided with his
permission.
print('optimization with the function 2D_num')
def func1(r1e,r1g,r2e,r2g,df,kex,pg, nb_pulse):
Rcalc_array=np.zeros(len(nb_pulse))
kge=(1-pg)*kex; keg=pg*kex
M0=np.matrix([[pg],[1-pg]])
tcp=Tc/(nb_pulse)
N = Tc/tcp
Mint=M0
Kb=keg;
Ka=kge;
Rinf=r2g;
dw=df #*math.pi*2
j=complex(0,1);
t3 = j;
t4 = t3*dw;
t5 = Kb*Kb;
t8 = 2*t3*Kb*dw;
t10 = 2*Kb*Ka;
t11 = dw*dw;
t14 = 2*t3*Ka*dw;
t15 = Ka*Ka;
t17 = cmath.sqrt(t5-t8+t10-t11+t14+t15);
t21 = np.exp((-Kb+t4-Ka+t17)*tcp/4);
t22 = 1./t17;
t28 = np.exp((-Kb+t4-Ka-t17)*tcp/4);
t31 = t21*t22*Ka-t28*t22*Ka;
t33 = cmath.sqrt(t5+t8+t10-t11-t14+t15);
t34 = Kb+t4-Ka+t33;
t37 = np.exp((-Kb-t4-Ka+t33)*tcp/2);
t39 = 1./t33;
t41 = Kb+t4-Ka-t33;
t44 = np.exp((-Kb-t4-Ka-t33)*tcp/2);
t47 = t34*t37*t39/2-t41*t44*t39/2;
t49 = Kb-t4-Ka-t17;
t51 = t21*t49*t22;
t52 = Kb-t4-Ka+t17;
t54 = t28*t52*t22;
t55 = -t51+t54;
t60 = t37*t39*Ka-t44*t39*Ka;
t62 = t31*t47+t55*t60/2;
t63 = 1./Ka;
t68 = -t52*t63*t51/2+t49*t63*t54/2;
t69 = t62*t68/2;
t72 = t37*t41*t39;
t76 = t44*t34*t39;
t78 = -t34*t63*t72/2+t41*t63*t76/2;
t80 = -t72+t76;
t82 = t31*t78/2+t55*t80/4;
t83 = t82*t55/2;
t88 = t52*t21*t22/2-t49*t28*t22/2;
t91 = t88*t47+t68*t60/2;
t92 = t91*t88;
t95 = t88*t78/2+t68*t80/4;
t96 = t95*t31;
t97 = t69+t83;
t98 = t97*t97;
t99 = t92+t96;
t102 = t99*t99;
t108 = t62*t88+t82*t31;
t112 = np.sqrt(t98-2*t99*t97+t102+4*(t91*t68/2+t95*t55/2)*t108);
t113 = t69+t83-t92-t96-t112;
t115 = N/2;
t116 = np.power(t69/2+t83/2+t92/2+t96/2+t112/2,t115)
t118 = 1./t112;
t120 = t69+t83-t92-t96+t112;
t122 = np.power(t69/2+t83/2+t92/2+t96/2-t112/2,t115)
t127 = 1./t108;
t139 =
1./(Ka+Kb)*((-t113*t116*t118/2+t120*t122*t118/2)*Kb+(-t113*t127*t116*t120*t118/2+t120*t127*t122*t113*t118/2)*Ka/2);
intensity0 = pg
intensity = np.real(t139)*np.exp(-Tc*Rinf)
Rcalc=-(1./Tc)*np.log(intensity/intensity0);
return Rcalc
if a==6:
# This function does the numerical solution of the BlochMcConnell
equations. It neglects the z component. The difference to the above "2D_slow
and 2D_fast" is that it uses complex matrices, this reduces the number of
matrix manipulations and thus speeds up the process.
print('optimization with the function 2D complex')
def func1(r1e,r1g,r2e,r2g,df,kex,pg, nb_pulse): # def
func1(R2E,R2G,fg,kge,keg, cpmgfreq,tcpmg):
kge=(1-pg)*kex; keg=pg*kex
Rcalc_array=np.zeros(len(nb_pulse))
Rr = -1.*np.matrix([[r2g, 0.],[0., r2e]])
Rex = -1.*np.matrix([[kge,-keg],[-kge, keg]])
RCS = complex(0.,-1.)*np.matrix([[0. ,0.],[0.,df]])
R = Rr + Rex + RCS
cR = np.conj(R)
M0=np.matrix([[keg/(keg+kge)],[kge/(keg+kge)]])
for k in range(len(nb_pulse)):
tcp=Tc/(4.*nb_pulse[k]) #tcp=Tc/(4.0*nb_pulse[k])
prop_2 =
np.dot(np.dot(expm(R*tcp),expm(cR*2.*tcp)),expm(R*tcp))
prop_total =mpower(prop_2,nb_pulse[k])
Moft = prop_total*M0
Mgx=Moft[0].real/M0[0]
Rcalc=-(1./Tc)*math.log(Mgx);
Rcalc_array[k]=Rcalc
return Rcalc_array
else: print ('ERROR! enter a number between 1 and 5 to chose the
optimization function')
return func1
print('optimization function number: '),int(sh.cell(8,1).value)
fp0=MakeFunction(int(sh.cell(8,1).value))
'''function which simulates a flat R2eff profile (non exchange model)
according to the number of pulses during the constant time:'''
def fpS0(v0,nb_pulse): #v0=R2inf
Rcalc_array=np.zeros(len(nb_pulse))
for k in range(len(nb_pulse)):
Rcalc_array[k]=v0
return Rcalc_array
'''khi2 calculation function between 2 same-size arrays'''
def chi(Rc, Rd, data_err):
KHI=(Rc-Rd)/data_err
return KHI
'''error function for the exchange model'''
def e(p, x, my_data, R1E, R1G, R2E, data_err):
KHI2=np.zeros(0) #len(my_data[0])=nb of R2eff
if p[1]<0: p[1]=-p[1]
if p[2]<0: p[2]=-p[2]
for kk in range(len(my_data)): #len(my_data)=nb of fields
KHI2=np.append(KHI2, chi(fp0(R1E,R1G,R2E,
p[3+kk],p[0]*my_fields[kk],p[1], p[2],x[kk]), my_data[kk], data_err[kk]))
return KHI2
'''error function for the non exchange model'''
def eS(p, x, my_data, R1E, R1G, R2E, data_err):
KHI2=np.zeros(0)
for kk in range(len(my_data)): #len(my_data)=nb of fields
KHI2=np.append(KHI2,chi(fpS0(p[kk], x[kk]), my_data[kk],
data_err[kk]))
return KHI2
KHI2={} #contains the KHI2 of the exchange model fit compared with
data, for each MonteCarlo simulation, for all residues
pSTAT=[] #contains the pValue for all residues
f=open(PathOutputFile+str('fitting.txt'), 'w')
f.close()
for kres in range(len(my_residues)):
'''non fitted parameters'''
R1E=float(sh.cell(21,1+res).value)
R1G=float(sh.cell(22,1+res).value)
R2E=float(sh.cell(23,1+res).value)
name=my_residues[kres]
print('------R2eff fitting calculation for residue: '), name
data_EXP={}
x={}
data_EXP_error={}
paramCPMG_0=paramCPMG0[kres]
index=np.zeros(nb_fields)
for nf in range(nb_fields):
data_EXP[nf]=np.delete(my_R2eff[nf, kres], 0)
#deletion of the first point, from the reference spectrum (no pulse applied)
data_EXP_error[nf]=np.delete(my_errors[nf, kres], 0)
#deletion of the first point, from the reference spectrum (no pulse applied)
x[nf]=np.delete(freq[nf], 0)
#deletion of the first point, from the reference spectrum (no pulse applied)
for kk in range(len(x[nf])):
#search for the repetition of the first point
if x[nf][kk]==0:
index[nf]=kk
if index[nf]!=0:
x[nf]=np.delete(x[nf], index[nf])
data_EXP[nf]=np.delete(data_EXP[nf], index[nf])
data_EXP_error[nf]=np.delete(data_EXP_error[nf],
index[nf])
'''param=[ df, kex, pg, R2G_f1, R2G_f2]'''
p, success =leastsq(e, paramCPMG_0, args=(x,data_EXP, R1E, R1G,
R2E,data_EXP_error), maxfev=nbCALLfunc)
pS, success =leastsq(eS, paramCPMG_0[3:], args=(x, data_EXP, R1E,
R1G, R2E, data_EXP_error), maxfev=nbCALLfunc)
print('initial parameters: '), paramCPMG_0
print('fitted parameters: '), p
yS={} #contains the fitted R2eff with
non-exchange model for all fields for one residue
y={} #contains the fitted R2eff with
exchange model for all fields for one residue
chi2EMF_n={} #contains the khi2 between R2eff data
and fitted R2eff for the exchange model
chi2SMF_n={} #contains the khi2 between R2eff data
and fitted R2eff for the non exchange model
for nf in range(nb_fields):
yS[nf]=fpS0(pS[nf], x[nf])
y[nf]=fp0(R1E,R1G,R2E, p[3+nf],p[0]*my_fields[nf],p[1],
p[2],x[nf])
chi2EMF_n[nf]=chi(y[nf], data_EXP[nf],
data_EXP_error[nf])*chi(y[nf], data_EXP[nf], data_EXP_error[nf])
chi2SMF_n[nf]=chi(yS[nf], data_EXP[nf],
data_EXP_error[nf])*chi(yS[nf], data_EXP[nf], data_EXP_error[nf])
'''Ftest'''
KHI2_emf=sum(chi2EMF_n[0])+sum(chi2EMF_n[1]) #sum of the khi2 for
only 2 fields, exchange model
KHI2_smf=sum(chi2SMF_n[0])+sum(chi2SMF_n[1]) #sum of the khi2 for
only 2 fields, non exchange model
print('KHI2 for the model with exchange : '), KHI2_emf
print('KHI2 for the model without exchange : '), KHI2_smf
P_SMF=nb_fields; P_EMF=nb_fields+3
###number of fitted parameters
Diff_SMF=len(x[0])-P_SMF; Diff_EMF=len(x[0])-P_EMF
Ftest=stats.f_value(KHI2_smf, KHI2_emf, Diff_SMF, Diff_EMF) # F
test
Pvalue=stats.f.sf(Ftest, Diff_SMF-Diff_EMF, Diff_EMF) # p
value calculation
print('Pvalue: '), Pvalue
pSTAT.append(Pvalue)
'''data plot and fit curves plot '''
parameters=[round(p[0],2), round(p[1],2), round(p[2],2)]
for nf in range(nb_fields):
parameters.append(round(p[3+nf],2))
parameters.append(round(KHI2_emf))
parameters.append(round(KHI2_smf))
parameters.append(Pvalue)
parametersNAME=['delta(ppm)', 'kex (s-1)', 'pG ']
for nf in range(nb_fields):
parametersNAME.append(str('R2G_')+str(int(my_fields[nf]))+str('(s-1)'))
parametersNAME.append('khi2(ex)')
parametersNAME.append('khi2(non ex)')
parametersNAME.append('Pvalue')
fig = plt.figure()
PlotData={}
NameFields=[]
for nf in range(nb_fields):
if index[nf]!=0:
plot(x[nf][:index[nf]], y[nf][:index[nf]]) #x[0][:22], y[:22]
else: plot(x[nf], y[nf])
plot(x[nf], yS[nf])
PlotData[nf]=errorbar(x[nf],data_EXP[nf],
yerr=data_EXP_error[nf], fmt='o')
NameFields.append(str(my_fields[nf])+str('MHz'))
plotdata=[]
for nf in range(nb_fields):
plotdata.append(PlotData[nf])
legend(plotdata, NameFields)
f=open(PathOutputFile+str('fitting.txt'), 'a')
if kres == 0:
f.write('param:')
f.write(' ')
for jj in range(len(parameters)):
f.write(parametersNAME[jj])
f.write(' ')
f.write('\n')
f.write(name)
f.write(' ')
for jj in range(len(parameters)):
f.write(str(parameters[jj]))
if jj != len(parameters)-1:
f.write(' ')
plt.title(name)
plt.savefig(PathOutputFile+name)
#plt.pause(0.5)
f.close()
if nbMC!=0: #test if MC calculation is expected
nbCall=100 #number of e function calls to fit the
pseudo-data, calculated from the fit curve + gaussian noise
print ('-----MonteCarlo calculations:-----'), nbMC
R_pdE={}
R_pdS={}
for nf in range(nb_fields):
R_pdE[nf]=fp0(R1E,R1G,R2E,p[3+nf],p[0]*my_fields[nf],p[1],p[2], x[nf])
R_pdS[nf]=fpS0(pS[0], x[nf])
KHI2[kres]=[]
pMC={}
for jj in range(nbMC):
#if nbMC%10==0:
# print ('MonteCarlo loading: '), jj
loading=0
data={}
for nf in range(nb_fields):
data[nf]=np.zeros(len(R_pdE[nf]))
for j in range(len(data[nf])):
data[nf][j]=gauss(R_pdE[nf][j],
data_EXP_error[nf][j])
pE, success=leastsq(e, p, args=(x,data, R1E, R1G,
R2E,data_EXP_error), maxfev=nbCall)
pSmc, success =leastsq(eS, pS, args=(x, data, R1E,
R1G, R2E, data_EXP_error), maxfev=nbCall)
numbercallofe=0
yE={}
k2={}
for nf in range(nb_fields):
yE[nf]=fp0(R1E,R1G,R2E,pE[3+nf],pE[0]*my_fields[nf],pE[1],pE[2], x[nf])
k2[nf]=chi(yE[nf],
data[nf],data_EXP_error[nf])*chi(yE[nf], data[nf], data_EXP_error[nf])
KHI2[kres].append(sum(k2[0])+sum(k2[1]))
pMC[jj]=pE
# print ('KHI2 for the residue '), name
# print('is: '), KHI2[kres]
fig = plt.figure()
pylab.hist(KHI2[kres], bins=nbMC)
DeviationStandard=np.std(KHI2[kres])
figtext(0.7, 0.75, 'Standard deviation:')
figtext(0.7, 0.72, round(DeviationStandard, 2))
figtext(0.7, 0.67, 'Khi2 exp:')
figtext(0.7, 0.64, round(KHI2_emf, 2))
legend([name])
plt.title(str('Khi2 for ')+str(nbMC)+str(' MonteCarlo
simulations'))
plt.savefig(PathOutputFile+name+str('_Khi2'))
# plt.pause(0.5)
pp={}
for NbPar in range(3+nb_fields):
pp[NbPar]=np.zeros(nbMC)
for jj in range(nbMC):
pp[NbPar][jj]=pMC[jj][NbPar]
fig = plt.figure()
pylab.hist(pp[NbPar], bins=nbMC)
DeviationStandard=np.std(pp[NbPar])
figtext(0.7, 0.75, 'Standard deviation:')
figtext(0.7, 0.72, round(DeviationStandard, 2))
figtext(0.7, 0.67, parametersNAME[NbPar]+str(':'))
figtext(0.7, 0.64,parameters[NbPar] )
legend([name])
plt.title(str(sh.cell(12+NbPar,0).value)+str('
(')+str(nbMC)+str(' MonteCarlo simulations)'))
plt.savefig(PathOutputFile+name+str('_')+str(NbPar))
# plt.pause(0.5)
else: print('-----No MonteCarlo calculation-----')
You will find in the attached program six different implementations, that can be chosen with the variable a set to 1-6. In practice, however, we are using only the "2D complex" numerical approach (a=6), or the Maple-derived approach (a=5). The latter is faster, and is most-often used. It has not been published, though. The numerical Bloch-McConnell treatment is described in many introductory NMR books, and also in e.g. McConnell, H. J Chem Phys 1958, 28, 430–431. (the original reference) Palmer, A. G. Chem Rev 2004, 104, 3623–3640. Palmer, A. G.; Massi, F. Chem Rev 2006, 106, 1700–1719. Baldwin, A. J.; Kay, L. E. J Biomol NMR 2013. The latter two are R1rho, but the formalism of Bloch-McConnell is the same. Paul Hi Paul, |