adding log Av grid spacing - explicit dust parameter grid weights now

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

@@ -24,6 +24,8 @@
 from beast.physicsmodel.grid import SpectralGrid, SEDGrid
 from beast.physicsmodel.priormodel import PriorDustModel
 
+from beast.physicsmodel.grid_weights import compute_grid_weights
+
 # from beast.external.eztables import Table
 from astropy.table import Table
 from beast.tools.helpers import generator
@@ -321,7 +323,7 @@ def make_extinguished_grid(
     # Create the sampling mesh
     # ========================
     # basically the dot product from all input 1d vectors
-    # setup interation over the full dust parameter grid
+    # setup iteration over the full dust parameter grid
 
     # setup the dust prior models
     av_prior = PriorDustModel(av_prior_model)
@@ -432,22 +434,22 @@ def make_extinguished_grid(
             #   moved here in 2023 to support distance based dust priors
             dists = g0.grid["distance"].data
             if av_prior_model["name"] == "step":
-                av_weights = av_prior(np.full((len(dists)), Av), y=dists)
+                av_prior_weights = av_prior(np.full((len(dists)), Av), y=dists)
             else:
-                av_weights = av_prior(Av)
+                av_prior_weights = av_prior(Av)
             if rv_prior_model["name"] == "step":
-                rv_weights = rv_prior(np.full((len(dists)), Rv), y=dists)
+                rv_prior_weights = rv_prior(np.full((len(dists)), Rv), y=dists)
             else:
-                rv_weights = rv_prior(Rv)
+                rv_prior_weights = rv_prior(Rv)
             if fA_prior_model["name"] == "step":
-                f_A_weights = fA_prior(np.full((len(dists)), f_A), y=dists)
+                f_A_prior_weights = fA_prior(np.full((len(dists)), f_A), y=dists)
             else:
                 if with_fA:
-                    f_A_weights = fA_prior(f_A)
+                    f_A_prior_weights = fA_prior(f_A)
                 else:
-                    f_A_weights = 1.0
+                    f_A_prior_weights = 1.0
 
-            dust_prior_weight = av_weights * rv_weights * f_A_weights
+            dust_prior_weight = av_prior_weights * rv_prior_weights * f_A_prior_weights
 
             # get new attributes if exist
             for key in list(temp_results.grid.keys()):
@@ -482,6 +484,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

@@ -14,10 +14,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,
@@ -84,7 +81,7 @@ def compute_distance_age_mass_metallicity_weights(
 
     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)
@@ -154,7 +151,7 @@ def compute_age_mass_metallicity_weights(
         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)
         age_prior = PriorAgeModel(age_prior_model)
         age_prior_weights = age_prior(uniq_ages)
 
@@ -169,7 +166,7 @@ def compute_age_mass_metallicity_weights(
             # compute the mass weights
             if len(aindxs) > 1:
                 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 = PriorMassModel(mass_prior_model)
                 mass_prior_weights = mass_prior(cur_masses)
             else:
@@ -194,7 +191,7 @@ def compute_age_mass_metallicity_weights(
     # 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 @@ def make_extinguished_sed_grid(
     """
 
     # 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])