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

@@ -23,7 +23,9 @@
 
 from beast.physicsmodel.stars import stellib
 from beast.physicsmodel.grid import SpectralGrid, SEDGrid
-from beast.physicsmodel.grid_and_prior_weights import compute_av_rv_fA_prior_weights
+from beast.physicsmodel.priormodel import PriorDustModel
+
+from beast.physicsmodel.grid_weights import compute_grid_weights
 
 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 integration over the full dust parameter grid
+    # setup interation over the full dust parameter grid
 
     if with_fA:
 
@@ -425,15 +427,25 @@ def make_extinguished_grid(
                 cols["Rv"][N0 * count : N0 * (count + 1)] = Rv
 
             # compute the dust weights
-            dust_prior_weight = compute_av_rv_fA_prior_weights(
-                Av,
-                Rv,
-                f_A,
-                g0.grid["distance"].data,
-                av_prior_model=av_prior_model,
-                rv_prior_model=rv_prior_model,
-                fA_prior_model=fA_prior_model,
-            )
+            #   moved here in 2023 to support distance based dust priors
+            dists = g0.grid["distance"].data
+            if av_prior_model["name"] == "step":
+                av_prior_weights = av_prior(np.full((len(dists)), Av), y=dists)
+            else:
+                av_prior_weights = av_prior(Av)
+            if rv_prior_model["name"] == "step":
+                rv_prior_weights = rv_prior(np.full((len(dists)), Rv), y=dists)
+            else:
+                rv_prior_weights = rv_prior(Rv)
+            if fA_prior_model["name"] == "step":
+                f_A_prior_weights = fA_prior(np.full((len(dists)), f_A), y=dists)
+            else:
+                if with_fA:
+                    f_A_prior_weights = fA_prior(f_A)
+                else:
+                    f_A_prior_weights = 1.0
+
+            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()):
@@ -468,6 +480,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 @@ 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)
@@ -157,7 +154,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)
         if isinstance(age_prior_model, dict):
             age_prior = PriorAgeModel(age_prior_model)
         else:
@@ -187,7 +184,7 @@ def compute_age_mass_metallicity_weights(
                 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 @@ def compute_age_mass_metallicity_weights(
                         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 @@ 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