fix tests; start on multispec tests.

tests/test_eline.py

@@ -59,16 +59,23 @@ def test_eline_parsing():
     assert model._fit_eline.sum() == (len(model._use_eline) - len(fix_lines))
 
 
-def test_nebline_phot_addition(get_sps):
-    fnames = [f"sdss_{b}0" for b in "ugriz"]
-    filts = observate.load_filters(fnames)
-
+def build_obs(filts):
     obs = dict(filters=filts,
                wavelength=np.linspace(3000, 9000, 1000),
                spectrum=np.ones(1000),
-               unc=np.ones(1000)*0.1)
+               unc=np.ones(1000)*0.1,
+               maggies=np.ones(len(filts))*1e-7,
+               maggies_unc=np.ones(len(filts))*1e-8)
     sdat, pdat = from_oldstyle(obs)
     obslist = [sdat, pdat]
+    [obs.rectify() for obs in obslist]
+    return obslist
+
+
+def test_nebline_phot_addition(get_sps):
+    fnames = [f"sdss_{b}0" for b in "ugriz"]
+    filts = observate.load_filters(fnames)
+    obslist = build_obs(filts)
 
     sps = get_sps
 
@@ -104,13 +111,8 @@ def test_filtersets(get_sps):
     """
     fnames = [f"sdss_{b}0" for b in "ugriz"]
     flist = observate.load_filters(fnames)
-
-    obs = dict(wavelength=np.linspace(3000, 9000, 1000),
-               spectrum=np.ones(1000),
-               unc=np.ones(1000)*0.1,
-               filters=fnames)
-    sdat, pdat = from_oldstyle(obs)
-    obslist = [sdat, pdat]
+    obslist = build_obs(flist)
+    sdat, pdat = obslist
 
     sps = get_sps
 
@@ -149,13 +151,7 @@ def test_eline_implementation(get_sps):
     test_eline_parsing()
 
     filters = observate.load_filters([f"sdss_{b}0" for b in "ugriz"])
-    obs = dict(filters=filters,
-               wavelength=np.linspace(3000, 9000, 1000),
-               spectrum=np.ones(1000),
-               unc=np.ones(1000)*0.1,
-               maggies=np.ones(len(filters))*1e-7,
-               maggies_unc=np.ones(len(filters))*1e-8)
-    obslist = from_oldstyle(obs)
+    obslist = build_obs(filters)
 
     model_pars = TemplateLibrary["parametric_sfh"]
     model_pars.update(TemplateLibrary["nebular"])

tests/test_lnlike.py

@@ -0,0 +1,90 @@
+#!/usr/bin/env python
+# -*- coding: utf-8 -*-
+
+import sys
+import numpy as np
+
+import pytest
+
+from prospect.sources import CSPSpecBasis
+from prospect.models import SpecModel, templates
+from prospect.observation import Spectrum, Photometry
+from prospect.likelihood import NoiseModel
+from prospect.likelihood.likelihood import compute_lnlike
+from prospect.fitting import lnprobfn
+
+
+@pytest.fixture
+def get_sps():
+    sps = CSPSpecBasis(zcontinuous=1)
+    return sps
+
+
+def build_model(add_neb=False, add_outlier=False):
+    model_params = templates.TemplateLibrary["parametric_sfh"]
+    if add_neb:
+        model_params.update(templates.TemplateLibrary["nebular"])
+    if add_outlier:
+        model_params.update(templates.TemplateLibrary["outlier_model"])
+        model_params["f_outlier_phot"]["isfree"] = True
+        model_params["f_outlier_phot"]["init"] = 0.05
+    return SpecModel(model_params)
+
+
+def build_obs(multispec=True, add_outlier=True):
+    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 = [Spectrum(wavelength=np.linspace(4000, wsplit, N),
+                     flux=np.ones(N), uncertainty=np.ones(N) / 10,
+                     mask=slice(None))]
+
+    if add_outlier:
+        phot[0].noise = NoiseModel(frac_out_name='f_outlier_phot',
+                                   nsigma_out_name='nsigma_outlier_phot')
+
+    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_lnlike_shape(get_sps):
+    # testing lnprobfn
+    sps = get_sps
+
+    for add_out in [True, False]:
+        observations = build_obs(add_outlier=add_out)
+        model = build_model(add_neb=add_out, add_outlier=add_out)
+
+        model.set_parameters(model.theta)
+        [obs.noise.update(**model.params) for obs in observations
+        if obs.noise is not None]
+        predictions, x = model.predict(model.theta, observations, sps=sps)
+
+        # check you get a scalar lnp for each observation
+        lnp_data = [compute_lnlike(pred, obs, vectors={}) for pred, obs
+                    in zip(predictions, observations)]
+        assert np.all([np.isscalar(p) for p in lnp_data]), f"failed for add_outlier={add_out}"
+        assert len(lnp_data) == len(observations), f"failed for add_outlier={add_out}"
+
+        # check lnprobfn returns scalar
+        lnp = lnprobfn(model.theta, model=model, observations=observations, sps=sps)
+
+        assert np.isscalar(lnp), f"failed for add_outlier={add_out}"
+
+    # %timeit model.prior_product(model.theta)
+    # arr = np.zeros(model.ndim)
+    # arr[-1] = 1
+    # theta = model.theta.copy()
+    # %timeit predictions, x = model.predict(theta + np.random.uniform(-0.1, 0.1) * arr, observations=observations, sps=sps)
+    # %timeit lnp_data = [compute_lnlike(pred, obs, vectors={}) for pred, obs in zip(predictions, observations)]
+    # %timeit lnp = lnprobfn(theta + np.random.uniform(0, 3) * arr, model=model, observations=observations, sps=sps)

tests/test_multispec.py

@@ -0,0 +1,102 @@
+#!/usr/bin/env python
+# -*- coding: utf-8 -*-
+
+import sys
+import numpy as np
+
+import pytest
+
+from sedpy.observate import load_filters
+from prospect.sources import CSPSpecBasis
+from prospect.models import SpecModel, templates
+from prospect.observation import Spectrum, Photometry
+
+
+@pytest.fixture(scope="module")
+def build_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=True):
+    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 = [Spectrum(wavelength=np.linspace(4000, wsplit, N),
+                     flux=np.ones(N), uncertainty=np.ones(N) / 10,
+                     mask=slice(None))]
+
+    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_prediction_nodata(build_sps):
+    sps = build_sps
+    model = build_model(add_neb=True)
+    sobs, pobs = build_obs(multispec=False)
+    pobs.flux = None
+    pobs.uncertainty = None
+    sobs.wavelength = None
+    sobs.flux = None
+    sobs.uncertainty = None
+    pred, mfrac = model.predict(model.theta, observations=[sobs, pobs], sps=sps)
+    assert len(pred[0]) == len(sps.wavelengths)
+    assert len(pred[1]) == len(pobs.filterset)
+
+
+def test_multispec(build_sps):
+    sps = build_sps
+
+    obslist_single = build_obs(multispec=False)
+    obslist_multi = build_obs(multispec=True)
+    model = build_model(add_neb=True)
+
+    preds_single, mfrac = model.predict(model.theta, observations=obslist_single, sps=sps)
+    preds_multi, mfrac = model.predict(model.theta, observations=obslist_multi, sps=sps)
+
+    assert len(preds_single) == 2
+    assert len(preds_multi) == 3
+    assert np.allclose(preds_single[-1], preds_multi[-1])
+
+    # TODO: turn this plot into an actual test
+    #import matplotlib.pyplot as pl
+    #fig, ax = pl.subplots()
+    #ax.plot(obslist_single[0].wavelength, predictions_single[0])
+    #for p, o in zip(predictions, obslist):
+    #    if o.kind == "photometry":
+    #        ax.plot(o.wavelength, p, "o")
+    #    else:
+    #        ax.plot(o.wavelength, p)
+
+
+def test_multires():
+    # Test the smoothing of multiple spectra to different resolutions
+    # - give the same wavelength array different instrumental resolutions, assert similar but different, and that smoothing by the difference gives the right answer
+    # Test the use of two differernt smoothings, physical and instrumental
+    # - give an obs with no instrument smoothing and one with, make sure they are different
+    pass
+
+
+def test_multinoise():
+    pass
+
+
+def test_multical():
+    pass