Merge pull request #124 from lcjohnso/adjust_priors_lcj

Adjust Av and Age priors, allow verify_params failures

Failures are understood and have to do with the development version of astropy

beast/physicsmodel/grid_and_prior_weights.py

@@ -28,7 +28,7 @@
 __all__ = ['compute_age_mass_metallicity_weights',
            'compute_bin_boundaries']
 
-def compute_age_mass_metallicity_weights(_tgrid):
+def compute_age_mass_metallicity_weights(_tgrid, **kwargs):
     """ 
     Computes the age-mass-metallicity grid and prior weights 
     on the BEAST model spectra grid
@@ -62,7 +62,7 @@ def compute_age_mass_metallicity_weights(_tgrid):
         uniq_ages = np.unique(_tgrid[zindxs]['logA']) 
 
         # compute the age weights
-        age_grid_weights = compute_age_grid_weights(uniq_ages)  
+        age_grid_weights = compute_age_grid_weights(uniq_ages, **kwargs)  
         age_prior_weights = compute_age_prior_weights(uniq_ages)  
 
         for ak, age_val in enumerate(uniq_ages):

beast/physicsmodel/grid_weights.py

@@ -47,31 +47,39 @@ def compute_bin_boundaries(tab):
     tab2[1:-1] = temp
     return tab2
 
-def compute_age_grid_weights(logages):
-    """ Computes the age weights to provide constant star formation rate
-    (in linear age)
+def compute_age_grid_weights(logages, constantSFR=True):
+    """ Computes the age weights to set prior on star formation history
 
     Keywords
     --------
     logages : numpy vector
        log(ages)
 
+    constantSFR : boolean
+       Sets assumption of star formation history: flat in log or linear age
+       Default = True (constant in linear age)
+
     Returns
     -------
     age_weights : numpy vector
        total masses at each age for a constant SFR in linear age
     """
-    # ages need to be monotonically increasing
-    aindxs = np.argsort(logages)   
-
-    # Computes the bin boundaries in log
-    logages_bounds = compute_bin_boundaries(logages[aindxs])    
-
-    # initialize the age weights
-    age_weights = np.full(len(aindxs),0.0)
-
-    # Returns the age weight as a numpy array
-    age_weights[aindxs] = np.diff(10**(logages_bounds))
+    if constantSFR:
+        # ages need to be monotonically increasing
+        aindxs = np.argsort(logages)   
+
+        # Computes the bin boundaries in log
+        logages_bounds = compute_bin_boundaries(logages[aindxs])    
+
+        # initialize the age weights
+        age_weights = np.full(len(aindxs),0.0)
+
+        # Returns the age weight as a numpy array
+        age_weights[aindxs] = np.diff(10**(logages_bounds))
+
+    else:
+        # Returns the age weight as a numpy array
+        age_weights = np.ones(len(logages))
 
     # return in the order that logages was passed
     return age_weights    

beast/physicsmodel/model_grid.py

@@ -185,7 +185,7 @@ def add_stellar_priors(project, specgrid, verbose=True, **kwargs):
         if verbose:
             print('Make Prior Weights')
 
-        compute_age_mass_metallicity_weights(specgrid.grid)
+        compute_age_mass_metallicity_weights(specgrid.grid, **kwargs)
 
         #write to disk
         if hasattr(specgrid, 'writeHDF'):

beast/physicsmodel/prior_weights_dust.py

@@ -85,6 +85,24 @@ def _two_lognorm(xs,
     normalization = np.trapz(pointwise, x=xs)
     return pointwise / normalization
 
+def _exponential(x, a=2.0, N=1.):
+    """
+    Exponential distribution
+    Parameters
+    ----------
+    x: vector
+       x values
+    a: float
+       Decay Rate parameter in exp: N*e^-ax
+       Distribution Mean = 1/a
+    N: float
+       Multiplicative factor
+    Returns
+    -------
+    exponential computed on the x grid
+    """
+    return N * np.exp(-1. * a * x)
+
 class PriorWeightsDust():
     """
     Compute the priors as weights given the input grid
@@ -195,6 +213,7 @@ def set_av_weights(self, model={'name': 'flat'}):
           flat = flat prior on linear A(V)
           lognormal = lognormal prior on linear A(V)
           two_lognormal = two lognormal prior on linear A(V)
+          exponential = exponential prior on linear A(V)
         """
         if model['name'] == 'flat':
             self.av_priors = np.full(len(self.av_vals),1.0)
@@ -211,6 +230,10 @@ def set_av_weights(self, model={'name': 'flat'}):
                                           sigma2=model['sigma2'],
                                           N1=model['N1'],
                                           N2=model['N2'])
+        elif model['name'] == 'exponential':
+            self.av_priors = _exponential(self.av_vals,
+                                          a=model['a'],
+                                          N=model['N'])
         else:
             print('**error in setting the A(V) dust prior weights!**')
             print('**model ' + model['name'] + ' not supported**')

beast/tools/verify_params.py

@@ -13,39 +13,43 @@
 from sys import exit
 
 
-def verify_range(param, param_name, param_lim):
+def verify_range(param, param_name, param_lim, noexit=False):
     # check if input param limits make sense
 
     if any(p < param_lim[0] for p in param):
         print((param_name + ' min value not physical.'))
-        exit()
+        if not noexit: exit()
 
     if any(p > param_lim[1] for p in param):
         print((param_name + ' max value not physical.'))
-        exit()
+        if not noexit: exit()
 
-def check_grid(param, param_name, param_lim):
+def check_grid(param, param_name, param_lim, noexit=False):
     # check if input param limits and grid initialisation make sense
 
     param_min, param_max, param_step = param
 
     if param_min < param_lim[0]:
         print((param_name + ' min value not physical.'))
-        exit()
+        if not noexit: exit()
 
     if param_max > param_lim[1]:
         print((param_name + ' max value not physical.'))
-        exit()
+        if not noexit: exit()
 
     if param_min > param_max:
         print((param_name + ' min value greater than max'))
-        exit()
+        if not noexit: exit()
 
     if (param_max-param_min) < param_step:
-        print((param_name + ' step value greater than (max-min)'))
-        exit()
-
-def verify_one_input_format(param, param_name, param_format, param_lim):
+        if param_max-param_min == 0.:
+            print('Note: '+param_name+' grid is single-valued.')
+        else:
+            print((param_name + ' step value greater than (max-min)'))
+            if not noexit: exit()
+
+def verify_one_input_format(param, param_name, param_format, param_lim,
+                            noexit=False):
     """
     Test for an input parameter correctness of format and limits.
 
@@ -64,45 +68,51 @@ def verify_one_input_format(param, param_name, param_format, param_lim):
        [-inf, inf] when param is not constraint at all;
        None when concerns a str parameter
 
+    noexit: boolean
+       Override exit() commands if there is an error.
+       Default = False (exit on error)
     """
 
     if 'list' in param_format:
         if type(param) is not list:
-            print((param_name + ' is not in the right format - a list.'))
-            exit()
+            if param is None:
+                print('Warning: '+ param_name + ' is not defined.')
+            else:
+                print((param_name + ' is not in the right format - a list.'))
+                if not noexit: exit()
         elif 'float' in param_format:
             is_list_of_floats = all(type(item) is float for item in param)
             if not is_list_of_floats:
                 print((param_name +
                       ' is not in the right format - list of floats.'))
-                exit()
+                if not noexit: exit()
             elif 'grid' in param_format:
                 # when param is aranged from given [min, max, step],
                 # instead of a specific list of values
-                check_grid(param, param_name, param_lim)
+                check_grid(param, param_name, param_lim, noexit=noexit)
             else:
-                verify_range(param, param_name, param_lim)
+                verify_range(param, param_name, param_lim, noexit=noexit)
 
         if 'str' in param_format:
             if type(param) is not str:
                 print((param_name + ' is not in the right format - a string.')) 
-                exit()
+                if not noexit: exit()
             elif 'file' in param_format:
                 if not exists(param):
                     print((param_name +
                           ' does not exist. Please provide the file path.')) 
-                    exit()
+                    if not noexit: exit()
 
         if 'version' in param_format:
             if type(param) is not float:
                 print((param_name + ' is not in the right format - a float'))
-                exit()
+                if not noexit: exit()
             elif param not in param_lim:
                 print((param_name +
                     ' is an invalid number, leading to version of the isochrone.'))
-                exit()
+                if not noexit: exit()
 
-def verify_input_format(datamodel):
+def verify_input_format(datamodel, noexit=False):
 
     """
     Import BEAST input parameters from datamodel.
@@ -115,6 +125,10 @@ def verify_input_format(datamodel):
     datamodel: module
         Input parameters are initialized in datamodel
 
+    noexit: boolean
+        Override exit() commands if there is an error.
+        Default = False (exit on error)
+
     """
 
     parameters = [datamodel.z, datamodel.obsfile, \
@@ -130,7 +144,8 @@ def verify_input_format(datamodel):
 
     for i, param_ in enumerate(parameters):
         verify_one_input_format(param_, parameters_names[i],
-                                param_format[i], parameters_limits[i])
+                                param_format[i], parameters_limits[i],
+                                noexit=noexit)
 
 if __name__ == '__main__':
 

docs/beast_setup.rst

@@ -0,0 +1,73 @@
+Running the BEAST
+=================
+
+1) Define project and grid input parameters in datamodel.py
+
+2) Execute BEAST Run using ``python run_beast.py`` with appropriate task flags
+
+   * Default Full Stack Run: ``python run_beast.py -p -o -t -f``
+
+BEAST Data Model
+================
+
+Primary Parameters
+------------------
+
+Project Details
+
+* ``project``: sets name of working subdirectory
+* ``filters``: names of photometric filter passbands (matching library names)
+* ``basefilters``: short versions of passband names
+* ``obsfile``: filename for input flux data
+* ``obs_colnames``: column names in ``obsfile`` for observed fluxes
+* ``astfile``: filename for AST results file
+* ``distanceModulus``: distance modulus to target in mags
+
+Grid Definition Parameters
+
+* ``logt``: age grid range parameters (min, max, step)
+* ``z``: metallicity grid points
+* ``oiso``: stellar model definition
+* ``osl``: stellar library definition
+* ``extLaw``: extinction law definition
+* ``avs``: A_V grid range parameters (min, max, step)
+* ``rvs``: R_V grid range parameters (min, max, step)
+* ``fAs``: f_A grid range parameters (min, max, step)
+* ``*_prior_model``: prior model definitions for dust parameters (default: flat prior)
+
+Optional Features
+-----------------
+
+Add additional filters to grid
+^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
+Define list of filternames as ``additional_filters`` and alter ``add_spectral_properties`` call:
+
+``add_spectral_properties_kwargs = dict(filternames=filters + additional_filters)``
+
+Skip verify_params exit
+^^^^^^^^^^^^^^^^^^^^^^^
+Add ``noexit=True`` keyword to ``verify_input_format()`` call in run_beast.py:
+
+``verify_params.verify_input_format(datamodel, noexit=True)``
+
+Remove constant SFH prior
+^^^^^^^^^^^^^^^^^^^^^^^^^
+Add ``prior_kwargs`` to datamodel.py:
+
+``prior_kwargs = dict(constantSFR=False)``
+
+Add kwargs defining code block before ``add_stellar_priors()`` call in run_beast.py:
+
+.. code-block:: python
+
+  if hasattr(datamodel, 'prior_kwargs'):
+    prior_kwargs = datamodel.prior_kwargs
+  else:
+    prior_kwargs = {}
+
+Enable Exponential Av Prior
+^^^^^^^^^^^^^^^^^^^^^^^^^^^
+
+Set ``av_prior_model`` in datamodel.py:
+
+``av_prior_model = {'name': 'exponential', 'a': 2.0, 'N': 4.0}``

docs/index.rst

@@ -28,7 +28,10 @@ Installing the BEAST
 Setup the BEAST
 ===============
 
-Need documention on setting up the datamodel.py file for BEAST runs
+.. toctree::
+   :maxdepth: 1
+
+   beast_setup.rst
 
 Format of BEAST Files
 =====================