from __future__ import division
from __future__ import print_function
from __future__ import unicode_literals
import difflib
import glob
import inspect
import os
from astropy.io.fits.diff import FITSDiff
import prodimopy.read_slab as pread_slab
import numpy as np
[docs]
class CompareAbs(object):
"""
An "abstract" class for comparing some kind of |prodimo| model.
A subclass of this class needs to implement the necessary compare routine(s).
(example see :class:`Compare`)
.. todo::
* make tolerances configurable via config file (see punit checkModels).
"""
[docs]
def diffArrayRel(self, a, aref, diff):
"""
Checks the relative difference between two arrays.
If the arrays do not have the same shape they are considered
as inequal (return `False`).
Parameters
----------
a : array_like(float,ndim=whatever)
an array.
aref: array_like(float,ndim=same as a)
the reference array for comparison.
diff : float
if the values of the arrays differe only by <diff they are considered
as equal.
"""
if a.shape != aref.shape:
return False, None
da = np.absolute(a / aref - 1.0)
if da.max() >= diff:
return False, da
return True, da
[docs]
def diffArray(self, a, aref, rtol=1.0e-3, atol=0.0, nlimfrac=0.0, mask=None):
"""
Checks the difference between two arrays using the np.isclose routine but
additionally allows to set a limit for how many elements can be different
in the arrays to still be considered as equal (or similar).
If the arrays do not have the same shape they are considered
as inequal (return `False`).
Parameters
----------
a : array_like(float,ndim=whatever)
an array.
aref: array_like(float,ndim=same as a)
the reference array for comparison.
rtol : float
Relative tolerance as used in :func:`numpy.iclose`.
atol : float
Absolute tolerance as used in :func:`numpy.iclose`.
nlimfrac : float
Relative limit of elements that are allowed to be not equal.
If less than `nlim * len(aref)` elements are different, the arrays
are still considered as equal.
mask : mask for the array
Mask for the elements that should not be included in the evaluation.
Those elements are set to True in any case.
"""
if a is None and aref is None:
return True, None, None
dbool = np.isclose(a, aref, rtol=rtol, atol=atol, equal_nan=True)
# Calculate the difference in percent, only for the values that are not masked
# deal with zero values, in the reference array
maskZ = aref != 0
da = aref * 0.0
da[maskZ] = np.absolute(a[maskZ] / (aref[maskZ]) - 1.0) * 100
da[~maskZ] = np.absolute(a[~maskZ] + 1.0e-100 / (aref[~maskZ] + 1.0e-100) - 1.0) * 100
# set the value that are masked to zero (not considered in for e.g. max error)
if mask is not None:
dbool[mask] = True
da[mask] = 0.0
# set the values that are considered as similar to 0 ... because of atol
# da[dbool]=0.0
# check if the number of inequal elements is lower than allowed.
if nlimfrac > 0.0:
eq = (np.size(dbool) - np.count_nonzero(dbool)) < int(nlimfrac * np.size(dbool))
else:
eq = dbool.all()
# if da.max()>=diff:
# return False,da
# return True,da
return eq, da, dbool
[docs]
def diff(self, val, valref, rtol=1.0e-10, atol=0.0):
"""
Checks the relative difference between two values.
FIXME: is not inconsisten with diffArr
Parameters
----------
a : float
a value.
aref: float
the reference value
diff : float
if the two values differe only by <diff the are considered
as equal.
"""
dbool = np.isclose(val, valref)
d = abs(val / valref - 1.0) * 100
return dbool, d
[docs]
def diffFile(self, fname):
"""
Compares the file with the name given by `fname` (a file produced by |prodimo|)
for the the two models.
Parameters
----------
fname : str
the Filename to compare, (must exist for both models)
"""
# check if file exist
fe = os.path.isfile(self.m.directory + "/" + fname)
feref = os.path.isfile(self.mref.directory + "/" + fname)
# File does not exist in any model -> that is okay
if (not fe) and (not feref):
return True, None
# if the file exists in only one model, the routine returns false!
if fe ^ feref:
return False, None
# now compare the files
f = open(self.m.directory + "/" + fname)
fref = open(self.mref.directory + "/" + fname)
fcont = f.readlines()
fcontRef = fref.readlines()
diff = list(difflib.unified_diff(fcont, fcontRef, n=0))
if len(diff) == 0:
return True, None
else:
return False, diff
[docs]
def diffFitsFile(self, fname, atol=0.0, rtol=2.0e-7):
"""
Compares two fits Files using the astropy routines.
Parameters
----------
fname : str
the Filename to compare, (must exist for both models)
"""
# check if file exist
fe = os.path.isfile(self.m.directory + "/" + fname)
feref = os.path.isfile(self.mref.directory + "/" + fname)
# File does not exist in any model -> that is okay
if (not fe) and (not feref):
return True, None
# if the file exists in only one model, the routine returns false!
if fe ^ feref:
return False, None
# FIXME: in restart.dat this CONVERGENCE keywordk makes problems,
# I think the KEYWORD name is too long (see the fitsheader of reastart.fits.gz)
# do not test for MAINIT, so that also models that were rerun can be compared
try:
diff = FITSDiff(
self.m.directory + "/" + fname,
self.mref.directory + "/" + fname,
ignore_keywords=["CONVERGENCE", "MAINIT", "TIMESTAMP", "REVISION"],
numdiffs=2,
atol=atol,
rtol=rtol,
)
except Exception as e:
print("Exception in FITSDiff: ", e)
return False, None
if diff.identical:
return True, None
else:
return False, diff.report()
[docs]
def doAll(self):
"""
Utility function to call all `compare*` methods.
Prints out what function failed and the errors.
Assume a naming convention. The method needs to start
with `compare`. If the compare method uses the :func:`diffFile`,
the method name should start with `compareFile`.
"""
# defines the compare* methods
def iscompare(object) -> bool:
if inspect.ismethod(object):
return object.__name__.startswith("compare")
else:
return False
# print(inspect.getmembers(self))
for name, cfun in inspect.getmembers(self, predicate=iscompare):
print("{:30s}".format(name + ": "), end="")
ret = cfun()
ok = ret[0]
val = ret[1]
valb = None
if len(ret) > 2:
valb = ret[2]
if ok:
print("{:8s}".format("OK"), end="")
else:
print("{:8s}".format("FAILED"), end="")
# Special treatment for the case of comparing Files
if name.startswith("compareFile"):
if val is None and not ok:
print(" File does only exist in one model.")
elif val is not None:
for line in val:
if line.startswith("@@"):
print("")
print(line, end="")
else:
print("")
elif name.startswith("compareFitsFile"):
if val is None and not ok:
print(" File does only exist in one model.")
elif val is not None:
print(val)
else:
print("")
elif name.startswith("compareChemnet"):
if not ok:
if val is not None:
print(" Chemical networks differ: ")
print(" Number of different reactions (and or rates): ", val[0])
print(" Number of reactions only in reference model: ", val[1])
print(" Number of reactions only in other model: ", val[2])
else:
print(" Chemical networks are different.")
else:
print("")
else:
# TODO: also add the number of wrong points in the array
if val is not None:
if valb is not None:
nnoteqp = (np.size(valb) - np.count_nonzero(valb)) / np.size(valb) * 100
if np.abs(np.nanmax(val)) > 0.0:
maxidx = [int(value) for value in np.unravel_index(np.nanargmax(val), val.shape)]
else:
maxidx = "[-]"
print(
" different/Avg Err/Max Err/ Max Err Index: ",
"{:8.2e}%".format(nnoteqp),
"/",
"{:8.2e}%".format(np.nanmean(val)),
"/",
"{:8.2e}%".format(np.nanmax(val)),
"/",maxidx
)
# print the first 10 species (lines) that show the largest differences
if name.startswith("compareLineEstimates") and not ok:
sortidx = np.flip(np.argsort(val), axis=0)
lastIdents = list()
iprint = 0
for i in range(len(val)):
# if self.m.lineEstimates[sortidx[i]].flux>1.e-24 or self.mref.lineEstimates[sortidx[i]].flux>1.e-24:
# check also if the value is not masked (then it is true in valb anyway)
if (
not valb[sortidx[i]]
and self.m.lineEstimates[sortidx[i]].ident not in lastIdents
):
print(
"{:40s}".format(str(self.m.lineEstimates[sortidx[i]])),
"{:40s}".format(str(self.mref.lineEstimates[sortidx[i]])),
)
lastIdents.append(self.m.lineEstimates[sortidx[i]].ident)
iprint += 1
if iprint > 5:
break
print("")
# FIXME: currently only works for CompareMC
if (
name.startswith("compareAbundances")
and not ok
and isinstance(self, CompareMc)
):
# find out the ages with maximum error
# the age index for each species with the maximum error
maxa = np.argmax(val, axis=0)
# print(maxa[self.m.species.index("Cl+")])
# print(self.m.abundances[:,self.m.species.index("Cl+")])
# print(self.mref.abundances[:,self.mref.species.index("Cl+")])
# no make a list for those values and sort it (find the maximum species index)
diffspec = np.array([val[idx, i] for i, idx in enumerate(maxa)])
# print(diffspec,np.max(diffspec))
sortidx = np.flip(np.argsort(diffspec))
lastIdents = list()
iprint = 0
for i in range(len(sortidx)):
# if self.m.lineEstimates[sortidx[i]].flux>1.e-24 or self.mref.lineEstimates[sortidx[i]].flux>1.e-24:
if (
not valb[maxa[sortidx[i]], sortidx[i]]
and self.m.species[sortidx[i]] not in lastIdents
):
print(
"{:20s}".format(self.m.species[sortidx[i]]),
"{:10.3f}".format(val[maxa[sortidx[i]], sortidx[i]]),
# need +1 here because the 0 age index is not considerd in the comparison, see compare_abundances
"{:40s}".format(
str(self.m.abundances[maxa[sortidx[i]] + 1, sortidx[i]])
),
"{:40s}".format(
str(self.mref.abundances[maxa[sortidx[i]] + 1, sortidx[i]])
),
)
lastIdents.append(self.m.species[sortidx[i]])
iprint += 1
if iprint > 5:
break
elif val is None and not ok:
print(" Value/array is None in only one model or have different shapes.")
else:
print("")
[docs]
class Compare(CompareAbs):
"""
Class for comparing to |prodimo| models of type :class:`~prodimopy.read.Data_ProDiMo`
Every compare Function returns true or false, and the relative differences
(in case of arrays these are arrays).
Can be used in e.g. automatic testing routines or in simple command line
tools to compare |prodimo| model results.
"""
def __init__(self, model, modelref):
"""
Parameters
----------
model : :class:`~prodimopy.read_mc.Data_mc`
The model one wants to compare.
modelref : :class:`~prodimopy.read_mc.Data_mc`
The reference model to compare with.
Attributes
----------
"""
if model is None or modelref is None:
raise TypeError("Either model or modelref is None, no comparison possible.")
self.m = model
""" :class:`~prodimopy.read.Data_ProDiMo` :
The model one wants to compare.
"""
self.mref = modelref
""" :class:`~prodimopy.read.Data_ProDiMo` :
The reference model agains which model is compared..
"""
self.tolerances = {}
""" dict :
A dictionary of tolerances. The naming convention has to follow the naming
of the compare Function. e.g. If one wants to provide tolerances for
`compareAbundances` the key names are `Abundances.rtol` and `Abundances.atol`.
FIXME: This is still experimental and not implemented for all compare routines.
FIXME: This should be something that is in CompareAbs
"""
# the default allowed difference
# self.d=1.e-3
# self.dcdnmol=0.5 # the allowed differences for the column densities (radial and vertical)
# # chemistry is difficult and uncertain :) FIXME: this has to become better, there
# # a simply some columns which usually fail an require time-dependent chemistry and
# # the outcome seems to depend on numerical uncertainties ...
# # (e.g. 7% point in the column density goes crazy ... still oka
# self.fcdnmol=0.1
# self.lcdnmol=1.e-10 # lower limit for the column density to avoid problems with "empty" columns
# self.dTgas=5e-1 # FIXME: 50% Tgas is also quite sensitive, but lower would be better
# self.dZetaCR=self.d
# self.dZetaX=5.e-1 # FIXME: 50% , should not be that high
# self.dHX=3.e-1 # FIXME: 50% , should not be that high
# self.lZetaX=1.e-25 # lower limit for the check
# # self.specCompare=("e-","H2","CO","H2O","Ne+","Ne++","H+")
# # self.specCompare=("N2","N2#","CO#","H2O#","H3","H3+","HCO+","HN2+","SO2","SiO")
# # self.specCompare=("CO","CN")
self.specCompare = (
"e-",
"H",
"H2",
"CO",
"H2O",
"N2",
"N2#",
"CO#",
"H2O#",
"H3+",
"HCO+",
"HN2+",
"SO2",
"SiO",
"Ne+",
"Ne++",
"H+",
"OH",
"C+",
"S+",
"Si+",
"N+",
"CN",
"HCN",
"NH3",
)
# switch of some logging
self.m._log = False
self.mref._log = False
[docs]
def getTols(self, atol=None, rtol=None, nlimfrac=None):
"""
Get the tolerance values from the `self.tolerances` dictionary
Returns
-------
dict: containing {"atol": value, "rtol": value } if set. Can also be empty
"""
# This gives me the name of the caller function. Check if it is a compare
# function
callername = inspect.currentframe().f_back.f_code.co_name
if not callername.startswith("compare"):
raise ValueError("getTols can only be called from a compare function!")
# strip the compare
valname = callername.replace("compare", "")
# The dictionary that contains the return values, if any
retdict = {}
for key, defval in zip(["atol", "rtol", "nlimfrac"], [atol, rtol, nlimfrac]):
if valname + "." + key in self.tolerances:
retdict[key] = self.tolerances[valname + "." + key]
elif defval is not None:
retdict[key] = defval
return retdict
[docs]
def compareX(self):
"""
Compares the x grid points.
"""
return self.diffArray(self.m.x, self.mref.x, **self.getTols(rtol=1.0e-5))
[docs]
def compareZ(self):
"""
Compares the z grid points.
"""
ret = self.diffArray(self.m.z, self.mref.z, **self.getTols(rtol=1.0e-5))
return ret
[docs]
def compareLineFluxes(self):
"""
Compares the line fluxes
Currently assumes that both models include the same lines in the same order.
"""
if self.m.lines is None and self.mref.lines is None:
return True, None
if self.m.lines is not None and self.mref.lines is None:
return False, None
if self.m.lines is None and self.mref.lines is not None:
return False, None
mFluxes = np.array([x.flux for x in self.m.lines])
mrefFluxes = np.array([x.flux for x in self.mref.lines])
# if they have different shapes it will not work: catch this here
# FIXME: not nice but should not happen to much, and this way it at least does not crash
if mFluxes.shape != mrefFluxes.shape:
return False, None
return self.diffArray(mFluxes, mrefFluxes, **self.getTols(rtol=0.05))
[docs]
def compareLineEstimates(self):
"""
Compares the fluxes from the line estimates
Currently assumes that both models include the same lines in the same order.
.. todo::
* it might be worth to make this smarter and compare the individual line estimates.
currently the comparison fails if e.g. there are more level for one line ... but the rest is fine
"""
if self.m.lineEstimates is None and self.mref.lineEstimates is None:
return True, None
if self.m.lineEstimates is not None and self.mref.lineEstimates is None:
return False, None
if self.m.lineEstimates is None and self.mref.lineEstimates is not None:
return False, None
mFluxes = np.array([x.flux for x in self.m.lineEstimates])
mrefFluxes = np.array([x.flux for x in self.mref.lineEstimates])
# if they have different shapes it will not work: catch this here
# FIXME: not nice but should not happen to much, and this way it at least does not crash
if mFluxes.shape != mrefFluxes.shape:
return False, None
mask = np.logical_and(mFluxes < 1.0e-24, mrefFluxes < 1.0e-24)
return self.diffArray(
mFluxes, mrefFluxes, mask=mask, **self.getTols(rtol=3.0e-1, nlimfrac=0.01)
)
[docs]
def compareSED(self):
"""
Compares the SEDs
"""
if self.m.sed is None and self.mref.sed is None:
return True, None
if self.m.sed is not None and self.mref.sed is None:
return False, None
if self.m.sed is None and self.mref.sed is not None:
return False, None
return self.diffArray(self.m.sed.fnuErg, self.mref.sed.fnuErg, **self.getTols(rtol=1.0e-3))
[docs]
def compareContinuumImages(self):
"""
Compares some of the continuum images.
"""
if self.m.contImages is None and self.mref.contImages is None:
return True, None
if self.m.contImages is not None and self.mref.contImages is None:
return False, None
if self.m.contImages is None and self.mref.contImages is not None:
return False, None
for wl in [1, 10, 100, 1000]:
x, y, imm, immwl = self.m.contImages.getImage(wl)
xref, yref, imref, imrefwl = self.mref.contImages.getImage(wl)
f, d = self.diff(immwl, imrefwl)
if not f:
return False, d
f, d, b = self.diffArray(imm, imref)
if not f:
return False, d, b
f, d, b = self.diffArray(x, xref)
if not f:
return False, d, b
f, d, b = self.diffArray(y, yref)
if not f:
return False, d, b
return True, d, b
[docs]
def compareCdnmol(self):
"""
checks the vertical column densities
only the outermost points (e.g. total column density) are checked
"""
# specidxM=[self.m.spnames[spname] for spname in self.specCompare if spname in self.m.spnames]
# specidxMref=[self.mref.spnames[spname] for spname in self.specCompare if spname in self.mref.spnames]
# get rid of very small column densities
# for i in specidxM:
# self.m.cdnmol[self.m.cdnmol[:,0,i]<self.lcdnmol,0,i]=self.lcdnmol
#
# for i in specidxMref:
# self.mref.cdnmol[self.mref.cdnmol[:,0,i]<self.lcdnmol,0,i]=self.lcdnmol
return self.diffArray(
self.m.cdnmol[:, 0, :],
self.mref.cdnmol[:, 0, :],
**self.getTols(rtol=1.0e-1, nlimfrac=0.1),
)
# ok,d,b=self.diffArray(self.m.cdnmol[:,0,specidxM],self.mref.cdnmol[:,0,specidxMref],1.e-1,1.0,0.1)
# if false check if it is only a certain fraction of the columns, it than can
# be still okay
# is not really elegant I would say
# TODO: also somehow return the number of failed columns
# TODO: maybe merge rcdnmol
# if ok is False and d is not None:
# ok=True # and check if any columns faild
# for i in range(len(specidxMref)):
# faildcolumns=(d[:,i]>self.dcdnmol).sum()
# if ((float(faildcolumns)/float(len(d[:,i])))>self.fcdnmol):
# ok=False
# specs=list(self.mref.spnames.items())
# print(faildcolumns,specs[specidxM[i]])
# return ok,d,b
[docs]
def compareRcdnmol(self):
"""
Checks the radial column densities.
Only the outermost points are checked (i.e. the total radial column density)
"""
# specidxM=[self.m.spnames[spname] for spname in self.specCompare if spname in self.m.spnames]
# print(specidxM)
# specidxMref=[self.mref.spnames[spname] for spname in self.specCompare if spname in self.mref.spnames]
# for i in specidxM:
# self.m.rcdnmol[-1,self.m.rcdnmol[-1,:,i]<self.lcdnmol,i]=self.lcdnmol
#
# for i in specidxMref:
# self.mref.rcdnmol[-1,self.mref.rcdnmol[-1,:,i]<self.lcdnmol,i]=self.lcdnmol
return self.diffArray(
self.m.rcdnmol[-1, :, :],
self.mref.rcdnmol[-1, :, :],
**self.getTols(rtol=1.0e-1, nlimfrac=0.1),
)
# ok,d,b=self.diffArray(self.m.rcdnmol[-1,:,specidxM],self.mref.rcdnmol[-1,:,specidxMref],1.e-1,1.0,0.1)
# if false check if it is only a certain fraction of the columns, it than can
# be still okay
# is not really elegant I would say
# TODO: also somehow return the number of failed columns
# TODO: maybe merge Cdnmol
# if ok is False:
# ok=True # and check if any columns faild
# for i in range(len(specidxMref)):
# faildcolumns=(d[:,i]>self.dcdnmol).sum()
# if ((float(faildcolumns)/float(len(d[:,i])))>self.fcdnmol):
# print(faildcolumns,len(d[:,i]),list(self.mref.spnames)[i])
# ok=False
# return ok,d,b
[docs]
def compareSolved_chem(self):
"""
Compares the solved_chem field. The flag for the quality of the chemical solution.
"""
return self.diffArray(
self.m.solved_chem, self.mref.solved_chem, **self.getTols(rtol=1.0e-10, nlimfrac=0.1)
)
[docs]
def compareAbundances(self):
"""
Compare the abundances.
"""
if self.m.nspec != self.mref.nspec:
return False, None
# set low values to zero
# self.m.nmol[self.m.nmol<self.lAbundances]=self.lAbundances
# self.mref.nmol[self.mref.nmol<self.lAbundances]=self.lAbundances
# print(self.m.nmol[:,0,0])
# print(self.mref.nmol[:,0,0])
# only check a small selection of species
# spec=("H2","CO","H2O","CO#","H2O#","H3+")
specs = self.mref.spnames
# want to compare the abundances
amref = self.mref.nmol[:, :, :] * 0.0
am = self.m.nmol[:, :, :] * 0.0
for i in range(len(specs)):
amref[:, :, i] = self.mref.nmol[:, :, i] / self.mref.nHtot
am[:, :, i] = self.m.nmol[:, :, i] / self.m.nHtot
# specidxM=[self.m.spnames[idx] for idx in specs if idx in self.m.spnames]
# print(specidxM)
# specidxMref=[self.mref.spnames[idx] for idx in specs if idx in self.mref.spnames]
# print(specidxMref)
# self.mref.solved_chem!=1 and
mask = np.logical_and(am < 1.0e-20, amref < 1.0e-20)
masknotsolved = np.logical_or(self.m.solved_chem != 1, self.mref.solved_chem != 1)
# do not know a better ways
for i in range(len(specs)):
mask[:, :, i] = np.logical_or(mask[:, :, i], masknotsolved)
return self.diffArray(am, amref, mask=mask, **self.getTols(rtol=1.0e-3, nlimfrac=0.05))
[docs]
def compareElements(self):
"""
Compares the numbers from the Elemenst.out file.
"""
ok, val = self.diff(self.m.elements.muHamu, self.mref.elements.muHamu, rtol=1.0e-5)
if not ok:
return ok, val
ok, val, b = self.diffArray(
np.array(list(self.m.elements.amu.values())),
np.array(list(self.mref.elements.amu.values())),
rtol=1.0e-5,
)
if not ok:
return ok, val, b
ok, val, b = self.diffArray(
np.array(list(self.m.elements.massRatio.values())),
np.array(list(self.mref.elements.massRatio.values())),
rtol=1.0e-5,
)
if not ok:
return ok, val, b
return self.diffArray(
np.array(list(self.m.elements.abun12.values())),
np.array(list(self.mref.elements.abun12.values())),
rtol=1.0e-5,
)
[docs]
def compareSpecies(self):
"""
Compares the numbers from the Species.out file.
"""
ok, val, b = self.diffArray(
np.array(list(self.m.species.mass.values())),
np.array(list(self.mref.species.mass.values())),
rtol=1.0e-5,
)
if not ok:
return ok, val, b
ok, val, b = self.diffArray(
np.array(list(self.m.species.charge.values())),
np.array(list(self.mref.species.charge.values())),
rtol=1.0e-5,
)
if not ok:
return ok, val, b
return self.diffArray(
np.array(list(self.m.species.chemPot.values())),
np.array(list(self.mref.species.chemPot.values())),
rtol=1.0e-5,
)
[docs]
def compareDustOpacEnv(self):
"""
Compares the dust opacities for an envelope model (optional output).
"""
if self.m.env_dust is None and self.mref.env_dust is None:
return True, None
if self.m.env_dust is not None and self.mref.env_dust is None:
return False, None
if self.m.env_dust is None and self.mref.env_dust is not None:
return False, None
return self.diffArray(
self.m.env_dust.kext, self.mref.env_dust.kext, **self.getTols(rtol=1.0e-3)
)
[docs]
def compareDustOpac(self):
"""
Compares the dust opacities.
"""
if self.m.dust.kext.shape != self.mref.dust.kext.shape:
return False, None
# self.tolerances["DustOpac.rtol"]=0.06
return self.diffArray(self.m.dust.kext, self.mref.dust.kext, **self.getTols(rtol=1.0e-3))
[docs]
def compareDust_charge_mean(self):
"""
Compares the dust charge.
"""
if self.m.dust_charge_mean is None and self.mref.dust_charge_mean is None:
return True, None
if self.m.dust_charge_mean is None and self.mref.dust_charge_mean is not None:
return False, None
if self.m.dust_charge_mean is not None and self.mref.dust_charge_mean is None:
return False, None
if self.m.dust_charge_mean.shape != self.mref.dust_charge_mean.shape:
return False, None
return self.diffArray(self.m.dust_charge_mean, self.mref.dust_charge_mean)
[docs]
def compareStarSpec(self):
"""
Compares the input Stellar spectrum, from X-rays to mm
"""
return self.diffArray(self.m.star.Inu, self.mref.star.Inu)
[docs]
def compareNHtot(self):
"""
checks the total hydrogen number density
"""
return self.diffArray(self.m.nHtot, self.mref.nHtot)
[docs]
def compareNd(self):
"""
Checks the (total) dust number density
"""
return self.diffArray(self.m.nd, self.mref.nd)
[docs]
def compareRhog(self):
"""
checks the gas density
"""
return self.diffArray(self.m.rhog, self.mref.rhog)
[docs]
def compareVelocity(self):
"""
checks the velocity structure.
If in one of the models velocity is None it is assumed it is an old model
and the comparison returns Tru
"""
if (self.m.velocity is None) ^ (self.mref.velocity is None):
print(
"WARN: Ignore velocity comparison as it does not exist in one of the models. ",
end="",
)
return True, None, None
# if both are None it is also fine
if (self.m.velocity is None) and (self.mref.velocity is None):
return True, None, None
return self.diffArray(self.m.velocity, self.mref.velocity, rtol=1.0e-5)
[docs]
def compareRhod(self):
"""
checks the dust density
"""
return self.diffArray(self.m.rhod, self.mref.rhod)
[docs]
def compareD2g(self):
"""
checks the resulting dust to gass mass ratio
"""
return self.diffArray(self.m.d2g, self.mref.d2g)
[docs]
def compareDamean(self):
"""
checks the resulting mean grain size
"""
return self.diffArray(self.m.damean, self.mref.damean)
[docs]
def compareFsize_Rho(self):
"""
checks the dust grain size distribution, if it exists.
"""
return self.diffArray(self.m.dust.fsize_rho, self.mref.dust.fsize_rho)
[docs]
def compareTg(self):
"""
checks the gas temperature
"""
return self.diffArray(self.m.tg, self.mref.tg, **self.getTols(rtol=2.0e-2, nlimfrac=0.01))
[docs]
def compareTd(self):
"""
checks the dust temperature
"""
return self.diffArray(self.m.td, self.mref.td, **self.getTols(rtol=2.0e-3, nlimfrac=0.01))
[docs]
def compareRadFields(self):
"""
checks the radiation fields
"""
return self.diffArray(
self.m.radFields, self.mref.radFields, **self.getTols(rtol=2.0e-2, nlimfrac=0.01)
)
[docs]
def compareHeat(self):
"""
checks the heating rates
"""
return self.diffArray(
self.m.heat, self.mref.heat, **self.getTols(rtol=1.0e-3, nlimfrac=0.01)
)
[docs]
def compareCool(self):
"""
checks the cooling rates
"""
return self.diffArray(
self.m.cool, self.mref.cool, **self.getTols(rtol=1.0e-3, nlimfrac=0.01)
)
[docs]
def compareZetaCR(self):
"""
checks the cosmic ray ionisation rate
"""
return self.diffArray(self.m.zetaCR, self.mref.zetaCR, **self.getTols(rtol=1.0e-5))
[docs]
def compareZetaX(self):
"""
checks the Xray ionisation rate
"""
mask = np.logical_and(self.m.zetaX < 1.0e-25, self.mref.zetaX < 1.0e-25)
return self.diffArray(
self.m.zetaX, self.mref.zetaX, mask=mask, **self.getTols(rtol=1.0e-1, nlimfrac=0.05)
)
[docs]
def compareHX(self):
"""
checks the Xray energy deposition rate
"""
# set low values to zero
# self.m.Hx[self.m.zetaX < self.lZetaX]=self.lZetaX
# self.mref.Hx[self.mref.zetaX < self.lZetaX]=self.lZetaX
return self.diffArray(self.m.Hx, self.mref.Hx, **self.getTols(rtol=1.0e-1, nlimfrac=0.05))
[docs]
def compareFLiTsSpec(self):
"""
Compares the FLits spectra if they are available
"""
# As a flits spec is optional cover the None cases
if self.m.FLiTsSpec is None and self.mref.FLiTsSpec is None:
return True, None
if self.m.FLiTsSpec is not None and self.mref.FLiTsSpec is None:
return False, None
if self.m.FLiTsSpec is None and self.mref.FLiTsSpec is not None:
return False, None
# check if the number of wl points is the same, otherwise we cannot compare
if self.m.FLiTsSpec.wl.shape != self.mref.FLiTsSpec.wl.shape:
return False, None
# first check the wavelength, if this differs, comparing fluxes makes no sense
eq, da, dbool = self.diffArray(self.m.FLiTsSpec.wl, self.mref.FLiTsSpec.wl)
if not eq:
return eq, da, dbool
# now compare the fluxes, first the continuum
eq, da, dbool = self.diffArray(
self.m.FLiTsSpec.flux_cont, self.mref.FLiTsSpec.flux_cont, **self.getTols(rtol=1.0e-5)
)
if not eq:
return eq, da, dbool
# now compare the fluxes
return self.diffArray(
self.m.FLiTsSpec.flux, self.mref.FLiTsSpec.flux, **self.getTols(rtol=1.0e-5)
)
[docs]
def compareChemnet(self):
"""
Compare the chemical network. Uses the compare routine from :class:`~prodimopy.chemistry.network.ReactionNetwork`
"""
ret = self.mref.chemnet.compare(self.m.chemnet, printresults=False, log=False)
return (ret[0], [len(ret[1]), len(ret[2]), len(ret[3])], None)
[docs]
def compareFileCheckNetworkLog(self):
"""
Makes a file comparison with CheckNetwork.log
"""
return self.diffFile("CheckNetwork.log")
[docs]
def compareFileCheckChemLog(self):
"""
Makes a file comparison with CheckChem.log
"""
return self.diffFile("CheckChem.log")
[docs]
def compareFitsFileRestart(self):
"""
Makes a fits file comparison with restart.fits.gz
"""
return self.diffFitsFile("restart.fits.gz", rtol=1.0e-3, atol=1.0e-40)
[docs]
def compareFitsFileMie(self):
"""
Makes a fits file comparison with Mie.fits.gz
"""
return self.diffFitsFile("Mie.fits.gz", rtol=1.0e-5)
[docs]
def compareFitsFileFLiTs(self):
"""
Makes a fits file comparison with ProDiMoForFLits.fits
"""
return self.diffFitsFile("ProDiMoForFLiTs.fits", rtol=1.0e-3, atol=1.0e-40)
[docs]
def compareFitsFileLineCubes(self):
"""
Makes a fits file comparison with for all line cubes (LINE_3D_*)
"""
fcubes = glob.glob(self.m.directory + "/LINE_3D_???.fits")
# FIXME: rtol is pretty high
# FIXME: this deos not work if the other model actually has linecubes
# e.g. it would say that it is okay
if fcubes is not None:
for fname in fcubes:
ok, val = self.diffFitsFile(os.path.basename(fname), rtol=1.0e-3)
if not ok: # stop the whole thing if one cube went wrong
return ok, val
return True, None
[docs]
class CompareMc(CompareAbs):
"""
Class for comparing two |prodimo| models of type :class:`~prodimopy.read_mc.Data_mc`
Every compare Function returns true or false, and the relative differences
(in case of arrays these are arrays).
Can be used in e.g. automatic testing routines or in simple command line
tools to compare |prodimo| model results.
"""
def __init__(self, model, modelref):
"""
Parameters
----------
model : :class:`~prodimopy.read_mc.Data_mc`
The model one wants to compare.
modelref : :class:`~prodimopy.read_mc.Data_mc`
The reference model to compare with.
Attributes
----------
"""
if model is None or modelref is None:
raise TypeError("Either model or modelref is None, no comparison possible.")
self.m = model
""" :class:`~prodimopy.read_mc.Data_mc` :
The model one wants to compare.
"""
self.mref = modelref
""" :class:`~prodimopy.read_mc.Data_mc` :
The reference model to compare with.
"""
[docs]
def compareAbundances(self):
"""
Compares the abundances of two molecular cloud (0D chemistry) models.
Assumes that both models used the same number of ages and species in the same order.
"""
# set low values to equal numbers, which means they are kind of ignored
# in the comparison
# self.m.abundances[self.m.abundances<self.lAbundances]=self.lAbundances
# self.mref.abundances[self.mref.abundances<self.lAbundances]=self.lAbundances
# Do not consider the first age entry at it is the initial conditions
# that can vary from model two model and are not really a result
mask = np.logical_and(
self.m.abundances[1:, :] < 1.0e-20, self.mref.abundances[1:, :] < 1.0e-20
)
return self.diffArray(
self.m.abundances[1:, :], self.mref.abundances[1:, :], rtol=5.0e-2, mask=mask
)
[docs]
def compareRatecoefficients(self):
"""
Compares the rate coefficients of two molecular cloud (0D chemistry) models.
Assumes that both models have exactly the same chemical reactions in the same order.
"""
if len(self.m.ratecoefficients) != len(self.mref.ratecoefficients):
return False, None, None
mask = np.logical_and(
self.m.ratecoefficients < 1.0e-50, self.mref.ratecoefficients < 1.0e-50
)
return self.diffArray(
self.m.ratecoefficients, self.mref.ratecoefficients, rtol=2.0e-7, mask=mask
)
[docs]
class CompareSlab1D(CompareAbs):
"""
Class for comparing two |prodimo| SLAB models of type :class:`~prodimopy.read_slab.slab1D`
Every compare Function returns true or false, and the relative differences
(in case of arrays these are arrays).
Can be used in e.g. automatic testing routines or in simple command line
tools to compare |prodimo| model results.
"""
def __init__(self, model: pread_slab.slab1D, modelref: pread_slab.slab1D):
"""
Parameters
----------
model :
The model one wants to compare.
modelref :
The reference model to compare with.
Attributes
----------
"""
if model is None or modelref is None:
raise TypeError("Either model or modelref is None, no comparison possible.")
self.m = model
""" The model one wants to compare. """
self.mref = modelref
""" The reference model to compare with. """
[docs]
def compareFlux(self):
"""
Compares the fluxes of the slab1D models.
"""
if self.m.flux is None and self.mref.flux is None:
return True, None
if self.m.flux is not None and self.mref.flux is None:
return False, None
if self.m.flux is None and self.mref.flux is not None:
return False, None
return self.diffArray(self.m.flux, self.mref.flux, rtol=1.0e-5)
[docs]
def compareSource_function(self):
"""
Compares the source functions of the slab1D models.
"""
if self.m.source_function is None and self.mref.source_function is None:
return True, None
if self.m.source_function is not None and self.mref.source_function is None:
return False, None
if self.m.source_function is None and self.mref.source_function is not None:
return False, None
return self.diffArray(self.m.source_function, self.mref.source_function, rtol=1.0e-5)
[docs]
class CompareSlab(CompareAbs):
"""
Class for comparing two |prodimo| SLAB models of type :class:`~prodimopy.read_slab.slab_data`
Every compare Function returns true or false, and the relative differences
(in case of arrays these are arrays).
Can be used in e.g. automatic testing routines or in simple command line
tools to compare |prodimo| model results.
"""
def __init__(self, model, modelref):
"""
Parameters
----------
model : :class:`~prodimopy.read_slab.slab_data`
The model one wants to compare.
modelref : :class:`~prodimopy.read_slab.slab_data`
The reference model to compare with.
Attributes
----------
"""
if model is None or modelref is None:
raise TypeError("Either model or modelref is None, no comparison possible.")
self.m = model
""" :class:`~prodimopy.read_slab.slab_data` :
The model one wants to compare.
"""
self.mref = modelref
""" :class:`~prodimopy.read_slab.slab_data` :
The reference model to compare with.
"""
[docs]
def compareLbLFluxes(self):
"""
Compares the FNLTE and FLTE of the linedata for all the models in line-by-line mode
If number of models is not the same return false.
Assumes that both models used the same number of species in the same order.
"""
if len(self.mref.models) != len(self.m.models):
return False, None, None
for m, mref in zip(self.m.models, self.mref.models):
arrm = np.array(m.linedata[["FNLTE", "FLTE"]])
arrmref = np.array(mref.linedata[["FNLTE", "FLTE"]])
if (arrm.shape != arrmref.shape):
print("")
print(m)
print("")
print("Reference Model: ")
print(mref)
raise ValueError("Shapes of the line fluxes differ, cannot compare.")
eq, da, dbool = self.diffArray(arrm, arrmref, rtol=1.0e-4)
# Current stop at the first error not great
if not eq:
return eq, da, dbool
# That means all good
return eq, da, dbool
[docs]
def compareLbLTau(self):
"""
Compares the FNLTE and FLTE of the linedata for all the models in line-by-line mode
If number of models is not the same return false.
Assumes that both models used the same number of species in the same order.
"""
if len(self.mref.models) != len(self.m.models):
return False, None, None
for m, mref in zip(self.m.models, self.mref.models):
arrm = np.array(m.linedata[["tauNLTE", "tauLTE"]])
arrmref = np.array(mref.linedata[["tauNLTE", "tauLTE"]])
eq, da, dbool = self.diffArray(arrm, arrmref, rtol=1.0e-4)
# Current stop at the first error not great
if not eq:
return eq, da, dbool
# That means all got
return eq, da, dbool
[docs]
def compareLbLPop(self):
"""
Compares the LTE and NLTE level populations of the linedata for all the models in line-by-line mode
If number of models is not the same return false.
Assumes that both models used the same number of species in the same order.
"""
if len(self.mref.models) != len(self.m.models):
return False, None, None
for m, mref in zip(self.m.models, self.mref.models):
arrm = np.array(m.linedata[["pop", "ltepop"]])
arrmref = np.array(mref.linedata[["pop", "ltepop"]])
eq, da, dbool = self.diffArray(arrm, arrmref, rtol=1.0e-4)
# Current stop at the first error not great
if not eq:
return eq, da, dbool
for m, mref in zip(self.m.models, self.mref.models):
arrm = np.array(m.leveldata[["pop", "ltepop"]])
arrmref = np.array(mref.leveldata[["pop", "ltepop"]])
eq, da, dbool = self.diffArray(arrm, arrmref, rtol=1.0e-4)
# Current stop at the first error not great
if not eq:
return eq, da, dbool
# That means all got
return eq, da, dbool
# def compareOverlapFluxes(self):
# """
# Compares the FNLTE and FLTE of the linedata for all the models in line overlap mode
# If number of models is not the same return false.
#
#
# Assumes that both models used the same number of species in the same order.
# """
#
# if len(self.mref.models)!=len(self.m.models):
# return False,None,None
#
# # for m,mref in zip(self.m.models,self.mref.models):
# # arrm=np.array(m.overlapLTE)
# # arrmref=np.array(mref.overlapLTE)
# # print(arrm)
# # print(arrmref)
# #
# # eq,da,dbool=self.diffArray(arrm,arrmref,rtol=1.e-4)
# #
# # # Current stop at the first error not great
# # if not eq: return eq,da,dbool
#
# for m,mref in zip(self.m.models,self.mref.models):
# arrm=np.array(m.overlapNLTE)
# arrmref=np.array(mref.overlapNLTE)
#
# eq,da,dbool=self.diffArray(arrm,arrmref,rtol=1.e-4)
#
# # Current stop at the first error not great
# if not eq: return eq,da,dbool
#
# # That means all got
# return eq,da,dbool
# def compareOverlapTau(self):
# """
# Compares the LTE and NLTE optical depths for all the models in line overlap mode
# If number of models is not the same return false.
#
#
# Assumes that both models used the same number of species in the same order.
# """
#
# if len(self.mref.models)!=len(self.m.models):
# return False,None,None
#
# for m,mref in zip(self.m.models,self.mref.models):
# arrm=np.array(m.overlapTauLTE)
# arrmref=np.array(mref.overlapTauLTE)
#
# eq,da,dbool=self.diffArray(arrm,arrmref,rtol=1.e-4)
#
# # Current stop at the first error not great
# if not eq: return eq,da,dbool
#
# for m,mref in zip(self.m.models,self.mref.models):
# arrm=np.array(m.overlapTauNLTE)
# arrmref=np.array(mref.overlapTauNLTE)
#
# eq,da,dbool=self.diffArray(arrm,arrmref,rtol=1.e-4)
#
# # Current stop at the first error not great
# if not eq: return eq,da,dbool
#
# # That means all got
# return eq,da,dbool
# def compareOverlapFreq(self):
# """
# Compares the frequency grid of all the models in line overlap mode
# If number of models is not the same return false.
#
#
# Assumes that both models used the same number of species in the same order.
# """
#
# if len(self.mref.models)!=len(self.m.models):
# return False,None,None
#
# for m,mref in zip(self.m.models,self.mref.models):
# arrm=np.array(m.overlapFreq)
# arrmref=np.array(mref.overlapFreq)
#
# eq,da,dbool=self.diffArray(arrm,arrmref,rtol=1.e-4)
#
# # Current stop at the first error not great
# if not eq: return eq,da,dbool
#
# # That means all got
# return eq,da,dbool
[docs]
def eval_model_type(modelDir):
"""
Try to guess the model type (e.g. mc, full prodimo etc.).
Default is always prodimo.
Possible types:
`prodimo` .... full prodimo model (:class:`prodimopy.read.Data_ProDiMo`)
`mc` ......... molecular cloud model (:class:`prodimopy.read_mc.Data_mc`)
Returns
-------
str either `prodimo`, `prodimoTD`, `slab`, `slab1D` or `mc`
FIXME: this is maybe not the best place for this routine
FIXME: provide constants for the values (what's the best way to do this in pyhton?)
FIXME: this is just a quick hack, it would be better to use the parameters in Parameter.in
"""
if (
os.path.isfile(modelDir + "/ProDiMo_01.out")
or os.path.isfile(modelDir + "/ProDiMo_001.out")
or os.path.isfile(modelDir + "/ProDiMo_0001.out")
or os.path.isfile(modelDir + "/ProDiMo_00001.out")
):
mtype = "prodimoTD"
elif os.path.isfile(modelDir + "/Molecular_Cloud_Input.in"):
mtype = "mc"
elif os.path.isfile(os.path.join(modelDir, "SlabInput.in")):
mtype = "slab"
# maybe a bit of an overkill,
f=open(os.path.join(modelDir, "SlabInput.in"),"r")
for line in f:
if "! slab_1D" in line.strip():
val=line.split("! slab_1D")[0].strip().lower()
if (val==".true."or val=="true"):
mtype="slab1D"
break
else:
mtype = "prodimo"
return mtype
# def compareFlineEstimates(self):
# '''
# Compares the FlineEstimaes
# '''
#
# if len(self.m.lineEstimates) !=len(self.mref.lineEstimates):
# return False,None
#
#
# self.diff(self.m.lineEstimate[i].flux, self.mref.lineEstimate[i].flux, self.dLineFluxes)
#
# # Compare fluxes
# for i in range(len(self.m.lines)):
# f,d=
# if f == False:
# return False,d
#
# return True,None
#