diff --git a/beast/physicsmodel/grid_and_prior_weights.py b/beast/physicsmodel/grid_and_prior_weights.py index d5e079f2e..82bc1bb69 100644 --- a/beast/physicsmodel/grid_and_prior_weights.py +++ b/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): diff --git a/beast/physicsmodel/grid_weights.py b/beast/physicsmodel/grid_weights.py index f2d11def2..5f9895547 100644 --- a/beast/physicsmodel/grid_weights.py +++ b/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 diff --git a/beast/physicsmodel/model_grid.py b/beast/physicsmodel/model_grid.py index 745c74a14..6d6e615bb 100644 --- a/beast/physicsmodel/model_grid.py +++ b/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'): diff --git a/beast/physicsmodel/prior_weights_dust.py b/beast/physicsmodel/prior_weights_dust.py index f3e24a588..0d5815bca 100644 --- a/beast/physicsmodel/prior_weights_dust.py +++ b/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**') diff --git a/beast/tools/verify_params.py b/beast/tools/verify_params.py index 35bc30d41..58079a4ef 100644 --- a/beast/tools/verify_params.py +++ b/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__': diff --git a/docs/beast_setup.rst b/docs/beast_setup.rst new file mode 100644 index 000000000..8a7858ac3 --- /dev/null +++ b/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}`` diff --git a/docs/index.rst b/docs/index.rst index 0a9d584bc..27424cea9 100644 --- a/docs/index.rst +++ b/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 =====================