diff --git a/doc/dataformat.rst b/doc/dataformat.rst index 36ea4ac2..29051163 100644 --- a/doc/dataformat.rst +++ b/doc/dataformat.rst @@ -79,7 +79,7 @@ For a single observation, you might do something like: .. code-block:: python def build_obs(N): - from prospect.data import Spectrum + from prospect.observation import Spectrum N = 1000 # number of wavelength points spec = Spectrum(wavelength=np.linspace(3000, 5000, N), flux=np.zeros(N), uncertainty=np.ones(N)) # ensure that this is a valid observation for fitting @@ -95,7 +95,7 @@ For photometry this might look like: .. code-block:: python def build_obs(N): - from prospect.data import Photometry + from prospect.observation import Photometry # valid sedpy filter names fnames = list([f"sdss_{b}0" for b in "ugriz"]) Nf = len(fnames) @@ -113,7 +113,7 @@ A tool exists to convert old combined observation dictionaries to a list of .. code-block:: python - from prospect.data import from_oldstyle + from prospect.observation import from_oldstyle # dummy observation dictionary with just a spectrum N = 1000 obs = dict(wavelength=np.linspace(3000, 5000, N), spectrum=np.zeros(N), unc=np.ones(N), @@ -129,7 +129,7 @@ The :py:meth:`build_obs` function The :py:meth:`build_obs` function in the parameter file is written by the user. It should take a dictionary of command line arguments as keyword arguments. It -should return a list of :py:class:`prospect.data.Observation` instances, +should return a list of :py:class:`prospect.observation.Observation` instances, described above. Other than that, the contents can be anything. Within this function you might @@ -142,7 +142,7 @@ The point of this function is that you don't have to *externally* convert your data format to be what |Codename| expects and keep another version of files lying around: the conversion happens *within* the code itself. Again, the only requirement is that the function can take a ``run_params`` dictionary as keyword -arguments and that it return :py:class:`prospect.data.Observation` instances, as +arguments and that it return :py:class:`prospect.observation.Observation` instances, as described above. Each observation instance should correspond to a particular dataset (e.g. a broadband photomtric SED, the spectrum from a particular instrument, or the spectrum from a particular night) that shares instrumental diff --git a/doc/faq.rst b/doc/faq.rst index 3dbbfad8..1e3728d0 100644 --- a/doc/faq.rst +++ b/doc/faq.rst @@ -40,49 +40,6 @@ sampling algorithm being used, and the kind of data that you have. Hours or even days per fit is not uncommon for more complex models. -Can I fit my spectrum too? --------------------------- -There are several extra considerations that come up when fitting spectroscopy - - 1) Wavelength range and resolution. - Prospector is based on FSPS, which uses the MILES spectral library. These - have a resolution of ~2.5A FWHM from 3750AA - 7200AA restframe, and much - lower (R~200 or so, but not actually well defined) outside this range. - Higher resolution data (after including both velocity dispersion and - instrumental resolution) or spectral points outside this range cannot yet - be fit. - - 2) Relatedly, line spread function. - Prospector includes methods for FFT based smoothing of the spectra, - assuming a gaussian LSF (in either wavelength or velocity space). There is - also the possibility of FFT based smoothing for wavelength dependent - Gaussian dispersion (i.e. sigma_lambda = f(lambda) with f possibly a - polynomial of lambda). In practice the smoothed spectra will be a - combination of the library resolution plus whatever FFT smoothing is - applied. Hopefully this can be made to match your actual data resolution, - which is a combination of the physical velocity dispersion and the - instrumental resolution. The smoothing is controlled by the parameters - `sigma_smooth`, `smooth_type`, and `fftsmooth` - - 3) Nebular emission. - While prospector/FSPS include self-consistent nebular emission, the - treatment is probably not flexible enough at the moment to fit high S/N, - high resolution data including nebular emission (e.g. due to deviations of - line ratios from Cloudy predictions or to complicated gas kinematics that - are different than stellar kinematics). Thus fitting nebular lines should - take adavantage of the nebular line amplitude optimization/marginalization - capabilities. For very low resolution data this is less of an issue. - - 4) Spectrophotometric calibration. - There are various options for dealing with the spectral continuum shape - depending on how well you think your spectra are calibrated and if you - also fit photometry to tie down the continuum shape. You can optimize out - a polynomial "calibration" vector, or simultaneously fit for a polynomial - and marginalize over the polynomial coefficients (this allows you to place - priors on the accuracy of the spectrophotometric calibration). Or you can - just take the spectrum as perfectly calibrated. - - How do I fit for redshift as well as other parameters? ------------------------------------------------------ The simplest way is to just let the parameter specification for ``"zred"`` diff --git a/doc/index.rst b/doc/index.rst index a6817275..d82c2b4e 100644 --- a/doc/index.rst +++ b/doc/index.rst @@ -30,6 +30,7 @@ Prospector allows you to: installation usage dataformat + spectra models sfhs nebular diff --git a/doc/quickstart.rst b/doc/quickstart.rst index 6ca448de..eabc38bf 100644 --- a/doc/quickstart.rst +++ b/doc/quickstart.rst @@ -45,7 +45,7 @@ include an empty Spectrum data set to force a prediction of the full spectrum. .. code:: python from sedpy.observate import load_filters - from prospect.data import Photometry, Spectrum + from prospect.observation import Photometry, Spectrum filters = load_filters([f"sdss_{b}0" for b in bands]) maggies = np.array([10**(-0.4 * cat[0][f"cModelMag_{b}"]) for b in bands]) @@ -59,6 +59,7 @@ include an empty Spectrum data set to force a prediction of the full spectrum. for obs in observations: obs.redshift = shdus[2].data[0]["z"] + In principle we could also add noise models for the spectral and photometric data (e.g. to fit for the photometric noise floor), but we'll make the default assumption of iid Gaussian noise for the moment. @@ -69,18 +70,25 @@ Build a Model Here we will get a default parameter set for a simple parametric SFH, and add a set of parameters describing nebular emission. We'll also fix the redshift to -the value given by SDSS. This model has 5 free parameters, each of which has a -an associated prior distribution. These default prior distributions can and -should be replaced or adjusted depending on your particular science question. +the value given by SDSS. This model has 5 free parameters, each of which has an +associated prior distribution. These default prior distributions can and should +be replaced or adjusted depending on your particular science question. Here +we'll just change the prior distribution for stellar mass, as an example. .. code:: python + # Get a baseline set of model parameters from prospect.models.templates import TemplateLibrary from prospect.models import SpecModel model_params = TemplateLibrary["parametric_sfh"] model_params.update(TemplateLibrary["nebular"]) model_params["zred"]["init"] = obs["redshift"] + # Adjust the prior for mass, giving a wider range + from prospect.models import priors + model_params["mass"]["prior"] = priors.LogUniform(mini=1e6, maxi=1e13) + + # Instantiate the model using this parameter dictionary model = SpecModel(model_params) assert len(model.free_params) == 5 print(model) @@ -100,13 +108,15 @@ spectral and isochrone libraries are being used. sps = CSPSpecBasis(zcontinuous=1) print(sps.ssp.libraries) +For piecewise constant and other flexible SFHs use `FastStepBasis` instead of +`CSPSpecBasis`. Make a prediction ----------------- We can now predict our data for any set of parameters. This will take a little time because fsps is building and caching the SSPs. Subsequent calls to predict -will be faster. Here we'll just make the predicition for the current value of +will be faster. Here we'll just make the prediction for the current value of the free parameters. .. code:: python @@ -114,23 +124,36 @@ the free parameters. current_parameters = ",".join([f"{p}={v}" for p, v in zip(model.free_params, model.theta)]) print(current_parameters) (spec, phot), mfrac = model.predict(model.theta, observations, sps=sps) - print(phot / obs["maggies"]) + print("filter,observed,predicted") + for i, f in enumerate(obs["filters"]): + print(f"{f.name},{obs['maggies'][i]},{phot[i]}") Run a fit --------- -Since we can make predictions and we have data and uncertainties, we should be -able to construct a likelihood function. Here we'll use the pre-defined default +Since we can make predictions and we have (photometric) data and uncertainties, +we should be able to construct a likelihood function, and then combine with the +priors to sample the posterior. Here we'll use the pre-defined default posterior probability function. We also set some sampling related keywords to -make the fit go a little faster, though it should still take of order tens of -minutes. +make the fit go a little faster (but give rougher posterior estimates), though +it should still take of order tens of minutes. .. code:: python from prospect.fitting import lnprobfn, fit_model + + # just the photometry + obs = [observations[1]] + + # posterior probability of the initial parameters given the photometry + lnp = lnprobfn(model.theta, model, observations=obs, sps=sps) + + # now do the posterior sampling fitting_kwargs = dict(nlive_init=400, nested_method="rwalk", nested_target_n_effective=10000, nested_dlogz_init=0.05) - output = fit_model(obs, model, sps, optimize=False, dynesty=True, lnprobfn=lnprobfn, **fitting_kwargs) + output = fit_model(obs, model, sps, lnprobfn=lnprobfn, + optimize=False, dynesty=True, + **fitting_kwargs) result, duration = output["sampling"] The ``result`` is a dictionary with keys giving the Monte Carlo samples of @@ -146,7 +169,7 @@ as follows: from prospect.io import write_results as writer hfile = "./quickstart_dynesty_mcmc.h5" - writer.write_hdf5(hfile, {}, model, obs, + writer.write_hdf5(hfile, {}, model, observations, output["sampling"][0], None, sps=sps, tsample=output["sampling"][1], @@ -179,7 +202,7 @@ We'll mark the highest probablity posterior sample as well. from prospect.plotting import corner nsamples, ndim = out["chain"].shape cfig, axes = pl.subplots(ndim, ndim, figsize=(10,9)) - axes = corner.allcorner(out["chain"].T, out["theta_labels"], axes, weights=out["weights"], color="royalblue", show_titles=True) + #axes = corner.allcorner(out["chain"].T, out["theta_labels"], axes, weights=out["weights"], color="royalblue", show_titles=True) from prospect.plotting.utils import best_sample pbest = best_sample(out) diff --git a/doc/spectra.rst b/doc/spectra.rst new file mode 100644 index 00000000..ef0a96a3 --- /dev/null +++ b/doc/spectra.rst @@ -0,0 +1,64 @@ +Fitting Spectra +================ + +There are several extra considerations that come up when fitting spectroscopy. + +Wavelength range, resolution, and linespread function +----------------------------------------------------- + +Prospector is based on FSPS, which uses stellar spectral libraries with given +resolution. The empirical MILES library has a resolution of ~2.5A FWHM from +3750AA - 7200AA restframe, and much lower (R~200 or so, but not actually well +defined) outside this range. Higher resolution data (after including both +velocity dispersion and instrumental resolution) or spectral points outside this +range cannot yet be fit. + +Prospector includes methods for FFT based smoothing of the spectra, assuming a +Gaussian LSF (in either wavelength or velocity space). There is also the +possibility of FFT based smoothing for wavelength dependent Gaussian dispersion +(i.e. sigma_lambda = f(lambda) with f possibly a polynomial of lambda). In +practice the smoothed model spectra will be a combination of the library resolution +plus whatever FFT smoothing is applied. Hopefully this can be made to match your +actual data resolution, which is a combination of the physical velocity +dispersion and the instrumental resolution. The smoothing is controlled by the +parameters `sigma_smooth`, `smooth_type`, and `fftsmooth` + +For undersampled spectra, a special :py:class:`UnderSampledSpectrum` class +exists that will integrate the model (smoothed by velocity dispersion and +intrumental resolution) over the supplied pixels. + + +Instrumental Response & Spectrophotometric Calibration +--------------------- + +There are various options for dealing with the spectral continuum shape +depending on how well you think your spectra are calibrated and if you also fit +photometry to tie down the continuum shape. You can optimize out a polynomial +"calibration" vector, or simultaneously fit for a polynomial and marginalize +over the polynomial coefficients (this allows you to place priors on the +accuracy of the spectrophotometric calibration). Or you can just take the +spectrum as perfectly calibrated. + +Particular treatments can be implemented using different mix-in classes, e.g. + +.. code-block:: python + + from prospect.observation import Spectrum, PolyOptCal + class PolySpectrum(PolyOptCal, Spectrum): + pass + spec = PolySpectrum(wavelength=np.linspace(3000, 5000, N), + flux=np.zeros(N), + uncertainty=np.ones(N), + polynomial_order=5) + + +Nebular emission +---------------- + +While prospector/FSPS include self-consistent nebular emission, the treatment is +probably not flexible enough at the moment to fit high S/N, high resolution data +including nebular emission (e.g. due to deviations of line ratios from Cloudy +predictions or to complicated gas kinematics that are different than stellar +kinematics). Thus fitting nebular lines should take adavantage of the nebular +line amplitude optimization/marginalization capabilities. For very low +resolution data this is less of an issue. \ No newline at end of file diff --git a/prospect/io/read_results.py b/prospect/io/read_results.py index cff75423..0e9dac33 100644 --- a/prospect/io/read_results.py +++ b/prospect/io/read_results.py @@ -240,9 +240,9 @@ def get_sps(res): "same as the FSPS libraries that you are using now ({})".format(flib, rlib)) # If fitting and reading in are happening in different python versions, # ensure string comparison doesn't throw error: - if type(flib[0]) == 'bytes': + if isinstance(flib[0], bytes): flib = [i.decode() for i in flib] - if type(rlib[0]) == 'bytes': + if isinstance(rlib[0], bytes): rlib = [i.decode() for i in rlib] assert (flib[0] == rlib[0]) and (flib[1] == rlib[1]), liberr diff --git a/prospect/io/write_results.py b/prospect/io/write_results.py index 0ce30d8b..cf3d312c 100644 --- a/prospect/io/write_results.py +++ b/prospect/io/write_results.py @@ -279,7 +279,7 @@ def chain_to_struct(chain, model=None, names=None, **extras): struct : A structured ndarray of parameter values. """ - indict = type(chain) == dict + indict = isinstance(chain, dict) if indict: return dict_to_struct(chain) else: diff --git a/prospect/models/__init__.py b/prospect/models/__init__.py index 1fe98241..728647e8 100644 --- a/prospect/models/__init__.py +++ b/prospect/models/__init__.py @@ -6,14 +6,13 @@ """ -from .sedmodel import ProspectorParams, SpecModel -from .sedmodel import PolySpecModel, SplineSpecModel -from .sedmodel import AGNSpecModel +from .parameters import ProspectorParams +from .sedmodel import SpecModel, HyperSpecModel, AGNSpecModel __all__ = ["ProspectorParams", "SpecModel", - "PolySpecModel", "SplineSpecModel", - "LineSpecModel", "AGNSpecModel" + "HyperSpecModel", + "AGNSpecModel" ] diff --git a/prospect/models/sedmodel.py b/prospect/models/sedmodel.py index 41abbbba..b7e73477 100644 --- a/prospect/models/sedmodel.py +++ b/prospect/models/sedmodel.py @@ -22,10 +22,8 @@ __all__ = ["SpecModel", - "PolySpecModel", "SplineSpecModel", - "HyperSpecModel", "HyperPolySpecModel", - "AGNSpecModel", - "PolyFitModel"] + "HyperSpecModel", + "AGNSpecModel"] class SpecModel(ProspectorParams): @@ -103,7 +101,28 @@ def predict(self, theta, observations=None, sps=None, **extras): will be `mfrac` the ratio of the surviving stellar mass to the stellar mass formed. """ + self.predict_init(theta, sps) + # generate predictions for likelihood + # this assumes all spectral datasets (if present) occur first + # because they can change the line strengths during marginalization. + predictions = [self.predict_obs(obs) for obs in observations] + + return predictions, self._mfrac + + def predict_init(self, theta, sps): + """Generate the physical model on the model wavelength grid, and cache + many quantities used in common for all kinds of predictions. + + Parameters + ---------- + theta : ndarray of shape ``(ndim,)`` + Vector of free model parameter values. + + sps : + An `sps` object to be used in the model generation. It must have + the :py:func:`get_galaxy_spectrum` method defined. + """ # generate and cache intrinsic model spectrum and info self.set_parameters(theta) self._wave, self._spec, self._mfrac = sps.get_galaxy_spectrum(**self.params) @@ -131,13 +150,6 @@ def predict(self, theta, observations=None, sps=None, **extras): self._smooth_spec = self.add_dla(self._wave, self._smooth_spec) self._smooth_spec = self.add_damping_wing(self._wave, self._smooth_spec) - # generate predictions for likelihood - # this assumes all spectral datasets (if present) occur first - # because they can change the line strengths during marginalization. - predictions = [self.predict_obs(obs) for obs in observations] - - return predictions, self._mfrac - def predict_obs(self, obs): if obs.kind == "spectrum": prediction = self.predict_spec(obs) @@ -205,7 +217,7 @@ def predict_intrinsic(self, obs_dummy, continuum_only=True, **extras): return inst_spec - def predict_spec(self, obs, **extras): + def predict_spec(self, obs): """Generate a prediction for the observed spectrum. This method assumes that the parameters have been set and that the following attributes are present and correct @@ -253,6 +265,7 @@ def predict_spec(self, obs, **extras): obs_wave = self.observed_wave(self._wave, do_wavecal=False) # get output wavelength vector + # TODO: remove this and require all Spectrum instances to have a wavelength array self._outwave = obs.wavelength if self._outwave is None: self._outwave = obs_wave @@ -285,13 +298,19 @@ def predict_spec(self, obs, **extras): inst_spec[emask] += self._fix_eline_spec.sum(axis=1) # --- (de-) apply calibration --- - self._speccal = self.spec_calibration(obs=obs, spec=inst_spec, **extras) - inst_spec = inst_spec * self._speccal + extra_mask = self._fit_eline_pixelmask + if not extra_mask.any(): + extra_mask = True # all pixels are ok + response = obs.compute_response(spec=inst_spec, + extra_mask=extra_mask, + **self.params) + inst_spec = inst_spec * response # --- fit and add lines if necessary --- emask = self._fit_eline_pixelmask if emask.any(): - # We need the spectroscopic covariance matrix to do emission line optimization and marginalization + # We need the spectroscopic covariance matrix to do emission line + # optimization and marginalization spec_unc = None # FIXME: do this only if the noise model is non-trivial, and make sure masking is consistent #vectors = obs.noise.populate_vectors(obs) @@ -300,7 +319,8 @@ def predict_spec(self, obs, **extras): inst_spec[emask] += self._fit_eline_spec.sum(axis=1) # --- cache intrinsic spectrum for this observation --- - self._sed = inst_spec / self._speccal + self._sed = inst_spec / response + self._speccal = response return inst_spec @@ -633,6 +653,7 @@ def fit_mle_elines(self, obs, calibrated_spec, sigma_spec=None): # generate line amplitudes in observed flux units units_factor = self.flux_norm() / (1 + self._zred) + # FIXME: use obs.response instead of _speccal, remove all references to speccal calib_factor = np.interp(self._ewave_obs[idx], nebwave, self._speccal[emask]) linecal = units_factor * calib_factor alpha_breve = self._eline_lum[idx] * linecal @@ -804,9 +825,6 @@ def wave_to_x(self, wavelength=None, mask=slice(None), **extras): x = 2.0 * (x / (x[mask]).max()) - 1.0 return x - def spec_calibration(self, **kwargs): - return np.ones_like(self._outwave) - def absolute_rest_maggies(self, filterset): """Return absolute rest-frame maggies (=10**(-0.4*M)) of the last computed spectrum. @@ -842,160 +860,6 @@ def absolute_rest_maggies(self, filterset): return abs_rest_maggies -class PolySpecModel(SpecModel): - - """This is a subclass of *SpecModel* that generates the multiplicative - calibration vector at each model `predict` call as the maximum likelihood - chebyshev polynomial describing the ratio between the observed and the model - spectrum. - """ - - def _available_parameters(self): - # These should both be attached to the Observation instance as attributes - pars = [("polynomial_order", "order of the polynomial to fit"), - ("polynomial_regularization", "vector of length `polyorder` providing regularization for each polynomial term"), - ("median_polynomial", "if > 0, median smooth with a kernel of width order/range/median_polynomial before fitting") - ] - - return pars - - def spec_calibration(self, theta=None, obs=None, spec=None, **kwargs): - """Implements a Chebyshev polynomial calibration model. This uses - least-squares to find the maximum-likelihood Chebyshev polynomial of a - certain order describing the ratio of the observed spectrum to the model - spectrum, conditional on all other parameters. If emission lines are - being marginalized out, they are excluded from the least-squares fit. - - Parameters - ---------- - obs : Instance of `Spectrum` - - spec : ndarray of shape (nwave,) - The model spectrum. - - Returns - ------- - cal : ndarray of shape (nwave,) - A polynomial given by :math:`\sum_{m=0}^M a_{m} * T_m(x)`. - """ - if theta is not None: - self.set_parameters(theta) - - # polynomial order and regularization, from the obs object or overridden from the model params - order = np.squeeze(getattr(obs, "polynomial_order", 0)) - order = np.squeeze(self.params.get('polyorder', order)) - reg = np.squeeze(getattr(obs, "polynomial_regularization", 0.)) - reg = self.params.get('poly_regularization', reg) - - polyopt = ((order > 0) & - (obs.get('spectrum', None) is not None)) - if polyopt: - # generate mask - # remove region around emission lines if doing analytical marginalization - mask = obs.get('mask', np.ones_like(obs['wavelength'], dtype=bool)).copy() - if self.params.get('marginalize_elines', False): - mask[self._fit_eline_pixelmask] = 0 - - # map unmasked wavelengths to the interval -1, 1 - # masked wavelengths may have x>1, x<-1 - x = self.wave_to_x(obs["wavelength"], mask) - y = (obs['spectrum'] / spec)[mask] - 1.0 - - if self.params.get('median_polynomial', 0) > 0: - kernel_factor = self.params["median_polynomial"] - knl = int((x.max() - x.min()) / order / kernel_factor) - knl += int((knl % 2) == 0) - y = medfilt(y, knl) - yerr = (obs['unc'] / spec)[mask] - yvar = yerr**2 - A = chebvander(x[mask], order) - ATA = np.dot(A.T, A / yvar[:, None]) - if np.any(reg > 0): - ATA += reg**2 * np.eye(order+1) - ATAinv = np.linalg.inv(ATA) - c = np.dot(ATAinv, np.dot(A.T, y / yvar)) - Afull = chebvander(x, order) - poly = np.dot(Afull, c) - self._poly_coeffs = c - else: - poly = np.zeros_like(self._outwave) - - return (1.0 + poly) - - -class SplineSpecModel(SpecModel): - - """This is a subclass of *SpecModel* that generates the multiplicative - calibration vector at each model `predict` call as the maximum likelihood - cubic spline describing the ratio between the observed and the model - spectrum. - """ - - def _available_parameters(self): - pars = [("spline_knot_wave", "vector of wavelengths for the location of the spline knots"), - ("spline_knot_spacing", "spacing between knots, in units of wavelength"), - ("spline_knot_n", "number of interior knots between minimum and maximum unmasked wavelength") - ] - - return pars - - def spec_calibration(self, theta=None, obs=None, spec=None, **kwargs): - """Implements a spline calibration model. This fits a cubic spline with - determined knot locations to the ratio of the observed spectrum to the - model spectrum. If emission lines are being marginalized out, they are - excluded from the least-squares fit. - - The knot locations must be specified as model parameters, either - explicitly or via a number of knots or knot spacing (in angstroms) - """ - if theta is not None: - self.set_parameters(theta) - - splineopt = True - if splineopt: - mask, (wlo, whi) = self.obs_to_mask(obs) - y = (obs['spectrum'] / spec)[mask] - 1.0 - yerr = (obs['unc'] / spec)[mask] # HACK - knots_x = self.make_knots(wlo, whi, as_x=True, **self.params) - x = 2.0 * (obs["wavelength"] - wlo) / (whi - wlo) - 1.0 - tck = splrep(x[mask], y[mask], w=1/yerr[mask], k=3, task=-1, t=knots_x) - self._spline_coeffs = tck - spl = BSpline(*tck) - spline = spl(x) - else: - spline = np.zeros_like(self._outwave) - return (1.0 + spline) - - def make_knots(self, wlo, whi, as_x=True, **params): - """ - """ - if "spline_knot_wave" in params: - knots = np.squeeze(params["spline_knot_wave"]) - elif "spline_knot_spacing" in params: - s = np.squeeze(params["spline_knot_spacing"]) - # we drop the start so knots are interior - knots = np.arange(wlo, whi, s)[1:] - elif "spline_knot_n" in params: - n = np.squeeze(params["spline_knot_n"]) - # we need to drop the endpoints so knots are interior - knots = np.linspace(wlo, whi, n)[1:-1] - else: - raise KeyError("No valid spline knot specification in self.params") - - if as_x: - knots = 2.0 * (knots - wlo) / (whi - wlo) - 1.0 - - return knots - - def obs_to_mask(self, obs): - mask = obs.get('mask', np.ones_like(obs['wavelength'], dtype=bool)).copy() - if self.params.get('marginalize_elines', False): - mask[self._fit_eline_pixelmask] = 0 - w = obs["wavelength"] - wrange = w[mask].min(), w[mask].max() - return mask, wrange - - class AGNSpecModel(SpecModel): # TODO: simplify this to use SpecModel methods @@ -1033,7 +897,7 @@ def init_aline_info(self): assert np.abs(self.emline_info["wave"][59] - 4863) < 2 self._aline_lum[ainds] = afluxes - def predict_spec(self, obs, **extras): + def predict_spec(self, obs): """Generate a prediction for the observed spectrum. This method assumes that the parameters have been set and that the following attributes are present and correct @@ -1109,13 +973,14 @@ def predict_spec(self, obs, **extras): smooth_spec[emask] += self._agn_eline_spec # --- calibration --- - self._speccal = self.spec_calibration(obs=obs, spec=smooth_spec, **extras) - calibrated_spec = smooth_spec * self._speccal + response = obs.compute_response(spec=smooth_spec, **self.params) + inst_spec = smooth_spec * response # --- cache intrinsic spectrum --- - self._sed = calibrated_spec / self._speccal + self._sed = inst_spec / response + self._speccal = response - return calibrated_spec + return inst_spec def predict_lines(self, obs, **extras): """Generate a prediction for the observed nebular line fluxes, including @@ -1199,54 +1064,9 @@ def predict_aline_spec(self, line_indices, wave): return aline_spec -class PolyFitModel(SpecModel): - - """This is a subclass of :py:class:`SpecModel` that generates the - multiplicative calibration vector as a Chebyshev polynomial described by the - ``'poly_coeffs'`` parameter of the model, which may be free (fittable) - """ - - def spec_calibration(self, theta=None, obs=None, **kwargs): - """Implements a Chebyshev polynomial calibration model. This only - occurs if ``"poly_coeffs"`` is present in the :py:attr:`params` - dictionary, otherwise the value of ``params["spec_norm"]`` is returned. - - :param theta: (optional) - If given, set :py:attr:`params` using this vector before - calculating the calibration polynomial. ndarray of shape - ``(ndim,)`` - - :param obs: - A dictionary of observational data, must contain the key - ``"wavelength"`` - - :returns cal: - If ``params["cal_type"]`` is ``"poly"``, a polynomial given by - :math:`\times (\Sum_{m=0}^M```'poly_coeffs'[m]``:math:` \times T_n(x))`. - """ - if theta is not None: - self.set_parameters(theta) - - if ('poly_coeffs' in self.params): - mask = obs.get('mask', slice(None)) - # map unmasked wavelengths to the interval -1, 1 - # masked wavelengths may have x>1, x<-1 - x = self.wave_to_x(obs["wavelength"], mask) - # get coefficients. - c = self.params['poly_coeffs'] - poly = chebval(x, c) - return poly - else: - return 1.0 - - class HyperSpecModel(ProspectorHyperParams, SpecModel): pass -class HyperPolySpecModel(ProspectorHyperParams, PolySpecModel): - pass - - def ln_mvn(x, mean=None, cov=None): """Calculates the natural logarithm of the multivariate normal PDF diff --git a/prospect/observation/__init__.py b/prospect/observation/__init__.py index 282e93c6..46bf1ef4 100644 --- a/prospect/observation/__init__.py +++ b/prospect/observation/__init__.py @@ -1,10 +1,13 @@ # -*- coding: utf-8 -*- from .observation import Observation -from .observation import Photometry, Spectrum, Lines, UndersampledSpectrum, IntrinsicSpectrum +from .observation import Photometry, Spectrum, Lines +from .observation import UndersampledSpectrum, IntrinsicSpectrum +from .observation import PolyOptCal, SplineOptCal from .observation import from_oldstyle, from_serial __all__ = ["Observation", "Photometry", "Spectrum", "Lines", "UndersampledSpectrum", "InstrinsicSpectrum", + "PolyOptCal", "SplineOptCal", "from_oldstyle", "from_serial"] diff --git a/prospect/observation/observation.py b/prospect/observation/observation.py index 1d0cd5e5..60f916c3 100644 --- a/prospect/observation/observation.py +++ b/prospect/observation/observation.py @@ -3,6 +3,10 @@ import json import numpy as np +from numpy.polynomial.chebyshev import chebval, chebvander +from scipy.interpolate import splrep, BSpline +from scipy.signal import medfilt + from sedpy.observate import FilterSet from sedpy.smoothing import smooth_fft from sedpy.observate import rebin @@ -13,6 +17,7 @@ __all__ = ["Observation", "Spectrum", "Photometry", "Lines", "UndersampledSpectrum", "IntrinsicSpectrum", + "SplineOptCal", "PolyOptCal", "from_oldstyle", "from_serial", "obstypes"] @@ -41,7 +46,7 @@ class Observation: noise : """ - kind = "observation" + _kind = "observation" logify_spectrum = False alias = {} _meta = ("kind", "name") @@ -288,12 +293,12 @@ class Spectrum(Observation): _meta = ("kind", "name", "lambda_pad") _data = ("wavelength", "flux", "uncertainty", "mask", - "resolution", "calibration") + "resolution", "response") def __init__(self, wavelength=None, resolution=None, - calibration=None, + response=None, name=None, lambda_pad=100, polynomial_order=0., @@ -322,11 +327,9 @@ def __init__(self, super(Spectrum, self).__init__(name=name, **kwargs) self.lambda_pad = lambda_pad self.resolution = resolution - self.calibration = calibration + self.response = response self.instrument_smoothing_parameters = dict(smoothtype="vel", fftsmooth=True) self.wavelength = wavelength - self.polynomial_order = polynomial_order - self.polynomial_regularization = 0. @property def wavelength(self): @@ -437,25 +440,11 @@ def instrumental_smoothing(self, model_wave_obsframe, model_flux, zred=0, libres return outspec_padded[self._unpadded_inds] -class UndersampledSpectrum(Spectrum): - """ - As for Spectrum, but account for pixelization effects when pixels - undersamplesthe instrumental LSF. - """ - - def _pixelize(self, outwave, inwave, influx): - return rebin(outwave, inwave, influx) - - -class IntrinsicSpectrum(Spectrum): - - """ - As for Spectrum, but predictions for this object type will not include - polynomial fitting or fitting of the emission line strengths (previously fit - and cached emission luminosities will be used.) - """ - - _kind = "intrinsic" + def compute_response(self, **extras): + if self.response is not None: + return self.response + else: + return 1.0 class Lines(Spectrum): @@ -473,6 +462,7 @@ class Lines(Spectrum): def __init__(self, line_ind=None, + line_names=None, name=None, **kwargs): @@ -496,12 +486,256 @@ def __init__(self, dispersion (:math:`= c \, \sigma_\lambda / \lambda = c \, \FWHM_\lambda / 2.355 / \lambda = c / (2.355 \, R_\lambda)` where :math:`c=2.998e5 {\rm km}/{\rm s}` + line_ind : iterable of string, optional + Names of the lines. + :param calibration: not sure yet .... """ super(Lines, self).__init__(name=name, resolution=None, **kwargs) assert (line_ind is not None), "You must identify the lines by their index in the FSPS emission line array" self.line_ind = np.array(line_ind).as_type(int) + self.line_names = line_names + + +class UndersampledSpectrum(Spectrum): + """ + As for Spectrum, but account for pixelization effects when pixels + undersamplesthe instrumental LSF. + """ + #TODO: Implement as a convolution with a square kernel (or sinc in frequency space) + + def _pixelize(self, outwave, inwave, influx): + return rebin(outwave, inwave, influx) + + +class IntrinsicSpectrum(Spectrum): + + """ + As for Spectrum, but predictions for this object type will not include + polynomial fitting or fitting of the emission line strengths (previously fit + and cached emission luminosities will be used.) + """ + + _kind = "intrinsic" + + +class PolyOptCal: + + """A mixin class that allows for optimization of a Chebyshev response + function given a model spectrum. + """ + + def __init__(self, *args, + polynomial_order=0, + polynomial_regularization=0, + median_polynomial=0, + **kwargs): + super(PolyOptCal, self).__init__(*args, **kwargs) + self.polynomial_order = polynomial_order + self.polynomial_regularization = polynomial_regularization + self.median_polynomial = median_polynomial + + def _available_parameters(self): + # These should both be attached to the Observation instance as attributes + pars = [("polynomial_order", "order of the polynomial to fit"), + ("polynomial_regularization", "vector of length `polyorder` providing regularization for each polynomial term"), + ("median_polynomial", "if > 0, median smooth with a kernel of width order/range/median_polynomial before fitting") + ] + + return pars + + def compute_response(self, spec=None, extra_mask=True, **kwargs): + """Implements a Chebyshev polynomial response function model. This uses + least-squares to find the maximum-likelihood Chebyshev polynomial of a + certain order describing the ratio of the observed spectrum to the model + spectrum. If emission lines are being marginalized out, they should be + excluded from the least-squares fit using the ``extra_mask`` keyword + + Parameters + ---------- + spec : ndarray of shape (Spectrum().ndata,) + The model spectrum on the observed wavelength grid + + extra_mask : ndarray of Booleans of shape (Spectrum().ndata,) + The extra mask to be applied. True=use data, False=mask data + + Returns + ------- + response : ndarray of shape (nwave,) + A polynomial given by :math:`\sum_{m=0}^M a_{m} * T_m(x)`. + """ + + order = self.polynomial_order + reg = self.polynomial_regularization + mask = self.mask & extra_mask + assert (self.mask.sum() > order), f"Not enough points to constrain polynomial of order {order}" + + polyopt = (order > 0) + if (not polyopt): + print("no polynomial") + self.response = np.ones_like(self.wavelength) + return self.response + + x = wave_to_x(self.wavelength, mask) + y = (self.flux / spec - 1.0)[mask] + yerr = (self.uncertainty / spec)[mask] + yvar = yerr**2 + + if self.median_polynomial > 0: + kernel_factor = self.median_polynomial + knl = int((x.max() - x.min()) / order / kernel_factor) + knl += int((knl % 2) == 0) + y = medfilt(y, knl) + + Afull = chebvander(x, order) + A = Afull[mask, :] + ATA = np.dot(A.T, A / yvar[:, None]) + if np.any(reg > 0): + ATA += reg**2 * np.eye(order+1) + c = np.linalg.solve(ATA, np.dot(A.T, y / yvar)) + + poly = np.dot(Afull, c) + + self._chebyshev_coefficients = c + self.response = poly + 1.0 + return self.response + + +class SplineOptCal: + + """A mixin class that allows for optimization of a Chebyshev response + function given a model spectrum. + """ + + + def __init__(self, *args, + spline_knot_wave=None, + spline_knot_spacing=None, + spline_knot_n=None, + **kwargs): + super(SplineOptCal, self).__init__(*args, **kwargs) + + self.params = {} + if spline_knot_wave is not None: + self.params["spline_knot_wave"] = spline_knot_wave + elif spline_knot_spacing is not None: + self.params["spline_knot_spacing"] = spline_knot_spacing + elif spline_knot_n is not None: + self.params["spline_knot_n"] = spline_knot_n + + # build and cache the knots + w = self.wavelength[self.mask] + (wlo, whi) = w.min(), w.min() + self.wave_x = 2.0 * (self.wavelength - wlo) / (whi - wlo) - 1.0 + self.knots_x = self.make_knots(wlo, whi, as_x=True, **self.params) + + + def _available_parameters(self): + pars = [("spline_knot_wave", "vector of wavelengths for the location of the spline knots"), + ("spline_knot_spacing", "spacing between knots, in units of wavelength"), + ("spline_knot_n", "number of interior knots between minimum and maximum unmasked wavelength") + ] + + return pars + + def compute_response(self, spec=None, extra_mask=True, **extras): + """Implements a spline response function model. This fits a cubic + spline with determined knot locations to the ratio of the observed + spectrum to the model spectrum. If emission lines are being + marginalized out, they are excluded from the least-squares fit. + + The knot locations must be specified as model parameters, either + explicitly or via a number of knots or knot spacing (in angstroms) + """ + mask = self.mask & extra_mask + + + splineopt = True + if ~splineopt: + self.response = np.ones_like(self.wavelength) + return self.response + + y = (self.flux / spec - 1.0)[mask] + yerr = (self.uncertainty / spec)[mask] # HACK + tck = splrep(self.wave_x[mask], y[mask], w=1/yerr[mask], k=3, task=-1, t=self.knots_x) + self._spline_coeffs = tck + spl = BSpline(*tck) + spline = spl(self.wave_x) + + self.response = (1.0 + spline) + return self.response + + def make_knots(self, wlo, whi, as_x=True, **params): + """Can we move this to instantiation? + """ + if "spline_knot_wave" in params: + knots = np.squeeze(params["spline_knot_wave"]) + elif "spline_knot_spacing" in params: + s = np.squeeze(params["spline_knot_spacing"]) + # we drop the start so knots are interior + knots = np.arange(wlo, whi, s)[1:] + elif "spline_knot_n" in params: + n = np.squeeze(params["spline_knot_n"]) + # we need to drop the endpoints so knots are interior + knots = np.linspace(wlo, whi, n)[1:-1] + else: + raise KeyError("No valid spline knot specification in self.params") + + if as_x: + knots = 2.0 * (knots - wlo) / (whi - wlo) - 1.0 + + return knots + + +class PolyFitCal: + + """This is a mixin class that generates the + multiplicative response vector as a Chebyshev polynomial described by the + ``poly_param_name`` parameter of the model, which may be free (fittable) + """ + + def __init__(self, *args, poly_param_name=None, **kwargs): + super(SplineOptCal, self).__init(*args, **kwargs) + self.poly_param_name = poly_param_name + + def _available_parameters(self): + pars = [(self.poly_param_name, "vector of polynomial chabyshev coefficients")] + + return pars + + def compute_response(self, **kwargs): + """Implements a Chebyshev polynomial calibration model. This only + occurs if ``"poly_coeffs"`` is present in the :py:attr:`params` + dictionary, otherwise the value of ``params["spec_norm"]`` is returned. + + :param theta: (optional) + If given, set :py:attr:`params` using this vector before + calculating the calibration polynomial. ndarray of shape + ``(ndim,)`` + + :param obs: + A dictionary of observational data, must contain the key + ``"wavelength"`` + + :returns cal: + If ``params["cal_type"]`` is ``"poly"``, a polynomial given by + :math:`\times (\Sum_{m=0}^M```'poly_coeffs'[m]``:math:` \times T_n(x))`. + """ + + if self.poly_param_name in kwargs: + mask = self.get('mask', slice(None)) + # map unmasked wavelengths to the interval -1, 1 + # masked wavelengths may have x>1, x<-1 + x = wave_to_x(self.wavelength, mask) + # get coefficients. + c = kwargs[self.poly_param_name] + poly = chebval(x, c) + else: + poly = 1.0 + + self.response = poly + return self.response obstypes = dict(photometry=Photometry, @@ -536,3 +770,11 @@ def from_serial(arr, meta): return obs + +def wave_to_x(wavelength=None, mask=slice(None), **extras): + """Map unmasked wavelengths to the interval -1, 1 + masked wavelengths may have x>1, x<-1 + """ + x = wavelength - (wavelength[mask]).min() + x = 2.0 * (x / (x[mask]).max()) - 1.0 + return x \ No newline at end of file diff --git a/tests/test_polycal.py b/tests/test_polycal.py new file mode 100644 index 00000000..b4ea4488 --- /dev/null +++ b/tests/test_polycal.py @@ -0,0 +1,80 @@ +#!/usr/bin/env python +# -*- coding: utf-8 -*- + +import numpy as np +import pytest + +from prospect.sources import CSPSpecBasis +from prospect.models import SpecModel, templates +from prospect.observation import Spectrum, Photometry, PolyOptCal + + +class PolySpectrum(PolyOptCal, Spectrum): + pass + + + +@pytest.fixture +def get_sps(): + sps = CSPSpecBasis(zcontinuous=1) + return sps + + +def build_model(add_neb=False): + model_params = templates.TemplateLibrary["parametric_sfh"] + if add_neb: + model_params.update(templates.TemplateLibrary["nebular"]) + return SpecModel(model_params) + + +def build_obs(multispec=False): + N = 1500 * (2 - multispec) + wmax = 7000 + wsplit = wmax - N * multispec + + fnames = list([f"sdss_{b}0" for b in "ugriz"]) + Nf = len(fnames) + phot = [Photometry(filters=fnames, + flux=np.ones(Nf), + uncertainty=np.ones(Nf)/10)] + spec = [PolySpectrum(wavelength=np.linspace(4000, wsplit, N), + flux=np.ones(N), + uncertainty=np.ones(N) / 10, + mask=slice(None), + polynomial_order=5) + ] + + if multispec: + spec += [Spectrum(wavelength=np.linspace(wsplit+1, wmax, N), + flux=np.ones(N), uncertainty=np.ones(N) / 10, + mask=slice(None))] + + obslist = spec + phot + [obs.rectify() for obs in obslist] + return obslist + + +def test_polycal(get_sps, plot=False): + """Make sure the polynomial optimization works + """ + sps = get_sps + observations = build_obs() + model = build_model() + + preds, extra = model.predict(model.theta, observations=observations, sps=sps) + obs = observations[0] + + assert np.any(obs.response != 0) + + if plot: + import matplotlib.pyplot as pl + fig, axes = pl.subplots(3, 1, sharex=True) + ax = axes[0] + ax.plot(obs.wavelength, obs.flux, label="obseved flux (ones)") + ax.plot(obs.wavelength, preds[0], label="model flux (times response)") + ax = axes[1] + ax.plot(obs.wavelength, obs.response, label="instrumental response (polynomial)") + ax = axes[2] + ax.plot(obs.wavelength, preds[0]/ obs.response, label="intrinsic model spectrum") + ax.set_xlabel("wavelength") + [ax.legend() for ax in axes] \ No newline at end of file