Commit d198ca29 authored by Jorn Bruggeman's avatar Jorn Bruggeman
Browse files

compute C_e, B_e and fractionation internally; allow organelle-specific CA

parent dd424697
......@@ -23,3 +23,28 @@ class DataFile(object):
def grabColumn(self,name,dtype=float):
icol = self.labels.index(name)
return numpy.array([d[icol] for d in self.data],dtype=dtype)
def grabColumnMasked(self,name):
icol = self.labels.index(name)
coldata = numpy.ma.array(numpy.empty((len(self.data),)),mask=True)
for i,rowdata in enumerate(self.data):
if rowdata[icol]!='NA': coldata[i] = float(rowdata[icol])
return coldata
def renameColumn(self,oldname,newname):
icol = self.labels.index(oldname)
self.labels[icol] = newname
def addColumn(self,label,data):
assert len(data)==len(self.data)
self.labels.append(label)
for row,newvalue in zip(self.data,data): row.append(newvalue)
def write(self,path):
with open(path,'w') as f:
f.write('%s\n' % ('\t'.join(self.labels),))
for row in self.data:
values = [str(value) for value in row]
for i,value in enumerate(row):
if getattr(value,'mask',False): values[i] = 'NA'
f.write('%s\n' % ('\t'.join(values)))
import argparse
parser = argparse.ArgumentParser(description='Fit model to data (data+model defined in problem.py).')
parser.add_argument('--skipfit', action='store_true',help='Skip fitting and show result of current parameterization in model.py')
parser.add_argument('--save', help='Path to save parameter estimates to.')
args = parser.parse_args()
import sys
......@@ -9,30 +10,33 @@ sys.path.append('..')
import numpy
import problem
# Standard settings (minimum, mximum, log scale) for different parameters.
# Standard settings (minimum, maximum, log scale) for different parameters.
pardata = {'alpha':(1e-3,1e5,True),
'eps_FIX':(-30.,-10.,False),
'eps_FIX':(-30.,-5.,False),
'FB_vu_per_F_CAL':(0.,2.,False),
'k_f':(1e-1,1e1,True),
'c':(1e-3,1e1,True),
'n':(0.,5.,False),
'k_ei':(1e-5,1e-3,True),
'k_ei':(1e-5,2e-3,True),
'k_eiB':(1e-7,1e-3,True),
'k_plus_i_CA':(1e-6,3000.,True),
'Km':(0.01,0.2,False),
'f_chl':(0.0001,1.,True),
'f_cv':(0.0001,1.,True),
'T_a':(1e-5,1.,True),
'T_b':(1e-5,1.,True),
'k_plus_CA':(.1,1e8,True),
'k_plus_CA':(.01,1e6,True),
'k_plus_i_CA':(.01,1e6,True),
'k_plus_x_CA':(.01,1e6,True),
'k_plus_v_CA':(.01,1e6,True),
'pH_i':(7.,8.2,False),
'pH_x':(7.,8.2,False),
'pH_v':(7.,8.2,False),
}
# Create optimization problem (this will also create the model itself, as member "model")
pars = 'k_ei','T_a','T_b','eps_FIX','k_plus_CA' #,'pH_i','pH_x','pH_v'
pars = 'k_ei','T_a','T_b','eps_FIX','k_plus_i_CA','k_plus_x_CA','k_plus_v_CA','pH_i','pH_x','pH_v'
prob = problem.Problem(pars)
problem.model.single_CA_activity = False
prob.fit_FB_vu = False
# Create arrays with parameter bounds, initial estimates [taken from model.py], whether to log-transform.
......@@ -53,6 +57,7 @@ if prob.fit_FB_vu:
if args.skipfit:
p = p_ini
prob.verbose = True
else:
# Run the optimization
import optimize
......@@ -65,6 +70,11 @@ print 'Optimization result:'
for par,value in zip(pars,p):
print ' %s = %s' % (par,value)
if args.save:
with open(args.save,'w') as f:
for par,value in zip(pars,p):
f.write('%s: %s\n' % (par,value))
# Evaluate model at the optimum to ensure all parameters are set to their optimal value.
fitness,results = prob.evaluateFitness(p,return_results=True)
......
......@@ -16,15 +16,20 @@ class Model(object):
#self.k_plus_x_CA = 2700 # s-1, Hopkinson et al. Table S3
#self.k_plus_v_CA = 0. # s-1, no CA
self.k_plus_CA = 4e7
self.k_plus_CA = 1e3
self.k_plus_i_CA = self.k_plus_CA
self.k_plus_x_CA = self.k_plus_CA
self.k_plus_v_CA = self.k_plus_CA
self.single_CA_activity = True
self.pH_i = 7.
self.pH_x = 7.9
self.pH_v = 7.1
self.T_a = 0.0007
self.T_b = 0.012
self.T_b = 0.03
# Spheroid parameters for chloroplast and vesicle
self.f_x = 1.1
self.f_v = 0.6
......@@ -39,8 +44,8 @@ class Model(object):
self.eps_cb_noCA = -13.0 # per mille, Zeebe Table 3.2.5, p 184, 24 degrees Celsius!
self.eps_bc_noCA = -22.0 # per mille, Zeebe Table 3.2.5, p 184, 24 degrees Celsius!
self.eps_cb4 = -11.
self.eps_bc4 = -20.
self.eps_cb4 = -11. # per mille, Zeebe Table 3.2.5, p 184, 24 degrees Celsius!
self.eps_bc4 = -20. # per mille, Zeebe Table 3.2.5, p 184, 24 degrees Celsius!
# To ensure difference is equal to equilibrium fractionation factor 10.3,
# which is used below to compute deltaC_e from deltaB_e.
......@@ -53,66 +58,106 @@ class Model(object):
for k,v in sorted(info.items(),cmp=lambda x,y: cmp(x[0],y[0])): print '%s = %s' % (k,v)
# Catch common unit errors
assert info['C_e']<1. and info['C_e']>0.001, 'Invalid C_e=%s. CO2 [aq] should lie between 0.001 and 1 mol m-3' % info['C_e']
assert info['DIC']>=0.5 and info['DIC']<30, 'Invalid DIC=%s. Expected value between 0.5 and 30 mol m-3' % info['DIC']
assert info['pH']>6. and info['pH']<9., 'Invalid pH=%s. Expected value between 6 and 9' % info['pH']
assert info['r']<0.0001, 'Invalid r=%s. radius in m should be less than 0.0001' % info['r']
self.k_plus_i_CA = self.k_plus_CA
self.k_plus_x_CA = self.k_plus_CA
self.k_plus_v_CA = self.k_plus_CA
if self.single_CA_activity:
self.k_plus_i_CA = self.k_plus_CA
self.k_plus_x_CA = self.k_plus_CA
self.k_plus_v_CA = self.k_plus_CA
self.k_ex = self.k_ei
self.k_ev = self.k_ei
# Shortcuts to commonly used arguments
F_FIX = info['F_FIX']*2
F_CAL = info['F_CAL']*2
C_e = info['C_e']
B_e = info['B_e']
F_FIX = info['F_FIX']*2 # NB factor 2 because only active during day
F_CAL = info['F_CAL']*2 # NB factor 2 because only active during day
r = info['r']
# H+ concentrations in mol L-1 from prescribed pH
H_e = 10.**-info['pH']
H_i = 10.**-self.pH_i
H_v = 10.**-self.pH_v
H_x = 10.**-self.pH_x
# Compute equilibrium constants accroding to Lueker et al 2000 [as per EPOCA recommendation]:
# K1_star = [HCO3-]*[H+]/[CO2] = k_plus/k_minus
# K2_star = [CO3--]*[H+]/[HCO3-]
S = 35.
TK = info['temp'] + 273.15
#pK1_star = 3670.7/TK - 62.008 + 9.7944*numpy.log(TK) - 0.0118*S + 0.000116*S*S # Millero p.664 (1995) using Mehrbach et al. data on seawater scale
pK1_star = 3633.86/TK - 61.2172 + 9.67770*numpy.log(TK) - 0.011555*S + 0.0001152*S*S # Lueker et al 2000 as per EPOCA recommendation
pK2_star = 471.78/TK + 25.9290 - 3.16967*numpy.log(TK) - 0.01781 *S + 0.0001122*S*S # Lueker et al 2000 as per EPOCA recommendation
K1_star = 10.**(-pK1_star)
K2_star = 10.**(-pK2_star)
if verbose: print 'K1_star = %s, K2_star = %s' % (K1_star,K2_star)
# Compute k+ using temperature-dependent relation of Zeebe book, p 110, table 2.3.1
k_plus_noCA = numpy.exp(1246.98 - 6.19e4/TK - 183.0*numpy.log(TK))
# Hydroxylation rate constant
R_igc = 8.3144621
k4_plus = 4.7e7 * numpy.exp(-23200/(R_igc*TK))
# Ionic product of water ([H+]*[OH-]) in M*M
# Source: EPOCA/OA recommendations Eq 1.10
K_W_star = numpy.exp(-13847.26/TK + 148.9652 - 23.6521*numpy.log(TK) + (118.67/TK - 5.977 + 1.0495*numpy.log(TK))*numpy.sqrt(S) - 0.01615*S)
OH_i = K_W_star/H_i
OH_x = K_W_star/H_x
OH_v = K_W_star/H_v
if verbose: print 'OH concentrations: cytosol = %s, chloroplast = %s, vesicle = %s' % (OH_i,OH_x,OH_v)
# Compute B_e and C_e from equilibrium constants, pH and total DIC.
CO3_per_HCO3_e = K2_star/H_e
HCO3_per_CO2_e = K1_star/H_e
B_e = info['DIC']/(1./HCO3_per_CO2_e + 1. + CO3_per_HCO3_e)
C_e = B_e/HCO3_per_CO2_e
if verbose: print 'B_e = %.4f mol m-3, C_e = %.4f mol m-3' % (B_e,C_e)
# Carbon fractionation factors from Zhang et al 1995
eps_CO2_g = 0.0049*info['temp'] - 1.31
eps_HCO3_g = -0.114 *info['temp'] + 10.78
eps_CO3_g = -0.052 *info['temp'] + 7.22
alpha_CO2_HCO3 = (1.+eps_CO2_g/1000)/(1.+eps_HCO3_g/1000)
alpha_CO3_HCO3 = (1.+eps_CO3_g/1000)/(1.+eps_HCO3_g/1000)
alpha_HCO3 = (1.+info['deltaC_DIC']/1000.)*info['DIC']/(alpha_CO2_HCO3*C_e + B_e + alpha_CO3_HCO3*(info['DIC']-C_e-B_e))
deltaB_e = (alpha_HCO3-1.)*1000
deltaC_e = (alpha_HCO3*alpha_CO2_HCO3-1.)*1000
if verbose: print 'deltaB_e = %.3f, deltaC_e = %.3f' % (deltaB_e,deltaC_e)
# Scale kinetic fractionation factors so that equilibrium fractionation
# is identical to reference value from Zhang et al 1995
eps_CO2_HCO3 = (alpha_CO2_HCO3-1)*1000
fac = eps_CO2_HCO3/(self.eps_bc_CA-self.eps_cb_CA)
self.eps_bc_CA *= fac
self.eps_cb_CA *= fac
fac = eps_CO2_HCO3/(self.eps_bc_noCA-self.eps_cb_noCA)
self.eps_cb_noCA *= fac
self.eps_bc_noCA *= fac
fac = eps_CO2_HCO3/(self.eps_bc4-self.eps_cb4)
self.eps_cb4 *= fac
self.eps_bc4 *= fac
SA = 4*numpy.pi*r**2 # cell surface area in m2, assuming spherical shape
V_cell = 4./3.*numpy.pi*r**3 # cell volume in m3, assuming spherical shape
if verbose: print 'division rate based on carbon estimated from volume = %.3f d-1' % (F_FIX/(20e3*V_cell)*86400,)
if verbose: print 'PIC:POC = %.3f' % (F_CAL/F_FIX,)
# Compute surface area and volume of organelles, and remaining volume of cytoplasm
SA_x,V_x = oblateSpheroid(self.f_x*r,self.f_x*r/self.S_x)
SA_v,V_v = oblateSpheroid(self.f_v*r,self.f_v*r/self.S_v)
V_i = V_cell - V_v - V_x # cytosol volume (m3)
if verbose: print 'relative SA of chloroplast = %.3f, vesicle = %.3f' % (SA_x/SA,SA_v/SA)
if verbose: print 'relative volume of chloroplast = %.3f, vesicle = %.3f, cytosol = %.3f' % (V_x/V_cell,V_v/V_cell,V_i/V_cell)
# Compute utilization without HCO3 permeability regulation, and from that compute final permeability
Ut = (F_FIX+F_CAL)/(self.k_ei*C_e*SA + self.T_a*self.k_ei*B_e*SA)
P_Bi = self.k_ei*(self.T_a + self.T_b*Ut)
P_Bx = P_Bi
P_Bv = P_Bi
if verbose: print 'membrane permeability for HCO3 relative to CO2 = %.4f (utilization = %.3f)' % (P_Bi/self.k_ei,Ut)
# Compute k+ using temperature-dependent relation of Zeebe book, p 110, table 2.3.1
TK = info['temp'] + 273.15
k_plus_noCA = numpy.exp(1246.98 - 6.19e4/TK - 183.0*numpy.log(TK))
# Hydroxylation rate constant
R_igc = 8.3144621
k4_plus = 4.7e7 * numpy.exp(-23200/(R_igc*TK))
# Compute K1_star = [HCO3-]*[H+]/[CO2] = k_plus/k_minus
S = 35.
#pK1_star = 3670.7/TK - 62.008 + 9.7944*numpy.log(TK) - 0.0118*S + 0.000116*S*S # Millero p.664 (1995) using Mehrbach et al. data on seawater scale
pK1_star = 3633.86/TK - 61.2172 + 9.67770*numpy.log(TK) - 0.011555*S + 0.0001152*S*S # Lueker et al 2000 as per EPOCA recommendation
K1_star = 10.**(-pK1_star)
if verbose: print 'K1_star = %s' % K1_star
# Ionic product of water ([H+]*[OH-]) in M*M
# Source: EPOCA/OA recommendations Eq 1.10
K_W_star = numpy.exp(-13847.26/TK + 148.9652 - 23.6521*numpy.log(TK) + (118.67/TK - 5.977 + 1.0495*numpy.log(TK))*numpy.sqrt(S) - 0.01615*S)
OH_i = K_W_star/H_i
OH_x = K_W_star/H_x
OH_v = K_W_star/H_v
# Compute CO2-specific hydration rate (combine hydration rate constant with and without CA, and hydroxilation rate)
k_specHyd_i = k_plus_noCA + self.k_plus_i_CA + k4_plus*OH_i
k_specHyd_x = k_plus_noCA + self.k_plus_x_CA + k4_plus*OH_x
......@@ -145,9 +190,9 @@ class Model(object):
A[5,:] = 0,P_Bv*SA_v,0,0,k_specHyd_v*V_v,-P_Bv*SA_v-k_specDeh_v*V_v
b[:] = -self.k_ei*SA*C_e,-P_Bi*SA*B_e,F_FIX,0,0,F_CAL
concs = numpy.linalg.solve(A,b)
if (concs<0).any(): return
C_i,B_i,C_x,B_x,C_v,B_v = concs
if verbose: print C_e,B_e,C_i,B_i,C_x,B_x,C_v,B_v
if verbose: print 'C_i = %.4f,B_i = %.4f,C_x = %.4f,B_x = %.4f,C_v = %.4f,B_v = %.4f' % (C_i,B_i,C_x,B_x,C_v,B_v)
if (concs<0).any(): return
if verbose: print 'Actual B_i/C_i = %.1f. Equilibrium value = %.1f' % (B_i/C_i,K1_star/H_i)
if verbose: print 'Actual B_x/C_x = %.1f. Equilibrium value = %.1f' % (B_x/C_x,K1_star/H_x)
if verbose: print 'Actual B_v/C_v = %.1f. Equilibrium value = %.1f' % (B_v/C_v,K1_star/H_v)
......@@ -203,8 +248,8 @@ class Model(object):
A[3,:] = (0.0, FB_xu, FHyd_x, -FB_xo-FDeh_x, 0.0, 0.0)
A[4,:] = (FC_vu, 0.0, 0.0, 0.0, -FC_vo-FHyd_v, FDeh_v)
A[5,:] = (0.0, FB_vu, 0.0, 0.0, FHyd_v, -FB_vo-FDeh_v-F_CAL)
b[:] = (-eps_bc_i*FDeh_i + eps_cb_i*FHyd_i - info['deltaC_e']*FC_iu,
eps_bc_i*FDeh_i - eps_cb_i*FHyd_i - info['deltaB_e']*FB_iu,
b[:] = (-eps_bc_i*FDeh_i + eps_cb_i*FHyd_i - deltaC_e*FC_iu,
eps_bc_i*FDeh_i - eps_cb_i*FHyd_i - deltaB_e*FB_iu,
-eps_bc_x*FDeh_x + eps_cb_x*FHyd_x + self.eps_FIX*F_FIX,
eps_bc_x*FDeh_x - eps_cb_x*FHyd_x,
-eps_bc_v*FDeh_v + eps_cb_v*FHyd_v,
......@@ -219,7 +264,7 @@ class Model(object):
# Compute delta13C of POC and PIC using isotopic composition of relevant pools and fractionation factors.
delta13C_org_mod = deltaC_x + self.eps_FIX
delta13C_cal_mod = deltaB_v + self.eps_CAL
if verbose: print delta13C_org_mod-info['deltaB_e']
if verbose: print delta13C_org_mod-deltaB_e
result = info.copy()
result.update({'FB_iu':FB_iu,'FB_xu':FB_xu,'FB_vu':FB_vu,
......
import argparse
import numpy
import model
parser = argparse.ArgumentParser(description='Plot figures to analyze model behaviour.')
parser.add_argument('--load', help='Path to load parameter estimates from.')
args = parser.parse_args()
m = model.Model()
# Model parameters (usually taken from best fit obtained earlier)
#m.k_ei = 0.004717235859
#m.T_a = 0.000814610301367
#m.T_b = 0.116871311168
#m.eps_FIX = -10.6106670595
#m.k_plus_CA = 2583.94260624
if args.load:
with open(args.load,'rU') as f:
for l in f:
par,value = l.split(': ')
print 'Setting %s to %s.' % (par,value)
setattr(m,par,float(value))
# Range of cell radius to try (micrometer)
rs = 1.,3.,9.
# Constant environment
temp = 15.
delta13C_HCO3 = 0. #0. #-3.5
delta13C_DIC = 0. #0. #-3.5
pH = 8.2
S = 35.
# Set division rate in per day.
div = 1.
eps_CO2_g = 0.0049*temp - 1.31 # From Zhang et al 1995
eps_HCO3_g = -0.114 *temp + 10.78 # From Zhang et al 1995
delta13C_CO2 = delta13C_HCO3 + eps_CO2_g - eps_HCO3_g
# Compute cytosolic K1_star = [HCO3-]*[H+]/[CO2] = k_plus_i/k_minus_i
TK = temp + 273.15
#pK1_star = 3670.7/TK - 62.008 + 9.7944*numpy.log(TK) - 0.0118*S + 0.000116*S*S # Millero p.664 (1995) using Mehrbach et al. data on seawater scale
......@@ -60,22 +62,19 @@ for ir,r in enumerate(rs):
PIC_POCs = numpy.linspace(0.,2.5,100)
# pCO2 range to evaluate
pCO2s = numpy.linspace(100.,1500.,101)
DICs = numpy.linspace(0.5,25.,101)
#pCO2s = 367.,
delta13C_cals = numpy.ma.array(numpy.empty((pCO2s.size,PIC_POCs.size)),mask=True)
delta13C_orgs = numpy.ma.array(numpy.empty((pCO2s.size,PIC_POCs.size)),mask=True)
for i,pCO2 in enumerate(pCO2s):
delta13C_cals = numpy.ma.array(numpy.empty((DICs.size,PIC_POCs.size)),mask=True)
delta13C_orgs = numpy.ma.array(numpy.empty((DICs.size,PIC_POCs.size)),mask=True)
for i,DIC in enumerate(DICs):
for j,PIC_POC in enumerate(PIC_POCs):
# Compute fixation and calcification rates (mol/s)
F_FIX = C_org*div/86400.
F_CAL = F_FIX*PIC_POC
# Have model estimate bicarbonate fluxes and resulting delta13C values.
# Note pCO2 must be provided as fraction, i.e., ppm/1e6
C_e = 37.*pCO2*1e-6 # NB 37. from Zeebe book p 67 Fig 1.5.20, verify 37*367*1e-6 = 0.0136 mol m-3 = 13.6 uM
B_e = B_e_per_C_e*C_e
result = m.getBicarbonateFluxes(r=r,F_FIX=F_FIX,F_CAL=F_CAL,div=div,C_e=C_e,B_e=B_e,temp=temp,deltaB_e=delta13C_HCO3,deltaC_e=delta13C_CO2)
result = m.getBicarbonateFluxes(r=r,F_FIX=F_FIX,F_CAL=F_CAL,div=div,DIC=DIC,temp=temp,deltaC_DIC=delta13C_DIC,pH=pH)
if result is not None:
delta13C_cals[i,j] = result['delta13C_cal']
delta13C_orgs[i,j] = result['delta13C_org']
......@@ -83,15 +82,15 @@ for ir,r in enumerate(rs):
pylab.figure(fig.number)
pylab.subplot(len(rs),2,ir*2+1)
cs = numpy.arange(numpy.floor(delta13C_cals.min()),numpy.ceil(delta13C_cals.max())+0.1,1.)
pc=pylab.contourf(PIC_POCs,pCO2s,delta13C_cals,cs,cmap=cmap)
pylab.contour(PIC_POCs,pCO2s,delta13C_cals,cs,colors='k',linewidths=.5,linestyles='-')
pc=pylab.contourf(PIC_POCs,DICs,delta13C_cals,cs,cmap=cmap)
pylab.contour(PIC_POCs,DICs,delta13C_cals,cs,colors='k',linewidths=.5,linestyles='-')
pylab.colorbar(pc)
if ir+1==len(rs): pylab.xlabel('PIC:POC (-)')
pylab.ylabel('pCO2 (ppm)')
pylab.ylabel('DIC (mol m-3)')
pylab.subplot(len(rs),2,ir*2+2)
cs = numpy.arange(numpy.floor(delta13C_orgs.min()),numpy.ceil(delta13C_orgs.max())+0.1,1.)
pc=pylab.contourf(PIC_POCs,pCO2s,delta13C_orgs,cs,cmap=cmap)
pylab.contour(PIC_POCs,pCO2s,delta13C_orgs,cs,colors='k',linewidths=.5,linestyles='-')
pc=pylab.contourf(PIC_POCs,DICs,delta13C_orgs,cs,cmap=cmap)
pylab.contour(PIC_POCs,DICs,delta13C_orgs,cs,colors='k',linewidths=.5,linestyles='-')
pylab.colorbar(pc)
if ir+1==len(rs): pylab.xlabel('PIC:POC (-)')
......
......@@ -10,7 +10,8 @@ import numpy
import datafile
import model
df = datafile.DataFile()
#df = datafile.DataFile()
df = datafile.DataFile('data20150430_new.txt')
# Filter dataset by requiring that the following columns do not have any missing data.
df.requireColumn('C_org_assim_rate')
......@@ -22,7 +23,8 @@ df.requireColumn('pCO2')
df.requireColumn('Strain')
df.requireColumn('Mol_Corg_per_cell')
df.requireColumn('d13CDIC')
#df.exclude('Strain','CBR') # Pelagicus
df.exclude('Strain','CBR') # Ros Pelagicus
df.exclude('Strain','GOC') # Ros G.oceanica
#df.exclude('Strain','1216') # Ehux
#df.exclude('Strain','1256') # Ehux
#df.exclude('Strain','1211') # G.oceanica
......@@ -38,20 +40,25 @@ rs = df.grabColumn('Rc')*1e-6 # from um to m
pCO2s = df.grabColumn('pCO2')*1e-6 # from ppm to fraction
strains = df.grabColumn('Strain',dtype=str)
Cpercell = df.grabColumn('Mol_Corg_per_cell')
d13C_HCO3s = df.grabColumn('dHCO3')
d13CDICs = df.grabColumn('d13CDIC')
d13C_CO2s = df.grabColumn('dCO2')
CO2s = df.grabColumn('CO2')*1e3 # from mol/L to mol/m3
HCO3s = df.grabColumn('HCO3')*1e3 # from mol/L to mol/m3
DICs = df.grabColumn('DICmM')
pHs = df.grabColumn('pH')
TEMPs = df.grabColumn('TEMP')
divs = df.grabColumn('Division_rate')
unique_strains = tuple(set(strains))
print 'Remaining strains:'
for strain in unique_strains: print ' %s' % strain
var_org = 0.2008
var_cal = 0.0262
class Problem(optimize.OptimizationProblem):
def __init__(self,parameternames):
self.model = model.Model()
self.parameternames = parameternames
self.fit_FB_vu = False
self.verbose = False
def evaluateFitness(self,parameters,return_results=False):
if self.fit_FB_vu:
......@@ -74,8 +81,8 @@ class Problem(optimize.OptimizationProblem):
results = []
for i in range(len(rs)):
# Let model estimate bicarbonate fluxes and compute resulting isotopic composition.
result = self.model.getBicarbonateFluxes(r=rs[i],F_FIX=F_FIXs[i],F_CAL=F_CALs[i],C_e=CO2s[i],deltaB_e=d13C_HCO3s[i],deltaC_e=d13C_CO2s[i],temp=TEMPs[i],div=divs[i],
deltaC_org=d13Corgs[i],deltaC_cal=d13Ccals[i],fixed_FB_vu=FB_vus[i],B_e=HCO3s[i],verbose=True)
result = self.model.getBicarbonateFluxes(r=rs[i],F_FIX=F_FIXs[i],F_CAL=F_CALs[i],DIC=DICs[i],temp=TEMPs[i],div=divs[i],
deltaC_DIC=d13CDICs[i], deltaC_org=d13Corgs[i],deltaC_cal=d13Ccals[i],fixed_FB_vu=FB_vus[i],pH=pHs[i],verbose=self.verbose)
# If the parameter set was invalid, return maximum penalty (fitness = -infinity)
if result is None: return -numpy.Inf
......@@ -83,13 +90,13 @@ class Problem(optimize.OptimizationProblem):
# Increment sum of squares for delta13C misfit of POC and PIC.
ssq_org += (d13Corgs[i]-result['delta13C_org'])**2
ssq_cal += (d13Ccals[i]-result['delta13C_cal'])**2
print 'RESULT: ',d13Corgs[i]-d13CDICs[i],result['delta13C_org']-d13CDICs[i],d13Ccals[i]-d13CDICs[i],result['delta13C_cal']-d13CDICs[i]
#print 'RESULT: ',d13Corgs[i]-d13CDICs[i],result['delta13C_org']-d13CDICs[i],d13Ccals[i]-d13CDICs[i],result['delta13C_cal']-d13CDICs[i]
results.append(result)
# Compute total sum of squares, combining misfit of delta13C_org and delta13C_cal.
# It would be conceivable to weigh these according to their expected variance.
ssq = ssq_org + ssq_cal
ssq = ssq_org/var_org + ssq_cal/var_cal
print parameters,numpy.sqrt(ssq_org/rs.size),numpy.sqrt(ssq_cal/rs.size),ssq
# Fitness = likelihood = proportional to -SSQ for normally distributed errors
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment