Add log Av grid spacing option

beast/observationmodel/observations.py

@@ -8,7 +8,7 @@
 
 from beast.observationmodel.vega import Vega
 from beast.physicsmodel.priormodel import PriorAgeModel, PriorMassModel
-from beast.physicsmodel.grid_weights_stars import compute_bin_boundaries
+from beast.physicsmodel.grid_weights import compute_bin_boundaries
 
 __all__ = ["Observations", "gen_SimObs_from_sedgrid"]
 

beast/physicsmodel/creategrid.py

@@ -25,6 +25,8 @@
 from beast.physicsmodel.grid import SpectralGrid, SEDGrid
 from beast.physicsmodel.grid_and_prior_weights import compute_av_rv_fA_prior_weights
 
+from beast.physicsmodel.grid_weights import compute_grid_weights
+
 from astropy.table import Table
 from beast.tools.helpers import generator
 from beast.tools import helpers
@@ -321,7 +323,7 @@ def make_extinguished_grid(
     # Create the sampling mesh
     # ========================
     # basically the dot product from all input 1d vectors
-    # setup integration over the full dust parameter grid
+    # setup interation over the full dust parameter grid
 
     if with_fA:
 
@@ -468,6 +470,25 @@ def make_extinguished_grid(
 
         _lamb = cols.pop("lamb")
 
+        # now add the grid weights
+        av_grid_weights = compute_grid_weights(avs)
+        for cav, cav_gweight in zip(avs, av_grid_weights):
+            gvals = cols["Av"] == cav
+            cols["weight"][gvals] *= cav_gweight
+            cols["grid_weight"][gvals] *= cav_gweight
+
+        rv_grid_weights = compute_grid_weights(rvs)
+        for rav, rav_gweight in zip(rvs, rv_grid_weights):
+            gvals = cols["Rv"] == rav
+            cols["weight"][gvals] *= rav_gweight
+            cols["grid_weight"][gvals] *= rav_gweight
+
+        fA_grid_weights = compute_grid_weights(fAs)
+        for cfA, cfA_gweight in zip(fAs, fA_grid_weights):
+            gvals = cols["f_A"] == cfA
+            cols["weight"][gvals] *= cfA_gweight
+            cols["grid_weight"][gvals] *= cfA_gweight
+
         # free the memory of temp_results
         # del temp_results
         # del tempgrid

beast/physicsmodel/grid_and_prior_weights.py

@@ -15,10 +15,7 @@
 
 import numpy as np
 
-from beast.physicsmodel.grid_weights_stars import compute_distance_grid_weights
-from beast.physicsmodel.grid_weights_stars import compute_age_grid_weights
-from beast.physicsmodel.grid_weights_stars import compute_mass_grid_weights
-from beast.physicsmodel.grid_weights_stars import compute_metallicity_grid_weights
+from beast.physicsmodel.grid_weights import compute_grid_weights
 
 from beast.physicsmodel.priormodel import (
     PriorAgeModel,
@@ -87,7 +84,7 @@
 
     if n_dist > 1:
         # get the distance weights
-        dist_grid_weights = compute_distance_grid_weights(uniq_dists)
+        dist_grid_weights = compute_grid_weights(uniq_dists)
         dist_grid_weights /= np.sum(dist_grid_weights)
         dist_prior = PriorDistanceModel(distance_prior_model)
         dist_prior_weights = dist_prior(uniq_dists)
@@ -157,7 +154,7 @@
         uniq_ages = np.unique(_tgrid[zindxs]["logA"])
 
         # compute the age weights
-        age_grid_weights = compute_age_grid_weights(uniq_ages)
+        age_grid_weights = compute_grid_weights(uniq_ages, log=True)
         if isinstance(age_prior_model, dict):
             age_prior = PriorAgeModel(age_prior_model)
         else:
@@ -187,7 +184,7 @@
                 cur_masses = np.unique(_tgrid_single_age["M_ini"])
                 n_masses = len(_tgrid_single_age["M_ini"])
                 if len(cur_masses) < n_masses:
-                    umass_grid_weights = compute_mass_grid_weights(cur_masses)
+                    umass_grid_weights = compute_grid_weights(cur_masses)
                     umass_prior_weights = mass_prior(cur_masses)
                     mass_grid_weights = np.zeros(n_masses, dtype=float)
                     mass_prior_weights = np.zeros(n_masses, dtype=float)
@@ -197,7 +194,7 @@
                         mass_prior_weights[gvals] = umass_prior_weights[k]
                 else:
                     cur_masses = _tgrid_single_age["M_ini"]
-                    mass_grid_weights = compute_mass_grid_weights(cur_masses)
+                    mass_grid_weights = compute_grid_weights(cur_masses)
                     mass_prior_weights = mass_prior(cur_masses)
 
             else:
@@ -222,7 +219,7 @@
     # ensure that the metallicity prior is uniform
     if len(uniq_Zs) > 1:
         # get the metallicity weights
-        met_grid_weights = compute_metallicity_grid_weights(uniq_Zs)
+        met_grid_weights = compute_grid_weights(uniq_Zs)
         met_grid_weights /= np.sum(met_grid_weights)
         met_prior = PriorMetallicityModel(met_prior_model)
         met_prior_weights = met_prior(uniq_Zs)

beast/physicsmodel/grid_weights.py

@@ -0,0 +1,77 @@
+import numpy as np
+
+__all__ = ["compute_grid_weights", "compute_bin_boundaries"]
+
+
+def compute_grid_weights(in_x, log=False):
+    """
+    Compute the grid weights.  Needed for marginalization (aka integration).  The
+    weights are the relative widths of of each x bin.
+
+    Parameters
+    ----------
+    x : numpy array
+       centers of each bin
+
+    log : boolean
+       set if values are in log units
+
+    Returns
+    -------
+    weights : numpy array
+        weights as bin widths divided by the average width
+    """
+    # ensure x values are monotonically increasing
+    sindxs = np.argsort(in_x)
+    x = in_x[sindxs]
+
+    n_x = len(x)
+    bin_hdiffs = np.diff(x) / 2.0
+
+    # define the bin min and max boundaries
+    # handling the two edge cases
+    bin_mins = np.zeros(n_x)
+    bin_mins[1:] = x[1:] - bin_hdiffs
+    bin_mins[0] = x[0] - bin_hdiffs[0]
+
+    bin_maxs = np.zeros(n_x)
+    bin_maxs[0:-1] = x[0:-1] + bin_hdiffs
+    bin_maxs[-1] = x[-1] + bin_hdiffs[-1]
+
+    if log:
+        weights = (10**bin_maxs) - (10**bin_mins)
+    else:
+        weights = bin_maxs - bin_mins
+
+    # put the weights in the same order as in_x
+    out_weights = np.zeros(n_x)
+    out_weights[sindxs] = weights
+
+    # return normalized weights to avoid numerical issues
+    return out_weights / np.average(out_weights)
+
+
+def compute_bin_boundaries(tab):
+    """
+    Computes the boundaries of bins
+
+    The bin boundaries are defined as the midpoint between each value in tab.
+    At the two edges, 1/2 of the bin width is subtracted/added to the
+    min/max of tab.
+
+    Parameters
+    ----------
+    tab : numpy array
+       centers of each bin
+
+    Returns
+    -------
+    tab2 : numpy array
+       boundaries of the bins
+    """
+    temp = tab[1:] - np.diff(tab) / 2.0
+    tab2 = np.zeros(len(tab) + 1)
+    tab2[0] = tab[0] - np.diff(tab)[0] / 2.0
+    tab2[-1] = tab[-1] + np.diff(tab)[-1] / 2.0
+    tab2[1:-1] = temp
+    return tab2

beast/physicsmodel/model_grid.py

@@ -435,7 +435,17 @@
     """
 
     # create the dust grid arrays
-    avs = np.arange(av[0], av[1] + 0.5 * av[2], av[2])
+    if len(av) > 3:
+        # check if a log grid is requested
+        if av[3] == "log":
+            print("generating a log av grid")
+            avs = 10 ** np.arange(np.log10(av[0]), np.log10(av[1]), av[2])
+        else:
+            print("generating a linear av grid")
+            avs = np.arange(av[0], av[1] + 0.5 * av[2], av[2])
+    else:
+        print("generating a linear av grid")
+        avs = np.arange(av[0], av[1] + 0.5 * av[2], av[2])
     rvs = np.arange(rv[0], rv[1] + 0.5 * rv[2], rv[2])
     if fA is not None:
         fAs = np.arange(fA[0], fA[1] + 0.5 * fA[2], fA[2])

beast/physicsmodel/priormodel.py

@@ -3,7 +3,7 @@
 from scipy.integrate import quad
 import astropy.units as u
 
-from beast.physicsmodel.grid_weights_stars import compute_bin_boundaries
+from beast.physicsmodel.grid_weights import compute_bin_boundaries
 import beast.physicsmodel.priormodel_functions as pmfuncs