https://github.com/steven-murray/pycamb
Tip revision: 6db4aea707d6594d166ced8ce1202245b0bc53f7 authored by Steven Murray on 16 November 2015, 02:44:22 UTC
more formatting
more formatting
Tip revision: 6db4aea
generatePyCamb.py
#!/usr/bin/env python
import numpy.f2py
import re
import warnings
# All the numerical parameters that we will let the user pass to camb.
# If the parameter name starts with @ that signifies it is *not* part of the camb params type but
# is defined elsewhere. For example, the dark energy parameters w_lam and cs2_lam are defined in the
# LambdaGeneral module in equations.F90
numericalParams = ['omegab', 'omegac', 'omegav', 'omegan', 'H0', 'TCMB',
'yhe', 'Num_Nu_massless', 'Num_Nu_massive', 'omegak',
'reion__redshift', 'reion__optical_depth', 'reion__fraction', "reion__delta_redshift",
'Scalar_initial_condition', 'scalar_index', 'scalar_amp',
'scalar_running', 'tensor_index', 'tensor_ratio', 'NonLinear',
'@lAccuracyBoost', '@lSampleBoost', '@w_lam', '@cs2_lam',
'@AccuracyBoost', 'transfer__kmax', 'transfer__k_per_logint',
'@ThreadNum',
]
defaultValues = {'@AccuracyBoost':1.0,
'@lAccuracyBoost':1.0,
'@lSampleBoost':1.0,
'@ThreadNum':0,
'scalar_amp':2.1e-9,
'@w_lam':-1.0,
'@w_perturb':False,
'@cs2_lam':1.0,
'transfer__kmax':2,
'transfer__k_per_logint':0,
'NonLinear':0}
# The boolean options that the user can pass to camb.
logicalParameters = ['WantScalars', 'WantTensors', 'reion__reionization', 'reion__use_optical_depth', '@w_perturb', 'DoLensing']
# These map the friendly names in the parameters above to the initial power vectors
alias = {
'scalar_index': 'InitPower__an(1)',
'scalar_amp': 'InitPower__ScalarPowerAmp(1)',
'scalar_running':'InitPower__n_run(1)',
'tensor_index':'InitPower__ant(1)',
'tensor_ratio': 'InitPower__rat(1)'
}
nparams = len(numericalParams) + len(logicalParameters)
endl = "\n"
from os import system
def template_file(infile, outfile, parameters):
text = open(infile).read()
for param, value in parameters.items():
pattern = r"\$%s\$" % param
count = 0
for i in xrange(1000):
text, c = re.subn(pattern, str(value), text)
count += c
if c == 0: break
else:
raise ValueError("Got into a long loop on template_file. Did your template values contain $$ symbols?")
if count == 0:
warnings.warn("Unused template parameter: %s" % param)
if "$" in text: raise ValueError("Template contains unfilled parameters")
open(outfile, "w").write(text)
#Generate the fortran subroutine called by the first one to set the CAMBparams instance's values from the vector we pass it.
def makeFortranParamSetter(numericalParameters, logicalParameters, defaultValues):
code = ""
nn = len(numericalParameters)
nm = len(numericalParameters) + len(logicalParameters)
for p, friendly_name in enumerate(numericalParameters):
if alias.has_key(friendly_name):
name = alias[friendly_name]
else:
name = friendly_name
fortranName = "P%%%s" % name.replace('__', '%')
if name.startswith("@"):
fortranName = "%s" % (name.lstrip("@").replace('__', '%'))
if name not in defaultValues:
raise ValueError("Global parameter %s requires default argument (or it is retained between runs). Edit the top of %s" % (name, __file__))
if friendly_name in defaultValues:
# print "%s has default" % name
code += """
if (paramVec(%d) .ne. -1.6375e30) then
%s = paramVec(%d)
else
%s = %e
endif
""" % (p + 1, fortranName, p + 1, fortranName, defaultValues[friendly_name])
else:
# print "%s has NO default" % name
code += "if (paramVec(%d) .ne. -1.6375e30) %s = paramVec(%d)\n" % (p + 1, fortranName, p + 1)
for p, friendly_name in enumerate(logicalParameters):
if alias.has_key(friendly_name):
name = alias[friendly_name]
else:
name = friendly_name
fortranName = "P%%%s" % name.replace('__', '%')
if name.startswith("@"):
fortranName = "%s" % (name.replace('__', '%').lstrip("@"))
if name not in defaultValues:
raise ValueError("Global parameter %s require default argument (or it is retained between runs). Edit the top of %s" % (name, __file__))
if friendly_name in defaultValues:
if defaultValues[friendly_name]:
defaultValue = '.true.'
else:
defaultValue = '.false.'
code += """
if (paramVec(%d) .ne. -1.6375e30) then
if(paramVec(%d) .ne. 0.0) then
%s=.true.
else
%s=.false.
endif
else
%s=%s
endif
""" % (nn + p + 1, nn + p + 1, fortranName, fortranName, fortranName, defaultValue) + endl
else: #No default value
code += """
if (paramVec(%d) .ne. -1.6375e30) then
if(paramVec(%d) .ne. 0.0) then
%s=.true.
else
%s=.false.
endif
endif
""" % (nn + p + 1, nn + p + 1, fortranName, fortranName)
return code
import pprint
#Generate the python code that will wrap the fortran.
def makePython(numericalParameters, logicalParameters, defaultValues):
param_string = (", ".join([name.lstrip("@") for name in numericalParameters + logicalParameters])) + "\n"
alias_string = "\n\t".join(["%s ---> %s" % (name.strip("@"), alias[name]) for name in alias])
defparam_string = "\n".join("%s (%f)" % (name.lstrip("@"), param) for (name, param) in defaultValues.items())
params = {"numericalParameters":numericalParameters,
"logicalParameters":logicalParameters,
"defaultValues":defaultValues,
"param_string":param_string,
"alias_string":alias_string,
"defparam_string":defparam_string,
}
template_file("templates/pycamb.py.template", "src/__init__.py", params)
def makeFortran(numericalParams, logicalParameters, defaultValues):
number_parameters = len(numericalParams) + len(logicalParameters)
param_caller_function = makeFortranParamSetter(numericalParams, logicalParameters, defaultValues)
params = {"number_parameters":number_parameters, "param_caller_function":param_caller_function}
template_file("templates/py_camb_wrap.f90.template", "src/py_camb_wrap.f90", params)
def main():
makeFortran(numericalParams, logicalParameters, defaultValues)
makePython(numericalParams, logicalParameters, defaultValues)
if __name__ == "__main__":
main()