import cProfile
import io
import pstats
import warnings
from concurrent.futures import ProcessPoolExecutor, as_completed
from copy import deepcopy
import astropy.stats as astats
# import astropy.wcs as wcs
import astropy.units as u
import numpy as np
import numpy.ma as ma
from astropy import log
# from astropy.io import fits
from astropy.io.fits.header import Header
from astropy.nddata import CCDData
from astropy.table import Column, Table
from lmfit import Minimizer, Parameters # , fit_report
from scipy.interpolate import interpn
from scipy.optimize import least_squares as _scipy_least_squares
from pdrtpy.pbar import get_progress_bar
from .. import utils
from .fitmap import FitMap
from .toolbase import ToolBase
log.setLevel("WARNING") # see issue 163
# ---------------------------------------------------------------------------
# Module-level helpers for parallel pixel fitting.
# Must be at module level so ProcessPoolExecutor can pickle them.
# ---------------------------------------------------------------------------
# Per-worker-process cache of RegularGridInterpolators. Set once by _init_worker.
_worker_model_interps = None
def _init_worker(model_points, model_values):
"""Initializer for ProcessPoolExecutor workers. Builds RegularGridInterpolators
from plain numpy arrays once per worker process so they are not rebuilt for
every pixel."""
from scipy.interpolate import RegularGridInterpolator
global _worker_model_interps
_worker_model_interps = [
RegularGridInterpolator(pts, vals, method="linear", bounds_error=True)
for pts, vals in zip(model_points, model_values, strict=False)
]
def _fit_pixel_worker(
j, obs_data_j, obs_err_j, init_density, init_rf, minn, maxn, minfuv, maxfuv, nan_policy, minimize_kwargs
):
"""Fit a single spatial pixel. Runs in a worker process spawned by
ProcessPoolExecutor. Uses the interpolators set up by _init_worker.
Returns (j, MinimizerResult)."""
import numpy as _np
from lmfit import Minimizer as _Minimizer, Parameters as _Parameters
def _residual(params):
parvals = params.valuesdict()
d = parvals["density"]
rf = parvals["radiation_field"]
mvalue = _np.array([float(interp((d, rf))) for interp in _worker_model_interps])
return (obs_data_j - mvalue) / obs_err_j
params = _Parameters()
params.add("density", min=minn, max=maxn, value=init_density)
params.add("radiation_field", min=minfuv, max=maxfuv, value=init_rf)
minimizer = _Minimizer(_residual, params=None, nan_policy=nan_policy)
return j, minimizer.minimize(params=params, **minimize_kwargs)
# ---------------------------------------------------------------------------
# Per-pixel proxy for joint fitting results.
# Presents the same interface as lmfit.MinimizerResult so that
# fit_result[j].params, fit_result[j].chisqr, etc. work unchanged.
# ---------------------------------------------------------------------------
class _PixelResult:
"""Thin per-pixel view into a joint MinimizerResult.
Returned by fit_result[j] when joint_fit='hybrid' or 'fast' was used. Exposes the
same attributes as a full lmfit.MinimizerResult for the two fitted
parameters (density, radiation_field) of a single pixel."""
def __init__(self, density_val, density_err, rf_val, rf_err, chisqr, redchi, ndata, success, residual):
self.params = Parameters()
self.params.add("density", value=density_val)
self.params["density"].stderr = density_err
self.params.add("radiation_field", value=rf_val)
self.params["radiation_field"].stderr = rf_err
self.chisqr = chisqr
self.redchi = redchi
self.success = success
self.residual = residual
self.method = "least_squares"
self.nvarys = 2
self.ndata = ndata
self.nfree = max(ndata - 2, 1)
self.nfev = None # not meaningful per-pixel in a joint fit
self.errorbars = density_err is not None and rf_err is not None
[docs]
class LineRatioFit(ToolBase):
"""Tool to fit observations of intensity ratios to a set of PDR models.
Takes as input a set of observations with errors represented as
:class:`~pdrtpy.measurement.Measurement` and a
:class:`~pdrtpy.modelset.ModelSet` for the models to fit. Observations
should be spectral line or continuum intensities, either spatial maps or
single pixel values at the same spatial resolution.
At least 3 observations are needed to make at least 2 ratios. Once the
fit is done, :class:`~pdrtpy.plot.LineRatioPlot` can be used to view the results.
Parameters
----------
modelset : :class:`~pdrtpy.modelset.ModelSet`
The set of PDR models to use for fitting.
measurements : list or dict of :class:`~pdrtpy.measurement.Measurement`, optional
Input measurements to be fit. If dict, keys should be Measurement identifiers.
"""
def __init__(self, modelset, measurements=None):
super().__init__() # needed?
if isinstance(modelset, str):
# may need to disable this
self._initialize_modelTable(modelset)
self._modelset = modelset
self._init_measurements(measurements)
self._set_measurementnaxis()
self._modelratios = None
self._modelnaxis = None
self._set_model_files_used()
self._observedratios = None
self._chisq = None
self._reduced_chisq = None
self._likelihood = None
self._radiation_field = None
self._density = None
self.radiation_field_unit = None
self.radiation_field_type = None
self.density_unit = None
self.density_type = None
self._fitresult = None
self._fitparam = None
self._minimizer = None
self._deltasq = None
self._ratiocount = None
# phase-space limits (issue #236): (lower, upper) in linear model units,
# populated by _normalize_phase_space_limits. Default None = unrestricted.
self._density_range = None
self._radiation_field_range = None
@property
def fit_result(self):
"""The result of the fitting procedure, including fit statistics, variable values and uncertainties, and correlations. One :class:`lmfit.minimizer.MinimizerResult` per pixel.
Returns
-------
:class:`~pdrtpy.tool.FitMap`
"""
return self._fitresult
@property
def modelset(self):
"""The underlying :class:`~pdrtpy.modelset.ModelSet`"""
return self._modelset
@property
def measurements(self):
"""The stored measurements as a dictionary with Measurement IDs as keys.
Returns
-------
dict of :class:`~pdrtpy.measurement.Measurement`
"""
return self._measurements
@property
def measurementIDs(self):
"""The stored measurement IDs, which are strings.
Returns
-------
:class:`dict_keys`
"""
if self._measurements is None:
return None
return self._measurements.keys()
@property
def observed_ratios(self):
"""The list of the observed line ratios that have been input so far.
Returns
-------
list of str
"""
return list(self._observedratios.keys())
@property
def ratiocount(self):
"""The number of ratios that match models available in the current :class:`~pdrtpy.modelset.ModelSet` given the current set of measurements.
Returns
-------
int
"""
# call to modelset._get_ratio_MA 01938elements is expensive. So set it once for each run.
if self._ratiocount is None:
self._ratiocount = self._modelset.ratiocount(self.measurementIDs)
return self._ratiocount
@property
def density(self):
"""The computed hydrogen nucleus density value(s).
Returns
-------
:class:`~pdrtpy.measurement.Measurement`
"""
return self._density
@property
def radiation_field(self):
"""The computed radiation field value(s).
Returns
-------
:class:`~pdrtpy.measurement.Measurement`
"""
return self._radiation_field
[docs]
def chisq(self, min=False):
r"""The computed chisquare value(s).
Parameters
----------
min : bool, optional
If True, return the minimum :math:`\chi^2`. For map inputs, returns a
spatial map of the minimum; for single-pixel inputs, returns a scalar.
If False, returns the full :math:`\chi^2` as a function of density and
radiation field. Default: False.
Returns
-------
:class:`~pdrtpy.measurement.Measurement`
"""
if min:
return self._chisq_min
else:
return self._chisq
[docs]
def reduced_chisq(self, min=False):
r"""The computed reduced chisquare value(s).
Parameters
----------
min : bool, optional
If True, return the minimum :math:`\chi_\nu^2`. For map inputs, returns
a spatial map of the minimum; for single-pixel inputs, returns a scalar.
If False, returns the full :math:`\chi_\nu^2` as a function of density and
radiation field. Default: False.
Returns
-------
:class:`~pdrtpy.measurement.Measurement`
"""
if min:
return self._reduced_chisq_min
else:
return self._reduced_chisq
def _init_measurements(self, m):
"""Initialize the measurements from an input list or dict.
Parameters
----------
m : list, tuple, or dict
The input Measurements. If dict, keys must be valid measurement identifiers.
"""
self._masks = dict() # need to save these so they can be reset later
if m is None:
self._measurements = None
elif isinstance(m, list) or isinstance(m, tuple):
self._measurements = dict()
for mm in m:
self._measurements[mm.id] = mm
self._masks[mm.id] = deepcopy(mm.mask)
elif isinstance(m, dict):
self._measurements = deepcopy(m)
for key in m:
self._masks[key] = deepcopy(m[key].mask)
else:
raise ValueError("Input measurements must be list, tuple, or dict")
def _set_model_files_used(self):
self._model_files_used = dict()
if self._measurements is None:
return
for x in self._modelset.find_files(self.measurementIDs):
self._model_files_used[x[0]] = x[1]
def _check_shapes(self, d):
# ugly
s1 = d[utils.firstkey(d)].shape
return np.all([m.shape == s1 for m in d.values()])
def _check_header(self, kw, value=NotImplemented):
"""Check to see if any of the given keyword values differ for the input measurements.
Parameters
----------
kw : str
The keyword to check.
value : any, optional
If given and not ``NotImplemented``, check that all *kw* values equal this.
The default is ``NotImplemented`` (not ``None``) so that ``None`` can be
checked against explicitly.
"""
d = self._measurements
fk = utils.firstkey(d)
in_wcs = False
try:
if value == NotImplemented:
h = d[fk].header
if kw not in h:
# CTYPES etc can also be in WCS so check that too
s1 = d[fk].wcs.to_header()[kw]
in_wcs = True
else:
s1 = h[kw]
else:
s1 = value
if in_wcs:
return np.all([m.wcs.to_header()[kw] == s1 for m in d.values()])
else:
return np.all([m.header[kw] == s1 for m in d.values()])
except KeyError:
print(f"WARNING: {kw} keyword not present in all Measurements")
return False
def _check_measurement_shapes(self):
if self._measurements is None:
return False
return self._check_shapes(self._measurements)
def _check_ratio_shapes(self):
if self._observedratios is None:
return False
return self._check_shapes(self._observedratios)
def _check_model_shapes(self):
if self._modelratios is None:
return False
return self._check_shapes(self._modelratios)
[docs]
def add_measurement(self, m):
r"""Add a Measurement to the internal dictionary used to compute ratios.
The measurement may be in intensity units (:math:`{\rm erg~s}^{-1}` :math:`{\rm cm}^{-2}`) or integrated intensity (K km/s).
Parameters
----------
m : :class:`~pdrtpy.measurement.Measurement`
A Measurement instance to be added to this tool.
"""
if self._measurements:
self._measurements[m.id] = m
else:
self._init_measurements(m)
self._set_model_files_used()
[docs]
def remove_measurement(self, id):
"""Delete a measurement from the internal dictionary used to compute ratios.
Parameters
----------
id : str
The measurement identifier.
Raises
------
KeyError
If id not in existing Measurements.
"""
del self._measurements[id]
self._set_model_files_used()
[docs]
def read_models(self, unit=u.dimensionless_unscaled):
"""Given a list of measurement IDs, find and open the FITS files with matching ratios and populate ``_modelratios``.
Uses :class:`pdrtpy.measurement.Measurement` as a storage mechanism.
Parameters
----------
unit : str or :class:`astropy.units.Unit`, optional
Units of the data.
"""
self._modelratios = self._modelset.get_models(self.measurementIDs, model_type="ratio")
if self.ratiocount < 2:
msg = (
"Not enough ratios. You need to provide at least 3 observations that can be used to compute 2 ratios"
f" that are covered by the ModelSet. From your observations, {self.ratiocount:d} ratio(s)"
)
if self.ratiocount > 0:
msg += f" {list(self._modelratios.keys())}"
msg += " can be computed."
raise Exception(msg)
k = utils.firstkey(self._modelratios)
self._modelnaxis = self._modelratios[k].wcs.naxis
self._modelshape = np.array(self._modelratios[k].data.shape)
if not self.density_unit:
self.density_unit = self._modelratios[k].wcs.wcs.cunit[0]
self.density_type = self._modelratios[k].wcs.wcs.ctype[0]
if not self.radiation_field_unit:
self.radiation_field_unit = self._modelratios[k].wcs.wcs.cunit[1]
self.radiation_field_type = self._modelratios[k].wcs.wcs.ctype[1]
# for wk2006 models, Habing units in the Y axis cause problems when trying to assign them
# to a WCS (see note in modelset.py). So get from header instead.
if self.radiation_field_unit.to_string() == "":
try:
self.radiation_field_unit = u.Unit(self._modelratios[k].header["CUNIT2"])
except KeyError as err:
raise Exception(
f"Keyword CUNIT2 is required in file {self._model_files_used[k]} FITS header to describe units"
" of interstellar radiation field"
) from err
if not self._check_model_shapes():
warnings.warn("Trimming all model grids to match H2 grid: log(n) = 1-5, log(G0) = 1-5", stacklevel=2)
utils._trim_all_to_H2(self._modelratios)
def _normalize_phase_space_limits(self, radiation_field_range, density_range):
"""Normalize and validate user-supplied phase-space limits (issue #236).
Converts ``radiation_field_range`` and ``density_range`` into
``(lower, upper)`` tuples of plain floats in the model's linear axis
units, substituting the model grid extremes for any ``None`` bound.
Results are stored in :attr:`_radiation_field_range` and
:attr:`_density_range` for use by the coarse grid search and the refine
step. Must be called after :meth:`read_models` so that
``density_unit`` / ``radiation_field_unit`` are set.
Parameters
----------
radiation_field_range : :class:`~astropy.units.Quantity`, sequence, or None
Allowed radiation field window; see :meth:`run`.
density_range : :class:`~astropy.units.Quantity`, sequence, or None
Allowed density window; see :meth:`run`.
Raises
------
ValueError
If a range is malformed, has incompatible units, or does not
overlap the model grid.
"""
fk = utils.firstkey(self._modelratios)
# linear (not log) physical axis values with units: x=density, y=radiation field
x, y = utils.get_xy_from_wcs(self._modelratios[fk], quantity=True, linear=True)
self._density_range = self._normalize_one_range(density_range, x, self.density_unit, "density_range")
self._radiation_field_range = self._normalize_one_range(
radiation_field_range, y, self.radiation_field_unit, "radiation_field_range"
)
@staticmethod
def _normalize_one_range(user_range, axis, unit, name):
"""Normalize a single phase-space limit to ``(lower, upper)`` floats.
Parameters
----------
user_range : :class:`~astropy.units.Quantity`, sequence, or None
The user-supplied range. A length-2 Quantity, or a length-2 sequence
whose elements are each ``None`` or a scalar Quantity. ``None`` (or a
``None`` element) means no limit on that side.
axis : :class:`~astropy.units.Quantity`
The model grid axis values (linear, with units) for this parameter.
unit : :class:`~astropy.units.Unit`
The model's linear axis unit to convert the bounds into.
name : str
Keyword name, used in error messages.
Returns
-------
tuple of float
``(lower, upper)`` in ``unit``, with grid extremes substituted for
unspecified bounds.
"""
grid_min = float(np.min(axis.to(unit).value))
grid_max = float(np.max(axis.to(unit).value))
if user_range is None:
return None # unrestricted
# Extract the two bounds, each a scalar Quantity or None.
if isinstance(user_range, u.Quantity):
if user_range.isscalar or len(user_range) != 2:
raise ValueError(f"{name} must have exactly two elements [lower, upper]")
bounds = [user_range[0], user_range[1]]
else:
try:
bounds = list(user_range)
except TypeError as err:
raise ValueError(f"{name} must be a length-2 Quantity or sequence [lower, upper]") from err
if len(bounds) != 2:
raise ValueError(f"{name} must have exactly two elements [lower, upper]")
limits = [grid_min, grid_max] # defaults substituted for None bounds
for i, b in enumerate(bounds):
if b is None:
continue
if not isinstance(b, u.Quantity):
raise ValueError(
f"{name} bounds must be astropy Quantities (with units); got {b!r}. For a one-sided limit put"
" the unit inside the brackets, e.g. [None, 1e4/u.cm**3]."
)
try:
limits[i] = float(b.to(unit).value)
except u.UnitConversionError as err:
raise ValueError(f"{name} bound {b} is not convertible to model units ({unit})") from err
lo, hi = limits
if lo > hi:
raise ValueError(f"{name} lower bound ({lo}) exceeds upper bound ({hi})")
if hi < grid_min or lo > grid_max:
raise ValueError(
f"{name} [{lo}, {hi}] {unit} does not overlap the model grid coverage [{grid_min}, {grid_max}] {unit}"
)
return (lo, hi)
def _check_compatibility(self):
"""Check that all Measurements are compatible (beams, coordinate systems, shapes) so that the computation can commence.
Raises
------
Exception
If headers and shapes don't match; warns if no beam present.
"""
if not self._check_measurement_shapes():
raise Exception("Your input Measurements have different dimensions")
# Check the beam sizes
# @Todo do the convolution ourselves if requested.
if not self._check_header("BMAJ"):
raise Exception(
"Beam major axis (BMAJ) of your input Measurements do not match. Please convolve all maps to the same"
" beam size"
)
if not self._check_header("BMIN"):
raise Exception(
"Beam minor axis (BMIN) of your input Measurements do not match. Please convolve all maps to the same"
" beam size"
)
if not self._check_header("BPA"):
raise Exception(
"Beam position angle (BPA) of your input Measurements do not match. Please convolve all maps to the"
" same beam size"
)
# Check the coordinate systems only if there is more than one pixel
m1 = self._measurements[utils.firstkey(self._measurements)]
if utils.is_image(m1):
if not self._check_header("CTYPE1"):
raise Exception(
"CTYPE1 of your input Measurements do not match. Please ensure coordinates of all Measurements are"
" the same."
)
if not self._check_header("CTYPE2"):
raise Exception(
"CTYPE2 of your input Measurements do not match. Please ensure coordinates of all Measurements are"
" the same."
)
# Only allow beam = None if single value measurements.
# if not utils.is_image(m1):
# if self._check_header("BMAJ", None) or self._check_header("BMIN", None) or self._check_header("BPA", None):
# utils.warn(self, "No beam parameters in Measurement headers, assuming they are all equal.")
# if not self._check_header("BUNIT") ...
[docs]
def run(self, **kwargs):
"""Run the full computation using all the added observations.
Checks compatibility of input observations (beam parameters, coordinate
types, axes lengths) and raises exceptions if they don’t match.
Parameters
----------
mask : list or None, optional
Indicate how to mask image observations before computing density and
radiation field. Possible values:
- ``[‘mad’, multiplier]`` — mask values between +/- multiplier*mad_std
- ``[‘data’, (low, high)]`` — mask data values between low and high
- ``[‘clip’, (low, high)]`` — mask data values outside [low, high]
- ``[‘error’, (low, high)]`` — mask where error pixel is below low or above high
- ``None`` — no masking (default)
radiation_field_range : :class:`~astropy.units.Quantity` or sequence, optional
Restrict the fit to radiation field values within
``[lower, upper]`` (inclusive), excluding unphysical regimes.
Accepts either a length-2 :class:`~astropy.units.Quantity`
(e.g. ``[100, 1000]*habing_unit``) or a length-2 sequence whose
elements are each ``None`` or a scalar
:class:`~astropy.units.Quantity` (e.g. ``[None, 1000*habing_unit]``
for an upper limit only). ``None`` on a side means no limit there.
Note the unit must be *inside* the brackets for a one-sided limit,
since ``None*unit`` is not allowed. Default: ``None`` (no restriction).
density_range : :class:`~astropy.units.Quantity` or sequence, optional
Restrict the fit to density values within ``[lower, upper]``
(inclusive). Same accepted forms as ``radiation_field_range``,
e.g. ``[1e3, 1e4]/u.cm**3`` or ``[None, 1e4/u.cm**3]``. ``None`` on
a side means no limit there. Default: ``None`` (no restriction).
method : str, optional
Fitting method. Default: ``’leastsq’`` (Levenberg-Marquardt). See
https://lmfit-py.readthedocs.io/en/latest/fitting.html#fit-methods-table.
nan_policy : str, optional
Action if fit returns NaN. One of ``’raise’`` (default), ``’propagate’``, ``’omit’``.
workers : int or None, optional
Worker processes for parallel pixel fitting. ``None`` uses serial fitting.
``-1`` uses all CPUs. Ignored for single-pixel fits, emcee, and when
``joint_fit`` is not None. Default: None.
joint_fit : str or None, optional
Controls joint pixel fitting:
- ``None`` (default): serial or parallel per-pixel fitting.
- ``’hybrid’``: joint scipy fit + single-pixel re-fits for boundary pixels. ~3× faster than ``workers=-1``.
- ``’fast’``: joint scipy fit only, no post-processing. ~11× faster but ~7–9% accuracy loss near boundaries.
All joint-fit modes use TRF. Ignored for single-pixel fits and emcee.
Raises
------
Exception
If no models match the input observations, observations are incompatible,
parameters are unrecognized, or NaN is encountered.
ValueError
If ``radiation_field_range`` or ``density_range`` is malformed, has
the wrong units, or specifies a window that does not overlap the
model grid.
"""
# @todo global masking for 'data', 'clip', 'error' not entirely useful unless all data/error have same ranges.
# need something like ['data',['key1':(low,hi), 'key2',(low,hi),...], which is very complicated.
# or data/error cut based on histogram
kwargs_opts = {
"mask": None,
"method": "leastsq",
"nan_policy": "raise",
"refine": True,
# phase-space limits (issue #236)
"radiation_field_range": None,
"density_range": None,
# for emcee
"burn": 0,
"steps": 1000,
# parallelism
"workers": None,
# joint fitting: None | 'hybrid' | 'fast'
"joint_fit": None,
# debugging
"test": False,
"profile": False,
}
kwargs_opts.update(kwargs)
profile = kwargs_opts.pop("profile")
self._stats = None
if profile:
pr = cProfile.Profile()
pr.enable()
self._check_compatibility()
self._ratiocount = None
self.read_models()
# Normalize and validate user phase-space limits now that model axis
# units/extents are known (read_models sets density_unit/radiation_field_unit).
# Raises ValueError up front if a range does not overlap the grid.
self._normalize_phase_space_limits(kwargs_opts.pop("radiation_field_range"), kwargs_opts.pop("density_range"))
self._reset_masks()
self._mask_measurements(kwargs_opts["mask"])
kwargs_opts.pop("mask")
self._compute_valid_ratios()
if self.ratiocount == 0:
raise Exception("No models were found that match your data. Check ModelSet.supported_ratios.")
# Pre-flatten observed ratio data and errors once so _residual_single_pixel
# does not call flatten() on every minimizer iteration.
self._observedratios_flat = {
k: (v.data.flatten(), v.uncertainty.array.flatten()) for k, v in self._observedratios.items()
}
# eventually need to check that the maps overlap in real space.
self._compute_residual()
self._nan_policy = kwargs_opts["nan_policy"]
self._minimizer = Minimizer(self._residual_single_pixel, params=None, nan_policy=self._nan_policy)
# need to pop nan_policy and test so that it does not get passed to Minimzer.minimize()
kwargs_opts.pop("nan_policy", None)
kwargs_opts.pop("test", None)
self._compute_chisq()
self._coarse_density_radiation_field()
if kwargs_opts["refine"]:
kwargs_opts.pop("refine")
joint_fit = kwargs_opts.pop("joint_fit")
self._refine_density_radiation_field(joint_fit=joint_fit, **kwargs_opts)
if profile:
pr.disable()
s = io.StringIO()
sortby = pstats.SortKey.CUMULATIVE
ps = pstats.Stats(pr, stream=s).sort_stats(sortby)
ps.print_stats()
self._stats = s
def _reset_masks(self):
for m in self._measurements:
self._measurements[m].mask = deepcopy(self._masks[m])
def _mask_measurements(self, mask):
"""Set the mask on the measurements based on noise characteristics.
Prevents computing garbage n,G0 where observed ratios are noise / noise.
Parameters
----------
mask : list or None
Indicate how to mask image observations. See :meth:`run` for valid values.
"""
if mask is None:
return
if self._measurementnaxis == 0:
utils.warn(self, "Ignoring 'mask' parameter for single pixel observations")
return
if mask[0] == "mad":
for k, v in self._measurements.items():
sigcut = mask[1] * astats.mad_std(v.data, ignore_nan=True)
print(f"Masking {k} data between [{-sigcut:.1e},{sigcut:.1e}]")
masked_data = ma.masked_inside(v.data, -sigcut, sigcut, copy=True)
# CCDData/NDData do not use MaskArrays underneath but two nddata.arrays. Why??
# Make a copy so we are keeping references to data copies lying around.
v.mask = masked_data.mask.copy()
elif mask[0] == "data":
for k, v in self._measurements.items():
masked_data = ma.masked_inside(v.data, mask[1][0], mask[1][1], copy=True)
v.mask = masked_data.mask.copy()
print(f"Masking {k} data between [{mask[1][0]:.1e},{mask[1][1]:.1e}]")
elif mask[0] == "clip":
for k, v in self._measurements.items():
masked_data = ma.masked_outside(v.data, mask[1][0], mask[1][1], copy=True)
v.mask = masked_data.mask.copy()
print(f"Masking {k} data outside [{mask[1][0]:.1e},{mask[1][1]:.1e}]")
elif mask[0] == "error":
for k, v in self._measurements.items():
# error is StdDevUncertainty so must use _array to get at raw values
indices = np.where((v.error <= mask[1][0]) | (v.error >= mask[1][1]))
if v.mask is not None:
v.mask[indices] = True
else:
v.mask = np.full(v.data.shape, False)
v.mask[indices] = True
print(f"Masking {k} data where error outside [{mask[1][0]:.1e},{mask[1][1]:.1e}]")
else:
raise ValueError(f"Unrecognized mask parameter {mask[0]}. Valid values are 'mad','data','error'")
def _compute_valid_ratios(self):
"""Compute the valid observed ratio maps for the available model data"""
if not self._check_measurement_shapes():
raise Exception("Measurement maps have different dimensions")
# Note _find_ratio_elements does not handle case of OI+CII/FIR so
# we have to deal with that separately below.
z = self._modelset._find_ratio_elements(self.measurementIDs)
self._observedratios = dict()
for p in z:
label = p["numerator"] + "/" + p["denominator"]
# deepcopy workaround for bug: https://github.com/astropy/astropy/issues/9006
num = utils.convert_if_necessary(self._measurements[p["numerator"]])
denom = utils.convert_if_necessary(self._measurements[p["denominator"]])
self._observedratios[label] = deepcopy(num / denom)
# @TODO create a meaningful header for the ratio map
self._ratioHeader(p["numerator"], p["denominator"], label)
self._observedshape = self._observedratios[label].data.shape
self._observedshape = np.array(self._observedshape)
self._add_oi_cii_fir()
def _add_oi_cii_fir(self):
"""add special case ([O I] 63 micron + [C II] 158 micron)/IFIR to observed ratios"""
m = self.measurementIDs
if "CII_158" in m and "FIR" in m:
if "OI_63" in m:
lab = "OI_63+CII_158/FIR"
oi = utils.convert_if_necessary(self._measurements["OI_63"])
cii = utils.convert_if_necessary(self._measurements["CII_158"])
a = deepcopy(oi + cii)
b = deepcopy(self._measurements["FIR"])
self._observedratios[lab] = a / b
self._observedratios[lab].meta = deepcopy(b.header)
self._ratioHeader("OI_63+CII_158", "FIR", lab)
if "OI_145" in m:
lab = "OI_145+CII_158/FIR"
oi = utils.convert_if_necessary(self._measurements["OI_145"])
cii = utils.convert_if_necessary(self._measurements["CII_158"])
aa = deepcopy(oi + cii)
bb = deepcopy(self._measurements["FIR"])
self._observedratios[lab] = aa / bb
self._observedratios[lab].meta = deepcopy(bb.header)
self._ratioHeader("OI_145+CII_158", "FIR", lab)
# function to minimize in single-pixel case
def _residual_single_pixel(self, params, index):
parvals = params.valuesdict()
mvalue = np.empty(self.ratiocount)
dvalue = np.empty(self.ratiocount)
evalue = np.empty(self.ratiocount)
i = 0
# @todo create temporary model, data, error as np arrays so the loop over k dictionary key
# isn't necessary.
# OR save residuals from computeDelta so you don't do this twice.
for k in self._modelratios:
mvalue[i] = self._modelratios[k].get(parvals["density"], parvals["radiation_field"])
dvalue[i] = self._observedratios_flat[k][0][index]
evalue[i] = self._observedratios_flat[k][1][index]
i = i + 1
return (dvalue - mvalue) / evalue
def _residual_multi_pixel(self, params, index):
# this is currently slower than the 'dumb' way of residual_single_pixel!
parvals = params.valuesdict()
return self._interp_resid(parvals["density"], parvals["radiation_field"], index)
def _compute_residual(self):
"""Compute the residual values from the observed ratios and models"""
if self.ratiocount < 2:
msg = (
"Not enough ratios. You need to provide at least 3 observations that can be used to compute 2 ratios"
f" that are covered by the ModelSet. From your observations, {self.ratiocount:d} ratio(s)"
)
if self.ratiocount > 0:
msg += f" {list(self._modelratios.keys())}"
msg += " can be computed."
raise Exception(msg)
if not self._check_ratio_shapes():
raise Exception("Observed ratio maps have different dimensions")
self._residual = dict()
for r in self._observedratios:
modelpix = self._modelratios[r].data.flatten()
mdata = ma.masked_invalid(self._observedratios[r].value)
merror = ma.masked_invalid(self._observedratios[r].error)
# Vectorized: broadcast modelpix over all spatial pixels at once.
# modelpix shape: (n_model_pix,); mdata shape: (ny, nx) or scalar.
# modelpix_exp shape: (n_model_pix, 1, ...) to broadcast against mdata.
# result shape: (n_model_pix, ny, nx) or (n_model_pix,).
modelpix_exp = modelpix.reshape((-1,) + (1,) * mdata.ndim)
residuals_arr = ma.masked_invalid((mdata[np.newaxis, ...] - modelpix_exp) / merror[np.newaxis, ...])
# result order is g0,n,y,x
# Catch the case of a single pixel
if self._observedratios[r].is_single_pixel():
newshape = np.hstack(self._modelratios[r].shape)
_meta = deepcopy(self._modelratios[r].meta)
_wcs = deepcopy(self._modelratios[r].wcs)
# clean potential crap
_meta.pop("", None)
_meta.pop("TITLE", None)
else:
newshape = np.hstack((self._modelratios[r].shape, self._observedratios[r].shape))
_meta = deepcopy(self._observedratios[r].meta)
_wcs = deepcopy(self._observedratios[r].wcs)
# result order is g0,n,y,x
_qq = np.squeeze(np.reshape(residuals_arr, newshape))
self._residual[r] = CCDData(_qq, unit="adu", wcs=_wcs, meta=_meta)
self._fancy_index_residual()
def _fancy_index_residual(self):
# create a custom numpy array for interpolating the residuals
# during fitting.
fk = utils.firstkey(self._modelratios)
resid_index = np.array(np.arange(len(self._residual)))
# use linear interpolation grid from model
rad_index = self._modelratios[fk]._world_axis_lin[1]
density_index = self._modelratios[fk]._world_axis_lin[0]
image_index = np.arange(self._observedratios[fk].size)
self._interpgrid = (resid_index, rad_index, density_index, image_index)
# put the residual data into a single numpy array
self._interpvalues = np.stack([self._residual[r].data for r in self._residual.keys()])
# the flatten the spatial indices, so indices are then [radiation_field,density,pixel]
self._interpvalues = self._interpvalues.reshape(*self._interpvalues.shape[:-2], -1)
def _interp_resid(self, density, radiation_field, pixel):
# density and radiation field must be linear not logarithmic values.
# pixel is 0-based pixel index into flattened map data array.
vals = list()
for i in range(len(self._residual)):
vals.append(interpn(self._interpgrid, self._interpvalues, (i, radiation_field, density, pixel)))
return np.array(vals)
def _compute_chisq(self):
"""Compute the chi-squared values from observed ratios and models"""
if self.ratiocount < 2:
raise Exception(f"Not enough ratios to compute chisq. Need 2, got {self.ratiocount:d}")
sumary = sum(self._residual[r]._data ** 2 for r in self._residual)
self._dof = len(self._residual) - 1
k = utils.firstkey(self._residual)
_wcs = deepcopy(self._residual[k].wcs)
_meta = deepcopy(self._residual[k].meta)
self._chisq = CCDData(sumary, unit="adu", wcs=_wcs, meta=_meta)
self._reduced_chisq = self._chisq.divide(self._dof)
# must make a copy here otherwise the header is an OrderDict
# instead of astropy.io.fits.header.Header
self._reduced_chisq.header = Header(deepcopy(self._chisq.header))
self._fixheader(self._chisq)
self._fixheader(self._reduced_chisq)
utils.comment("Chi-squared", self._chisq)
utils.comment(f"Reduced Chi-squared (DOF={self._dof:d})", self._reduced_chisq)
self._makehistory(self._chisq)
self._makehistory(self._reduced_chisq)
[docs]
def write_chisq(self, chi="chisq.fits", rchi="rchisq.fits", overwrite=True):
"""Write the chisq and reduced-chisq data to a file.
Parameters
----------
chi : str, optional
FITS file to write the chisq map to. Default: ``"chisq.fits"``.
rchi : str, optional
FITS file to write the reduced chisq map to. Default: ``"rchisq.fits"``.
"""
self._chisq.write(chi, overwrite=overwrite, hdu_mask="MASK", output_verify="silentfix")
self._reduced_chisq.write(rchi, overwrite=overwrite, hdu_mask="MASK", output_verify="silentfix")
@staticmethod
def _build_joint_params(valid_pixels, dflat, rflat, minn, maxn, minfuv, maxfuv):
"""Build an lmfit Parameters object with density_{i} and radiation_field_{i}
for each valid pixel i (0-indexed over valid_pixels). Initial values come
from the coarse fit (dflat[j], rflat[j])."""
params = Parameters()
for i, j in enumerate(valid_pixels):
params.add(f"density_{i}", min=minn, max=maxn, value=dflat[j])
params.add(f"radiation_field_{i}", min=minfuv, max=maxfuv, value=rflat[j])
return params
@staticmethod
def _build_jac_sparsity(n_valid, n_ratios):
"""Build a block-diagonal CSR sparsity matrix for joint pixel fitting.
Shape: (n_valid * n_ratios, 2 * n_valid). Pixel i's residual rows
(i*n_ratios : (i+1)*n_ratios) are non-zero only at columns 2*i
(density_i) and 2*i+1 (radiation_field_i)."""
from scipy.sparse import lil_matrix
sparsity = lil_matrix((n_valid * n_ratios, 2 * n_valid), dtype=np.int8)
for i in range(n_valid):
sparsity[i * n_ratios : (i + 1) * n_ratios, 2 * i : 2 * i + 2] = 1
return sparsity.tocsr()
def _refine_density_radiation_field(self, joint_fit=None, **kwargs):
workers = kwargs.pop("workers", None)
if kwargs["method"] != "emcee":
kwargs.pop("steps")
kwargs.pop("burn")
progress = kwargs.pop("progress", True) # progress bar
use_parallel = workers is not None and workers != 1
else:
progress = kwargs.get("progress", False) # keep the progress keyword for emcee, get vs pop
use_parallel = False # emcee manages its own parallelism
joint_fit = None # emcee manages its own parallelism
# First get the range of density n and radiation field FUV from the
# model space, in order to provide them to the Parameters object.
# Since the wk2006 H2 models have a smaller model space,
# we have to check if H2 is in one of the models used.
keys = list(self._modelratios.keys())
if utils._has_H2(keys):
# this will get the index for the first modelratio that has H2 in it
# https://stackoverflow.com/questions/2170900/get-first-list-index-containing-sub-string
i = next(idx for idx, s in enumerate(keys) if "H2" in s)
fk = keys[i]
else:
fk = keys[0]
# this will work regardless if Y axis is Habing, CGS, Draine, etc
# We do this in linear space because we want the fit report to be in linear space
# which is what observers care about.
x, y = utils.get_xy_from_wcs(self._modelratios[fk], linear=True)
minn = x[0]
maxn = x[-1]
minfuv = y[0]
maxfuv = y[-1]
# Clamp bounds to the user's phase-space window (issue #236). These four
# scalars feed the serial, parallel, and joint refine paths, so this
# single clamp restricts them all.
if self._density_range is not None:
minn = max(minn, self._density_range[0])
maxn = min(maxn, self._density_range[1])
if self._radiation_field_range is not None:
minfuv = max(minfuv, self._radiation_field_range[0])
maxfuv = min(maxfuv, self._radiation_field_range[1])
if self._radiation_field is None or self._density is None:
startn = x[int(len(x) / 2)]
startfuv = y[int(len(y) / 2)]
else:
startn = np.nanmean(self._density.value)
startfuv = np.nanmean(self._radiation_field.value)
self._fitparam = Parameters()
# @todo let user limit these ranges. either they pass in Parameters or
# limits = {'density':[low,hi], 'radiation_field':[low,hi]}
# in any unit and convert
self._fitparam.add("density", min=minn, max=maxn, value=startn)
self._fitparam.add("radiation_field", min=minfuv, max=maxfuv, value=startfuv)
# self._fitparam.pretty_print()
rf = np.empty(self._observedratios[fk].size)
den = np.empty(self._observedratios[fk].size)
rfe = np.empty(self._observedratios[fk].size)
dene = np.empty(self._observedratios[fk].size)
chi = np.empty(self._observedratios[fk].size)
rchi = np.empty(self._observedratios[fk].size)
dflat = self._density.value.flatten()
rflat = self._radiation_field.value.flatten()
fmdata = np.empty(self._observedratios[fk].size, dtype=object)
fm_mask = np.full(shape=self._observedratios[fk].data.shape, fill_value=False).flatten()
count = 0
excount = 0
# turn off progress bar for single pixel or emcee prints out multiple bars.
if self._observedratios[fk].size == 1:
progress = False
use_parallel = False
joint_fit = None
if joint_fit in ("hybrid", "fast"):
# ------------------------------------------------------------------
# Joint pixel fitting: single scipy.optimize.least_squares
# call for all pixels with block-diagonal jac_sparsity.
# Bypasses lmfit to avoid its sparse-Jacobian covariance bug
# (element-wise * instead of matrix @ when computing J^T J).
# ------------------------------------------------------------------
ratio_keys = list(self._modelratios.keys())
n_ratios = self.ratiocount
interps = [self._modelratios[k]._interp_lin for k in ratio_keys]
valid_mask = ~(np.isnan(dflat) | np.isnan(rflat))
valid_pixels = np.where(valid_mask)[0] # indices into flattened map
n_valid = len(valid_pixels)
obs_data = np.array(
[self._observedratios_flat[k][0][valid_pixels] for k in ratio_keys]
) # (n_ratios, n_valid)
obs_err = np.array(
[self._observedratios_flat[k][1][valid_pixels] for k in ratio_keys]
) # (n_ratios, n_valid)
# Residual takes a plain numpy array: x[::2]=density, x[1::2]=rf
def _joint_residual(x):
pts = np.column_stack([x[::2], x[1::2]]) # (n_valid, 2)
mvalues = np.array([interp(pts) for interp in interps]) # (n_ratios, n_valid)
return ((obs_data - mvalues) / obs_err).flatten() # (n_valid * n_ratios,)
# Build x0 (initial guess) and bounds from coarse-fit values
x0 = np.empty(2 * n_valid)
x0[::2] = dflat[valid_pixels]
x0[1::2] = rflat[valid_pixels]
lb = np.empty(2 * n_valid)
ub = np.empty(2 * n_valid)
lb[::2] = minn
ub[::2] = maxn
lb[1::2] = minfuv
ub[1::2] = maxfuv
sparsity = self._build_jac_sparsity(n_valid, n_ratios)
joint_result = _scipy_least_squares(
_joint_residual,
x0,
bounds=(lb, ub),
jac_sparsity=sparsity,
tr_solver="lsmr",
method="trf",
)
# Mark NaN/invalid pixels as masked
fm_mask[~valid_mask] = True
den[~valid_mask] = dene[~valid_mask] = np.nan
rf[~valid_mask] = rfe[~valid_mask] = np.nan
chi[~valid_mask] = rchi[~valid_mask] = np.nan
# Compute per-pixel stderr from block-diagonal Jacobian blocks.
# J_i = jac[i*n_ratios:(i+1)*n_ratios, 2*i:2*i+2] shape (n_ratios, 2)
# cov_i = inv(J_i.T @ J_i); stderr = sqrt(diag(cov_i))
jac = joint_result.jac # sparse or dense
all_resid = joint_result.fun.reshape(n_valid, n_ratios)
for i, j in enumerate(valid_pixels):
den[j] = joint_result.x[2 * i]
rf[j] = joint_result.x[2 * i + 1]
chi[j] = float(np.sum(all_resid[i] ** 2))
rchi[j] = chi[j] / max(n_ratios - 1, 1)
# Per-pixel stderr from local Jacobian block
r0, r1 = i * n_ratios, (i + 1) * n_ratios
c0, c1 = 2 * i, 2 * i + 2
J_i = jac[r0:r1, c0:c1]
if hasattr(J_i, "toarray"):
J_i = J_i.toarray()
try:
JTJ = J_i.T @ J_i
cov_i = np.linalg.inv(JTJ)
stderr_den = float(np.sqrt(max(cov_i[0, 0], 0.0)))
stderr_rf = float(np.sqrt(max(cov_i[1, 1], 0.0)))
except np.linalg.LinAlgError:
stderr_den = stderr_rf = None
dene[j] = stderr_den
rfe[j] = stderr_rf
fmdata[j] = _PixelResult(
den[j],
dene[j],
rf[j],
rfe[j],
chi[j],
rchi[j],
n_ratios,
joint_result.success,
all_resid[i],
)
count += 1
# ------------------------------------------------------------------
# Hybrid post-processing ('hybrid' only, skipped for 'fast'):
# Re-fit pixels that did not move from their coarse initial guess.
# These are typically pixels initialised at a model boundary where
# the global LSMR convergence criterion is satisfied before the
# local 2-parameter problem for that pixel is resolved.
# ------------------------------------------------------------------
if joint_fit == "hybrid":
stuck = np.isclose(joint_result.x[::2], x0[::2], rtol=1e-4, atol=0) & np.isclose(
joint_result.x[1::2], x0[1::2], rtol=1e-4, atol=0
)
print(f"Refitting {len(stuck)} stuck pixels...")
for i in np.where(stuck)[0]:
j = valid_pixels[i]
obs_j = obs_data[:, i]
err_j = obs_err[:, i]
def _sp_resid(x, _obs=obs_j, _err=err_j):
pts = x.reshape(1, 2)
mvals = np.array([interp(pts)[0] for interp in interps])
return (_obs - mvals) / _err
r = _scipy_least_squares(
_sp_resid,
[x0[2 * i], x0[2 * i + 1]],
bounds=([minn, minfuv], [maxn, maxfuv]),
method="trf",
)
den[j] = r.x[0]
rf[j] = r.x[1]
resid_i = r.fun
chi[j] = float(np.sum(resid_i**2))
rchi[j] = chi[j] / max(n_ratios - 1, 1)
J_i = r.jac
try:
cov_i = np.linalg.inv(J_i.T @ J_i)
stderr_den = float(np.sqrt(max(cov_i[0, 0], 0.0)))
stderr_rf = float(np.sqrt(max(cov_i[1, 1], 0.0)))
except np.linalg.LinAlgError:
stderr_den = stderr_rf = None
dene[j] = stderr_den
rfe[j] = stderr_rf
fmdata[j] = _PixelResult(
den[j],
dene[j],
rf[j],
rfe[j],
chi[j],
rchi[j],
n_ratios,
r.success,
resid_i,
)
elif use_parallel:
# ------------------------------------------------------------------
# Parallel pixel loop using ProcessPoolExecutor.
# Model interpolators are built once per worker process by
# _init_worker to avoid re-sending large model arrays per task.
# ------------------------------------------------------------------
ratio_keys = list(self._modelratios.keys())
model_points = [self._modelratios[k]._world_axis_lin for k in ratio_keys]
model_values = [self._modelratios[k].data.T for k in ratio_keys]
obs_data_arr = np.array([self._observedratios_flat[k][0] for k in ratio_keys])
obs_err_arr = np.array([self._observedratios_flat[k][1] for k in ratio_keys])
nan_policy = getattr(self, "_nan_policy", "raise")
max_workers = None if workers == -1 else workers
futures = {}
with ProcessPoolExecutor(
max_workers=max_workers,
initializer=_init_worker,
initargs=(model_points, model_values),
) as pool:
for j in range(self._observedratios[fk].size):
if np.isnan(dflat[j]) or np.isnan(rflat[j]):
fmdata[j] = None
fm_mask[j] = True
den[j] = dene[j] = rf[j] = rfe[j] = chi[j] = rchi[j] = np.nan
else:
futures[
pool.submit(
_fit_pixel_worker,
j,
obs_data_arr[:, j],
obs_err_arr[:, j],
dflat[j],
rflat[j],
minn,
maxn,
minfuv,
maxfuv,
nan_policy,
dict(kwargs),
)
] = j
with get_progress_bar(progress, len(futures), leave=True, position=0) as pbar:
for fut in as_completed(futures):
j = futures[fut]
try:
_, result = fut.result()
fmdata[j] = result
count += 1
den[j] = result.params["density"].value
dene[j] = result.params["density"].stderr
rf[j] = result.params["radiation_field"].value
rfe[j] = result.params["radiation_field"].stderr
chi[j] = result.chisqr
rchi[j] = result.redchi
except ValueError:
excount += 1
fmdata[j] = None
fm_mask[j] = True
den[j] = dene[j] = rf[j] = rfe[j] = chi[j] = rchi[j] = np.nan
pbar.update(1)
else:
# ------------------------------------------------------------------
# Serial pixel loop (default).
# ------------------------------------------------------------------
with get_progress_bar(progress, self._observedratios[fk].size, leave=True, position=0) as pbar:
for j in range(self._observedratios[fk].size):
# use previous coarse fit as first guess
self._fitparam["density"].value = dflat[j]
self._fitparam["radiation_field"].value = rflat[j]
if np.isnan(dflat[j]) or np.isnan(rflat[j]):
fmdata[j] = None
fm_mask[j] = True
den[j] = np.nan
dene[j] = np.nan
rfe[j] = np.nan
rf[j] = np.nan
chi[j] = np.nan
rchi[j] = np.nan
else:
try:
self._minimizer.userargs = (j,)
fmdata[j] = self._minimizer.minimize(params=self._fitparam, **kwargs)
# if hasattr(fmdata[j],"success") ugh. not guaranteed
# if fmdata[j].errorbars:
count = count + 1
den[j] = fmdata[j].params["density"].value
dene[j] = fmdata[j].params["density"].stderr
rf[j] = fmdata[j].params["radiation_field"].value
rfe[j] = fmdata[j].params["radiation_field"].stderr
chi[j] = fmdata[j].chisqr
rchi[j] = fmdata[j].redchi
# else:
# fmdata[j] = None
# fm_mask[j] = True
except ValueError:
# print("At pixel %d, got valuerror %s with fitparams %s" %(j, exc,self._fitparam))
excount = excount + 1
fmdata[j] = None
fm_mask[j] = True
den[j] = np.nan
dene[j] = np.nan
rfe[j] = np.nan
rf[j] = np.nan
chi[j] = np.nan
rchi[j] = np.nan
pbar.update(1)
fmdata = fmdata.reshape(self._observedratios[fk].data.shape)
fm_mask = fm_mask.reshape(self._observedratios[fk].data.shape)
# ff_mask = ff_mask | np.logical_not(np.isfinite(/*something*/))
self._fitresult = FitMap(fmdata, wcs=self._observedratios[fk].wcs, mask=fm_mask, name="result")
# print(f"fitted {count} of {self._observedratios[fk].size} pixels")
# print(f"got {excount} exceptions")
rshape = self._radiation_field.data.shape
dshape = self._density.data.shape
self._radiation_field.data = rf.reshape(rshape)
self._radiation_field.uncertainty.array = rfe.reshape(rshape)
self._density.data = den.reshape(dshape)
self._density.uncertainty.array = dene.reshape(dshape)
self._chisq_min.data = chi.reshape(dshape)
self._reduced_chisq_min.data = rchi.reshape(dshape)
def _coarse_density_radiation_field(self):
"""Compute the best-fit density and radiation field spatial maps
by searching for the minimum chi-squared at each spatial pixel."""
if self._chisq is None or self._reduced_chisq is None:
return
# get the chisq minima of each pixel along the g,n axes
fk = utils.firstkey(self._modelratios)
mshape = self._modelratios[fk].shape
# Wolfire 2006 models have NAXIS=2, while 2020+ have NAXIS=3.
# Deal with it.
# @see Measurement squeeze parameter. This should no longer be needed
if self._modelnaxis == 2:
firstindex = 0
secondindex = 1
thirdindex = 2
fourthindex = 3
elif self._modelnaxis == 3:
if mshape[0] != 1:
raise Exception(f"Unexpected NAXIS3 != 1 in model {fk}")
firstindex = 1
secondindex = 2
thirdindex = 3
fourthindex = 4
# Restrict the coarse grid search to the allowed phase-space window
# (issue #236). Out-of-range model cells are set to +inf on local copies
# (the stored full chisq cubes are left untouched) so the argmin below
# is confined to the physical region. Validation in
# _normalize_phase_space_limits guarantees >=1 in-range cell survives.
allowed = self._phase_space_grid_mask(self._reduced_chisq.data.ndim, firstindex, secondindex)
if allowed is None:
masked_rchi = self._reduced_chisq.data
masked_chi = self._chisq.data
else:
masked_rchi = np.where(allowed, self._reduced_chisq.data, np.inf)
masked_chi = np.where(allowed, self._chisq.data, np.inf)
rchi_min = np.amin(masked_rchi, (firstindex, secondindex))
chi_min = np.amin(masked_chi, (firstindex, secondindex))
gnxy = np.where(masked_rchi == rchi_min)
gi = gnxy[firstindex]
ni = gnxy[secondindex]
# spatial_idx is which indices are the spatial axes
if len(gnxy) >= 4:
# astronomical spatial indices
spatial_idx = (gnxy[thirdindex], gnxy[fourthindex])
else:
spatial_idx = 0
# model n,g0 indices
model_idx = np.transpose(np.array([ni, gi]))
if self._modelnaxis == 3:
# add 3rd axis to model_idx
model_idx = np.insert(model_idx, 0, [0], axis=1)
fk2 = utils.firstkey(self._observedratios)
newshape = self._observedratios[fk2].shape
g0 = 10 ** (self._modelratios[fk].wcs.wcs_pix2world(model_idx, 0))[:, 1]
n = 10 ** (self._modelratios[fk].wcs.wcs_pix2world(model_idx, 0))[:, 0]
self._radiation_field = deepcopy(self._observedratios[fk2])
if spatial_idx == 0 and newshape == (1,):
self._radiation_field.data = g0
self._radiation_field.uncertainty.array = np.array([np.nan])
else:
if self.has_vectors:
self._radiation_field.data = g0
else: # Measurement with image
# note this will reshape g0 in radiation_field for us!
self._radiation_field.data[spatial_idx] = g0
# We cannot mask nans because numpy does not support writing
# MaskedArrays to a file. Will get a not implemented error.
# Therefore just copy the nans over from the input observations.
self._radiation_field.data[np.isnan(self._observedratios[fk2])] = np.nan
# kluge because we dont know how to properly calcultate uncertainty on this.
# self._radiation_field.uncertainty.array=np.zeroes(self._radiation_field.uncertainty.array)
self._radiation_field.uncertainty.array[:] = np.nan
self._radiation_field.unit = self.radiation_field_unit
self._radiation_field.uncertainty.unit = self.radiation_field_unit
self._density = deepcopy(self._observedratios[fk2])
if spatial_idx == 0 and newshape == (1,):
self._density.data = n
self._density.uncertainty.array = np.array([np.nan])
else:
if self.has_vectors: # Measurement with data vector
self._density.data = n
else: # Measurement with image
# note this will reshape g0 in radiation_field for us!
self._density.data[spatial_idx] = n
self._density.data[np.isnan(self._observedratios[fk2])] = np.nan
# kluge because we dont know how to properly calcultate undertainty on this.
# self._density.uncertainty.array=np.zeroes(self._density.uncertainty.array)
self._density.uncertainty.array[:] = np.nan
self._density.unit = self.density_unit
self._density.uncertainty.unit = self.density_unit
# this raises exception, CCDData enforces both units the same
# self._density.uncertainty.unit = u.dimensionless_unscaled
# fix the headers
self._density_radiation_field_header()
# now save copies of the 2D min chisquares
self._chisq_min = deepcopy(self._observedratios[fk2])
if spatial_idx == 0 and newshape == (1,):
self._chisq_min.data = np.array([chi_min])
elif self._modelnaxis == 2:
self._chisq_min.data = chi_min
else:
self._chisq_min.data = chi_min[0, :, :]
self._chisq_min.data[np.isnan(self._observedratios[fk2])] = np.nan
self._chisq_min.unit = u.dimensionless_unscaled
self._chisq_min.uncertainty.array = [0.0]
self._chisq_min.uncertainty.unit = u.dimensionless_unscaled
self._reduced_chisq_min = deepcopy(self._observedratios[fk2])
if spatial_idx == 0 and newshape == (1,):
self._reduced_chisq_min.data = np.array([rchi_min])
else:
if self._modelnaxis == 2:
self._reduced_chisq_min.data = rchi_min
else:
self._reduced_chisq_min.data = rchi_min[0, :, :]
self._reduced_chisq_min.data[np.isnan(self._observedratios[fk2])] = np.nan
self._reduced_chisq_min.unit = u.dimensionless_unscaled
self._reduced_chisq_min.uncertainty.array = [0.0]
self._reduced_chisq_min.uncertainty.unit = u.dimensionless_unscaled
# update histories
utils.setkey("BUNIT", "Minimum Chi-squared", self._chisq_min)
utils.setkey("BUNIT", f"Minimum Reduced Chi-squared (DOF={self._dof:d})", self._reduced_chisq_min)
self._makehistory(self._reduced_chisq_min)
self._makehistory(self._chisq_min)
def _phase_space_grid_mask(self, ndim, firstindex, secondindex):
"""Boolean mask selecting model grid cells inside the allowed
phase-space window (issue #236), or ``None`` if unrestricted.
The returned array is broadcastable to the chi-squared cube: it is 1
along every axis except the radiation-field axis (``firstindex``) and
the density axis (``secondindex``).
Parameters
----------
ndim : int
Number of dimensions of the chi-squared cube.
firstindex : int
Cube axis index of the radiation field.
secondindex : int
Cube axis index of the density.
Returns
-------
:class:`numpy.ndarray` of bool or None
``None`` when no phase-space limits are active.
"""
if self._density_range is None and self._radiation_field_range is None:
return None
fk = utils.firstkey(self._modelratios)
# linear physical axis values with units: x=density, y=radiation field
x, y = utils.get_xy_from_wcs(self._modelratios[fk], quantity=True, linear=True)
xden = x.to(self.density_unit).value # values along secondindex
yrf = y.to(self.radiation_field_unit).value # values along firstindex
den_ok = np.ones(len(xden), dtype=bool)
rf_ok = np.ones(len(yrf), dtype=bool)
if self._density_range is not None:
dlo, dhi = self._density_range
den_ok = (xden >= dlo) & (xden <= dhi)
if self._radiation_field_range is not None:
rlo, rhi = self._radiation_field_range
rf_ok = (yrf >= rlo) & (yrf <= rhi)
den_shape = [1] * ndim
den_shape[secondindex] = len(xden)
rf_shape = [1] * ndim
rf_shape[firstindex] = len(yrf)
return rf_ok.reshape(rf_shape) & den_ok.reshape(den_shape)
def _makehistory(self, image):
"""Add HISTORY keyword indicating how density and radiation field were computed.
Parameters
----------
image : :class:`astropy.io.fits.ImageHDU`, :class:`astropy.nddata.CCDData`, or :class:`~pdrtpy.measurement.Measurement`
The image to which to add the history.
"""
s = "Measurements provided: " + str(list(self._measurements.keys()))
utils.history(s, image)
s = "Ratios used: " + str(list(self._residual.keys()))
utils.history(s, image)
utils.signature(image)
utils.dataminmax(image)
def _ratioHeader(self, numerator, denominator, label):
"""Add the RATIO identifier to the appropriate image.
Parameters
----------
numerator : str
Numerator key of the line ratio.
denominator : str
Denominator key of the line ratio.
label : str
Ratio key indicating which observation image (Measurement) to use.
"""
utils.addkey("RATIO", label, self._observedratios[label])
utils.dataminmax(self._observedratios[label])
utils.signature(self._observedratios[label])
def _fixheader(self, image):
"""Put additional axis and header values into an image.
Parameters
----------
image : :class:`astropy.io.fits.ImageHDU`, :class:`astropy.nddata.CCDData`, or :class:`~pdrtpy.measurement.Measurement`
The image to which to add the header values.
"""
if self._modelnaxis == 2:
naxis = len(image.shape)
else:
naxis = len(image.shape) - 1
ax1 = str(naxis - 1)
ax2 = str(naxis)
if "NAXIS" not in image.header:
utils.setkey("NAXIS", naxis, image)
utils.setkey("NAXIS" + ax1, image.shape[1], image)
utils.setkey("NAXIS" + ax2, image.shape[0], image)
utils.setkey("CTYPE" + ax1, self.density_type, image)
utils.setkey("CTYPE" + ax2, self.radiation_field_type, image)
utils.setkey("CUNIT" + ax1, str(self.density_unit), image)
utils.setkey("CUNIT" + ax2, str(self.radiation_field_unit), image)
fk = utils.firstkey(self._modelratios)
mod = self._modelratios[fk]
utils.setkey("CDELT" + ax1, mod.wcs.wcs.cdelt[0], image)
utils.setkey("CDELT" + ax2, mod.wcs.wcs.cdelt[1], image)
utils.setkey("CRVAL" + ax1, mod.wcs.wcs.crval[0], image)
utils.setkey("CRVAL" + ax2, mod.wcs.wcs.crval[1], image)
utils.setkey("CRPIX" + ax1, mod.wcs.wcs.crpix[0], image)
utils.setkey("CRPIX" + ax2, mod.wcs.wcs.crpix[1], image)
def _density_radiation_field_header(self):
"""Common header items in the density and radiation field FITS files"""
self._density.header.pop("RATIO")
self._radiation_field.header.pop("RATIO")
# note: must use to_string() here or astropy.io.fits.Card complains
# about the value being a Unit. Oddly it doesn't complain for the
# data units. Go figure.
utils.setkey("BUNIT", self.density_unit.to_string(), self._density)
utils.comment("Best-fit H2 volume density", self._density)
utils.setkey("BUNIT", self.radiation_field_unit.to_string(), self._radiation_field)
utils.comment("Best-fit interstellar radiation field", self._radiation_field)
self._makehistory(self._density)
self._makehistory(self._radiation_field)
# convert from OrderedDict to astropy.io.fits.header.Header
self._density.header = Header(self._density.header)
self._radiation_field.header = Header(self._radiation_field.header)
self._density._identifier = "H2 Volume Density"
self._radiation_field._identifier = "Radiation Field"
@property
def table(self):
# @TODO: make this work for map data ?
r"""Construct the table of input Measurements and, if the fit has been run, the density, radiation field, and :math:`\chi^2` values.
Returns
-------
:class:`astropy.table.Table`
"""
v = self._measurements.values()
# This only works for astropy version >= 4.1
# and for some reason requirements.txt did not install it
# for a test user.
# t = Table(self._measurements,
# units=[m.unit for m in v]
# )
t = Table()
cols = [Column(data=d, unit=d.unit) for d in v]
t.add_columns(cols=cols, names=[m.id for m in v])
if self._observedratios is not None:
v = self._observedratios.values()
cols = [Column(data=d, unit=d.unit) for d in v]
t.add_columns(cols=cols, names=[m.id for m in v])
if self.radiation_field is not None:
t.add_column(col=Column(self.radiation_field, unit=self.radiation_field_unit), name=self.radiation_field.id)
if self.density is not None:
t.add_column(col=Column(self.density, unit=self.density_unit), name=self.density.id)
if self._chisq_min is not None:
t.add_column(col=Column(self._chisq_min, unit=None), name="Chi-square")
for j in t.columns:
t[j].format = "3.2E"
return t
