Docs and test for simulate observations

beast/observationmodel/observations.py

@@ -8,6 +8,8 @@
 
 from astropy.table import Table, Column
 
+from beast.observationmodel.vega import Vega
+
 __all__ = ['Observations', 'gen_SimObs_from_sedgrid']
 
 
@@ -152,7 +154,8 @@ def enumobs(self):
 
 
 def gen_SimObs_from_sedgrid(sedgrid, sedgrid_noisemodel,
-                            nsim=100, compl_filter='F475W'):
+                            nsim=100, compl_filter='F475W',
+                            ranseed=None, vega_fname=None):
     """
     Generate simulated observations using the physics and observation grids.
     The priors are sampled as they give the ensemble model for the stellar
@@ -177,6 +180,14 @@ def gen_SimObs_from_sedgrid(sedgrid, sedgrid_noisemodel,
     compl_filter : str
         filter to use for completeness (required for toothpick model)
 
+    ranseed : int
+        used to set the seed to make the results reproducable
+        useful for testing
+
+    vega_fname : string
+        filename for the vega info
+        usefule for testing
+
     Outputs
     -------
     simtable : astropy Table
@@ -211,9 +222,17 @@ def gen_SimObs_from_sedgrid(sedgrid, sedgrid_noisemodel,
     # need to sum to 1
     gridweights = gridweights/np.sum(gridweights)
 
+    # set the random seed - mainly for testing
+    if not None:
+        np.random.seed(ranseed)
+
     # sample to get the indexes of the picked models
     indx = range(n_models)
     sim_indx = np.random.choice(indx, size=nsim, p=gridweights)
+    print(sim_indx)
+
+    # get the vega fluxes for the filters
+    _, vega_flux, _ = Vega(source=vega_fname).getFlux(sedgrid.filters)
 
     # setup the output table
     ot = Table()
@@ -224,7 +243,7 @@ def gen_SimObs_from_sedgrid(sedgrid, sedgrid_noisemodel,
         simflux_wbias = flux[sim_indx, k] + model_bias[sim_indx, k]
         simflux = np.random.normal(loc=simflux_wbias,
                                    scale=model_unc[sim_indx, k])
-        ot[colname] = Column(simflux)
+        ot[colname] = Column(simflux/vega_flux[k])
     # model parmaeters
     for qname in qnames:
         ot[qname] = Column(sedgrid[qname][sim_indx])

beast/observationmodel/tests/test_simobs.py

@@ -0,0 +1,60 @@
+import os.path
+
+import numpy as np
+
+from astropy.utils.data import download_file
+from astropy.tests.helper import remote_data
+from astropy.table import Table
+
+from ..noisemodel import generic_noisemodel as noisemodel
+from ...physicsmodel.grid import FileSEDGrid
+from beast.observationmodel.observations import gen_SimObs_from_sedgrid
+
+
+def _download_rename(filename):
+    """
+    Download a file and rename it to have the right extension
+
+    Otherwise, downloaded file will not have an extension at all
+    """
+    url_loc = 'http://www.stsci.edu/~kgordon/beast/'
+    fname_dld = download_file('%s%s' % (url_loc, filename))
+    extension = filename.split('.')[-1]
+    fname = '%s.%s' % (fname_dld, extension)
+    os.rename(fname_dld, fname)
+    return fname
+
+
+@remote_data
+def test_trim_grid():
+
+    # download the needed files
+    vega_fname = _download_rename('vega.hd5')
+    seds_fname = _download_rename('beast_example_phat_seds.grid.hd5')
+    noise_fname = _download_rename('beast_example_phat_noisemodel.hd5')
+
+    # download cached version of noisemodel on the sed grid
+    simobs_fname_cache = _download_rename('beast_example_phat_simobs.fits')
+
+    ################
+
+    # get the physics model grid - includes priors
+    modelsedgrid = FileSEDGrid(seds_fname)
+
+    # read in the noise model - includes bias, unc, and completeness
+    noisegrid = noisemodel.get_noisemodelcat(noise_fname)
+
+    table_new = gen_SimObs_from_sedgrid(modelsedgrid, noisegrid,
+                                        nsim=100,
+                                        compl_filter="f475w",
+                                        ranseed=1234,
+                                        vega_fname=vega_fname)
+
+    # check that the simobs files are exactly the same
+    table_cache = Table.read(simobs_fname_cache)
+
+    assert len(table_new) == len(table_cache)
+
+    for tcolname in table_new.colnames:
+        np.testing.assert_equal(table_new[tcolname], table_cache[tcolname],
+                                '%s columns not equal' % tcolname)

beast/tools/simulate_obs.py

@@ -1,7 +1,4 @@
 import argparse
-# import numpy as np
-
-# from astropy.table import Table, Column
 
 from beast.physicsmodel.grid import FileSEDGrid
 import beast.observationmodel.noisemodel.generic_noisemodel as noisemodel

docs/index.rst

@@ -23,6 +23,7 @@ User Documentation
    Example production run workflow <workflow.rst>
    Running in parallel by using subgrids <subgrid_parallelism.rst>
    Generating AST inputs <generating_asts.rst>
+   Simulating observations <simulations.rst>
    Format of BEAST grid files <beast_grid_format.rst>
 
 Installation

docs/simulations.rst

@@ -0,0 +1,48 @@
+###########
+Simulations
+###########
+
+Simulations of observations are useful for testing the sensitivity
+of a specific set of observations to BEAST model parameters.  This is
+useful for quantifying the BEAST performance for proposed and actual
+observations.
+
+This is done using the
+`beast.observationmodel.observations.gen_SimObs_from_sedgrid` function.
+The script
+`tools/simulate_obs.py` provides a commandline interface.  The module
+uses already created BEAST physics and observation model grids
+by sampling the full nD prior function that is part of the physics
+model grid.  The observation model grid provides the information on
+the flux uncertainties and biases as well as the completeness.
+
+*********
+Toothpick
+*********
+
+The files for the physicsgrid and obsgrid files are required inputs to
+the script.  The output filename is also required.  Note that the extension
+of this file will determine the type of file output (e.g. filebase.fits for
+a FITS file).
+The number of observations to simulate is given by the `--nsim` parameter.
+The filter to use for the completeness function is given by the
+`--compl_filter` parameter.
+
+.. code:: shell
+
+   $ python simulate.py physicsgrid obsgrid outfile \
+                --nsim 200 --compl_filter f475w
+
+The output file gives the simulated data in the observed data columns
+identified in the physicsgrid file along with all the model parameters
+from the physicsgrid file.  The names of the observed columns are
+`band_rate` and the units are normalized Vega fluxes (to match how
+the observed data are given).
+
+********
+Trunchen
+********
+
+The code does not handle the trunchen model at this point.  It is
+straightforward to extend the code for this model, but it has not
+been done yet.