Merge pull request #214 from benw1/master

AST input list generation taking source density and background density maps into account

beast/examples/phat_small/datamodel.py

@@ -47,6 +47,7 @@
 #   input data MUST be in fluxes, NOT in magnitudes 
 #   fluxes MUST be in normalized Vega units
 obs_colnames = [ f.lower() + '_rate' for f in basefilters ]
+#obs_colnames = [ f.upper() + '_RATE' for f in basefilters ]
 
 # obsfile : string 
 #   pathname of the observed catalog
@@ -69,7 +70,7 @@
 # ast_realization_per_model : integer
 # Number of Realizations of each included AST model
 # to be put into the list. (Default = 20)
-ast_realization_per_model = 20
+ast_realization_per_model = 20 
 
 
 # ast_maglimit : float (single value or array with one value per filter)
@@ -86,6 +87,21 @@
 # If True, the ast list is produced with X,Y positions.
 # If False, the ast list is produced with only magnitudes.
 ast_with_positions = True
+
+# ast_source_density_table :  (string,optional)
+# Name of source density table from tools/create_source_density_map.py
+# If supplied, the ASTs will be repeated for each source density bin in the table
+ast_source_density_table = None
+
+# ast_background_table :  (string,optional)
+# Name of background table from tools/create_background_density_map.py
+# If supplied, the ASTs will be repeated for each background level bin in the table
+ast_background_table = None
+
+#ast_N_bins : (int, optional)
+#Number of source or background bins that you want ASTs repeated over
+#ast_N_bins = 8
+
 
 # ast_pixel_distribution : float (optional)
 # (Used if ast_with_positions is True), minimum pixel separation between AST
@@ -124,8 +140,8 @@
 # Distance unit (any length or units.mag)
 distance_unit = units.mag
 
-# velocity of galaxy
-velocity =  -300 * units.km / units.s # M31 velocity from SIMBAD
+# velocity of galaxy	
+velocity =  -300 * units.km / units.s # M31 velocity from SIMBAD	
 
 ################
 

beast/examples/phat_small/run_beast.py

@@ -28,6 +28,7 @@
 import beast.observationmodel.noisemodel.generic_noisemodel as noisemodel
 from beast.observationmodel.ast.make_ast_input_list import pick_models
 from beast.observationmodel.ast.make_ast_xy_list import pick_positions
+from beast.observationmodel.ast.make_ast_xy_list import pick_positions_from_map
 from beast.fitting import fit
 from beast.fitting import trim_grid
 from beast.physicsmodel.grid import FileSEDGrid
@@ -116,8 +117,9 @@
 
     if args.ast:
         # get the modesedgrid on which to grab input AST
-        modelsedgridfile = datamodel.project + '/' + datamodel.project + \
-            '_seds.grid.hd5'
+        modelsedgridfile = './' + datamodel.project + '/' + datamodel.project + '_seds.grid.hd5'
+        modelsedgrid = FileSEDGrid(modelsedgridfile)
+        modelsedgridfile = datamodel.project + '/' + datamodel.project + '_seds.grid.hd5'
 
         N_models = datamodel.ast_models_selected_per_age
         Nfilters = datamodel.ast_bands_above_maglimit
@@ -138,19 +140,22 @@
 
             # max. mags from the gst observation cat.
             mag_cuts = min_mags + tmp_cuts
-
+        print(modelsedgrid)
         outfile = './' + datamodel.project + '/' + datamodel.project + '_inputAST.txt'
-        pick_models(modelsedgridfile, datamodel.filters, mag_cuts, Nfilter=Nfilters,
-                    N_stars=N_models, Nrealize=Nrealize, outfile=outfile)
-
+        chosen_seds = pick_models(modelsedgridfile, datamodel.filters, mag_cuts, Nfilter=Nfilters, N_stars=N_models, Nrealize=Nrealize, outfile=outfile)
         if datamodel.ast_with_positions == True:
             separation = datamodel.ast_pixel_distribution
             filename = datamodel.project + '/' + datamodel.project + '_inputAST.txt'
 
             if datamodel.ast_reference_image is not None:
                 pick_positions(obsdata, filename, separation,
-                               refimage=datamodel.ast_reference_image)
-            else:
+                               	refimage=datamodel.ast_reference_image)
+            if datamodel.ast_source_density_table is not None:
+                pick_positions_from_map(chosen_seds, datamodel.ast_source_density_table, 'sourcedens', datamodel.ast_N_bins, datamodel.ast_realization_per_model, outfile=filename, refimage=datamodel.ast_reference_image, Nrealize=1)
+
+            if datamodel.ast_background_table is not None:
+                pick_positions_from_map(chosen_seds, datamodel.ast_background_table, 'median_bg', datamodel.ast_N_bins, datamodel.ast_realization_per_model, outfile=filename, refimage=datamodel.ast_reference_image, Nrealize=1)
+            if datamodel.ast_reference_image is None and datamodel.ast_source_density_table is None and datamodel.ast_background_table is None:
                 pick_positions(obsdata, filename, separation)
 
     if args.observationmodel:

beast/observationmodel/ast/make_ast_xy_list.py

@@ -9,43 +9,46 @@
 
 from ...tools.pbar import Pbar
 
-
-def pick_positions_per_background(chosen_seds, bg_map, N_bg_bins,
+def pick_positions_from_map(chosen_seds, input_map, input_column, N_bins, Npermodel,
                                   outfile=None, refimage=None, Nrealize=1):
     """
-    Spreads a set of fake stars across regions of similar background
-    density, given a background density map file generated by 'create
-    background density map' in the tools directory.
+    Spreads a set of fake stars across regions of similar values, 
+    given a map file generated by 'create background density map' or
+    'create stellar density map' in the tools directory.
 
-    The tiles of the given background map are divided across a given
-    number of background intensity bins. Each background bin will then
-    have its own set of tiles, which constitute a region on the image.
+    The tiles of the given map are divided across a given
+    number of bins. Each bin will then have its own set of tiles, 
+    which constitute a region on the image.
 
-    Then, for each background bin, the given set of fake stars is
-    duplicated, and the stars are assigned random positions within this
-    region.
+    Then, for each bin, the given set of fake stars is duplicated, 
+    and the stars are assigned random positions within this region.
 
     This way, it can be ensured that enough ASTs are performed for each
-    regime of diffuse background emission, making it possible to have a
-    separate noise model for each of these regions.
+    regime of the map, making it possible to have a separate noise model 
+    for each of these regions.
 
     Parameters
     ----------
 
     chosen_seds: astropy Table
         Table containing fake stars to be duplicated and assigned positions
 
-    bg_map: str
-        Path to a fits file containing a background map. Each row in the
+    input_map: str
+        Path to a fits file containing a map. Each row in the
         fits table should represent a tile of the map. The table should
         have columns describing for each tile: the minimum and maximum
         RA, the minimum and maximum DEC,and a value which represents the
         background density.
 
-    N_bg_bins: int
+    input_column: str
+        Name of the column in the map file that contains of the values of 
+        interest (e.g. 'median_bg' for background or 'sourcdedens' for 
+        stellar density)
+
+    N_bins: int
         The number of bins for the range of background density values.
-        The bins will be picked on a linear grid, rangin from the
-        minimum to the maximum background value of the map. Then, each
+        The bins will be picked on a linear grid, ranging from the
+        minimum to the maximum value of the map. Then, each
         tile will be put in a bin, so that a set of tiles of the map is
         obtained for each range of background density values.
 
@@ -55,7 +58,7 @@ def pick_positions_per_background(chosen_seds, bg_map, N_bg_bins,
         instead.
 
     Nrealize: integer
-        The number of times each model shoud be repeated for each
+        The number of times each model should be repeated for each
         background regime. This is to sample the variance due to
         variations within each region, for each individual model.
 
@@ -66,30 +69,32 @@ def pick_positions_per_background(chosen_seds, bg_map, N_bg_bins,
     ascii file of this table, written to outfile
 
     """
-
+    
     # Load the background map
-    bg = Table.read(bg_map)
-    tile_bg_vals = bg['median_bg']
-    min_bg = np.amin(tile_bg_vals)
-    max_bg = np.amax(tile_bg_vals)
+    print(Npermodel,' repeats of each model in each map bin')
+    vals = Table.read(input_map)
+    tile_vals = vals[input_column]
+    min_val = np.amin(tile_vals)
+    max_val = np.amax(tile_vals)
 
-    # Create the background bins
+    # Create the bins
     # [min, ., ., ., max]
-    bg_bins = np.linspace(min_bg - 0.01 * abs(min_bg),
-                          max_bg + 0.01 * abs(max_bg), N_bg_bins + 1)
-
+    bins = np.linspace(min_val - 0.01 * abs(min_val),
+                          max_val + 0.01 * abs(max_val), N_bins + 1)
+    
     # Find which bin each tile belongs to
     # e.g. one of these numbers: 0 [1, 2, 3, 4, 5] 6
     # We have purposely chosen our bin boundaries so that no points fall
     # outside of the [1,5] range
-    bgbin_foreach_tile = np.digitize(tile_bg_vals, bg_bins)
+    bin_foreach_tile = np.digitize(tile_vals, bins)
     # Invert this (the [0] is to dereference the tuple (i,) returned by
     # nonzero)
-    tiles_foreach_bgbin = [np.nonzero(bgbin_foreach_tile == b + 1)[0]
-                           for b in range(N_bg_bins)]
+    tiles_foreach_bin = [np.nonzero(bin_foreach_tile == b + 1)[0]
+                           for b in range(N_bins)]
 
     # Remove empty bins
-    tile_sets = [tile_set for tile_set in tiles_foreach_bgbin if len(tile_set)]
+    tile_sets = [tile_set for tile_set in tiles_foreach_bin if len(tile_set)]
+    print(len(tile_sets), ' non-empty map bins found between ', min_val, 'and', max_val)
 
     # Repeat the seds Nrealize times (sample each on at Nrealize
     # different positions, in each region)
@@ -104,10 +109,10 @@ def pick_positions_per_background(chosen_seds, bg_map, N_bg_bins,
     ys = np.zeros(len(out_table))
     bin_indices = np.zeros(len(out_table))
 
-    tile_ra_min = bg['min_ra']
-    tile_dec_min = bg['min_dec']
-    tile_ra_delta = bg['max_ra'] - tile_ra_min
-    tile_dec_delta = bg['max_dec'] - tile_dec_min
+    tile_ra_min = vals['min_ra']
+    tile_dec_min = vals['min_dec']
+    tile_ra_delta = vals['max_ra'] - tile_ra_min
+    tile_dec_delta = vals['max_dec'] - tile_dec_min
 
     if refimage is None:
         wcs = None
@@ -116,7 +121,7 @@ def pick_positions_per_background(chosen_seds, bg_map, N_bg_bins,
         wcs = WCS(imagehdu.header)
 
     pbar = Pbar(len(tile_sets),
-                desc='{} models per background bin'.format(Nseds_per_region))
+                desc='{} models per map bin'.format(Nseds_per_region/Npermodel))
     for bin_index, tile_set in pbar.iterover(enumerate(tile_sets)):
         start = bin_index * Nseds_per_region
         stop = start + Nseds_per_region
@@ -168,6 +173,7 @@ def pick_positions_per_background(chosen_seds, bg_map, N_bg_bins,
     return out_table
 
 
+
 def pick_positions(catalog, filename, separation, refimage=None):
     """
     Assigns positions to fake star list generated by pick_models

beast/tools/create_source_density_map.py

@@ -85,7 +85,7 @@ def make_source_dens_map(catfile,
     ra_delt = dec_delt
     n_x = np.fix(np.round(math.cos(0.5*(max_dec+min_dec)*math.pi/180.)
                           *(max_ra-min_ra)/ra_delt))
-    ra_delt *= -1.
+    #ra_delt *= -1. #Not sure why the ra delta would want to be negative...
 
     n_x = int(np.max([n_x,1]))
     n_y = int(np.max([n_y,1]))
@@ -185,6 +185,17 @@ def make_source_dens_map(catfile,
     cat[nonzero_indxs].write(catfile.replace('.fits',
                                              '_with_sourceden.fits'),
                              overwrite=True)
+
+    bin_details = Table(names=['i_ra', 'i_dec', 'sourcedens', 'min_ra', 'max_ra', 'min_dec', 'max_dec'])
+
+    for i in range(n_x):
+
+        for j in range(n_y):
+
+            bin_details.add_row([i, j, npts_map[i, j], ra_limits[i], ra_limits[i + 1], dec_limits[j], dec_limits[j + 1]])
+
+    bin_details.write(catfile.replace('.fits', '_source_den_bins_table.fits'))
+
 
 if __name__ == '__main__':
 

docs/generating_asts.rst

@@ -8,8 +8,10 @@ Once the input lists have been run through the user's photometry program and eac
 
 Generating BEAST-friendly lists of artificial star tests
 
-1) Run "run_small.py -p" to produce the physics model grid file "project_name_seds.grid.hd5".
-2) Run "run_small.py -a"   This will use the datamodel to find everything it needs to make ASTs (filters, limits, SED grid, etc.).  It will produce a list of fake stars in all bands using the datamodel photometry catalog to trim the inputs at the proper magnitudes.  Currently, this script generates fake stars uniformly sampling log(age) space and randomly drawing from the metallicities in the model grid.
+1) Run "run_beast.py -p" to produce the physics model grid file "project_name_seds.grid.hd5".
+2) If you wish to repeat ASTs over multiple bins of stellar density, run "tools/create_source_density_map.py" on your observed star catalog.
+   If you wish to repeat ASTs over multiple bins of background brightness, run "tools/create_background_density_map.py" on your observed star catalog and reference image.
+3) Run "run_beast.py -a"   This will use the datamodel to find everything it needs to make ASTs (filters, limits, SED grid, etc.).  It will produce a list of fake stars in all bands using the datamodel photometry catalog to trim the inputs at the proper magnitudes.  Currently, this script generates fake stars uniformly sampling log(age) space and randomly drawing from the metallicities in the model grid.
 
 Functions
 =========
@@ -48,6 +50,14 @@ ast_with_positions :  (bool,optional)
 If True, the ast list is produced with X,Y positions.
 If False, the ast list is produced with only magnitudes.
 
+ast_multiple_source_density : string (optional, if defined the file named 
+will be used to determine spatial distribution of ASTs, and all ASTs will be repeated 
+within each source density region.
+
+ast_multiple_background brightness : string (optional, if defined the file named 
+will be used to determine spatial distribution of ASTs, and all ASTs will be repeated 
+within each background brightness region.
+
 ast_pixel_distribution : float (optional)
 (Used if ast_with_positions is True), minimum pixel separation between AST
 position and catalog star used to determine the AST spatial distribution.