import numpy as np
import math
import os
import shutil
import collections
import glob
from typing import Any, Callable
from numpy.typing import NDArray
GRIDINPUTFILE = "ParameterGrid.in"
""" The additional input file which includes all the parameters that are included in the grid. """
[docs]
def chgriddir(gridname: str):
"""
Changes the current working directory to the grid directory.
Parameters
----------
gridname :
The name of the grid (the directory with the models).
"""
if not os.getcwd().endswith("/" + gridname):
os.chdir(gridname)
[docs]
def genparamentry(name: str, value: Any)-> str:
"""
Generates one Parameter entry in the ProDiMo Parameter.in format.
If name is ``fcarbon`` then it it is assumed that the DIANA standard
opacities are used,and a proper ``! NDUST`` entry is created (with a fixed vacuum fraction of 0.25).
Parameters
----------
name :
The name of the ProDiMo parameter.
value : Any
The value of the Parameter, is converted to a string.
Returns
-------
str
The full str that can be written to a file in Parameter.in format.
"""
if name == "fcarbon":
fsil = 1.0 - (0.25 + value)
entry = "3 ! NDUST \n"
entry = entry + " " + str(fsil) + " Mg0.7Fe0.3SiO3[s] \n"
entry = entry + " " + str(value) + " amC-Zubko[s] \n"
entry = entry + " 0.25 vacuum[s] \n"
else:
entry = str(value) + " ! " + name + " \n"
return entry
[docs]
def genvalues(param: tuple) -> NDArray[Any]:
"""
Generates the values for the given parameter, depending on type.
Parameters
----------
param :
A tuple with x entries depending on type. see :func:`make_grid`.
Returns
-------
NDArray[Any]
A numpy array with the generated values for that `param`. Can be of any type (str,float,int,boolean).
"""
# boolean value
if len(param) == 1:
# don't use the python representation, this should work for both gfortran and ifort
return np.array([".true.", ".false."])
elif len(param) == 2: # fixed values
return np.array(param[1])
else: # normal Parameters
start = param[1]
end = param[2]
steps = param[3]
linlog = param[4]
if linlog == "lin":
return np.linspace(start, end, steps)
elif linlog == "log":
return np.logspace(math.log10(start), math.log10(end), steps, base=10.0)
else:
print(f"ERROR: Unknown spacing type {linlog} for {param[0]}")
return None
[docs]
def get_modeldirs(modeldirs: list | None = None, gridname: str | None = None):
"""
Returns all grid models of the current grid (current directory) or for the grid given by `gridname`.
Parameters
----------
modeldirs :
The list of model directories. If `None` all directories within the current directory starting with `model` will be selected.
gridname :
The name of the grid (it's directoryname to be precise). If gridname is given search is done with ``gridname+"/model*/"``
Returns
-------
list
If `modeldirs` is not `None` `modeldirs` is returned (so nothing is done).
Otherwise all directory names staring with `[gridname/]model*/` are returned.
"""
if modeldirs is None:
if gridname is not None:
modeldirs = glob.glob(os.path.join(gridname, "model*/"))
else:
modeldirs = glob.glob("model*/")
return modeldirs
[docs]
def make_grid(gridname: str, params: list, indir: str | None = None) -> tuple[list, tuple]:
"""
Generates a grid of prodimo models.
A directory with the name ``gridname`` is created. Within this directory for
each parameter combination a directory is created. This directory is initially
copied from the directory given by ``indir``, where certain files (*.out,*.log,*.fits) are ignored.
The generated parameters are added to a ``ParameterGrid.in`` file. If ``ParameterGrid.in`` does not exist
it will be created.
Parameters
----------
gridname :
The name of the grid (the directory with the models).
params :
The list of parameters. Depending of the type different entries per Parameter are required.
- numeric: Entry has to be like ``(paramname,startvalue,endvalue,numberofsteps,type)``. ``type`` can either be ``"lin"`` or ``"log"`` for linear or logarithmic spacing
- boolean: Entry is simply ``(paramname)``
Special parameters:
- fcarbon: Only works if DIANA opacities are used. Changes the carbon fraction and the silicate part accordingly.
indir :
an input or starting directory containing the usual |prodimo| input
files. This directory is copied into a new model.
If it is `None` only the ParameterGrid.in are created
Returns
-------
tuple[list,tuple]
First entry is the list of modelnames (director values). Second one is a tuple with an entry for each
parameter, giving the generated values for this parameter.
"""
# make the list for the meshgrid
values = list()
for param in params:
values.append(genvalues(param))
# this is just to return the generate parameter values for e.g. logging
gridvals = np.meshgrid(*values, sparse=True)
# this is the real grid
grid = np.meshgrid(*values)
# create the directory for the grid, if it exists this throws an exception
os.mkdir(gridname)
# if something happens in the block below we want to delete the newly created grid directory
try:
# create an iterator that loops over all indices of the first grid array (all combinations)
it = np.nditer(grid[0], flags=["multi_index"])
imodel = 0
modelnames = list()
while not it.finished:
# create the model directory
modelname = "model" + "{:07d}".format(imodel)
modeldir = os.path.join(gridname,modelname)
if indir is not None:
shutil.copytree(
indir, modeldir, ignore=shutil.ignore_patterns("*.out", "*.log", "*.fits")
)
else:
os.mkdir(modeldir)
fparam = open(os.path.join(modeldir, GRIDINPUTFILE), "a+")
for iparam in range(len(params)):
# fparam.write(str(grid[iparam][it.multi_index])+" ! "+params[iparam][0] +" \n")
fparam.write(genparamentry(params[iparam][0], grid[iparam][it.multi_index]))
fparam.close()
imodel = imodel + 1
modelnames.append(modelname)
it.iternext()
except Exception as e:
os.rmdir(gridname)
raise e
return modelnames, gridvals
[docs]
def run_grid(gridname: str, runProDiMo: str | Callable, modeldirs: list | None = None):
"""
Runs the grid.
Changes to the grid directory and runs each model by either calling
a passed python function or a system command string (see parameters).
Parameters
----------
gridname :
The name of the grid (the directory with the models).
runProDiMo :
If it is a ``str`` the string is interpreted as a system command.
Any occurrence of `$MODELNAME$` in the given string is replaced by the actual model name (model directory).
If the parameter is a python function (a `Callable`). This function is called with the `modeldir` as a parameter.
modeldirs :
A list of all models in the grid (the directory name of each model).
If ``None`` the list is created automatically using :func:`prodimop.grid.get_modeldirs`.
Only useful to run a selection of models from the grid (i.e. the failed ones).
"""
chgriddir(gridname)
modeldirs = get_modeldirs(modeldirs)
for modeldir in modeldirs:
if isinstance(runProDiMo, collections.abc.Callable):
print("run " + modeldir + ", exec. function: " + runProDiMo.__name__)
runProDiMo(modeldir)
else:
runProDiMoCMD=runProDiMo
if not runProDiMoCMD.endswith(GRIDINPUTFILE):
runProDiMoCMD=runProDiMoCMD+" "+GRIDINPUTFILE
runProDiMoCMD = runProDiMoCMD.replace("$MODELNAME$", modeldir)
os.chdir(modeldir)
print("run " + modeldir + ", exec. command: " + runProDiMoCMD)
os.system(runProDiMoCMD)
os.chdir("..")
[docs]
def check_grid(gridname, modeldirs=None):
"""
Checks if all models look okay.
This routine checks if finished.out exists for all models of the grid.
Parameters
----------
gridname : str
The name of the grid (the directory with the models).
modeldirs : list
a list of all models in the grid (directory name of each model).
If `modeldirs` is `None` all directories with names starting with `model` are
considered as potential grid models.
"""
chgriddir(gridname)
# guess the model directories
modeldirs = get_modeldirs(modeldirs)
for modeldir in modeldirs:
if not os.path.isfile(modeldir + "/finished.out"):
print("Model " + modeldir + " failed:")
[docs]
def sel_lowest_chisquare(
gridname: str, modeldirs: list | None = None, tolerance: float | None = None
) -> tuple[list, list]:
"""
Selects the model(s) with the lowest chi squared. If a tolerance is
given all models within this tolerance, measured relative to the minimum
chi square are selected.
Parameters
----------
gridname :
The name of the grid (the directory with the models).
modeldirs :
a list of all models in the grid (directory name of each model).
If `modeldirs` is `None` all directories with names starting with `model` are
considered as potential grid models.
"""
chgriddir(gridname)
modeldirs = get_modeldirs(modeldirs)
nmodels = len(modeldirs)
chisquares = np.ndarray(shape=(nmodels))
chisquares[:] = 1.0e100
# go through all models and read the chi square
for i in range(nmodels):
modeldir = modeldirs[i]
if os.path.isfile(modeldir + "/finished.out"):
ffin = open(modeldir + "/finished.out")
for line in ffin:
if "total chi" in line:
val = float(line.split("=")[1])
chisquares[i] = val
if all(chisquares == 1.0e100):
print("ERROR: It seems none of the models properly finished.")
return [], []
minchi = np.min(chisquares)
if tolerance is not None:
tol = minchi * tolerance
idx = np.where(np.abs(chisquares - minchi) < tol)[0]
else:
idx = [np.argmin(chisquares)]
return np.array(modeldirs)[idx], chisquares[idx]
[docs]
def sel_param_val(gridname: str, param: str, value: float | int | str, modeldirs: list[str] | None = None)-> list[str]:
"""
Very primitive and inefficient method to select models having a certain
value of a parameter. It only searches in ``grid.GRIDINPUTFILE``.
Parameters
----------
gridname :
The (directory) name of the grid.
param :
The parameter name to search for.
value :
The value to search for. In case of numbers the read value from the grid Parameter file
is converted to the type of `value`.
modeldirs :
List of model directories of the grid. If ``None`` :func:`~prodimopy.grid.get_modeldirs` is used.
"""
modeldirs = get_modeldirs(modeldirs,gridname=gridname)
nmodels = len(modeldirs)
selmodels = list()
for i in range(nmodels):
modeldir = modeldirs[i]
filep=os.path.join(modeldir, GRIDINPUTFILE)
if os.path.isfile(filep):
ffin = open(filep)
for line in ffin:
if param in line:
if isinstance(value,float):
val = float(line.split("!")[0])
elif isinstance(value,int):
val = int(line.split("!")[0])
else:
val = (line.split("!")[0]).strip()
if val == value:
selmodels.append(modeldir)
return selmodels