From 6b7e5d24b2b3183d4a48c3ee04701a6091d658b1 Mon Sep 17 00:00:00 2001
From: Joanne Bogart <jrb@slac.stanford.edu>
Date: Wed, 6 Dec 2023 17:32:50 -0800
Subject: [PATCH 1/6] Change options --flux-only to --no-main, main-only to
 --no-flux

---
 skycatalogs/scripts/create_sc.py | 10 +++++-----
 1 file changed, 5 insertions(+), 5 deletions(-)

diff --git a/skycatalogs/scripts/create_sc.py b/skycatalogs/scripts/create_sc.py
index 0503eac0..e48c9ec5 100644
--- a/skycatalogs/scripts/create_sc.py
+++ b/skycatalogs/scripts/create_sc.py
@@ -58,10 +58,10 @@
 
 parser.add_argument('--skip-done', action='store_true',
                     help='If supplied skip existing data files; else overwrite with message')
-parser.add_argument('--flux-only', action='store_true',
-                    help='If supplied only do flux files. Main files must already exist')
-parser.add_argument('--main-only', action='store_true',
-                    help='If supplied only do main files, not flux files')
+parser.add_argument('--no-main', action='store_true',
+                    help='If supplied do not create main files. Note main files must exist in order to create flux files.')
+parser.add_argument('--no-flux', action='store_true',
+                    help='If supplied do not create flux files.')
 parser.add_argument('--flux-parallel', default=16, type=int,
                     help='Number of processes to run in parallel when computing fluxes')
 parser.add_argument('--provenance', '--prov', choices=['yaml'], help='''
@@ -136,7 +136,7 @@
                          knots_mag_cut=args.knots_magnitude_cut,
                          knots=(not args.no_knots), logname=logname,
                          skip_done=args.skip_done,
-                         flux_only=args.flux_only, main_only=args.main_only,
+                         no_main=args.no_main, no_flux=args.no_flux,
                          flux_parallel=args.flux_parallel,
                          galaxy_nside=args.galaxy_nside,
                          galaxy_stride=args.galaxy_stride,

From 8e84b40ea1ae6d73c7b249350c657026bae451cf Mon Sep 17 00:00:00 2001
From: Joanne Bogart <jrb@slac.stanford.edu>
Date: Wed, 6 Dec 2023 17:36:58 -0800
Subject: [PATCH 2/6] Create diffsky SED file  if needed when computing fluxes.
 Also rework DiffskySedGenerator class. No longer subclasses CatalogCreator

---
 skycatalogs/catalog_creator.py |  36 +++--
 skycatalogs/diffsky_sedgen.py  | 240 ++++++++++++---------------------
 2 files changed, 109 insertions(+), 167 deletions(-)

diff --git a/skycatalogs/catalog_creator.py b/skycatalogs/catalog_creator.py
index 7218e700..77c2c9a4 100644
--- a/skycatalogs/catalog_creator.py
+++ b/skycatalogs/catalog_creator.py
@@ -191,8 +191,8 @@ def __init__(self, parts, area_partition=None, skycatalog_root=None,
                  output_type='parquet', mag_cut=None,
                  sed_subdir='galaxyTopHatSED', knots_mag_cut=27.0,
                  knots=True, logname='skyCatalogs.creator',
-                 pkg_root=None, skip_done=False, flux_only=False,
-                 main_only=False, flux_parallel=16, galaxy_nside=32,
+                 pkg_root=None, skip_done=False, no_main=False,
+                 no_flux=False, flux_parallel=16, galaxy_nside=32,
                  galaxy_stride=1000000, provenance=None,
                  dc2=False, sn_object_type='sncosmo', galaxy_type='cosmodc2',
                  include_roman_flux=False, star_input_fmt='sqlite'):
@@ -230,8 +230,8 @@ def __init__(self, parts, area_partition=None, skycatalog_root=None,
         skip_done       If True, skip over files which already exist. Otherwise
                         (by default) overwrite with new version.
                         Output info message in either case if file exists.
-        flux_only       Only create flux files, not main files
-        main_only       Only create main files, not flux files
+        no_main         Do not create main files
+        no_flux         Do not create flux files
         flux_parallel   Number of processes to divide work of computing fluxes
         galaxy_nside    Healpix configuration value "nside" for galaxy output
         galaxy_stride   Max number of rows per galaxy row group
@@ -316,8 +316,8 @@ def __init__(self, parts, area_partition=None, skycatalog_root=None,
         self._logname = logname
         self._logger = logging.getLogger(logname)
         self._skip_done = skip_done
-        self._flux_only = flux_only
-        self._main_only = main_only
+        self._no_main = no_main
+        self._no_flux = no_flux
         self._flux_parallel = flux_parallel
         self._galaxy_nside = galaxy_nside
         self._provenance = provenance
@@ -362,14 +362,14 @@ def create(self, catalog_type):
         None
         """
         if catalog_type == ('galaxy'):
-            if not self._flux_only:
+            if not self._no_main:
                 self.create_galaxy_catalog()
-            if not self._main_only:
+            if not self._no_flux:
                 self.create_galaxy_flux_catalog()
         elif catalog_type == ('pointsource'):
-            if not self._flux_only:
+            if not self._no_main:
                 self.create_pointsource_catalog()
-            if not self._main_only:
+            if not self._no_flux:
                 self.create_pointsource_flux_catalog()
         else:
             raise NotImplementedError(f'CatalogCreator.create: unsupported catalog type {catalog_type}')
@@ -647,6 +647,7 @@ def create_galaxy_flux_catalog(self, config_file=None):
         '''
 
         from .skyCatalogs import open_catalog
+        self._sed_gen = None
 
         self._gal_flux_schema = make_galaxy_flux_schema(self._logname,
                                                         self._galaxy_type,
@@ -666,6 +667,17 @@ def create_galaxy_flux_catalog(self, config_file=None):
         self.object_type = 'galaxy'
         if self._galaxy_type == 'diffsky':
             self.object_type = 'diffsky_galaxy'
+            from .diffsky_sedgen import DiffskySedGenerator
+            # Default values are ok for all the diffsky-specific
+            # parameters: include_nonLSST_flux, sed_parallel, auto_loop,
+            # wave_ang_min, wave_ang_max, rel_err, n_per
+            self._sed_gen = DiffskySedGenerator(logname=self._logname,
+                                                galaxy_truth=self._galaxy_truth,
+                                                output_dir=self._output_dir,
+                                                skip_done=True,
+                                                sky_cat=self._cat)
+
+
         self._flux_template = self._cat.raw_config['object_types'][self.object_type]['flux_file_template']
 
         self._logger.info('Creating galaxy flux files')
@@ -707,6 +719,10 @@ def _create_galaxy_flux_pixel(self, pixel):
                 self._logger.info(f'Skipping regeneration of {output_path}')
                 return
 
+        if self._galaxy_type == 'diffsky':
+            # Generate SEDs if necessary
+            self._sed_gen.generate_pixel(pixel)
+
         # If there are multiple row groups, each is stored in a separate
         # object collection. Need to loop over them
         object_list = self._cat.get_object_type_by_hp(pixel, self.object_type)
diff --git a/skycatalogs/diffsky_sedgen.py b/skycatalogs/diffsky_sedgen.py
index 6ec8f912..66aeb7f2 100644
--- a/skycatalogs/diffsky_sedgen.py
+++ b/skycatalogs/diffsky_sedgen.py
@@ -1,6 +1,7 @@
 import galsim
 import numpy as np
 import os
+import logging
 import h5py
 from lsstdesc_diffsky import read_diffskypop_params
 from lsstdesc_diffsky.io_utils import load_healpixel
@@ -16,10 +17,11 @@
 
 __all__ = ['DiffskySedGenerator']
 
-def _calculate_sed_multi(send_conn,_redshift,_mah_params,_ms_params,_q_params,
-                            _fbulge_params,_fknot,_ssp_data,galaxy_id,n_per):
-    def _calc_seds(_redshift,_mah_params,_ms_params,_q_params,
-                    _fbulge_params,_fknot,_ssp_data,l_bnd,u_bnd,n_per):
+def _calculate_sed_multi(send_conn, _redshift, _mah_params, _ms_params,
+                         _q_params, _fbulge_params, _fknot, _ssp_data,
+                         galaxy_id, n_per):
+    def _calc_seds(_redshift, _mah_params, _ms_params, _q_params,
+                    _fbulge_params, _fknot, _ssp_data, l_bnd, u_bnd, n_per):
         # Documentation for calculation available here:
         # https://lsstdesc-diffsky.readthedocs.io/en/latest/demo_roman_rubin_2023_seds_singlemet.html
         args = (_redshift[l_bnd:u_bnd],
@@ -65,53 +67,41 @@ def _calc_seds(_redshift,_mah_params,_ms_params,_q_params,
         return out_list
 
 
-class DiffskySedGenerator(CatalogCreator):
+class DiffskySedGenerator():
     '''
     Used for evaluating and storing diffsky galaxy SEDs, which are represented with
     an adaptively thinned spectrum and stored in an hdf5 file.
     Parameters
     ----------
-    (all the things in CatalogCreator)
-    sed_parallel    How many processes to split SED calculations over.
+    logname         Where to write log output
+    galaxy_truth    GCRCatalogs name for galaxy truth catalog
+    output_dir      Where diffsky parquet files are
+    sky_cat         To get access to skyCatalogs main files already written
+                    Must be a non-null SkyCatalog object.
     auto_loop       Iterate over all pixels?
     rel_err         Target relative tolerance for flux integral.
     wave_ang_min    Minimum wavelength to keep in SEDs.
     wave_ang_max    Maximum wavelength to keep in SEDs.
     n_per           Number of SEDs to batch calculate in diffsky.
-                    Memory footprint increases nonlinearly with larger n_per.
+                    Memory footprint increases nonlinearly with larger n_per
+    sed_out         If SEDs are to go somewhere other than usual output_dir
     '''
 
-    def __init__(self, parts, area_partition=None, skycatalog_root=None,
-                 catalog_dir='.', galaxy_truth=None,
-                 star_truth=None, sn_truth=None,
-                 config_path=None, catalog_name='skyCatalog',
-                 output_type='parquet', mag_cut=None,
-                 sed_subdir='galaxyTopHatSED', knots_mag_cut=27.0,
-                 knots=True, logname='skyCatalogs.creator',
-                 pkg_root=None, skip_done=False, flux_only=False,
-                 main_only=False, flux_parallel=16, galaxy_nside=32,
-                 galaxy_stride=1000000, provenance=None,
-                 dc2=False, sn_object_type='sncosmo', galaxy_type='cosmodc2',
-                 include_nonLSST_flux=False,sed_parallel=1,auto_loop=False,
-                 wave_ang_min=500,wave_ang_max=100000,rel_err=0.03,n_per=100,sed_out=None):
-        super().__init__(parts, area_partition=area_partition, skycatalog_root=skycatalog_root,
-                 catalog_dir=catalog_dir, galaxy_truth=galaxy_truth,
-                 star_truth=star_truth, sn_truth=sn_truth,
-                 config_path=config_path, catalog_name=catalog_name,
-                 output_type=output_type, mag_cut=mag_cut,
-                 sed_subdir=sed_subdir, knots_mag_cut=knots_mag_cut,
-                 knots=knots, logname=logname,
-                 pkg_root=pkg_root, skip_done=skip_done, flux_only=flux_only,
-                 main_only=main_only, flux_parallel=flux_parallel, galaxy_nside=galaxy_nside,
-                 galaxy_stride=galaxy_stride, provenance=provenance,
-                 dc2=dc2, sn_object_type=sn_object_type, galaxy_type=galaxy_type,
-                 include_nonLSST_flux=include_nonLSST_flux)
-
-        self._sed_parallel = sed_parallel
+    def __init__(self, logname='skyCatalogs.creator', galaxy_truth=None,
+                 output_dir=None, sky_cat=None, skip_done=True,
+                 include_roman_flux=False, auto_loop=False,
+                 wave_ang_min=500, wave_ang_max=100000, rel_err=0.03,
+                 n_per=10000, sed_out=None):       ### temp setting for n_per
+        ## self._galaxy_truth = galaxy_truth
+        self._output_dir = output_dir
+        self._cat = sky_cat
+        self._logger = logging.getLogger(logname)
+        self._skip_done = skip_done
         self._n_per = n_per
         self._sed_out = sed_out
+
         # Setup thinned SSP templates for evaluating SED over
-        # ############### Temporary ############
+        # ############### Maybe temporary ############
         from pathlib import Path
         PACKAGE_SRC_DIR = os.path.dirname(os.path.abspath(str(Path(__file__))))
         SKYCATALOGDATA_ROOT = os.path.join(PACKAGE_SRC_DIR, "data")
@@ -119,10 +109,10 @@ def __init__(self, parts, area_partition=None, skycatalog_root=None,
         # ##############
 
         self._get_thinned_ssp_data(rel_err,wave_ang_min,wave_ang_max,
-                                   SSP_file_name=SINGLE_MET)            # this arg. also temporary
+                                   SSP_file_name=SINGLE_MET)
         import GCRCatalogs
+        gal_cat = GCRCatalogs.load_catalog(galaxy_truth)
 
-        gal_cat = GCRCatalogs.load_catalog(self._galaxy_truth)
         self._hdf5_root_dir = gal_cat.get_catalog_info()['catalog_root_dir']
         self._hdf5_name_template = gal_cat.get_catalog_info()['catalog_filename_template']
 
@@ -132,7 +122,7 @@ def __init__(self, parts, area_partition=None, skycatalog_root=None,
             for p in self._parts:
                 self.iterate_pixel()
 
-    def iterate_pixel(self,pix=None):
+    def iterate_pixel(self, pix=None):
         """
         Iterates to next/target pixel.
         Parameters
@@ -147,7 +137,7 @@ def iterate_pixel(self,pix=None):
         else:
             p = self._parts[pix]
         self._logger.info(f'Starting on pixel {p}')
-        self._generate_pixel(p)
+        self.generate_pixel(p)
         self._logger.info(f'Completed pixel {p}')
         self.pix_iter+=1
 
@@ -207,41 +197,46 @@ def _combine_col(self,cnt,col1,col2,col3):
             tmp[len(col1)+len(col2):len(col1)+len(col2)+len(col3)] = col3
         return tmp
 
-    def _load_diffsky_data(self,pixel):
+    def _load_diffsky_data(self, pixel):
 
         hdf5_file_path = os.path.join(self._hdf5_root_dir,
-            self._hdf5_name_template.format(0,1,pixel))
+                                      self._hdf5_name_template.format(0, 1,
+                                                                      pixel))
         mock1, metadata = load_healpixel(hdf5_file_path)
         diffsky_param_data1 = load_diffsky_params(mock1)
         hdf5_file_path = os.path.join(self._hdf5_root_dir,
-            self._hdf5_name_template.format(1,2,pixel))
+            self._hdf5_name_template.format(1, 2, pixel))
         mock2, metadata = load_healpixel(hdf5_file_path)
         diffsky_param_data2 = load_diffsky_params(mock2)
         hdf5_file_path = os.path.join(self._hdf5_root_dir,
-            self._hdf5_name_template.format(2,3,pixel))
+            self._hdf5_name_template.format(2, 3, pixel))
         mock3, metadata = load_healpixel(hdf5_file_path)
         diffsky_param_data3 = load_diffsky_params(mock3)
         cnt = len(mock1['galaxy_id'])+len(mock2['galaxy_id'])+len(mock3['galaxy_id'])
-        galaxy_id = self._combine_col(cnt,mock1['galaxy_id'],mock2['galaxy_id'],mock3['galaxy_id'])
-        redshift = self._combine_col(cnt,mock1['redshift'],mock2['redshift'],mock3['redshift'])
-        mah_params = self._combine_col(cnt,diffsky_param_data1.mah_params,
-                                            diffsky_param_data2.mah_params,
-                                            diffsky_param_data3.mah_params)
-        ms_params = self._combine_col(cnt,diffsky_param_data1.ms_params,
-                                            diffsky_param_data2.ms_params,
-                                            diffsky_param_data3.ms_params)
-        q_params = self._combine_col(cnt,diffsky_param_data1.q_params,
-                                            diffsky_param_data2.q_params,
-                                            diffsky_param_data3.q_params)
-        fbulge_params = self._combine_col(cnt,diffsky_param_data1.fbulge_params,
-                                            diffsky_param_data2.fbulge_params,
-                                            diffsky_param_data3.fbulge_params)
-        fknot = self._combine_col(cnt,diffsky_param_data1.fknot,diffsky_param_data2.fknot,
-                                            diffsky_param_data3.fknot)
+        galaxy_id = self._combine_col(cnt, mock1['galaxy_id'],
+                                      mock2['galaxy_id'], mock3['galaxy_id'])
+        redshift = self._combine_col(cnt, mock1['redshift'],
+                                     mock2['redshift'], mock3['redshift'])
+        mah_params = self._combine_col(cnt, diffsky_param_data1.mah_params,
+                                       diffsky_param_data2.mah_params,
+                                       diffsky_param_data3.mah_params)
+        ms_params = self._combine_col(cnt, diffsky_param_data1.ms_params,
+                                      diffsky_param_data2.ms_params,
+                                      diffsky_param_data3.ms_params)
+        q_params = self._combine_col(cnt, diffsky_param_data1.q_params,
+                                     diffsky_param_data2.q_params,
+                                     diffsky_param_data3.q_params)
+        fbulge_params = self._combine_col(cnt,
+                                          diffsky_param_data1.fbulge_params,
+                                          diffsky_param_data2.fbulge_params,
+                                          diffsky_param_data3.fbulge_params)
+        fknot = self._combine_col(cnt, diffsky_param_data1.fknot,
+                                  diffsky_param_data2.fknot,
+                                  diffsky_param_data3.fknot)
 
         return galaxy_id,redshift,mah_params,ms_params,q_params,fbulge_params,fknot
 
-    def _generate_pixel(self,pixel):
+    def generate_pixel(self, pixel):
         """
         Iterates to next/target pixel.
         Parameters
@@ -283,10 +278,6 @@ def _generate_pixel(self,pixel):
 
         # If there are multiple row groups, each is stored in a separate
         # object collection. Need to loop over them
-        config_file = self.write_config(path_only=True)
-        if self._cat is None:
-            self._cat = open_catalog(config_file,
-                                 skycatalog_root=self._skycatalog_root)
         object_list = self._cat.get_object_type_by_hp(pixel, 'diffsky_galaxy')
 
         tmp_store = np.zeros((3,len(self.ssp_data.ssp_wave)),dtype='f4')
@@ -294,22 +285,8 @@ def _generate_pixel(self,pixel):
         for object_coll in object_list.get_collections():
             _galaxy_collection = object_coll
             galaxy_id = object_coll.get_native_attribute('galaxy_id')
-
-            l_bnd = 0
-            u_bnd = len(galaxy_id)
             rg_written = 0
-
-            self._logger.debug(f'Handling range {l_bnd} up to {u_bnd}')
-
-            # I couldn't get multiprocessing to work with the jit compilation
-            # in diffsky, but code is setup to use it
-            n_parallel = 1 # self._sed_parallel
-            if n_parallel == 1:
-                n_per = u_bnd - l_bnd
-            else:
-                n_per = int((u_bnd - l_bnd + n_parallel)/n_parallel)
-            lb = l_bnd
-            u = min(l_bnd + n_per, u_bnd)
+            self._logger.debug(f'Handling range 0 to {len(galaxy_id)}')
             readers = []
 
             # Limit objects to those in the matching skycatalog
@@ -317,84 +294,33 @@ def _generate_pixel(self,pixel):
 
             # Build output SED data chunks
             out_dict = []
-            if n_parallel == 1:
-                out_list = _calculate_sed_multi(None,
-                                                redshift[mask],
-                                                mah_params[mask],
-                                                ms_params[mask],
-                                                q_params[mask],
-                                                fbulge_params[mask],
-                                                fknot[mask],
-                                                self.ssp_data,
-                                                diffsky_galaxy_id[mask],
-                                                self._n_per)
-                ichunk=0
-                for chunk in out_list:
-                    print('chunk',ichunk)
-                    ichunk+=1
-                    for igid,gid in enumerate(chunk['galaxy_id']):
-                        tmp_store[0,:] = chunk['bulge'][igid,:]
-                        tmp_store[1,:] = chunk['disk'][igid,:]
-                        tmp_store[2,:] = chunk['knots'][igid,:]
-                        _ = h5_groups[gid//100000].create_dataset(str(gid),
-                                            maxshape=(3,len(self.ssp_data.ssp_wave),),
-                                            shape=(3,len(self.ssp_data.ssp_wave),),
-                                            dtype='f4', #compression="gzip",
-                                            #compression_opts=9,
-                                            data=tmp_store)
-            else:
-                tm = max(int(n_per), 5)  # Give ourselves a cushion
-                self._logger.info(f'Using timeout value {tm} for {n_per} sources')
-                p_list = []
-                for i in range(n_parallel):
-                    conn_rd, conn_wrt = Pipe(duplex=False)
-                    readers.append(conn_rd)
-
-                    # For debugging call directly
-                    proc = Process(target=_calculate_sed_multi,
-                                   name=f'proc_{i}',
-                                   args=(conn_wrt,
-                                        redshift[mask][lb:u],
-                                        mah_params[mask][lb:u],
-                                        ms_params[mask][lb:u],
-                                        q_params[mask][lb:u],
-                                        fbulge_params[mask][lb:u],
-                                        fknot[lb:u],
-                                        self.ssp_data,
-                                        diffsky_galaxy_id[mask][lb:u],
-                                        self._n_per))
-                    proc.start()
-                    p_list.append(proc)
-                    lb = u
-                    u = min(lb + n_per, u_bnd)
-
-                self._logger.debug('Processes started')
-                for i in range(n_parallel):
-                    ready = readers[i].poll(tm)
-                    if not ready:
-                        self._logger.error(f'Process {i} timed out after {tm} sec')
-                        sys.exit(1)
-                    dat = readers[i].recv()
-                    for chunk in dat:
-                        for igid,gid in enumerate(chunk['galaxy_id']):
-                            tmp_store[0,:] = chunk['bulge'][igid,:]
-                            tmp_store[1,:] = chunk['disk'][igid,:]
-                            tmp_store[2,:] = chunk['knots'][igid,:]
-                            _ = h5_groups[gid//100000].create_dataset(str(gid),
-                                                maxshape=(3,len(self.ssp_data.ssp_wave),),
-                                                shape=(3,len(self.ssp_data.ssp_wave),),
-                                                dtype='f4',# compression="gzip",
-                                                #compression_opts=9,
-                                                data=tmp_store)
-                for p in p_list:
-                    p.join()
+
+            out_list = _calculate_sed_multi(None,
+                                            redshift[mask],
+                                            mah_params[mask],
+                                            ms_params[mask],
+                                            q_params[mask],
+                                            fbulge_params[mask],
+                                            fknot[mask],
+                                            self.ssp_data,
+                                            diffsky_galaxy_id[mask],
+                                            self._n_per)
+            ichunk=0
+            for chunk in out_list:
+                print('chunk',ichunk)
+                ichunk+=1
+                for igid,gid in enumerate(chunk['galaxy_id']):
+                    tmp_store[0,:] = chunk['bulge'][igid,:]
+                    tmp_store[1,:] = chunk['disk'][igid,:]
+                    tmp_store[2,:] = chunk['knots'][igid,:]
+                    _ = h5_groups[gid//100000].create_dataset(str(gid),
+                                        maxshape=(3,len(self.ssp_data.ssp_wave),),
+                                        shape=(3,len(self.ssp_data.ssp_wave),),
+                                        dtype='f4', #compression="gzip",
+                                        #compression_opts=9,
+                                        data=tmp_store)
 
             rg_written += 1
 
         f.close()
-        # if self._sed_out is not None:
-        #     output_path_final = os.path.join(self._output_dir, output_filename)
-        #     os.rename(output_path, output_path_final)
-        self._logger.debug(f'# row groups written to flux file: {rg_written}')
-        # if self._provenance == 'yaml':
-        #     self.write_provenance_file(output_path)
+        self._logger.debug(f'SEDs for galaxies in {rg_written} row groups have been written')

From 3deb435901874a0f1f962e177102509737153481 Mon Sep 17 00:00:00 2001
From: Joanne Bogart <jrb@slac.stanford.edu>
Date: Fri, 15 Dec 2023 15:42:50 -0800
Subject: [PATCH 3/6] Add script to just generate diffsky SEDs. Also minor bug
 fix or two in diffsky_sedgen

---
 skycatalogs/catalog_creator.py            |  10 +-
 skycatalogs/diffsky_sedgen.py             |  49 ++++------
 skycatalogs/scripts/create_diffsky_sed.py | 109 ++++++++++++++++++++++
 3 files changed, 134 insertions(+), 34 deletions(-)
 create mode 100644 skycatalogs/scripts/create_diffsky_sed.py

diff --git a/skycatalogs/catalog_creator.py b/skycatalogs/catalog_creator.py
index 77c2c9a4..1c7dccf0 100644
--- a/skycatalogs/catalog_creator.py
+++ b/skycatalogs/catalog_creator.py
@@ -719,13 +719,17 @@ def _create_galaxy_flux_pixel(self, pixel):
                 self._logger.info(f'Skipping regeneration of {output_path}')
                 return
 
+        # If there are multiple row groups, each is stored in a separate
+        # object collection. Need to loop over them
+        object_list = self._cat.get_object_type_by_hp(pixel, self.object_type)
+        if len(object_list) == 0:
+            self._logger.warning(f'Cannot create flux file for pixel {pixel} because main file does not exist or is empty')
+            return
+
         if self._galaxy_type == 'diffsky':
             # Generate SEDs if necessary
             self._sed_gen.generate_pixel(pixel)
 
-        # If there are multiple row groups, each is stored in a separate
-        # object collection. Need to loop over them
-        object_list = self._cat.get_object_type_by_hp(pixel, self.object_type)
         writer = None
         global _galaxy_collection
         global _instrument_needed
diff --git a/skycatalogs/diffsky_sedgen.py b/skycatalogs/diffsky_sedgen.py
index 66aeb7f2..fbaded89 100644
--- a/skycatalogs/diffsky_sedgen.py
+++ b/skycatalogs/diffsky_sedgen.py
@@ -9,7 +9,6 @@
 from lsstdesc_diffsky.legacy.roman_rubin_2023.dsps.data_loaders.load_ssp_data import load_ssp_templates_singlemet
 from lsstdesc_diffsky.legacy.roman_rubin_2023.dsps.data_loaders.defaults import SSPDataSingleMet
 from lsstdesc_diffsky.defaults import OUTER_RIM_COSMO_PARAMS
-#from lsstdesc_diffsky.sed.disk_bulge_sed_kernels_singlemet import calc_rest_sed_disk_bulge_knot_singlegal
 from lsstdesc_diffsky.sed.disk_bulge_sed_kernels_singlemet import calc_rest_sed_disk_bulge_knot_galpop
 all_diffskypop_params = read_diffskypop_params("roman_rubin_2023")
 from .catalog_creator import CatalogCreator
@@ -78,21 +77,24 @@ class DiffskySedGenerator():
     output_dir      Where diffsky parquet files are
     sky_cat         To get access to skyCatalogs main files already written
                     Must be a non-null SkyCatalog object.
-    auto_loop       Iterate over all pixels?
+    skip_done       If false, overwrite existing files
+    auto_loop       If true, immediately loop through pixels in 'parts'
     rel_err         Target relative tolerance for flux integral.
     wave_ang_min    Minimum wavelength to keep in SEDs.
     wave_ang_max    Maximum wavelength to keep in SEDs.
     n_per           Number of SEDs to batch calculate in diffsky.
                     Memory footprint increases nonlinearly with larger n_per
     sed_out         If SEDs are to go somewhere other than usual output_dir
+    parts           Pixels for which SEDs are created (only used if auto_loop
+                    is True
     '''
 
     def __init__(self, logname='skyCatalogs.creator', galaxy_truth=None,
                  output_dir=None, sky_cat=None, skip_done=True,
-                 include_roman_flux=False, auto_loop=False,
-                 wave_ang_min=500, wave_ang_max=100000, rel_err=0.03,
-                 n_per=10000, sed_out=None):       ### temp setting for n_per
-        ## self._galaxy_truth = galaxy_truth
+                 auto_loop=False,
+                 wave_ang_min=500, wave_ang_max=100000,
+                 rel_err=0.03, n_per=100000,
+                 sed_out=None, parts=None):
         self._output_dir = output_dir
         self._cat = sky_cat
         self._logger = logging.getLogger(logname)
@@ -116,30 +118,10 @@ def __init__(self, logname='skyCatalogs.creator', galaxy_truth=None,
         self._hdf5_root_dir = gal_cat.get_catalog_info()['catalog_root_dir']
         self._hdf5_name_template = gal_cat.get_catalog_info()['catalog_filename_template']
 
-        self.pix_iter = 0
         if auto_loop:
             # Loop over necessary pixels
             for p in self._parts:
-                self.iterate_pixel()
-
-    def iterate_pixel(self, pix=None):
-        """
-        Iterates to next/target pixel.
-        Parameters
-        ----------
-        pix     Optional pixel list index.
-        """
-        if pix is None:
-            p = self._parts[self.pix_iter]
-        elif pix>=len(self._parts):
-            self._logger.error(f'Tried to calculate for pixel index {pix} outsid pixel list.')
-            return
-        else:
-            p = self._parts[pix]
-        self._logger.info(f'Starting on pixel {p}')
-        self.generate_pixel(p)
-        self._logger.info(f'Completed pixel {p}')
-        self.pix_iter+=1
+                self.generate_pixel(p)
 
     def _get_thinned_ssp_data(  self,
                                 rel_err,
@@ -238,12 +220,18 @@ def _load_diffsky_data(self, pixel):
 
     def generate_pixel(self, pixel):
         """
-        Iterates to next/target pixel.
         Parameters
         ----------
         pixel     Pixel to generate SEDs for.
         """
 
+        # Check that main file for this pixel exists and has some
+        # data for us  If not, issue warning and return.
+        object_list = self._cat.get_object_type_by_hp(pixel, 'diffsky_galaxy')
+        if len(object_list) == 0:
+            self._logger.warning(f'No sky catalog data for pixel {pixel}\nCannot create SEDs')
+            return
+
         # Load diffsky galaxy data
         diffsky_galaxy_id,redshift,mah_params, \
         ms_params,q_params,fbulge_params,fknot \
@@ -276,12 +264,11 @@ def generate_pixel(self, pixel):
         for g in np.unique(diffsky_galaxy_id//100000):
             h5_groups[g]=f.create_group('galaxy/'+str(g))
 
-        # If there are multiple row groups, each is stored in a separate
-        # object collection. Need to loop over them
-        object_list = self._cat.get_object_type_by_hp(pixel, 'diffsky_galaxy')
 
         tmp_store = np.zeros((3,len(self.ssp_data.ssp_wave)),dtype='f4')
         # Loop over object collections to do SED calculations
+        # If there are multiple row groups, each is stored in a separate
+        # object collection. Need to loop over them
         for object_coll in object_list.get_collections():
             _galaxy_collection = object_coll
             galaxy_id = object_coll.get_native_attribute('galaxy_id')
diff --git a/skycatalogs/scripts/create_diffsky_sed.py b/skycatalogs/scripts/create_diffsky_sed.py
new file mode 100644
index 00000000..2421c165
--- /dev/null
+++ b/skycatalogs/scripts/create_diffsky_sed.py
@@ -0,0 +1,109 @@
+# Use DiffskySedGenerator class
+# For each healpixel for which SEDs are to be created, main skyCatalogs
+# parquet file must already exist
+'''
+Create sky catalogs for one or more healpixels. Invoke with --help
+for details
+'''
+
+import os
+import numpy as np
+import argparse
+import logging
+import yaml
+from skycatalogs.utils.common_utils import print_date, log_callinfo
+
+parser = argparse.ArgumentParser(description='''
+Create SEDs for diffsky galaxies.''',
+                                 formatter_class=argparse.ArgumentDefaultsHelpFormatter)
+parser.add_argument('--pixels', type=int, nargs='*', default=[9556],
+                    help='healpix pixels for which catalogs will be created')
+parser.add_argument('--config-path', help='path to a skyCatalogs config file')
+parser.add_argument('--skycatalog-root',
+                    help='''Root directory for sky catalogs, typically
+                    site- or user-dependent. If not specified, use value of
+                    environment variable SKYCATALOG_ROOT or value from config
+                    file''')
+parser.add_argument('--catalog-dir', '--cat-dir',
+                    help='output file directory relative to skycatalog_root',
+                    default='.')
+parser.add_argument('--output-dir',
+                    help='''If specified output SEDs here. Else write to
+                    dir used for diffsky catalogs''')
+parser.add_argument('--log-level', help='controls logging output',
+                    default='INFO', choices=['DEBUG', 'INFO', 'WARNING',
+                                             'ERROR'])
+parser.add_argument('--galaxy-truth', help='''GCRCatalogs name for galaxy
+                    truth catalog.  If not specified, default will be used''')
+parser.add_argument('--rel-err', type=float, default=0.03, help='''
+                     target tolerance for flux error''')
+parser.add_argument('--wave-ang-min', type=int, default=500,
+                     help='''Min wavelength to keep in SEDs''')
+parser.add_argument('--wave-ang-max', type=int, default=100000,
+                     help='''Max wavelength to keep in SEDs''')
+parser.add_argument('--n-per', type=int, default=100000, help='''
+                    number of galaxies to be processed together''')
+parser.add_argument('--overwrite', action='store_true',
+                    help='''If supplied overwrite existing data files;
+                    else skip with message''')
+parser.add_argument('--options-file', default=None, help='''
+                    path to yaml file associating option names with values.
+                    Values for any options included will take precedence.''')
+
+args = parser.parse_args()
+
+if args.options_file:
+    with open(args.options_file) as f:
+        opt_dict = yaml.safe_load(f)
+        for k in opt_dict:
+            if k in args:
+                args.__setattr__(k, opt_dict[k])
+            else:
+                raise ValueError(f'Unknown attribute "{k}" in options file {args.options_file}')
+logname = 'diffsky_sed.creator'
+logger = logging.getLogger(logname)
+logger.setLevel(args.log_level)
+
+ch = logging.StreamHandler()
+ch.setLevel(args.log_level)
+formatter = logging.Formatter('%(asctime)s - %(levelname)s - %(message)s')
+ch.setFormatter(formatter)
+
+logger.addHandler(ch)
+
+log_callinfo('create_diffsky_sed', args, logname)
+
+skycatalog_root = args.skycatalog_root
+if not skycatalog_root:
+    skycatalog_root = os.getenv('SKYCATALOG_ROOT')
+
+from skycatalogs.skyCatalogs import open_catalog
+
+sky_cat = open_catalog(args.config_path, skycatalog_root=skycatalog_root)
+parts = args.pixels
+# if args.provenance:
+#     provenance = args.provenance
+# else:
+#     provenance = None
+
+from skycatalogs.diffsky_sedgen import DiffskySedGenerator
+
+# hard-code for now.  Should be able to retrieve galaxy truth from sky_cat
+galaxy_truth = args.galaxy_truth
+if galaxy_truth is None:
+    galaxy_truth = 'roman_rubin_2023_v1.1.1_elais'
+
+creator = DiffskySedGenerator(logname=logname, galaxy_truth=galaxy_truth,
+                              output_dir=args.output_dir, sky_cat=sky_cat,
+                              skip_done=(not args.overwrite),
+                              rel_err=args.rel_err,
+                              wave_ang_min=args.wave_ang_min,
+                              wave_ang_max=args.wave_ang_max,
+                              n_per=args.n_per, sed_out=args.output_dir)
+
+for p in args.pixels:
+    creator.generate_pixel(p)
+    logger.info(f'Done with pixel {p}')
+
+logger.info('All done')
+print_date()

From 57153df5b88464ac79140654b529385fe396ed33 Mon Sep 17 00:00:00 2001
From: Joanne Bogart <jrb@slac.stanford.edu>
Date: Mon, 18 Dec 2023 14:08:29 -0800
Subject: [PATCH 4/6] cosmetic changes for flake8

---
 skycatalogs/catalog_creator.py            |  31 +++---
 skycatalogs/diffsky_sedgen.py             | 128 +++++++++++-----------
 skycatalogs/scripts/create_diffsky_sed.py |   7 +-
 3 files changed, 84 insertions(+), 82 deletions(-)

diff --git a/skycatalogs/catalog_creator.py b/skycatalogs/catalog_creator.py
index 1c7dccf0..7ed77928 100644
--- a/skycatalogs/catalog_creator.py
+++ b/skycatalogs/catalog_creator.py
@@ -34,6 +34,7 @@
 _MW_rv_constant = 3.1
 _nside_allowed = 2**np.arange(15)
 
+
 def _get_tophat_info(columns):
     '''
     Parameters
@@ -150,7 +151,8 @@ def _generate_subpixel_masks(ra, dec, subpixels, nside=32):
 # Collection of galaxy objects for current row group, current pixel
 # Used while doing flux computation
 
-def _do_galaxy_flux_chunk(send_conn, galaxy_collection, instrument_needed, l_bnd, u_bnd):
+def _do_galaxy_flux_chunk(send_conn, galaxy_collection, instrument_needed,
+                          l_bnd, u_bnd):
     '''
     output connection
     l_bnd, u_bnd     demarcates slice to process
@@ -501,8 +503,8 @@ def create_galaxy_pixel(self, pixel, gal_cat, arrow_schema):
             # to_fetch = all columns of interest in gal_cat
             non_sed = ['galaxy_id', 'ra', 'dec', 'redshift', 'redshiftHubble',
                        'peculiarVelocity', 'shear_1', 'shear_2',
-                       'convergence',
-                       'size_bulge_true', 'size_minor_bulge_true', 'sersic_bulge',
+                       'convergence', 'size_bulge_true',
+                       'size_minor_bulge_true', 'sersic_bulge',
                        'size_disk_true', 'size_minor_disk_true', 'sersic_disk']
             if self._dc2:
                 non_sed += ['ellipticity_1_disk_true_dc2',
@@ -535,7 +537,7 @@ def create_galaxy_pixel(self, pixel, gal_cat, arrow_schema):
                         'convergence', 'diskEllipticity1', 'diskEllipticity2',
                         'spheroidEllipticity1', 'spheroidEllipticity2',
                         'spheroidHalfLightRadiusArcsec',
-                        'diskHalfLightRadiusArcsec','um_source_galaxy_obs_sm']
+                        'diskHalfLightRadiusArcsec', 'um_source_galaxy_obs_sm']
 
         # df is not a dataframe!  It's just a dict
         if not self._mag_cut:
@@ -562,12 +564,14 @@ def create_galaxy_pixel(self, pixel, gal_cat, arrow_schema):
                 to_rename['ellipticity_2_bulge_true_dc2'] = 'ellipticity_2_bulge_true'
 
             if self._sed_subdir:
-                #  Generate full paths for disk and bulge SED files, even though
+                #  Generate full paths for disk and bulge SED files even though
                 #  we don't actually write the files here
                 df['bulge_sed_file_path'] =\
-                    generate_sed_path(df['galaxy_id'], self._sed_subdir, 'bulge')
+                    generate_sed_path(df['galaxy_id'], self._sed_subdir,
+                                      'bulge')
                 df['disk_sed_file_path'] =\
-                    generate_sed_path(df['galaxy_id'], self._sed_subdir, 'disk')
+                    generate_sed_path(df['galaxy_id'], self._sed_subdir,
+                                      'disk')
 
             if self._knots:
                 # adjust disk sed; create knots sed
@@ -605,7 +609,8 @@ def create_galaxy_pixel(self, pixel, gal_cat, arrow_schema):
 
                 if self._galaxy_type == 'cosmodc2':
                     compressed = self._make_tophat_columns(compressed,
-                                                           sed_disk_names, 'disk')
+                                                           sed_disk_names,
+                                                           'disk')
                     compressed = self._make_tophat_columns(compressed,
                                                            sed_bulge_names,
                                                            'bulge')
@@ -649,9 +654,9 @@ def create_galaxy_flux_catalog(self, config_file=None):
         from .skyCatalogs import open_catalog
         self._sed_gen = None
 
-        self._gal_flux_schema = make_galaxy_flux_schema(self._logname,
-                                                        self._galaxy_type,
-                                                        include_roman_flux=self._include_roman_flux)
+        self._gal_flux_schema =\
+            make_galaxy_flux_schema(self._logname, self._galaxy_type,
+                                    include_roman_flux=self._include_roman_flux)
         self._gal_flux_needed = [field.name for field in self._gal_flux_schema]
 
         if not config_file:
@@ -677,7 +682,6 @@ def create_galaxy_flux_catalog(self, config_file=None):
                                                 skip_done=True,
                                                 sky_cat=self._cat)
 
-
         self._flux_template = self._cat.raw_config['object_types'][self.object_type]['flux_file_template']
 
         self._logger.info('Creating galaxy flux files')
@@ -780,7 +784,8 @@ def _create_galaxy_flux_pixel(self, pixel):
                     # For debugging call directly
                     proc = Process(target=_do_galaxy_flux_chunk,
                                    name=f'proc_{i}',
-                                   args=(conn_wrt, _galaxy_collection, _instrument_needed, lb, u))
+                                   args=(conn_wrt, _galaxy_collection,
+                                         _instrument_needed, lb, u))
                     proc.start()
                     p_list.append(proc)
                     lb = u
diff --git a/skycatalogs/diffsky_sedgen.py b/skycatalogs/diffsky_sedgen.py
index fbaded89..baa123f6 100644
--- a/skycatalogs/diffsky_sedgen.py
+++ b/skycatalogs/diffsky_sedgen.py
@@ -11,16 +11,15 @@
 from lsstdesc_diffsky.defaults import OUTER_RIM_COSMO_PARAMS
 from lsstdesc_diffsky.sed.disk_bulge_sed_kernels_singlemet import calc_rest_sed_disk_bulge_knot_galpop
 all_diffskypop_params = read_diffskypop_params("roman_rubin_2023")
-from .catalog_creator import CatalogCreator
-from .skyCatalogs import open_catalog
 
 __all__ = ['DiffskySedGenerator']
 
+
 def _calculate_sed_multi(send_conn, _redshift, _mah_params, _ms_params,
                          _q_params, _fbulge_params, _fknot, _ssp_data,
                          galaxy_id, n_per):
     def _calc_seds(_redshift, _mah_params, _ms_params, _q_params,
-                    _fbulge_params, _fknot, _ssp_data, l_bnd, u_bnd, n_per):
+                   _fbulge_params, _fknot, _ssp_data, l_bnd, u_bnd, n_per):
         # Documentation for calculation available here:
         # https://lsstdesc-diffsky.readthedocs.io/en/latest/demo_roman_rubin_2023_seds_singlemet.html
         args = (_redshift[l_bnd:u_bnd],
@@ -42,21 +41,21 @@ def _calc_seds(_redshift, _mah_params, _ms_params, _q_params,
     flux_factor = 4.0204145742268754e-16
 
     # Iterate over chunks of n_per
-    l_bnd=0
-    u_bnd=len(galaxy_id)
+    l_bnd = 0
+    u_bnd = len(galaxy_id)
     out_list = []
     lb = l_bnd
     u = min(l_bnd + n_per, u_bnd)
     for i in range((u_bnd-l_bnd)//n_per+1):
-        if lb==u:
+        if lb == u:
             break
-        sed_info = _calc_seds(_redshift,_mah_params,_ms_params,_q_params,
-                                _fbulge_params,_fknot,_ssp_data,lb,u,n_per)
+        sed_info = _calc_seds(_redshift, _mah_params, _ms_params, _q_params,
+                              _fbulge_params, _fknot, _ssp_data, lb, u, n_per)
         # Accumulate output chunks
-        out_list.append({'galaxy_id':galaxy_id[lb:u],
-                        'bulge':sed_info.rest_sed_bulge*flux_factor,
-                        'disk':sed_info.rest_sed_diffuse_disk*flux_factor,
-                        'knots':sed_info.rest_sed_knot*flux_factor})
+        out_list.append({'galaxy_id': galaxy_id[lb:u],
+                         'bulge': sed_info.rest_sed_bulge*flux_factor,
+                         'disk': sed_info.rest_sed_diffuse_disk*flux_factor,
+                         'knots': sed_info.rest_sed_knot*flux_factor})
         lb = u
         u = min(lb + n_per, u_bnd)
 
@@ -68,8 +67,8 @@ def _calc_seds(_redshift, _mah_params, _ms_params, _q_params,
 
 class DiffskySedGenerator():
     '''
-    Used for evaluating and storing diffsky galaxy SEDs, which are represented with
-    an adaptively thinned spectrum and stored in an hdf5 file.
+    Used for evaluating and storing diffsky galaxy SEDs, which are
+    represented with an adaptively thinned spectrum and stored in an hdf5 file.
     Parameters
     ----------
     logname         Where to write log output
@@ -107,10 +106,11 @@ def __init__(self, logname='skyCatalogs.creator', galaxy_truth=None,
         from pathlib import Path
         PACKAGE_SRC_DIR = os.path.dirname(os.path.abspath(str(Path(__file__))))
         SKYCATALOGDATA_ROOT = os.path.join(PACKAGE_SRC_DIR, "data")
-        SINGLE_MET = os.path.join(SKYCATALOGDATA_ROOT, "dsps_ssp_data_singlemet.h5")
+        SINGLE_MET = os.path.join(SKYCATALOGDATA_ROOT,
+                                  "dsps_ssp_data_singlemet.h5")
         # ##############
 
-        self._get_thinned_ssp_data(rel_err,wave_ang_min,wave_ang_max,
+        self._get_thinned_ssp_data(rel_err, wave_ang_min, wave_ang_max,
                                    SSP_file_name=SINGLE_MET)
         import GCRCatalogs
         gal_cat = GCRCatalogs.load_catalog(galaxy_truth)
@@ -123,11 +123,11 @@ def __init__(self, logname='skyCatalogs.creator', galaxy_truth=None,
             for p in self._parts:
                 self.generate_pixel(p)
 
-    def _get_thinned_ssp_data(  self,
-                                rel_err,
-                                wave_ang_min,
-                                wave_ang_max,
-                                SSP_file_name='dsps_ssp_data_singlemet.h5'):
+    def _get_thinned_ssp_data(self,
+                              rel_err,
+                              wave_ang_min,
+                              wave_ang_max,
+                              SSP_file_name='dsps_ssp_data_singlemet.h5'):
         """
         Return thinned SSP templates.
         Parameters
@@ -153,7 +153,7 @@ def _get_thinned_ssp_data(  self,
             x=ssp_wave_nm[mask], f=np.sum(ssp_data.ssp_flux[:, mask], axis=0)
         )
         sed = galsim.SED(sed_lut, wave_type="nm",
-                        flux_type="flambda", redshift=0.0)
+                         flux_type="flambda", redshift=0.0)
 
         # Thin template sum to target tolerance.
         sed2 = sed.thin(rel_err=rel_err, fast_search=False)
@@ -164,16 +164,16 @@ def _get_thinned_ssp_data(  self,
         thin_ssp_wave_ang = ssp_data.ssp_wave[mask][mask2]
         thin_ssp_flux = ssp_data.ssp_flux[:, mask][:, mask2]
         self.ssp_data = SSPDataSingleMet(ssp_data.ssp_lg_age_gyr,
-                                        thin_ssp_wave_ang, thin_ssp_flux)
-
-    def _combine_col(self,cnt,col1,col2,col3):
-        if len(np.shape(col1))>1:
-            tmp = np.zeros((cnt,np.shape(col1)[1]),dtype=col1.dtype)
-            tmp[0:len(col1),:] = col1
-            tmp[len(col1):len(col1)+len(col2),:] = col2
-            tmp[len(col1)+len(col2):len(col1)+len(col2)+len(col3),:] = col3
+                                         thin_ssp_wave_ang, thin_ssp_flux)
+
+    def _combine_col(self, cnt, col1, col2, col3):
+        if len(np.shape(col1)) > 1:
+            tmp = np.zeros((cnt, np.shape(col1)[1]), dtype=col1.dtype)
+            tmp[0:len(col1), :] = col1
+            tmp[len(col1):len(col1)+len(col2), :] = col2
+            tmp[len(col1)+len(col2):len(col1)+len(col2)+len(col3), :] = col3
         else:
-            tmp = np.zeros(cnt,dtype=col1.dtype)
+            tmp = np.zeros(cnt, dtype=col1.dtype)
             tmp[0:len(col1)] = col1
             tmp[len(col1):len(col1)+len(col2)] = col2
             tmp[len(col1)+len(col2):len(col1)+len(col2)+len(col3)] = col3
@@ -187,11 +187,13 @@ def _load_diffsky_data(self, pixel):
         mock1, metadata = load_healpixel(hdf5_file_path)
         diffsky_param_data1 = load_diffsky_params(mock1)
         hdf5_file_path = os.path.join(self._hdf5_root_dir,
-            self._hdf5_name_template.format(1, 2, pixel))
+                                      self._hdf5_name_template.format(1, 2,
+                                                                      pixel))
         mock2, metadata = load_healpixel(hdf5_file_path)
         diffsky_param_data2 = load_diffsky_params(mock2)
         hdf5_file_path = os.path.join(self._hdf5_root_dir,
-            self._hdf5_name_template.format(2, 3, pixel))
+                                      self._hdf5_name_template.format(2, 3,
+                                                                      pixel))
         mock3, metadata = load_healpixel(hdf5_file_path)
         diffsky_param_data3 = load_diffsky_params(mock3)
         cnt = len(mock1['galaxy_id'])+len(mock2['galaxy_id'])+len(mock3['galaxy_id'])
@@ -216,7 +218,7 @@ def _load_diffsky_data(self, pixel):
                                   diffsky_param_data2.fknot,
                                   diffsky_param_data3.fknot)
 
-        return galaxy_id,redshift,mah_params,ms_params,q_params,fbulge_params,fknot
+        return galaxy_id, redshift, mah_params, ms_params, q_params, fbulge_params, fknot
 
     def generate_pixel(self, pixel):
         """
@@ -233,8 +235,8 @@ def generate_pixel(self, pixel):
             return
 
         # Load diffsky galaxy data
-        diffsky_galaxy_id,redshift,mah_params, \
-        ms_params,q_params,fbulge_params,fknot \
+        diffsky_galaxy_id, redshift, mah_params, \
+            ms_params, q_params, fbulge_params, fknot \
             = self._load_diffsky_data(pixel)
 
         # Set up output file
@@ -252,36 +254,31 @@ def generate_pixel(self, pixel):
                 return
 
         # Set up h5 file and metadata on wavelengths
-        f = h5py.File(output_path,'w',libver='latest')
+        f = h5py.File(output_path, 'w', libver='latest')
         _ = f.create_dataset('meta/wave_list',
-                            maxshape=(len(self.ssp_data.ssp_wave),),
-                            shape=(len(self.ssp_data.ssp_wave),),
-                            dtype='f4', #compression="gzip",
-                            #compression_opts=9,
-                            data=self.ssp_data.ssp_wave)
+                             maxshape=(len(self.ssp_data.ssp_wave),),
+                             shape=(len(self.ssp_data.ssp_wave),),
+                             dtype='f4',  # compression="gzip",
+                             # compression_opts=9,
+                             data=self.ssp_data.ssp_wave)
         # Set up hdf5 groups
-        h5_groups={}
+        h5_groups = {}
         for g in np.unique(diffsky_galaxy_id//100000):
-            h5_groups[g]=f.create_group('galaxy/'+str(g))
+            h5_groups[g] = f.create_group('galaxy/'+str(g))
 
-
-        tmp_store = np.zeros((3,len(self.ssp_data.ssp_wave)),dtype='f4')
+        tmp_store = np.zeros((3, len(self.ssp_data.ssp_wave)), dtype='f4')
         # Loop over object collections to do SED calculations
         # If there are multiple row groups, each is stored in a separate
         # object collection. Need to loop over them
         for object_coll in object_list.get_collections():
-            _galaxy_collection = object_coll
             galaxy_id = object_coll.get_native_attribute('galaxy_id')
             rg_written = 0
             self._logger.debug(f'Handling range 0 to {len(galaxy_id)}')
-            readers = []
 
             # Limit objects to those in the matching skycatalog
-            mask              = np.in1d(diffsky_galaxy_id,galaxy_id)
+            mask = np.in1d(diffsky_galaxy_id, galaxy_id)
 
             # Build output SED data chunks
-            out_dict = []
-
             out_list = _calculate_sed_multi(None,
                                             redshift[mask],
                                             mah_params[mask],
@@ -292,21 +289,22 @@ def generate_pixel(self, pixel):
                                             self.ssp_data,
                                             diffsky_galaxy_id[mask],
                                             self._n_per)
-            ichunk=0
+            ichunk = 0
             for chunk in out_list:
-                print('chunk',ichunk)
-                ichunk+=1
-                for igid,gid in enumerate(chunk['galaxy_id']):
-                    tmp_store[0,:] = chunk['bulge'][igid,:]
-                    tmp_store[1,:] = chunk['disk'][igid,:]
-                    tmp_store[2,:] = chunk['knots'][igid,:]
-                    _ = h5_groups[gid//100000].create_dataset(str(gid),
-                                        maxshape=(3,len(self.ssp_data.ssp_wave),),
-                                        shape=(3,len(self.ssp_data.ssp_wave),),
-                                        dtype='f4', #compression="gzip",
-                                        #compression_opts=9,
-                                        data=tmp_store)
-
+                print('chunk', ichunk)
+                ichunk += 1
+                for igid, gid in enumerate(chunk['galaxy_id']):
+                    tmp_store[0, :] = chunk['bulge'][igid, :]
+                    tmp_store[1, :] = chunk['disk'][igid, :]
+                    tmp_store[2, :] = chunk['knots'][igid, :]
+                    _ = h5_groups[gid//100000].\
+                        create_dataset(str(gid),
+                                       maxshape=(3,
+                                                 len(self.ssp_data.ssp_wave),),
+                                       shape=(3, len(self.ssp_data.ssp_wave),),
+                                       dtype='f4',  # compression="gzip",
+                                       # compression_opts=9,
+                                       data=tmp_store)
             rg_written += 1
 
         f.close()
diff --git a/skycatalogs/scripts/create_diffsky_sed.py b/skycatalogs/scripts/create_diffsky_sed.py
index 2421c165..fa22d6be 100644
--- a/skycatalogs/scripts/create_diffsky_sed.py
+++ b/skycatalogs/scripts/create_diffsky_sed.py
@@ -7,7 +7,6 @@
 '''
 
 import os
-import numpy as np
 import argparse
 import logging
 import yaml
@@ -36,11 +35,11 @@
 parser.add_argument('--galaxy-truth', help='''GCRCatalogs name for galaxy
                     truth catalog.  If not specified, default will be used''')
 parser.add_argument('--rel-err', type=float, default=0.03, help='''
-                     target tolerance for flux error''')
+                    target tolerance for flux error''')
 parser.add_argument('--wave-ang-min', type=int, default=500,
-                     help='''Min wavelength to keep in SEDs''')
+                    help='''Min wavelength to keep in SEDs''')
 parser.add_argument('--wave-ang-max', type=int, default=100000,
-                     help='''Max wavelength to keep in SEDs''')
+                    help='''Max wavelength to keep in SEDs''')
 parser.add_argument('--n-per', type=int, default=100000, help='''
                     number of galaxies to be processed together''')
 parser.add_argument('--overwrite', action='store_true',

From ffe92364e322464ad6adf45eeec457a13ffdfde2 Mon Sep 17 00:00:00 2001
From: Joanne Bogart <jrb@slac.stanford.edu>
Date: Mon, 18 Dec 2023 18:00:50 -0800
Subject: [PATCH 5/6] use sensible output dir for diffsky SEDs when none is
 passed in as argument

---
 skycatalogs/catalog_creator.py            | 2 ++
 skycatalogs/diffsky_sedgen.py             | 2 ++
 skycatalogs/scripts/create_diffsky_sed.py | 9 +++------
 3 files changed, 7 insertions(+), 6 deletions(-)

diff --git a/skycatalogs/catalog_creator.py b/skycatalogs/catalog_creator.py
index 7ed77928..b45d41a9 100644
--- a/skycatalogs/catalog_creator.py
+++ b/skycatalogs/catalog_creator.py
@@ -251,6 +251,8 @@ def __init__(self, parts, area_partition=None, skycatalog_root=None,
 
         _cosmo_cat = 'cosmodc2_v1.1.4_image_addon_knots'
         _diffsky_cat = 'roman_rubin_2023_v1.1.1_elais'
+        # When it's officially available, default will be newer version:
+        # _diffsky_cat = 'roman_rubin_2023_v1.1.2_elais'
         _star_db = '/global/cfs/cdirs/lsst/groups/SSim/DC2/dc2_stellar_healpixel.db'
         _sn_db = '/global/cfs/cdirs/lsst/groups/SSim/DC2/cosmoDC2_v1.1.4/sne_cosmoDC2_v1.1.4_MS_DDF_healpix.db'
 
diff --git a/skycatalogs/diffsky_sedgen.py b/skycatalogs/diffsky_sedgen.py
index baa123f6..d8eae214 100644
--- a/skycatalogs/diffsky_sedgen.py
+++ b/skycatalogs/diffsky_sedgen.py
@@ -118,6 +118,8 @@ def __init__(self, logname='skyCatalogs.creator', galaxy_truth=None,
         self._hdf5_root_dir = gal_cat.get_catalog_info()['catalog_root_dir']
         self._hdf5_name_template = gal_cat.get_catalog_info()['catalog_filename_template']
 
+        if self._output_dir is None:
+            self._output_dir = sky_cat._cat_dir
         if auto_loop:
             # Loop over necessary pixels
             for p in self._parts:
diff --git a/skycatalogs/scripts/create_diffsky_sed.py b/skycatalogs/scripts/create_diffsky_sed.py
index fa22d6be..555d0f66 100644
--- a/skycatalogs/scripts/create_diffsky_sed.py
+++ b/skycatalogs/scripts/create_diffsky_sed.py
@@ -79,18 +79,15 @@
 from skycatalogs.skyCatalogs import open_catalog
 
 sky_cat = open_catalog(args.config_path, skycatalog_root=skycatalog_root)
-parts = args.pixels
-# if args.provenance:
-#     provenance = args.provenance
-# else:
-#     provenance = None
-
 from skycatalogs.diffsky_sedgen import DiffskySedGenerator
 
 # hard-code for now.  Should be able to retrieve galaxy truth from sky_cat
 galaxy_truth = args.galaxy_truth
 if galaxy_truth is None:
     galaxy_truth = 'roman_rubin_2023_v1.1.1_elais'
+    # Change to
+    #  galaxy_truth = 'roman_rubin_2023_v1.1.2_elais'
+    # when v1.1.2 is officially in GCRCatalogs
 
 creator = DiffskySedGenerator(logname=logname, galaxy_truth=galaxy_truth,
                               output_dir=args.output_dir, sky_cat=sky_cat,

From ba80f358e883aae120ea4f62a7fe12434c32e598 Mon Sep 17 00:00:00 2001
From: Joanne Bogart <jrb@slac.stanford.edu>
Date: Sat, 27 Jan 2024 18:15:18 -0800
Subject: [PATCH 6/6] change diffsky galaxy default catalog version to 1.1.2

---
 skycatalogs/catalog_creator.py            | 4 +---
 skycatalogs/scripts/create_diffsky_sed.py | 5 +----
 2 files changed, 2 insertions(+), 7 deletions(-)

diff --git a/skycatalogs/catalog_creator.py b/skycatalogs/catalog_creator.py
index b45d41a9..3cf36788 100644
--- a/skycatalogs/catalog_creator.py
+++ b/skycatalogs/catalog_creator.py
@@ -250,9 +250,7 @@ def __init__(self, parts, area_partition=None, skycatalog_root=None,
         """
 
         _cosmo_cat = 'cosmodc2_v1.1.4_image_addon_knots'
-        _diffsky_cat = 'roman_rubin_2023_v1.1.1_elais'
-        # When it's officially available, default will be newer version:
-        # _diffsky_cat = 'roman_rubin_2023_v1.1.2_elais'
+        _diffsky_cat = 'roman_rubin_2023_v1.1.2_elais'
         _star_db = '/global/cfs/cdirs/lsst/groups/SSim/DC2/dc2_stellar_healpixel.db'
         _sn_db = '/global/cfs/cdirs/lsst/groups/SSim/DC2/cosmoDC2_v1.1.4/sne_cosmoDC2_v1.1.4_MS_DDF_healpix.db'
 
diff --git a/skycatalogs/scripts/create_diffsky_sed.py b/skycatalogs/scripts/create_diffsky_sed.py
index 555d0f66..844111ce 100644
--- a/skycatalogs/scripts/create_diffsky_sed.py
+++ b/skycatalogs/scripts/create_diffsky_sed.py
@@ -84,10 +84,7 @@
 # hard-code for now.  Should be able to retrieve galaxy truth from sky_cat
 galaxy_truth = args.galaxy_truth
 if galaxy_truth is None:
-    galaxy_truth = 'roman_rubin_2023_v1.1.1_elais'
-    # Change to
-    #  galaxy_truth = 'roman_rubin_2023_v1.1.2_elais'
-    # when v1.1.2 is officially in GCRCatalogs
+    galaxy_truth = 'roman_rubin_2023_v1.1.2_elais'
 
 creator = DiffskySedGenerator(logname=logname, galaxy_truth=galaxy_truth,
                               output_dir=args.output_dir, sky_cat=sky_cat,