Skip to content

Instantly share code, notes, and snippets.

@barronh
Last active January 13, 2021 20:13
Show Gist options
  • Save barronh/c6311e69d457c31013db42e8136c723c to your computer and use it in GitHub Desktop.
Save barronh/c6311e69d457c31013db42e8136c723c to your computer and use it in GitHub Desktop.
Patch necessary to update BETR-Global from Python2 to Python3
This patch can be used to update the folder dowloaded from
https://sourceforge.net/projects/betrs/files/BETR-Research-1.0.tar.gz/download.
The original file is Python2 compliant. After the patch, the tutorial works on
Google Colab using Python3.
diff -ru BETR-Global/BETRS.py BETR-Global3/BETRS.py
--- BETR-Global/BETRS.py 2010-09-28 03:03:24.000000000 +0000
+++ BETR-Global3/BETRS.py 2021-01-13 18:33:01.000000000 +0000
@@ -23,8 +23,9 @@
from numpy import *
import sys
import pdb
-import cPickle
+import _pickle as cPickle
import os
+import imp
## these imports reload modules in each run, for making development
@@ -49,19 +50,19 @@
import solver
import output
import modifyparams
-reload(readinput)
-reload(globalz)
-reload(zvalues)
-reload(tempcorrect)
-reload(volumes)
-reload(processes)
-reload(mkflowD)
-reload(mkbigD)
-reload(emissions)
-reload(helpers)
-reload(solver)
-reload(output)
-reload(modifyparams)
+imp.reload(readinput)
+imp.reload(globalz)
+imp.reload(zvalues)
+imp.reload(tempcorrect)
+imp.reload(volumes)
+imp.reload(processes)
+imp.reload(mkflowD)
+imp.reload(mkbigD)
+imp.reload(emissions)
+imp.reload(helpers)
+imp.reload(solver)
+imp.reload(output)
+imp.reload(modifyparams)
from readinput import *
from globalz import *
from zvalues import *
@@ -169,7 +170,7 @@
def update_compartments(self, compfile):
self.compdict = readCompartments(compfile)
def update_flows(self, flowdir, mean=False):
- self.flowdict = readFlows(self.compdict.keys(),flowdir)
+ self.flowdict = readFlows(list(self.compdict.keys()),flowdir)
def update_processlist(self, processfile):
self.proclist = readProcesses(self.compdict, processfile)
def update_control(self, controlfile):
@@ -230,7 +231,7 @@
#Calculating system-matrix indices indicating
#zero-volume compartments; setting self.killidx
killidx=[]
- for c in self.compdict.keys():
+ for c in list(self.compdict.keys()):
zerovols=where(sum(self.vdict[c]['bulk'], axis=1)==0)[0]
killidx.extend(tocell(zerovols+1,c,self))
self.killidx=sort(array(killidx))
@@ -260,14 +261,14 @@
os.mkdir(os.path.dirname(fn))
else:
if os.path.exists(fn):
- print('Attention: overwriting %s\n') % (fn)
+ print(('Attention: overwriting %s\n') % (fn))
of=open(fn,'w')
cPickle.dump(self.bigDlist,of,protocol=2)
of.close()
- print('wrote D-value - matrices to %s\n') % (fn)
+ print(('wrote D-value - matrices to %s\n') % (fn))
## consistency check for killidx
- print('Consistency check of zero-volume compartments: '),
+ print(('Consistency check of zero-volume compartments: '), end=' ')
kidx=self.mkkillidx_fromD()
if not all(kidx == self.killidx):
sys.exit('Detected inconsistencies between\
@@ -314,13 +315,18 @@
Averages the all mass rate matrices for each season (month) to
get a time invariant system. Solves for the steady-state.
'''
+ from scipy.sparse import csr_matrix
## average D-value matrices
- self.Davg = mean(self.bigDlist)
+ dls = 0
+ for dl in self.bigDlist:
+ dls = dls + dl
+ dlm = dls / len(self.bigDlist)
+ self.Davg = dlm
try:
q=self.emission.get_emission(0,self,type='array')
except AttributeError:
- print("update_emissions(<filename>) has to be called before\n"\
- +"solving for steady-state. Did nothing!\n")
+ print(("update_emissions(<filename>) has to be called before\n"\
+ +"solving for steady-state. Did nothing!\n"))
return()
ss_res =dot(-inv(self.Davg.toarray()), q)
self.ss_res=expandvec(ss_res, self.killidx)
@@ -374,8 +380,8 @@
try:
avgmatrix=expandmatrix(self.Davg, self.killidx).toarray()
except AttributeError:
- print("solve_ss() has to be called before acessing the averaged"\
- +" system matrix.")
+ print(("solve_ss() has to be called before acessing the averaged"\
+ +" system matrix."))
return()
print("OK")
return(avgmatrix)
@@ -388,8 +394,8 @@
try:
ssem=self.emission.get_emission(0,self,'array')
except AttributeError:
- print("update_emissions(<filename>) has to be called before"\
- +" acessing the emission vector. Did nothing!")
+ print(("update_emissions(<filename>) has to be called before"\
+ +" acessing the emission vector. Did nothing!"))
return()
print("Expanding emission vector ...")
ssem=expandvec(ssem, self.killidx)
@@ -403,7 +409,7 @@
self.parmean[field]=mean(self.par[field], axis=1)
def average_flows(self):
self.flowdictmean={}
- for item in self.flowdict.items():
+ for item in list(self.flowdict.items()):
meanmat=c_[item[1][:,0:2], mean(item[1][:,2:], axis=1)]
self.flowdictmean[item[0]]=meanmat
diff -ru BETR-Global/emissions.py BETR-Global3/emissions.py
--- BETR-Global/emissions.py 2010-09-27 12:13:32.000000000 +0000
+++ BETR-Global3/emissions.py 2021-01-13 17:35:59.000000000 +0000
@@ -21,7 +21,7 @@
try:
f=open(fn, 'r')
except IOError:
- print("emissions.py: emission file %s not found. Aborting!") % (fn)
+ print(("emissions.py: emission file %s not found. Aborting!") % (fn))
sys.exit(1)
lines=f.readlines()
f.close
diff -ru BETR-Global/Flows/default/flow55.txt BETR-Global3/Flows/default/flow55.txt
--- BETR-Global/Flows/default/flow55.txt 2010-09-24 02:18:36.000000000 +0000
+++ BETR-Global3/Flows/default/flow55.txt 2021-01-13 19:15:00.000000000 +0000
@@ -801,7 +801,7 @@
15 15 123000000 123000000 123000000 123000000 123000000 123000000 123000000 123000000 123000000 123000000 123000000 123000000
16 16 121000000 121000000 121000000 121000000 121000000 121000000 121000000 121000000 121000000 121000000 121000000 121000000
17 17 120000000 120000000 120000000 120000000 120000000 120000000 120000000 120000000 120000000 120000000 120000000 120000000
-18 18 130000000 130000000 130000000 130000000 130000000 130000000 130000000 130000000 130000000 130000000 130000000 130000000 `
+18 18 130000000 130000000 130000000 130000000 130000000 130000000 130000000 130000000 130000000 130000000 130000000 130000000
19 19 120000000 120000000 120000000 120000000 120000000 120000000 120000000 120000000 120000000 120000000 120000000 120000000
20 20 147000000 147000000 147000000 147000000 147000000 147000000 147000000 147000000 147000000 147000000 147000000 147000000
21 21 166000000 166000000 166000000 166000000 166000000 166000000 166000000 166000000 166000000 166000000 166000000 166000000
Only in BETR-Global/Flows/default: .svn
Only in BETR-Global/Flows: .svn
diff -ru BETR-Global/mkbigD.py BETR-Global3/mkbigD.py
--- BETR-Global/mkbigD.py 2010-09-24 02:18:54.000000000 +0000
+++ BETR-Global3/mkbigD.py 2021-01-13 17:35:59.000000000 +0000
@@ -27,16 +27,16 @@
'''takes a model object (m) as argument and returns a list of sparse
matrices which contain D-values from inter-cell flows for each timestep.'''
fdict=copy.deepcopy(m.Dflow)
- for tp in fdict.keys(): # check whether we have all compartments
+ for tp in list(fdict.keys()): # check whether we have all compartments
for c in tp:
- if c not in m.compdict.keys():
- print('''flows.py: Compartment %i not in compartment-list.
- Aborting !''') % (c)
+ if c not in list(m.compdict.keys()):
+ print(('''flows.py: Compartment %i not in compartment-list.
+ Aborting !''') % (c))
sys.exit(1)
# remove intra-cell flows
# this is for example oceanic sinking flux, which is dealt with
# by "betr_ocean_sinkflux" in "processes.py".
- for [k,v] in fdict.items():
+ for [k,v] in list(fdict.items()):
intercell=list(where(v[:,0]!=v[:,1])[0])
intracell=list(where(v[:,0]==v[:,1])[0])
fdict[k]=v[intercell,:]
@@ -45,7 +45,7 @@
for t in arange(0,m.nots):
matlist.append(sp.coo_matrix(matrix(zeros((m.matdim,m.matdim),
dtype=float64))))
- for [k,v] in fdict.items():
+ for [k,v] in list(fdict.items()):
if len(v)==0: continue
toidx=array([(x-1)*m.nocomp+k[1]-1 for x in v[:,1]])
fromidx=array([(x-1)*m.nocomp+k[0]-1 for x in v[:,0]])
@@ -60,19 +60,19 @@
def mkprocmat(m):
''' returns a list of sparse matrices, one for each timestep, that
contain D-values for intra-cell processes'''
- for tp in m.Dproc.keys(): # check whether we have all compartments
+ for tp in list(m.Dproc.keys()): # check whether we have all compartments
tp=(tp[0],tp[1])
for c in tp:
- if c not in m.compdict.keys():
- print('''mkprocmat.py: Compartment %i not in compartment-list.
- Aborting !''') % (c)
+ if c not in list(m.compdict.keys()):
+ print(('''mkprocmat.py: Compartment %i not in compartment-list.
+ Aborting !''') % (c))
sys.exit(1)
# construct list of large sparse matrices
matlist=[]
for t in arange(0,m.nots):
matlist.append(sp.coo_matrix(matrix(zeros((m.matdim,m.matdim),
dtype=float64))))
- for [k,v] in m.Dproc.items():
+ for [k,v] in list(m.Dproc.items()):
toidx=array([x*m.nocomp+k[1]-1 for x in range(0,m.nocells)])
fromidx=array([x*m.nocomp+k[0]-1 for x in range(0,m.nocells)])
ij=(toidx.astype(int),fromidx.astype(int))
@@ -118,7 +118,7 @@
# construct list of large sparse matrices
mlist=[]
diags=zeros((m.nots,m.matdim))
- for c in m.compdict.keys():
+ for c in list(m.compdict.keys()):
idx=tocell(arange(1, m.nocells+1),c,m)
d=m.zdict[c]['bulk']*m.vdict[c]['bulk']
for t in arange(0,m.nots):
diff -ru BETR-Global/mkflowD.py BETR-Global3/mkflowD.py
--- BETR-Global/mkflowD.py 2010-09-24 02:18:54.000000000 +0000
+++ BETR-Global3/mkflowD.py 2021-01-13 17:35:59.000000000 +0000
@@ -14,7 +14,7 @@
def mkflowD(model):
fdict=copy.deepcopy(model.flowdict)
- for f in fdict.keys():
+ for f in list(fdict.keys()):
fromcells=fdict[f][:,0].astype(int)
zvals=model.zdict[f[0]]['bulk'][fromcells-1,:]
fdict[f][:,2:]=fdict[f][:,2:]*zvals
Only in BETR-Global3/Output/default: ss_kg.nc
Only in BETR-Global3/Output/default: ss_mol_per_m3.nc
Only in BETR-Global3/Output/default: ss_out.cpk
diff -ru BETR-Global/output.py BETR-Global3/output.py
--- BETR-Global/output.py 2010-09-27 22:40:35.000000000 +0000
+++ BETR-Global3/output.py 2021-01-13 20:04:16.983520770 +0000
@@ -8,17 +8,18 @@
from numpy import *
import sys
import os
-import cPickle
+import pickle
import write_ncfile
-reload(write_ncfile)
+import imp
+imp.reload(write_ncfile)
from write_ncfile import *
def write_output_ss(m, fn, units, netcdf):
''' Output for steady-state results'''
out={}
- for c in m.compdict.keys():
+ for c in list(m.compdict.keys()):
out[c]={}
- idx=arange(c-1,m.matdim,len(m.compdict.keys()))
+ idx=arange(c-1,m.matdim,len(list(m.compdict.keys())))
## deal with zero-volumes / Z-values
varr=mean(m.vdict[c]['bulk'],axis=1)
zv_arr=mean(m.zdict[c]['bulk'],axis=1)*varr
@@ -46,17 +47,18 @@
os.mkdir(os.path.dirname(fncpk))
else:
if os.path.exists(fncpk):
- print('Attention: overwriting %s\n') % (fncpk)
- f=open(fncpk,'w')
- cPickle.dump(out,f)
+ print(('Attention: overwriting %s\n') % (fncpk))
+ f=open(fncpk,'wb')
+
+ pickle.dump(out,f)
f.close()
if netcdf:
for u in units:
fnnc=fn+'_'+u+'.nc'
if os.path.exists(fnnc):
- print('Attention: overwriting %s\n') % (fnnc)
+ print(('Attention: overwriting %s\n') % (fnnc))
outarray=zeros((m.nocomp,12,24))
- for c in m.compdict.keys():
+ for c in list(m.compdict.keys()):
outarray[c-1,:,:]=out[c][u].reshape(12,24)
writenc(outarray,fnnc,False,True,1,varname='V',unit=u)
@@ -67,9 +69,9 @@
if periods != (timesteps-1)/float(m.nots):
sys.exit('timesteps of output {0:d} not multiple of '
+'seasons ({1:d}): Aborting !'.format(timesteps,m.nots))
- for c in m.compdict.keys():
+ for c in list(m.compdict.keys()):
out[c]={}
- idx=arange(c-1,m.matdim,len(m.compdict.keys()))
+ idx=arange(c-1,m.matdim,len(list(m.compdict.keys())))
## deal with zero volumes / Z-values
varr=hstack((ones((idx.shape[0],1)),
tile(m.vdict[c]['bulk'], (1,periods))))
@@ -102,17 +104,17 @@
os.mkdir(os.path.dirname(fncpk))
else:
if os.path.exists(fncpk):
- print('Attention: overwriting %s\n') % (fncpk)
- f=open(fncpk,'w')
- cPickle.dump(out,f)
+ print(('Attention: overwriting %s\n') % (fncpk))
+ f=open(fncpk,'wb')
+ pickle.dump(out,f)
f.close()
if netcdf:
for u in units:
fnnc=fn+'_'+u+'.nc'
if os.path.exists(fnnc):
- print('Attention: overwriting %s\n') % (fnnc)
+ print(('Attention: overwriting %s\n') % (fnnc))
outarray=zeros((timesteps,m.nocomp,12,24))
- for c in m.compdict.keys():
+ for c in list(m.compdict.keys()):
for t in arange(0,timesteps):
outarray[t,c-1,:,:]=out[c][u][:,t].reshape(12,24)
writenc(outarray, fnnc, True, True, 1, varname='V',unit=u)
Only in BETR-Global/PY: modifyparams.pyc
diff -ru BETR-Global/PY/processes.py BETR-Global3/PY/processes.py
--- BETR-Global/PY/processes.py 2010-09-24 02:18:35.000000000 +0000
+++ BETR-Global3/PY/processes.py 2021-01-13 17:35:59.000000000 +0000
@@ -38,8 +38,8 @@
try:
getattr(self,p)
except AttributeError:
- print("processes.py: "
- +"Method %s not implemented !\n Aborting!") % (p)
+ print(("processes.py: "
+ +"Method %s not implemented !\n Aborting!") % (p))
sys.exit(1)
def getD(self):
@@ -54,7 +54,7 @@
# don't change this
D={}
procname=inspect.getframeinfo(inspect.currentframe())[2]
- for c in self.m.compdict.keys():
+ for c in list(self.m.compdict.keys()):
D[(c,c,procname)]=self.m.chempardict[c]['k_reac']\
*self.m.vdict[c]['bulk']\
*self.m.zdict[c]['bulk']
Only in BETR-Global/PY: processes.pyc
Only in BETR-Global3/PY: __pycache__
Only in BETR-Global/PY: .svn
diff -ru BETR-Global/PY/volumes.py BETR-Global3/PY/volumes.py
--- BETR-Global/PY/volumes.py 2010-09-24 02:18:35.000000000 +0000
+++ BETR-Global3/PY/volumes.py 2021-01-13 17:35:59.000000000 +0000
@@ -27,8 +27,8 @@
empty=[]
for i in arange(0,len(compdict)):
empty.append({})
- self.V=dict(zip(self.compdict.keys(), empty))
- for c in self.compdict.keys():
+ self.V=dict(list(zip(list(self.compdict.keys()), empty)))
+ for c in list(self.compdict.keys()):
try:
self.V[c]=getattr(self,'V'+str(c))()
except AttributeError:
Only in BETR-Global/PY: volumes.pyc
diff -ru BETR-Global/PY/zvalues.py BETR-Global3/PY/zvalues.py
--- BETR-Global/PY/zvalues.py 2010-09-24 02:18:35.000000000 +0000
+++ BETR-Global3/PY/zvalues.py 2021-01-13 17:36:00.000000000 +0000
@@ -30,8 +30,8 @@
empty=[]
for i in arange(0,len(compdict)):
empty.append({})
- self.Z=dict(zip(self.compdict.keys(), empty))
- for c in self.compdict.keys():
+ self.Z=dict(list(zip(list(self.compdict.keys()), empty)))
+ for c in list(self.compdict.keys()):
try:
self.Z[c]=getattr(self,'Z'+str(c))()
except AttributeError:
Only in BETR-Global/PY: zvalues.pyc
Only in BETR-Global3: __pycache__
diff -ru BETR-Global/readinput.py BETR-Global3/readinput.py
--- BETR-Global/readinput.py 2010-09-24 02:18:54.000000000 +0000
+++ BETR-Global3/readinput.py 2021-01-13 17:49:15.000000000 +0000
@@ -13,7 +13,7 @@
import csv
from globalz import *
import sys
-
+string = str
def _readEnvironment(partype, fn):
""" Reads a file containing environmental variables.
<partype> : 'const' or 'seasonal'. Determines whether the input file
@@ -57,7 +57,7 @@
## 1) increasing cell number and 2) increasing timestep
pars=sort(pars, axis=0, order=['CELL','TS'])
fields=list(pars.dtype.names)
- fields=filter(lambda x: x!='CELL' and x!='TS', fields)
+ fields=[x for x in fields if x!='CELL' and x!='TS']
pars=pars[fields]
pars=reshape(pars,(cells, ts))
return(pars)
@@ -101,10 +101,10 @@
line=f.readline()
if line == '': break
if line[0] =='#' or line == '\n': continue
- line=string.split(line)
+ line=str.split(line)
num=line[0]
del line[0]
- compdict[int(num)] = dict(zip(header,line))
+ compdict[int(num)] = dict(list(zip(header,line)))
return(compdict)
def readFlows(complist, flowdir):
@@ -128,6 +128,7 @@
if not (fro in complist and to in complist):
continue
else:
+ print(fn)
fmatrixdict[(fro,to)]=loadtxt(fn, comments='#')
return(fmatrixdict)
@@ -145,7 +146,7 @@
line=f.readline()
c=line[0]
count+=1
- header=string.split(line_old)
+ header=str.split(line_old)
header[0]=header[0].lstrip("#")
f.seek(0)
# jump over comments
@@ -175,7 +176,7 @@
if line[0] == '#' or line == '\n': continue
line=line.split()
comps=[int(x) for x in line[1:]]
- if all([x in compdict.keys() for x in comps]):
+ if all([x in list(compdict.keys()) for x in comps]):
pl=[line[0]]
pl.extend(comps)
proclist.append(tuple(pl))
diff -ru BETR-Global/solveBETR.py BETR-Global3/solveBETR.py
--- BETR-Global/solveBETR.py 2010-09-24 02:18:55.000000000 +0000
+++ BETR-Global3/solveBETR.py 2021-01-13 17:35:59.000000000 +0000
@@ -18,10 +18,10 @@
sys.exit(mmmdir + ' doesn not exist; aborting !')
print('Initialization of Model')
MMMdata.mmmdir = mmmdir
- print('\tModel directory: '+MMMdata.mmmdir)
+ print(('\tModel directory: '+MMMdata.mmmdir))
MMMdata.area,\
MMMdata.volume = self._readSize()
- print('\tRead '+str(MMMdata.area.nnz)+' compartment dimensions.')
+ print(('\tRead '+str(MMMdata.area.nnz)+' compartment dimensions.'))
MMMdata.stepdef,\
MMMdata.period,\
MMMdata.steps_per_period = self._readSteps(run)
@@ -30,10 +30,10 @@
## print('\t'+str(MMMdata.stepdef))
MMMdata.grid,\
MMMdata.comps = self._readDims()
- print('\tThe grid has '+str(len(MMMdata.grid[0]))+' x '+str(len(MMMdata.grid[1]))+\
- ' regions')
- print('\tEach region has '+str(len(MMMdata.comps))+' compartments:')
- print('\t'+str(MMMdata.comps))
+ print(('\tThe grid has '+str(len(MMMdata.grid[0]))+' x '+str(len(MMMdata.grid[1]))+\
+ ' regions'))
+ print(('\tEach region has '+str(len(MMMdata.comps))+' compartments:'))
+ print(('\t'+str(MMMdata.comps)))
MMMdata.cells = len(MMMdata.comps) * len(MMMdata.grid[0]) * len(MMMdata.grid[1])
MMMdata.matrices, MMMdata.killidx = self._getMatrices(ss)
## print('\t'+str(len(MMMdata.matrices))+\
@@ -66,7 +66,7 @@
Area = csr_matrix((array(Adata),array(ij)))
Vol = csr_matrix((array(Vdata),array(ij)))
- print(dir(Area))
+ print((dir(Area)))
if (Area.nnz != Vol.nnz):
sys.exit('Inconsistent Area / Volume information in SIZE.txt; aborting !')
return (Area, Vol)
@@ -150,7 +150,7 @@
A = zeros((cells,cells))
try:
f=open(filename, 'r')
- print("\treading file: %s") %(filename)
+ print(("\treading file: %s") %(filename))
except IOError:
sys.exit('Could not read ' + filename + ' ; aborting !')
lines=f.readlines()
@@ -214,7 +214,7 @@
[m,r,c,val]=[x.lstrip() for x in [x.rstrip() for x in lin.split(',')]]
self.em.setdefault(int(float(m)),[]).append((int(float(r)), int(float(c)) ,float(val)))
print('Initialization of emissions')
- print('\tRead emission data for '+str(len(self.em))+' months')
+ print(('\tRead emission data for '+str(len(self.em))+' months'))
print('--- end of initialization')
def get(self, monthidx, killidx):
@@ -322,7 +322,7 @@
meanblock = divz(meanblock, newVolBlock)
Acoarse[newReg[0]*7:(newReg[0]+1)*7,newReg[0]*7:(newReg[0]+1)*7]=meanblock
# average inter - regional flows
- for border in B.keys():
+ for border in list(B.keys()):
meanblock = zeros((7,7))
newVolBlock = tile(Vnew[border[0]-1,:], (7,1))
for fromReg, toReg in B[border]:
@@ -400,7 +400,7 @@
ts=time.time()
for m in range(0,self.no_months):
self.current_month = m
- print('Integrating month # %s') % (m)
+ print(('Integrating month # %s') % (m))
ig.set_initial_value(ig.y, ig.t)
self.current_season = mod(m, self.model.steps_per_period)
self.current_matrix_csr = self.model.matrices[self.current_season]
@@ -408,17 +408,17 @@
no_sub_steps = int(ceil(self.model.stepdef[self.current_season] / \
self.maxsubsteplength))
substeplength = self.model.stepdef[self.current_season] / no_sub_steps
- print('\tSeason %d\t %d substeps of length %.5f hours')\
- % (self.current_season,no_sub_steps, substeplength)
+ print(('\tSeason %d\t %d substeps of length %.5f hours')\
+ % (self.current_season,no_sub_steps, substeplength))
for i in range(0, no_sub_steps):
- print('\tIntegrating substep %d') % (i)
+ print(('\tIntegrating substep %d') % (i))
ig.integrate(ig.t + substeplength)
if (not ig.successful()):
sys.exit("Integration failed; aborting !")
- print('done with month %d') % (m)
+ print(('done with month %d') % (m))
res=hstack((res, ig.y))
te=time.time()
- print "time elapsed: %f" %(te-ts)
+ print("time elapsed: %f" %(te-ts))
return res
@@ -429,14 +429,14 @@
model.stepdef,\
model.period,\
model.steps_per_period = model._readSteps(emobj.run)
- print('\tDuration of a period (hrs): '+str(model.period))
- print('\t'+str(model.steps_per_period)+' steps per period:')
- print('\t'+str(model.stepdef))
+ print(('\tDuration of a period (hrs): '+str(model.period)))
+ print(('\t'+str(model.steps_per_period)+' steps per period:'))
+ print(('\t'+str(model.stepdef)))
model.matrices, model.killidx = model._getMatrices(ss=True)
- print('\t'+str(len(model.matrices))+\
+ print(('\t'+str(len(model.matrices))+\
' ('+str(model.cells)+' X '+str(model.cells)+\
- ') matrices were successfully read.')
- print('\tThere are '+str(len(model.killidx))+' zero volume compartments')
+ ') matrices were successfully read.'))
+ print(('\tThere are '+str(len(model.killidx))+' zero volume compartments'))
A = model.matrices[0]
Ainv = inv(A.toarray())
emission = emobj.get(0,model.killidx)
@@ -462,11 +462,11 @@
if not os.access(respath, os.W_OK):
try:
os.mkdir(respath)
- print('Created directory: '+respath)
+ print(('Created directory: '+respath))
except:
sys.exit('Unable to access '+respath+'; aborting')
else:
- print('Overwriting '+respath+' !')
+ print(('Overwriting '+respath+' !'))
resex = expandvec(result[:,0], model.killidx).reshape(model.cells,1)
for t in arange(1,result.shape[1]):
exvec = expandvec(result[:,t], model.killidx).reshape(model.cells,1)
@@ -490,18 +490,18 @@
rundir = os.path.join(subsdir,run)
if not os.access(rundir,os.W_OK):
try:
- print('Creating directory: '+rundir)
+ print(('Creating directory: '+rundir))
os.mkdir(rundir)
- print('Copying TEMPLATE/%s/* to %s') % (run, rundir)
+ print(('Copying TEMPLATE/%s/* to %s') % (run, rundir))
print(' You still have to provide an EMISSIONS.txt file')
commandstring = 'cp '+os.path.join(templatedir,'SS','*')+' '+rundir
os.system(commandstring)
except:
sys.exit('Unable to access '+rundir+'; aborting')
else:
- print('Using existing '+rundir+' !')
+ print(('Using existing '+rundir+' !'))
commandstring = 'cp '+os.path.join(templatedir,'*')+' '+subsdir
- print('Copying %s/* to %s') % (templatedir, subsdir)
+ print(('Copying %s/* to %s') % (templatedir, subsdir))
os.system(commandstring)
diff -ru BETR-Global/solver.py BETR-Global3/solver.py
--- BETR-Global/solver.py 2010-09-27 12:12:56.000000000 +0000
+++ BETR-Global3/solver.py 2021-01-13 17:40:36.000000000 +0000
@@ -14,10 +14,18 @@
import os.path
import time
from scipy.sparse import *
-from scipy.integrate.ode import *
+# from scipy.integrate.ode import *
from scipy.linalg import *
from helpers import *
+from scipy.integrate import ode
+get_return_code = ode.get_return_code
+set_f_params = ode.set_f_params
+set_integrator = ode.set_integrator
+set_jac_params = ode.set_jac_params
+set_solout = ode.set_solout
+
+
class Solver:
def __init__(self, model, stepfile, solvparamfile, initfile):
self.model=model
@@ -103,7 +111,7 @@
ig.set_initial_value(self.y0,0)
ts0=time.time()
for m in range(0, ts):
- print('Integrating month # %s') % (m)
+ print(('Integrating month # %s') % (m))
#ig.set_initial_value(ig.y, ig.t)
current_season = mod(m, self.steps_per_period)
@@ -115,20 +123,20 @@
no_sub_steps = int(ceil(self.stepdef[current_season]
/self.solvparms['maxsubsteplength']))
substeplength = self.stepdef[current_season] / no_sub_steps
- print('\tSeason %d\t %d substeps of length %.5f hours')\
- % (current_season,no_sub_steps, substeplength)
+ print(('\tSeason %d\t %d substeps of length %.5f hours')\
+ % (current_season,no_sub_steps, substeplength))
for i in range(0, no_sub_steps):
ts1=time.time()
- print('\tIntegrating substep {0:d} '.format(i)),
+ print(('\tIntegrating substep {0:d} '.format(i)), end=' ')
ig.integrate(ig.t + substeplength)
if (not ig.successful()):
sys.exit("Integration failed; aborting !")
te1=time.time()
- print(' [%.3f seconds]') % (te1-ts1)
- print('done with month %d') % (m)
+ print((' [%.3f seconds]') % (te1-ts1))
+ print(('done with month %d') % (m))
res=vstack((res, ig.y))
te0=time.time()
- print "time elapsed: %f minutes" % ((te0-ts0)/60.0)
+ print("time elapsed: %f minutes" % ((te0-ts0)/60.0))
res=res.T
return res
Only in BETR-Global: .svn
diff -ru BETR-Global/tempcorrect.py BETR-Global3/tempcorrect.py
--- BETR-Global/tempcorrect.py 2010-09-24 02:18:55.000000000 +0000
+++ BETR-Global3/tempcorrect.py 2021-01-13 17:35:59.000000000 +0000
@@ -30,7 +30,7 @@
## then perfect persistence would be much easier to implement
chem_spatial_dict={}
- for k in compdict.keys():
+ for k in list(compdict.keys()):
## initialize record-array
pa=zeros(par.shape, dtype=dtype([('k_reac','f8'), ('Kaw','f8'),
('Koa','f8'),('Kow','f8')]))
@@ -54,8 +54,8 @@
## weight degradation in air according to OH-radical concentration
## ATT : This special treatment of air1 and air2 should be generalized
## somehow.
- if 1 in chem_spatial_dict.keys():
+ if 1 in list(chem_spatial_dict.keys()):
chem_spatial_dict[1]['k_reac']*=par['OHair1'] # upper air
- if 2 in chem_spatial_dict.keys():
+ if 2 in list(chem_spatial_dict.keys()):
chem_spatial_dict[2]['k_reac']*=par['OHair2'] # lower air
return(chem_spatial_dict)
diff -ru BETR-Global/testrun.py BETR-Global3/testrun.py
--- BETR-Global/testrun.py 2010-09-28 03:30:26.000000000 +0000
+++ BETR-Global3/testrun.py 2021-01-13 17:35:59.000000000 +0000
@@ -1,88 +1 @@
-## This is a test-run for BETR-Reserach ####
-import BETRS
-reload(BETRS)
-from BETRS import *
-
-## constructs the 'default' - model, equivalent to BETR-Global for
-## chemical nr. 11 = PCB-118
-m=Model(11)
-
-## loads the emission-scenario derived from night-light emissionsy
-## (the same as in the demo-run of BETR-Global
-m.update_emissions('emissions_default.txt')
-
-## solves for steady-state
-m.solve_ss()
-m.output_ss(units=['mol_per_m3', 'kg','Pa'], netcdf=True)
-
-## look at results
-result=load('Output/default/ss_out.cpk')
-c=result[2]['mol_per_m3']
-cgrid=c.reshape(12,24)
-clat=mean(cgrid, axis=1)
-from pylab import *
-plot(clat); show()
-
-## get full steady-state system matrix
-ssmat=m.get_avgMatrix()
-
-### overall persistence
-## get emission vector
-ssem=m.get_ssEmission()
-Pov = sum(m.ss_res)/sum(ssem.toarray())/24
-print("Overall persistence of %s:") % (m.chemdict['Name'])
-print("%s days") % (Pov)
-
-## get full steady-state mass vector
-## mass flows from and to African soil (region 158, compartment 6)
-soilflow_to=ssmat[tocell(158,6,m),:]*m.ss_res
-soilflow_from=ssmat[:,tocell(158,6,m)]*m.ss_res[tocell(158,6,m)]
-
-## flows from soil
-happens=cell2regcomp(where(soilflow_to > 0 )[0],m)
-print("Flows happens to soil to compartments %s") %(happens[1])
-print("in regions %s:") % (happens[0])
-balance_to=soilflow_to[where(soilflow_to > 0 )[0]]
-print(balance_to)
-
-## flows to soil
-happens=cell2regcomp(where(soilflow_from > 0 )[0],m)
-print("Flows happens from soil to compartments %s") %(happens[1])
-print("in regions %s:") % (happens[0])
-balance_from=soilflow_from[where(soilflow_from > 0 )[0]]
-print(balance_from)
-
-## degradation
-deg=sum(soilflow_from)
-
-print(" ***** net balance ****")
-balance=(soilflow_to-soilflow_from)[where(soilflow_from > 0 )[0]]
-print("compatment: 2 3 4")
-print(" "+str(balance))
-print("degradation loss: %s") % (deg)
-
-## plot it
-import pylab
-pos=array([0,1,2,3])
-posc=pos+0.4
-fig=barh(pos,height=0.8,width=insert(balance,0,deg))
-yticks(posc,('degradation','air','vegetation','freshwater'))
-show()
-
-## solve dynamically
-m.update_emissions('emissions_firstyear.txt')
-m.update_solver()
-m.solve_dyn(5)
-m.output_dyn(units=['mol','mol_per_m3','kg','Pa'],netcdf=True)
-
-## plot timeseries in Northern European lower air
-dynresult=load('Output/default/dyn_out.cpk')
-plot(dynresult[2]['mol_per_m3'][37,:])
-fig=figure()
-
-log time-series
-fig2=plot(log(dynresult[2]['mol_per_m3'][37,:]))
-show()
-
-###############################################################################
-
+Non
\ No newline at end of file
diff -ru BETR-Global/validate.py BETR-Global3/validate.py
--- BETR-Global/validate.py 2010-09-27 08:49:55.000000000 +0000
+++ BETR-Global3/validate.py 2021-01-13 17:35:59.000000000 +0000
@@ -21,13 +21,14 @@
from numpy import *
import sys
from scipy import sparse
-import cPickle
+import pickle
import subprocess
import os
import shutil
from pylab import *
import BETRS
-reload(BETRS)
+import imp
+imp.reload(BETRS)
from BETRS import *
import helpers
@@ -43,11 +44,11 @@
def validate(level, chemicals):
print('Validating BETR-Research for chemicals')
print(chemicals)
- print('at level {0:.3e}'.format(level))
+ print(('at level {0:.3e}'.format(level)))
valid=True
difflist=[]
for chem in chemicals:
- print('Calculating system matrices for chemical %d') % (chem)
+ print(('Calculating system matrices for chemical %d') % (chem))
m=Model(chem,controlfile='control_validation.txt')
shutil.copy(os.path.join('Output','default','matrixdump.cpk'),
os.path.join('Validation'))
@@ -55,7 +56,7 @@
## load BETRS-output
fn=os.path.join('Output','default','matrixdump.cpk')
f=open(fn,'r')
- Dmatpy=cPickle.load(f)
+ Dmatpy=pickle.load(f)
f.close()
## load BETR-VBA output
@@ -66,7 +67,7 @@
sys.exit('Did not find the file {0:s}. Aborting !'.format(fnvba))
for mo in arange(0,12):
- print('comparing substance %s, month %i') % (chem,mo+1)
+ print(('comparing substance %s, month %i') % (chem,mo+1))
Dres=Dmatpy[mo]
## get empty volume elements
killidx=findemptycells(Dres)
@@ -101,10 +102,10 @@
diff1=where(logical_and(Dres==0, Dvba!=0))
diff2=where(logical_and(Dres!=0, Dvba==0))
if diff1[0].size > 0 or diff2[0].size > 0:
- print('Validation failed ! There are zeros in one matrix that'\
- +' are non-zero in the other:')
- print('research == 0, vba != 0: {0:s}'.format(diff1))
- print('research != 0, vba == 0: {0:s}'.format(diff2))
+ print(('Validation failed ! There are zeros in one matrix that'\
+ +' are non-zero in the other:'))
+ print(('research == 0, vba != 0: {0:s}'.format(diff1)))
+ print(('research != 0, vba == 0: {0:s}'.format(diff2)))
valid=False
else:
print('Test for equality of sparsity pattern: passed.')
@@ -116,23 +117,23 @@
fromc=list(cell2regcomp(s1[1],m)[1])
toc=list(cell2regcomp(s1[0],m)[1])
cells=list(cell2regcomp(s1[0],m)[0])
- ft=zip(fromc,toc)
+ ft=list(zip(fromc,toc))
diffdict={}
for i in enumerate(ft):
- if diffdict.has_key(i[1]):
+ if i[1] in diffdict:
diffdict[i[1]].append(cells[i[0]])
else:
diffdict[i[1]]=[cells[i[0]]]
if bool(diffdict):
- print('Validation not passed, relative differences > {0:.3e}'\
- +' occured :'.format(level))
+ print(('Validation not passed, relative differences > {0:.3e}'\
+ +' occured :'.format(level)))
print(diffdict)
print('erroneous from-to:')
- print(unique(ft))
+ print((unique(ft)))
valid=False
else:
- print('Matices equal within rtol={0:.3e}'.format(level))
- print(60*'*')
+ print(('Matices equal within rtol={0:.3e}'.format(level)))
+ print((60*'*'))
if valid:
print('Overall validation: PASSED')
else:
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment