From 0f8131d1333c4b1ead35677ca36a26a958b02b9d Mon Sep 17 00:00:00 2001
From: Benjamin Williams <ben@D-10-19-217-41.dhcp4.washington.edu>
Date: Thu, 3 May 2018 10:41:20 -0700
Subject: [PATCH 01/10] added fits table output to creas_source_density_map

---
 beast/tools/create_source_density_map.py | 11 +++++++++++
 1 file changed, 11 insertions(+)

diff --git a/beast/tools/create_source_density_map.py b/beast/tools/create_source_density_map.py
index 500871e1e..194d1ee23 100755
--- a/beast/tools/create_source_density_map.py
+++ b/beast/tools/create_source_density_map.py
@@ -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__':
 

From 806636ce0474da165e81ddefc2cb7f41011efdba Mon Sep 17 00:00:00 2001
From: Benjamin Williams <ben@D-10-19-217-41.dhcp4.washington.edu>
Date: Thu, 3 May 2018 11:38:44 -0700
Subject: [PATCH 02/10] added options to AST generation for background and
 source density to datamodel

---
 beast/examples/phat_small/datamodel.py | 10 ++++++++++
 1 file changed, 10 insertions(+)

diff --git a/beast/examples/phat_small/datamodel.py b/beast/examples/phat_small/datamodel.py
index 138db3bc2..9035fde94 100755
--- a/beast/examples/phat_small/datamodel.py
+++ b/beast/examples/phat_small/datamodel.py
@@ -86,6 +86,16 @@
 # 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_with_source_density :  (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_with_source_density = None
+
+# ast_with_background :  (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_with_background = None
                          
 # ast_pixel_distribution : float (optional)
 # (Used if ast_with_positions is True), minimum pixel separation between AST

From ba2b434e73df2951aa2a3550b3c42a314e1d9562 Mon Sep 17 00:00:00 2001
From: Benjamin Williams <ben@D-10-19-217-41.dhcp4.washington.edu>
Date: Thu, 3 May 2018 13:29:34 -0700
Subject: [PATCH 03/10] adding repeating ASTs over different background or
 density bins

---
 beast/examples/phat_small/datamodel.py        |   5 +
 beast/examples/phat_small/run_beast.py        |   7 +
 .../observationmodel/ast/make_ast_xy_list.py  | 139 ++++++++++++++++++
 docs/generating_asts.rst                      |  14 +-
 4 files changed, 163 insertions(+), 2 deletions(-)

diff --git a/beast/examples/phat_small/datamodel.py b/beast/examples/phat_small/datamodel.py
index 9035fde94..5fa24b706 100755
--- a/beast/examples/phat_small/datamodel.py
+++ b/beast/examples/phat_small/datamodel.py
@@ -96,6 +96,11 @@
 # 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_with_background = None
+
+#ast_N_dens_bins : (int, optional)
+#Number of source or background bins that you want ASTs repeated over
+ast_N_dens_bins = 10
+
                          
 # ast_pixel_distribution : float (optional)
 # (Used if ast_with_positions is True), minimum pixel separation between AST
diff --git a/beast/examples/phat_small/run_beast.py b/beast/examples/phat_small/run_beast.py
index ce6614eb4..b235a2704 100755
--- a/beast/examples/phat_small/run_beast.py
+++ b/beast/examples/phat_small/run_beast.py
@@ -145,6 +145,13 @@
             if datamodel.ast_reference_image is not None:
                 pick_positions(obsdata, filename, separation,
                                refimage=datamodel.ast_reference_image)
+            if datamodel.ast_source_density_table is not None:
+		pick_positions_per_density(datamodel.chosen_seds, datamodel.ast_source_density_table, datamodel.ast_N_dens_bins,
+                                  outfile=filename, refimage=datamodel.ast_reference_image, Nrealize=1):
+
+            if datamodel.ast_background_table is not None:
+		pick_positions_per_density(datamodel.chosen_seds, datamodel.ast_background_table, datamodel.ast_N_dens_bins,
+                                  outfile=filename, refimage=datamodel.ast_reference_image, Nrealize=1):
             else:
                 pick_positions(obsdata, filename, separation)
 
diff --git a/beast/observationmodel/ast/make_ast_xy_list.py b/beast/observationmodel/ast/make_ast_xy_list.py
index ae607e119..90475ac91 100644
--- a/beast/observationmodel/ast/make_ast_xy_list.py
+++ b/beast/observationmodel/ast/make_ast_xy_list.py
@@ -241,3 +241,142 @@ def pick_positions(catalog, filename, separation, refimage=None):
     
     ascii.write(astmags,filename,overwrite=True)
     
+def pick_positions_per_density(chosen_seds, dens_map, N_dens_bins,
+                                  outfile=None, refimage=None, Nrealize=1):
+    """
+    Spreads a set of fake stars across regions of similar source
+    density, given a source density map file generated by 'create
+    source density map' in the tools directory.
+
+    The tiles of the given density map are divided across a given
+    number of source density bins. Each density bin will then
+    have its own set of tiles, which constitute a region on the image.
+
+    Then, for each density 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 stellar crowding, 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
+
+    dens_map: str
+        Path to a fits file containing a density 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
+        source density.
+
+    N_dens_bins: int
+        The number of bins for the range of source density values.
+        The bins will be picked on a linear grid, rangin from the
+        minimum to the maximum density 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.
+
+    refimage: str
+        Path to fits image that is used for the positions. If none is
+        given, the ra and dec will be put in the x and y output columns
+        instead.
+
+    Nrealize: integer
+        The number of times each model shoud be repeated for each
+        background regime. This is to sample the variance due to
+        variations within each region, for each individual model.
+
+    Returns
+    -------
+    astropy Table: List of fake stars, with magnitudes and positions
+    - optionally -
+    ascii file of this table, written to outfile
+
+    """
+
+    # Load the density map
+    dens = Table.read(dens_map)
+    tile_dens_vals = dens['median_dens']
+    min_dens = np.amin(tile_dens_vals)
+    max_dens = np.amax(tile_dens_vals)
+
+    # Create the background bins
+    # [min, ., ., ., max]
+    dens_bins = np.linspace(min_dens - 0.01 * abs(min_dens),
+                          max_dens + 0.01 * abs(max_dens), N_dens_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
+    densbin_foreach_tile = np.digitize(tile_dens_vals, dens_bins)
+    # Invert this (the [0] is to dereference the tuple (i,) returned by
+    # nonzero)
+    tiles_foreach_densbin = [np.nonzero(densbin_foreach_tile == b + 1)[0]
+                           for b in range(N_dens_bins)]
+
+    # Remove empty bins
+    tile_sets = [tile_set for tile_set in tiles_foreach_densbin if len(tile_set)]
+
+    # Repeat the seds Nrealize times (sample each on at Nrealize
+    # different positions, in each region)
+    repeated_seds = np.repeat(chosen_seds, Nrealize)
+    Nseds_per_region = len(repeated_seds)
+    # For each set of tiles, repeat the seds and spread them evenly over
+    # the tiles
+    repeated_seds = np.repeat(repeated_seds, len(tile_sets))
+
+    out_table = Table(repeated_seds, names=chosen_seds.colnames)
+    ras = np.zeros(len(out_table))
+    decs = np.zeros(len(out_table))
+    bin_indices = np.zeros(len(out_table))
+
+    tile_ra_min = dens['min_ra']
+    tile_dec_min = dens['min_dec']
+    tile_ra_delta = dens['max_ra'] - tile_ra_min
+    tile_dec_delta = dens['max_dec'] - tile_dec_min
+
+    pbar = Pbar(len(tile_sets),
+                desc='{} models per background bin'.format(Nseds_per_region))
+    for bin_index, tile_set in pbar.iterover(enumerate(tile_sets)):
+        start = bin_index * Nseds_per_region
+        stop = start + Nseds_per_region
+        bin_indices[start:stop] = bin_index
+        for i in range(Nseds_per_region):
+            j = bin_index * Nseds_per_region + i
+            # Pick a random tile
+            tile = np.random.choice(tile_set)
+            # Within this tile, pick a random ra and dec
+            ras[j] = tile_ra_min[tile] + \
+                np.random.random_sample() * tile_ra_delta[tile]
+            decs[j] = tile_dec_min[tile] + \
+                np.random.random_sample() * tile_dec_delta[tile]
+
+    # I'm just mimicking the format that is produced by the examples
+    cs = []
+    cs.append(Column(np.zeros(len(out_table)), name='zeros'))
+    cs.append(Column(np.ones(len(out_table)), name='ones'))
+
+    if refimage is None:
+        cs.append(Column(ras, name='RA'))
+        cs.append(Column(decs, name='DEC'))
+    else:
+        imagehdu = fits.open(refimage)[1]
+        wcs = WCS(imagehdu.header)
+        xs, ys = wcs.all_world2pix(ras, decs, 0)
+        cs.append(Column(xs, name='X'))
+        cs.append(Column(ys, name='Y'))
+
+    for i, c in enumerate(cs):
+        out_table.add_column(c, index=i)  # insert these columns from the left
+
+    # Write out the table in ascii
+    if outfile:
+        formats = {k: '%.5f' for k in out_table.colnames}
+        ascii.write(out_table, outfile, overwrite=True, formats=formats)
+
+    return out_table
diff --git a/docs/generating_asts.rst b/docs/generating_asts.rst
index 31310a2d0..80c520a42 100644
--- a/docs/generating_asts.rst
+++ b/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.

From 95fac043aa3322a3e5325c51c512d3e1d1aa1ac2 Mon Sep 17 00:00:00 2001
From: Benjamin Williams <ben@D-10-19-48-105.dhcp4.washington.edu>
Date: Mon, 14 May 2018 14:24:13 -0700
Subject: [PATCH 04/10] working on ASTs in run_beast

---
 beast/examples/phat_small/datamodel.py | 18 +++++++++++-------
 beast/examples/phat_small/run_beast.py | 12 +++++++-----
 2 files changed, 18 insertions(+), 12 deletions(-)

diff --git a/beast/examples/phat_small/datamodel.py b/beast/examples/phat_small/datamodel.py
index 5fa24b706..a559d8d06 100755
--- a/beast/examples/phat_small/datamodel.py
+++ b/beast/examples/phat_small/datamodel.py
@@ -46,11 +46,13 @@
 #   need to match column names in the observed catalog,
 #   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.lower() + '_rate' for f in basefilters ]
+obs_colnames = [ f.upper() + '_RATE' for f in basefilters ]
 
 # obsfile : string 
 #   pathname of the observed catalog
-obsfile = 'data/b15_4band_det_27_A.fits'
+#obsfile = 'data/b15_4band_det_27_A.fits'
+obsfile = 'data/14675_LMC-5665ne-12232.st.fits'
 
 #------------------------------------------------------
 # Artificial Star Test Input File Generation Parameters
@@ -87,15 +89,16 @@
 # If False, the ast list is produced with only magnitudes.
 ast_with_positions = True
 
-# ast_with_source_density :  (string,optional)
+# 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_with_source_density = None
+ast_source_density_table = None
 
-# ast_with_background :  (string,optional)
+# 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_with_background = None
+#ast_background_table = None
+ast_background_table = "data/14675_LMC-5665ne-12232_F475W_drz_background_map.fits"
 
 #ast_N_dens_bins : (int, optional)
 #Number of source or background bins that you want ASTs repeated over
@@ -111,7 +114,8 @@
 # is True and no X and Y information  is present in the photometry catalog)	
 # Name of the reference image used by DOLPHOT when running the measured 
 # photometry.	            
-ast_reference_image = None
+#ast_reference_image = None
+ast_reference_image = "data/14675_LMC-5665ne-12232_F475W_drz.fits"
 
 
 #-------------------------------------------
diff --git a/beast/examples/phat_small/run_beast.py b/beast/examples/phat_small/run_beast.py
index b235a2704..cdfad1238 100755
--- a/beast/examples/phat_small/run_beast.py
+++ b/beast/examples/phat_small/run_beast.py
@@ -27,6 +27,8 @@
 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_per_density
+from beast.observationmodel.ast.make_ast_xy_list import pick_positions_per_background
 from beast.fitting import fit
 from beast.fitting import trim_grid
 from beast.physicsmodel.grid import FileSEDGrid
@@ -144,14 +146,14 @@
 
             if datamodel.ast_reference_image is not None:
                 pick_positions(obsdata, filename, separation,
-                               refimage=datamodel.ast_reference_image)
+                               	refimage=datamodel.ast_reference_image)
             if datamodel.ast_source_density_table is not None:
-		pick_positions_per_density(datamodel.chosen_seds, datamodel.ast_source_density_table, datamodel.ast_N_dens_bins,
-                                  outfile=filename, refimage=datamodel.ast_reference_image, Nrealize=1):
+                pick_positions_per_density(filename, datamodel.ast_source_density_table, datamodel.ast_N_dens_bins,
+                               	outfile=filename, refimage=datamodel.ast_reference_image, Nrealize=1)
 
             if datamodel.ast_background_table is not None:
-		pick_positions_per_density(datamodel.chosen_seds, datamodel.ast_background_table, datamodel.ast_N_dens_bins,
-                                  outfile=filename, refimage=datamodel.ast_reference_image, Nrealize=1):
+                pick_positions_per_background(filename, datamodel.ast_background_table, datamodel.ast_N_dens_bins,
+                               	outfile=filename, refimage=datamodel.ast_reference_image, Nrealize=1)
             else:
                 pick_positions(obsdata, filename, separation)
 

From cfbbf50cb6dbd65c50dd556c02a806a2630c4297 Mon Sep 17 00:00:00 2001
From: Benjamin Williams <ben@D-10-19-241-235.dhcp4.washington.edu>
Date: Fri, 18 May 2018 11:32:03 -0700
Subject: [PATCH 05/10] working on ASTs in run_beast

---
 beast/examples/phat_small/run_beast.py | 17 +++++++----------
 1 file changed, 7 insertions(+), 10 deletions(-)

diff --git a/beast/examples/phat_small/run_beast.py b/beast/examples/phat_small/run_beast.py
index 099f9352e..21792e447 100755
--- a/beast/examples/phat_small/run_beast.py
+++ b/beast/examples/phat_small/run_beast.py
@@ -111,10 +111,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)
-
+        print(modelsedgridfile)
         N_models = datamodel.ast_models_selected_per_age
         Nfilters = datamodel.ast_bands_above_maglimit
         Nrealize = datamodel.ast_realization_per_model
@@ -134,11 +133,9 @@
 
             # max. mags from the gst observation cat.
             mag_cuts = min_mags + tmp_cuts
-
+        print(modelsedgrid)
         outfile = './' + datamodel.project + '/' + datamodel.project + '_inputAST.txt'
-        pick_models(modelsedgrid, 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'
@@ -147,13 +144,13 @@
                 pick_positions(obsdata, filename, separation,
                                	refimage=datamodel.ast_reference_image)
             if datamodel.ast_source_density_table is not None:
-                pick_positions_per_density(filename, datamodel.ast_source_density_table, datamodel.ast_N_dens_bins,
+                pick_positions_per_density(chosen_seds, datamodel.ast_source_density_table, datamodel.ast_N_dens_bins,
                                	outfile=filename, refimage=datamodel.ast_reference_image, Nrealize=1)
 
             if datamodel.ast_background_table is not None:
-                pick_positions_per_background(filename, datamodel.ast_background_table, datamodel.ast_N_dens_bins,
+                pick_positions_per_background(chosen_seds, datamodel.ast_background_table, datamodel.ast_N_dens_bins,
                                	outfile=filename, refimage=datamodel.ast_reference_image, Nrealize=1)
-            else:
+            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:

From f12df7b5a44f906cec78fd420be28b8a674749c0 Mon Sep 17 00:00:00 2001
From: Benjamin Williams <ben@MacBook-Pro-6.domain>
Date: Thu, 24 May 2018 21:10:16 -0700
Subject: [PATCH 06/10] source density bins now have max_ra > min_ra and ast
 lists can now be made based on source density

---
 .../observationmodel/ast/make_ast_xy_list.py  | 58 ++++++++++++-------
 beast/tools/create_source_density_map.py      |  2 +-
 2 files changed, 39 insertions(+), 21 deletions(-)

diff --git a/beast/observationmodel/ast/make_ast_xy_list.py b/beast/observationmodel/ast/make_ast_xy_list.py
index 4d042e58e..50ea7e117 100644
--- a/beast/observationmodel/ast/make_ast_xy_list.py
+++ b/beast/observationmodel/ast/make_ast_xy_list.py
@@ -321,11 +321,10 @@ def pick_positions_per_density(chosen_seds, dens_map, N_dens_bins,
     min_dens = np.amin(tile_dens_vals)
     max_dens = np.amax(tile_dens_vals)
 
-    # Create the background bins
+    # Create the sourcedens bins
     # [min, ., ., ., max]
     dens_bins = np.linspace(min_dens - 0.01 * abs(min_dens),
                           max_dens + 0.01 * abs(max_dens), N_dens_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
@@ -348,8 +347,8 @@ def pick_positions_per_density(chosen_seds, dens_map, N_dens_bins,
     repeated_seds = np.repeat(repeated_seds, len(tile_sets))
 
     out_table = Table(repeated_seds, names=chosen_seds.colnames)
-    ras = np.zeros(len(out_table))
-    decs = np.zeros(len(out_table))
+    xs = np.zeros(len(out_table))
+    ys = np.zeros(len(out_table))
     bin_indices = np.zeros(len(out_table))
 
     tile_ra_min = dens['min_ra']
@@ -357,6 +356,14 @@ def pick_positions_per_density(chosen_seds, dens_map, N_dens_bins,
     tile_ra_delta = dens['max_ra'] - tile_ra_min
     tile_dec_delta = dens['max_dec'] - tile_dec_min
 
+    if refimage is None:
+        wcs = None
+    else:
+        imagehdu = fits.open(refimage)[1]
+        wcs = WCS(imagehdu.header)
+
+
+
     pbar = Pbar(len(tile_sets),
                 desc='{} models per background bin'.format(Nseds_per_region))
     for bin_index, tile_set in pbar.iterover(enumerate(tile_sets)):
@@ -364,27 +371,37 @@ def pick_positions_per_density(chosen_seds, dens_map, N_dens_bins,
         stop = start + Nseds_per_region
         bin_indices[start:stop] = bin_index
         for i in range(Nseds_per_region):
+            x = -1
+            y = -1
+            # Convert each ra,dec to x,y. If there are negative values, try again
+            while x < 0 or y < 0:
+                # Pick a random tile
+                tile = np.random.choice(tile_set)
+                # Within this tile, pick a random ra and dec
+                ra = tile_ra_min[tile] + \
+                     np.random.random_sample() * tile_ra_delta[tile]
+                dec = tile_dec_min[tile] + \
+                      np.random.random_sample() * tile_dec_delta[tile]
+
+                if wcs is None:
+                    x, y = ra, dec
+                    break
+                else:
+                    [x], [y] = wcs.all_world2pix(np.array([ra]), np.array([dec]), 0)
             j = bin_index * Nseds_per_region + i
-            # Pick a random tile
-            tile = np.random.choice(tile_set)
-            # Within this tile, pick a random ra and dec
-            ras[j] = tile_ra_min[tile] + \
-                np.random.random_sample() * tile_ra_delta[tile]
-            decs[j] = tile_dec_min[tile] + \
-                np.random.random_sample() * tile_dec_delta[tile]
+            xs[j] = x
+            ys[j] = y
+
 
     # I'm just mimicking the format that is produced by the examples
     cs = []
-    cs.append(Column(np.zeros(len(out_table)), name='zeros'))
-    cs.append(Column(np.ones(len(out_table)), name='ones'))
+    cs.append(Column(np.zeros(len(out_table), dtype=int), name='zeros'))
+    cs.append(Column(np.ones(len(out_table), dtype=int), name='ones'))
 
-    if refimage is None:
-        cs.append(Column(ras, name='RA'))
-        cs.append(Column(decs, name='DEC'))
+    if wcs is None:
+        cs.append(Column(xs, name='RA'))
+        cs.append(Column(ys, name='DEC'))
     else:
-        imagehdu = fits.open(refimage)[1]
-        wcs = WCS(imagehdu.header)
-        xs, ys = wcs.all_world2pix(ras, decs, 0)
         cs.append(Column(xs, name='X'))
         cs.append(Column(ys, name='Y'))
 
@@ -393,7 +410,8 @@ def pick_positions_per_density(chosen_seds, dens_map, N_dens_bins,
 
     # Write out the table in ascii
     if outfile:
-        formats = {k: '%.5f' for k in out_table.colnames}
+        formats = {k: '%.5f' for k in out_table.colnames[2:]}
         ascii.write(out_table, outfile, overwrite=True, formats=formats)
 
     return out_table
+
diff --git a/beast/tools/create_source_density_map.py b/beast/tools/create_source_density_map.py
index 4e980de24..e6e3159f3 100755
--- a/beast/tools/create_source_density_map.py
+++ b/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]))

From 0aa7e497d1175f982c6b382b9ab572b2dccdb8dc Mon Sep 17 00:00:00 2001
From: Benjamin Williams <ben@D-10-19-28-116.dhcp4.washington.edu>
Date: Fri, 25 May 2018 10:18:27 -0700
Subject: [PATCH 07/10] putting datamodel.py into a correct state for the
 standard example

---
 beast/examples/phat_small/datamodel.py | 11 +++--------
 1 file changed, 3 insertions(+), 8 deletions(-)

diff --git a/beast/examples/phat_small/datamodel.py b/beast/examples/phat_small/datamodel.py
index 0ba024cb6..ca3ad863f 100755
--- a/beast/examples/phat_small/datamodel.py
+++ b/beast/examples/phat_small/datamodel.py
@@ -51,8 +51,7 @@
 
 # obsfile : string 
 #   pathname of the observed catalog
-#obsfile = 'data/b15_4band_det_27_A.fits'
-obsfile = 'data/14675_LMC-5665ne-12232.st.fits'
+obsfile = 'data/b15_4band_det_27_A.fits'
 
 #------------------------------------------------------
 # Artificial Star Test Input File Generation Parameters
@@ -92,14 +91,12 @@
 # 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_source_density_table = "data/14675_LMC-5665ne-12232.st_source_den_bins_table.fits"
+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_background_table = "data/14675_LMC-5665ne-12232_F475W_drz_background_map.fits"
 
 #ast_N_dens_bins : (int, optional)
 #Number of source or background bins that you want ASTs repeated over
@@ -115,9 +112,7 @@
 # is True and no X and Y information  is present in the photometry catalog)	
 # Name of the reference image used by DOLPHOT when running the measured 
 # photometry.	            
-#ast_reference_image = None
-ast_reference_image = "data/14675_LMC-5665ne-12232_F475W_drz.fits"
-
+ast_reference_image = None
 
 #-------------------------------------------
 #Noise Model Artificial Star Test Parameters

From 965c8407aee5b8ae3ad6c67354bf193cfa0c8f26 Mon Sep 17 00:00:00 2001
From: Benjamin Williams <ben@D-10-18-246-137.dhcp4.washington.edu>
Date: Fri, 1 Jun 2018 12:57:22 -0700
Subject: [PATCH 08/10] merging functions

---
 beast/examples/phat_small/run_beast.py        |   4 +-
 .../observationmodel/ast/make_ast_xy_list.py  | 164 +++++++++++++++++-
 2 files changed, 165 insertions(+), 3 deletions(-)

diff --git a/beast/examples/phat_small/run_beast.py b/beast/examples/phat_small/run_beast.py
index 70349e37b..9eeb74a03 100755
--- a/beast/examples/phat_small/run_beast.py
+++ b/beast/examples/phat_small/run_beast.py
@@ -145,11 +145,11 @@
                 pick_positions(obsdata, filename, separation,
                                	refimage=datamodel.ast_reference_image)
             if datamodel.ast_source_density_table is not None:
-                pick_positions_per_density(chosen_seds, datamodel.ast_source_density_table, datamodel.ast_N_dens_bins,
+                pick_positions_from_map(chosen_seds, datamodel.ast_source_density_table, 'sourcedens', datamodel.ast_N_dens_bins,
                                	outfile=filename, refimage=datamodel.ast_reference_image, Nrealize=1)
 
             if datamodel.ast_background_table is not None:
-                pick_positions_per_background(chosen_seds, datamodel.ast_background_table, datamodel.ast_N_dens_bins,
+                pick_positions_from_map(chosen_seds, datamodel.ast_background_table, 'median_bg', datamodel.ast_N_dens_bins,
                                	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)
diff --git a/beast/observationmodel/ast/make_ast_xy_list.py b/beast/observationmodel/ast/make_ast_xy_list.py
index 50ea7e117..846ae665f 100644
--- a/beast/observationmodel/ast/make_ast_xy_list.py
+++ b/beast/observationmodel/ast/make_ast_xy_list.py
@@ -9,6 +9,168 @@
 
 from ...tools.pbar import Pbar
 
+def pick_positions_from_map(chosen_seds, input_map, input_column, N_bins,
+                                  outfile=None, refimage=None, Nrealize=1):
+    """
+    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 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 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 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
+
+    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.
+
+    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, 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.
+
+    refimage: str
+        Path to fits image that is used for the positions. If none is
+        given, the ra and dec will be put in the x and y output columns
+        instead.
+
+    Nrealize: integer
+        The number of times each model shoud be repeated for each
+        background regime. This is to sample the variance due to
+        variations within each region, for each individual model.
+
+    Returns
+    -------
+    astropy Table: List of fake stars, with magnitudes and positions
+    - optionally -
+    ascii file of this table, written to outfile
+
+    """
+    
+    # Load the background map
+    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
+    # [min, ., ., ., max]
+    bins = np.linspace(min_val - 0.01 * abs(min_al),
+                          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
+    bin_foreach_tile = np.digitize(tile_vals, bins)
+    # Invert this (the [0] is to dereference the tuple (i,) returned by
+    # nonzero)
+    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_bin if len(tile_set)]
+
+    # Repeat the seds Nrealize times (sample each on at Nrealize
+    # different positions, in each region)
+    repeated_seds = np.repeat(chosen_seds, Nrealize)
+    Nseds_per_region = len(repeated_seds)
+    # For each set of tiles, repeat the seds and spread them evenly over
+    # the tiles
+    repeated_seds = np.repeat(repeated_seds, len(tile_sets))
+
+    out_table = Table(repeated_seds, names=chosen_seds.colnames)
+    xs = np.zeros(len(out_table))
+    ys = np.zeros(len(out_table))
+    bin_indices = np.zeros(len(out_table))
+
+    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
+    else:
+        imagehdu = fits.open(refimage)[1]
+        wcs = WCS(imagehdu.header)
+
+    pbar = Pbar(len(tile_sets),
+                desc='{} models per background bin'.format(Nseds_per_region))
+    for bin_index, tile_set in pbar.iterover(enumerate(tile_sets)):
+        start = bin_index * Nseds_per_region
+        stop = start + Nseds_per_region
+        bin_indices[start:stop] = bin_index
+        for i in range(Nseds_per_region):
+            x = -1
+            y = -1
+            # Convert each ra,dec to x,y. If there are negative values, try again
+            while x < 0 or y < 0:
+                # Pick a random tile
+                tile = np.random.choice(tile_set)
+                # Within this tile, pick a random ra and dec
+                ra = tile_ra_min[tile] + \
+                     np.random.random_sample() * tile_ra_delta[tile]
+                dec = tile_dec_min[tile] + \
+                      np.random.random_sample() * tile_dec_delta[tile]
+
+                if wcs is None:
+                    x, y = ra, dec
+                    break
+                else:
+                    [x], [y] = wcs.all_world2pix(np.array([ra]), np.array([dec]), 0)
+
+            j = bin_index * Nseds_per_region + i
+            xs[j] = x
+            ys[j] = y
+
+
+    # I'm just mimicking the format that is produced by the examples
+    cs = []
+    cs.append(Column(np.zeros(len(out_table), dtype=int), name='zeros'))
+    cs.append(Column(np.ones(len(out_table), dtype=int), name='ones'))
+
+    if wcs is None:
+        cs.append(Column(xs, name='RA'))
+        cs.append(Column(ys, name='DEC'))
+    else:
+        cs.append(Column(xs, name='X'))
+        cs.append(Column(ys, name='Y'))
+
+    for i, c in enumerate(cs):
+        out_table.add_column(c, index=i)  # insert these columns from the left
+
+    # Write out the table in ascii
+    if outfile:
+        formats = {k: '%.5f' for k in out_table.colnames[2:]}
+        ascii.write(out_table, outfile, overwrite=True, formats=formats)
+
+    return out_table
+
+
 
 def pick_positions_per_background(chosen_seds, bg_map, N_bg_bins,
                                   outfile=None, refimage=None, Nrealize=1):
@@ -292,7 +454,7 @@ def pick_positions_per_density(chosen_seds, dens_map, N_dens_bins,
 
     N_dens_bins: int
         The number of bins for the range of source density values.
-        The bins will be picked on a linear grid, rangin from the
+        The bins will be picked on a linear grid, ranging from the
         minimum to the maximum density 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.

From 2da48fb4e7506e924699400a9064e846d121d650 Mon Sep 17 00:00:00 2001
From: Benjamin Williams <ben@D-10-18-246-137.dhcp4.washington.edu>
Date: Fri, 1 Jun 2018 14:10:25 -0700
Subject: [PATCH 09/10] AST functions merged

---
 beast/examples/phat_small/datamodel.py        |  14 +-
 beast/examples/phat_small/run_beast.py        |   9 +-
 .../observationmodel/ast/make_ast_xy_list.py  | 329 +-----------------
 3 files changed, 17 insertions(+), 335 deletions(-)

diff --git a/beast/examples/phat_small/datamodel.py b/beast/examples/phat_small/datamodel.py
index 0fa63d773..ab551107a 100755
--- a/beast/examples/phat_small/datamodel.py
+++ b/beast/examples/phat_small/datamodel.py
@@ -46,8 +46,8 @@
 #   need to match column names in the observed catalog,
 #   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 ]
+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
@@ -70,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)
@@ -98,9 +98,9 @@
 # If supplied, the ASTs will be repeated for each background level bin in the table
 ast_background_table = None
 
-#ast_N_dens_bins : (int, optional)
+#ast_N_bins : (int, optional)
 #Number of source or background bins that you want ASTs repeated over
-ast_N_dens_bins = 10
+#ast_N_bins = 8
 
                          
 # ast_pixel_distribution : float (optional)
@@ -114,6 +114,7 @@
 # photometry.	            
 ast_reference_image = None
 
+
 #-------------------------------------------
 #Noise Model Artificial Star Test Parameters
 #-------------------------------------------
@@ -139,9 +140,6 @@
 # Distance unit (any length or units.mag)
 distance_unit = units.mag
 
-# velocity of galaxy
-velocity =  -300 * units.km / units.s # M31 velocity from SIMBAD
-
 ################
 
 ### Stellar grid definition
diff --git a/beast/examples/phat_small/run_beast.py b/beast/examples/phat_small/run_beast.py
index 463ad8ee0..fb6b91b5f 100755
--- a/beast/examples/phat_small/run_beast.py
+++ b/beast/examples/phat_small/run_beast.py
@@ -28,8 +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_per_density
-from beast.observationmodel.ast.make_ast_xy_list import pick_positions_per_background
+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
@@ -152,12 +151,10 @@
                 pick_positions(obsdata, filename, separation,
                                	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_dens_bins,
-                               	outfile=filename, refimage=datamodel.ast_reference_image, Nrealize=1)
+                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_dens_bins,
-                               	outfile=filename, refimage=datamodel.ast_reference_image, Nrealize=1)
+                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)
 
diff --git a/beast/observationmodel/ast/make_ast_xy_list.py b/beast/observationmodel/ast/make_ast_xy_list.py
index 846ae665f..247b73507 100644
--- a/beast/observationmodel/ast/make_ast_xy_list.py
+++ b/beast/observationmodel/ast/make_ast_xy_list.py
@@ -9,7 +9,7 @@
 
 from ...tools.pbar import Pbar
 
-def pick_positions_from_map(chosen_seds, input_map, input_column, N_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 values, 
@@ -58,7 +58,7 @@ def pick_positions_from_map(chosen_seds, input_map, input_column, N_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.
 
@@ -71,16 +71,17 @@ def pick_positions_from_map(chosen_seds, input_map, input_column, N_bins,
     """
     
     # Load the background map
+    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]
-    bins = np.linspace(min_val - 0.01 * abs(min_al),
+    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
@@ -93,6 +94,7 @@ def pick_positions_from_map(chosen_seds, input_map, input_column, N_bins,
 
     # Remove empty bins
     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)
@@ -119,7 +121,7 @@ def pick_positions_from_map(chosen_seds, input_map, input_column, N_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
@@ -172,164 +174,6 @@ def pick_positions_from_map(chosen_seds, input_map, input_column, N_bins,
 
 
 
-def pick_positions_per_background(chosen_seds, bg_map, N_bg_bins,
-                                  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.
-
-    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.
-
-    Then, for each background 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.
-
-    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
-        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
-        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
-        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.
-
-    refimage: str
-        Path to fits image that is used for the positions. If none is
-        given, the ra and dec will be put in the x and y output columns
-        instead.
-
-    Nrealize: integer
-        The number of times each model shoud be repeated for each
-        background regime. This is to sample the variance due to
-        variations within each region, for each individual model.
-
-    Returns
-    -------
-    astropy Table: List of fake stars, with magnitudes and positions
-    - optionally -
-    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)
-
-    # Create the background 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)
-
-    # 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)
-    # 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)]
-
-    # Remove empty bins
-    tile_sets = [tile_set for tile_set in tiles_foreach_bgbin if len(tile_set)]
-
-    # Repeat the seds Nrealize times (sample each on at Nrealize
-    # different positions, in each region)
-    repeated_seds = np.repeat(chosen_seds, Nrealize)
-    Nseds_per_region = len(repeated_seds)
-    # For each set of tiles, repeat the seds and spread them evenly over
-    # the tiles
-    repeated_seds = np.repeat(repeated_seds, len(tile_sets))
-
-    out_table = Table(repeated_seds, names=chosen_seds.colnames)
-    xs = np.zeros(len(out_table))
-    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
-
-    if refimage is None:
-        wcs = None
-    else:
-        imagehdu = fits.open(refimage)[1]
-        wcs = WCS(imagehdu.header)
-
-    pbar = Pbar(len(tile_sets),
-                desc='{} models per background bin'.format(Nseds_per_region))
-    for bin_index, tile_set in pbar.iterover(enumerate(tile_sets)):
-        start = bin_index * Nseds_per_region
-        stop = start + Nseds_per_region
-        bin_indices[start:stop] = bin_index
-        for i in range(Nseds_per_region):
-            x = -1
-            y = -1
-            # Convert each ra,dec to x,y. If there are negative values, try again
-            while x < 0 or y < 0:
-                # Pick a random tile
-                tile = np.random.choice(tile_set)
-                # Within this tile, pick a random ra and dec
-                ra = tile_ra_min[tile] + \
-                     np.random.random_sample() * tile_ra_delta[tile]
-                dec = tile_dec_min[tile] + \
-                      np.random.random_sample() * tile_dec_delta[tile]
-
-                if wcs is None:
-                    x, y = ra, dec
-                    break
-                else:
-                    [x], [y] = wcs.all_world2pix(np.array([ra]), np.array([dec]), 0)
-
-            j = bin_index * Nseds_per_region + i
-            xs[j] = x
-            ys[j] = y
-
-
-    # I'm just mimicking the format that is produced by the examples
-    cs = []
-    cs.append(Column(np.zeros(len(out_table), dtype=int), name='zeros'))
-    cs.append(Column(np.ones(len(out_table), dtype=int), name='ones'))
-
-    if wcs is None:
-        cs.append(Column(xs, name='RA'))
-        cs.append(Column(ys, name='DEC'))
-    else:
-        cs.append(Column(xs, name='X'))
-        cs.append(Column(ys, name='Y'))
-
-    for i, c in enumerate(cs):
-        out_table.add_column(c, index=i)  # insert these columns from the left
-
-    # Write out the table in ascii
-    if outfile:
-        formats = {k: '%.5f' for k in out_table.colnames[2:]}
-        ascii.write(out_table, outfile, overwrite=True, formats=formats)
-
-    return out_table
-
-
 def pick_positions(catalog, filename, separation, refimage=None):
     """
     Assigns positions to fake star list generated by pick_models
@@ -420,160 +264,3 @@ def pick_positions(catalog, filename, separation, refimage=None):
     
     ascii.write(astmags,filename,overwrite=True)
     
-def pick_positions_per_density(chosen_seds, dens_map, N_dens_bins,
-                                  outfile=None, refimage=None, Nrealize=1):
-    """
-    Spreads a set of fake stars across regions of similar source
-    density, given a source density map file generated by 'create
-    source density map' in the tools directory.
-
-    The tiles of the given density map are divided across a given
-    number of source density bins. Each density bin will then
-    have its own set of tiles, which constitute a region on the image.
-
-    Then, for each density 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 stellar crowding, 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
-
-    dens_map: str
-        Path to a fits file containing a density 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
-        source density.
-
-    N_dens_bins: int
-        The number of bins for the range of source density values.
-        The bins will be picked on a linear grid, ranging from the
-        minimum to the maximum density 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.
-
-    refimage: str
-        Path to fits image that is used for the positions. If none is
-        given, the ra and dec will be put in the x and y output columns
-        instead.
-
-    Nrealize: integer
-        The number of times each model shoud be repeated for each
-        background regime. This is to sample the variance due to
-        variations within each region, for each individual model.
-
-    Returns
-    -------
-    astropy Table: List of fake stars, with magnitudes and positions
-    - optionally -
-    ascii file of this table, written to outfile
-
-    """
-
-    # Load the density map
-    dens = Table.read(dens_map)
-    tile_dens_vals = dens['sourcedens']
-    min_dens = np.amin(tile_dens_vals)
-    max_dens = np.amax(tile_dens_vals)
-
-    # Create the sourcedens bins
-    # [min, ., ., ., max]
-    dens_bins = np.linspace(min_dens - 0.01 * abs(min_dens),
-                          max_dens + 0.01 * abs(max_dens), N_dens_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
-    densbin_foreach_tile = np.digitize(tile_dens_vals, dens_bins)
-    # Invert this (the [0] is to dereference the tuple (i,) returned by
-    # nonzero)
-    tiles_foreach_densbin = [np.nonzero(densbin_foreach_tile == b + 1)[0]
-                           for b in range(N_dens_bins)]
-
-    # Remove empty bins
-    tile_sets = [tile_set for tile_set in tiles_foreach_densbin if len(tile_set)]
-
-    # Repeat the seds Nrealize times (sample each on at Nrealize
-    # different positions, in each region)
-    repeated_seds = np.repeat(chosen_seds, Nrealize)
-    Nseds_per_region = len(repeated_seds)
-    # For each set of tiles, repeat the seds and spread them evenly over
-    # the tiles
-    repeated_seds = np.repeat(repeated_seds, len(tile_sets))
-
-    out_table = Table(repeated_seds, names=chosen_seds.colnames)
-    xs = np.zeros(len(out_table))
-    ys = np.zeros(len(out_table))
-    bin_indices = np.zeros(len(out_table))
-
-    tile_ra_min = dens['min_ra']
-    tile_dec_min = dens['min_dec']
-    tile_ra_delta = dens['max_ra'] - tile_ra_min
-    tile_dec_delta = dens['max_dec'] - tile_dec_min
-
-    if refimage is None:
-        wcs = None
-    else:
-        imagehdu = fits.open(refimage)[1]
-        wcs = WCS(imagehdu.header)
-
-
-
-    pbar = Pbar(len(tile_sets),
-                desc='{} models per background bin'.format(Nseds_per_region))
-    for bin_index, tile_set in pbar.iterover(enumerate(tile_sets)):
-        start = bin_index * Nseds_per_region
-        stop = start + Nseds_per_region
-        bin_indices[start:stop] = bin_index
-        for i in range(Nseds_per_region):
-            x = -1
-            y = -1
-            # Convert each ra,dec to x,y. If there are negative values, try again
-            while x < 0 or y < 0:
-                # Pick a random tile
-                tile = np.random.choice(tile_set)
-                # Within this tile, pick a random ra and dec
-                ra = tile_ra_min[tile] + \
-                     np.random.random_sample() * tile_ra_delta[tile]
-                dec = tile_dec_min[tile] + \
-                      np.random.random_sample() * tile_dec_delta[tile]
-
-                if wcs is None:
-                    x, y = ra, dec
-                    break
-                else:
-                    [x], [y] = wcs.all_world2pix(np.array([ra]), np.array([dec]), 0)
-            j = bin_index * Nseds_per_region + i
-            xs[j] = x
-            ys[j] = y
-
-
-    # I'm just mimicking the format that is produced by the examples
-    cs = []
-    cs.append(Column(np.zeros(len(out_table), dtype=int), name='zeros'))
-    cs.append(Column(np.ones(len(out_table), dtype=int), name='ones'))
-
-    if wcs is None:
-        cs.append(Column(xs, name='RA'))
-        cs.append(Column(ys, name='DEC'))
-    else:
-        cs.append(Column(xs, name='X'))
-        cs.append(Column(ys, name='Y'))
-
-    for i, c in enumerate(cs):
-        out_table.add_column(c, index=i)  # insert these columns from the left
-
-    # Write out the table in ascii
-    if outfile:
-        formats = {k: '%.5f' for k in out_table.colnames[2:]}
-        ascii.write(out_table, outfile, overwrite=True, formats=formats)
-
-    return out_table
-

From 697b9b44fa9925e4c6d1028697db128f15310d9d Mon Sep 17 00:00:00 2001
From: Karl Gordon <kgordon@stsci.edu>
Date: Fri, 1 Jun 2018 18:26:48 -0400
Subject: [PATCH 10/10] Keeping the velocity of M31

---
 beast/examples/phat_small/datamodel.py | 3 +++
 1 file changed, 3 insertions(+)

diff --git a/beast/examples/phat_small/datamodel.py b/beast/examples/phat_small/datamodel.py
index ab551107a..c1e655841 100755
--- a/beast/examples/phat_small/datamodel.py
+++ b/beast/examples/phat_small/datamodel.py
@@ -140,6 +140,9 @@
 # Distance unit (any length or units.mag)
 distance_unit = units.mag
 
+# velocity of galaxy	
+velocity =  -300 * units.km / units.s # M31 velocity from SIMBAD	
+
 ################
 
 ### Stellar grid definition