From 7c6b76e7e2d7285cb8e261aa1fc2f72b69912cb1 Mon Sep 17 00:00:00 2001 From: Alex Drlica-Wagner Date: Fri, 6 Sep 2024 14:38:22 -0500 Subject: [PATCH 01/15] comments for bayes factor --- ugali/analysis/results.py | 11 +++++++++++ 1 file changed, 11 insertions(+) diff --git a/ugali/analysis/results.py b/ugali/analysis/results.py index 1f4425e..082ac42 100644 --- a/ugali/analysis/results.py +++ b/ugali/analysis/results.py @@ -131,6 +131,17 @@ def estimate_position_angle(self,param='position_angle',burn=None,clip=10.0,alph return ret def bayes_factor(self,param,burn=None,clip=10.0,bins=50): + """This is a very simplistic calculation of the Bayes Factor from the + posterior samples of a single parameter at the specific + parameter value of x=0 using the Savage-Dickey method for + nested hypothesis. It is currently assumping a flat prior + since the prior is not being passed in explicitly. + + More information can be found here: + https://www.sciencedirect.com/science/article/abs/pii/S0022249611000666 + + ADW 2024-09-06: This should be moved to stats. + """ # CAREFUL: Assumes a flat prior... try: data = self.samples.get(param,burn=burn,clip=clip) From 6e2dc9ebff100eed72ec30e522a757cd1e45e8ca Mon Sep 17 00:00:00 2001 From: Alex Drlica-Wagner Date: Tue, 10 Sep 2024 10:44:35 -0500 Subject: [PATCH 02/15] update submit --- ugali/analysis/farm.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/ugali/analysis/farm.py b/ugali/analysis/farm.py index 3073b9b..eb63050 100755 --- a/ugali/analysis/farm.py +++ b/ugali/analysis/farm.py @@ -128,7 +128,7 @@ def submit(self, pixels, queue=None, debug=False, configfile=None): logger.info('=== Submit Likelihood ===') for ii,pix in enumerate(pixels): msg = ' (%i/%i) pixel=%i nside=%i; (lon, lat) = (%.2f, %.2f)' - msg = msg%(ii+1,len(pixels),pix, self.nside_likelihood,lon[ii],lat[ii]) + msg = msg%(ii,len(pixels),pix, self.nside_likelihood,lon[ii],lat[ii]) logger.info(msg) # Create outfile name From 7442ae563e920bcaa724a8f8d23860403a218a09 Mon Sep 17 00:00:00 2001 From: Alex Drlica-Wagner Date: Tue, 10 Sep 2024 10:45:43 -0500 Subject: [PATCH 03/15] typo --- ugali/analysis/search.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/ugali/analysis/search.py b/ugali/analysis/search.py index 85b12a0..f2347b2 100644 --- a/ugali/analysis/search.py +++ b/ugali/analysis/search.py @@ -382,11 +382,11 @@ def finalizeObjects(self, objects): coordsys = self.config['coords']['coordsys'] if coordsys.lower() == 'gal': - print("GAL coordintes") + print("GAL coordinates") objs['GLON'],objs['GLAT'] = lon,lat objs['RA'],objs['DEC'] = gal2cel(lon,lat) else: - print("CEL coordintes") + print("CEL coordinates") objs['RA'],objs['DEC'] = lon,lat objs['GLON'],objs['GLAT'] = cel2gal(lon,lat) From 58e6f5b939b0dfb8d27fc141ffdb4ee18a2dc3e6 Mon Sep 17 00:00:00 2001 From: Alex Drlica-Wagner Date: Tue, 10 Sep 2024 10:55:52 -0500 Subject: [PATCH 04/15] typo --- ugali/observation/roi.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/ugali/observation/roi.py b/ugali/observation/roi.py index f79c49f..61c57ce 100644 --- a/ugali/observation/roi.py +++ b/ugali/observation/roi.py @@ -152,7 +152,7 @@ def plot(self, value=None, pixel=None): # ADW: Maybe these should be associated with the PixelRegion objects def inPixels(self,lon,lat,pixels): - """ Function for testing if coordintes in set of ROI pixels. """ + """ Function for testing if coordinates in set of ROI pixels. """ nside = self.config.params['coords']['nside_pixel'] return ugali.utils.healpix.in_pixels(lon,lat,pixels,nside) From e97d291526e1a4ea16f0c21a8b5f78112d16a7eb Mon Sep 17 00:00:00 2001 From: Alex Drlica-Wagner Date: Tue, 10 Sep 2024 10:59:29 -0500 Subject: [PATCH 05/15] condor batch --- ugali/utils/batch.py | 33 ++++++++++++++++++++++++++++----- 1 file changed, 28 insertions(+), 5 deletions(-) diff --git a/ugali/utils/batch.py b/ugali/utils/batch.py index 371c240..f98550a 100644 --- a/ugali/utils/batch.py +++ b/ugali/utils/batch.py @@ -28,7 +28,7 @@ ('lsf' ,['express','short','medium','long','xlong','xxl', 'kipac-ibq','bulletmpi']), ('slurm' ,[]), - ('condor',['local','vanilla','universe','grid']), + ('condor',['local','vanilla']), ]) # https://confluence.slac.stanford.edu/x/OaUlCw @@ -80,9 +80,7 @@ def factory(queue,**kwargs): elif name in CLUSTERS['slurm']+QUEUES['slurm']: batch = Slurm(**kwargs) elif name in CLUSTERS['condor']+QUEUES['condor']: - # Need to learn how to use condor first... batch = Condor(**kwargs) - raise Exception("Condor cluster not implemented") else: raise TypeError('Unexpected queue name: %s'%name) @@ -105,9 +103,17 @@ def __init__(self, **kwargs): self.submit_cmd = "submit %(opts)s %(command)s" self.jobs_cmd = "jobs" + def parse_options(self, **opts): + """ Parse command line options. """ + + # Default options for the cluster + options = odict(self.default_opts) + options.update(opts) + return ''.join('--%s %s '%(k,v) for k,v in options.items()) + def jobs(self): out = self.popen(self.jobs_cmd) - stdout = out.communicate()[0] + stdout = out.communicate()[0].decode() return stdout def njobs(self): @@ -312,13 +318,30 @@ def parse_options(self, **opts): class Condor(Batch): """ Not implemented yet... """ + + _defaults = odict([ + ('n', '50'), + ]) + + _mapping = odict([ + ('jobname','J'), + ('logfile','o'), + ('njobs','n'), + ]) + def __init__(self, **kwargs): super(Condor,self).__init__(**kwargs) logger.warning('Condor cluster is untested') - self.jobs_cmd = 'condor_q -u %s'%self.username + self.jobs_cmd = "cjobs -u %s"%(self.username) self.submit_cmd = "csub %(opts)s %(command)s" + def parse_options(self, **opts): + # Default options for the cluster + options = odict(self.default_opts) + options.update(opts) + return ''.join('-%s %s '%(k,v) for k,v in options.items()) + if __name__ == "__main__": import argparse description = __doc__ From 9ab1a6a8518fa0415cbb077a087d4e045967b82c Mon Sep 17 00:00:00 2001 From: Alex Drlica-Wagner Date: Tue, 10 Sep 2024 11:00:29 -0500 Subject: [PATCH 06/15] output keys --- ugali/utils/fileio.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/ugali/utils/fileio.py b/ugali/utils/fileio.py index b62bd6b..c8b848a 100644 --- a/ugali/utils/fileio.py +++ b/ugali/utils/fileio.py @@ -141,7 +141,7 @@ def load_headers(filenames,multiproc=False,**kwargs): out = [load_header(kw) for kw in kwargs] logger.debug('Concatenating headers...') - return np.rec.fromrecords([o.values() for o in out],names=out[0].keys()) + return np.rec.fromrecords([list(o.values()) for o in out], names=list(out[0].keys())) def load_file(kwargs): From 2ab090ed62004f452378235e249e47e9767876d4 Mon Sep 17 00:00:00 2001 From: Alex Drlica-Wagner Date: Tue, 10 Sep 2024 11:01:41 -0500 Subject: [PATCH 07/15] asscalar --- ugali/utils/mlab.py | 5 +++++ 1 file changed, 5 insertions(+) diff --git a/ugali/utils/mlab.py b/ugali/utils/mlab.py index ddf581c..eeb1407 100644 --- a/ugali/utils/mlab.py +++ b/ugali/utils/mlab.py @@ -17,6 +17,11 @@ def isstring(obj): """Python 2/3 compatible string check""" return isinstance(obj, six.string_types) +def asscalar(a): + """Python 2/3 compatible scalar conversion""" + try: return np.asscalar(a) + except AttributeError: return a.item() + def rec_append_fields(rec, names, arrs, dtypes=None): """ Return a new record array with field names populated with data From 7bfc9ce546f3936e40c25696d2b724660eaa83e2 Mon Sep 17 00:00:00 2001 From: Alex Drlica-Wagner Date: Tue, 10 Sep 2024 11:02:55 -0500 Subject: [PATCH 08/15] asscalar --- ugali/utils/projector.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/ugali/utils/projector.py b/ugali/utils/projector.py index 9de825a..fbdb5a2 100644 --- a/ugali/utils/projector.py +++ b/ugali/utils/projector.py @@ -8,7 +8,7 @@ import numpy as np from ugali.utils.logger import logger -from ugali.utils.mlab import isstring +from ugali.utils.mlab import isstring, asscalar ############################################################ @@ -117,7 +117,7 @@ def sphereToImage(self, lon, lat): lon_rotated, lat_rotated = self.rotator.rotate(lon.flat, lat.flat) x, y = self.sphere_to_image_func(lon_rotated, lat_rotated) - if scalar: return np.asscalar(x), np.asscalar(y) + if scalar: return asscalar(x), asscalar(y) else: return x.reshape(lon.shape), y.reshape(lat.shape) sphere2image = sphereToImage @@ -129,7 +129,7 @@ def imageToSphere(self, x, y): lon_rotated, lat_rotated = self.image_to_sphere_func(x.flat, y.flat) lon, lat = self.rotator.rotate(lon_rotated, lat_rotated, invert = True) - if scalar: return np.asscalar(lon), np.asscalar(lat) + if scalar: return asscalar(lon), asscalar(lat) else: return lon.reshape(x.shape), lat.reshape(y.shape) image2sphere = imageToSphere From e39fe58e78489ef5c16ef99944ef8349c03033d1 Mon Sep 17 00:00:00 2001 From: Alex Drlica-Wagner Date: Tue, 10 Sep 2024 11:03:53 -0500 Subject: [PATCH 09/15] hb_spread iterable and mass interpolation --- ugali/isochrone/model.py | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/ugali/isochrone/model.py b/ugali/isochrone/model.py index d673934..fa73f00 100644 --- a/ugali/isochrone/model.py +++ b/ugali/isochrone/model.py @@ -219,11 +219,11 @@ def sample(self, mode='data', mass_steps=1000, mass_min=0.1, full_data_range=Fal mass_init_min = self.mass_init[self.stage==self.hb_stage].min() mass_init_max = self.mass_init[self.stage==self.hb_stage].max() cut = (mass_init_array>mass_init_min)&(mass_init_array Date: Tue, 10 Sep 2024 11:06:17 -0500 Subject: [PATCH 10/15] merged movie --- ugali/pipeline/run_03.0_likelihood.py | 33 ++++++++++++++++++++------- 1 file changed, 25 insertions(+), 8 deletions(-) diff --git a/ugali/pipeline/run_03.0_likelihood.py b/ugali/pipeline/run_03.0_likelihood.py index 895963e..9947c25 100755 --- a/ugali/pipeline/run_03.0_likelihood.py +++ b/ugali/pipeline/run_03.0_likelihood.py @@ -24,6 +24,8 @@ def run(self): farm.submit_all(coords=self.opts.coords,queue=self.opts.queue,debug=self.opts.debug) if 'merge' in self.opts.run: + import numpy as np + logger.info("Running 'merge'...") mergefile = self.config.mergefile roifile = self.config.roifile @@ -42,13 +44,13 @@ def run(self): for pix in np.unique(superpixel): outfile = mergefile%pix if exists(outfile) and not self.opts.force: - logger.warn(" Found %s; skipping..."%outfile) + logger.warning(" Found %s; skipping..."%outfile) else: healpix.merge_partial_maps(infiles[superpixel == pix], outfile,multiproc=8) if exists(roifile) and not self.opts.force: - logger.warn(" Found %s; skipping..."%roifile) + logger.warning(" Found %s; skipping..."%roifile) else: logger.info(" Merging likelihood headers...") healpix.merge_likelihood_headers(infiles,roifile) @@ -64,7 +66,7 @@ def run(self): logfile = os.path.join(logdir,'scan_tar.log') cmd = 'tar --remove-files -cvzf %s %s'%(tarfile,scanfile) if exists(tarfile) and not self.opts.force: - logger.warn(" Found %s; skipping..."%tarfile) + logger.warning(" Found %s; skipping..."%tarfile) else: logger.info(" Tarring likelihood files...") logger.info(cmd) @@ -73,17 +75,32 @@ def run(self): if 'plot' in self.opts.run: # WARNING: Loading the full 3D healpix map is memory intensive. logger.info("Running 'plot'...") + import numpy as np # Should do this in environment variable import matplotlib matplotlib.use('Agg') import pylab as plt import ugali.utils.plotting as plotting - skymap = ugali.utils.skymap.readSparseHealpixMap(self.config.mergefile,'LOG_LIKELIHOOD')[1] - plotting.plotSkymap(skymap) + from skymap import DESSkymap + + filenames = self.config.mergefile.split('_%')[0]+'_*.fits' + infiles = np.array(sorted(glob.glob(filenames))) + hpxmaps = ugali.utils.healpix.read_partial_map(infiles,'LOG_LIKELIHOOD',fullsky=True)[2] + outdir = mkdir(self.config['output']['plotdir']) - basename = os.path.basename(self.config.mergefile.replace('.fits','.png')) - outfile = os.path.join(outdir,basename) - plt.savefig(outfile) + basename = os.path.basename(self.config.mergefile.split('_%')[0]) + for i,hpxmap in enumerate(hpxmaps): + #plotting.plotSkymap(hpxmap) + smap = DESSkymap() + smap.draw_hpxmap(hpxmap, cmap='gray_r', vmax=3.5) + smap.draw_inset_colorbar() + outfile = os.path.join(outdir,basename+'_%02d.png'%i) + print("Writing %s..."%outfile) + plt.savefig(outfile) + # Make the movie + cmd = "convert -delay 30 -loop 0 %s/*.png %s.gif"%(outdir,basename) + print(cmd) + subprocess.check_call(cmd,shell=True) if 'check' in self.opts.run: # Check the completion fraction From 6bb1c66d0751d155b7a46db798b980c0dfd55c72 Mon Sep 17 00:00:00 2001 From: Alex Drlica-Wagner Date: Tue, 10 Sep 2024 11:06:57 -0500 Subject: [PATCH 11/15] white space --- ugali/utils/healpix.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/ugali/utils/healpix.py b/ugali/utils/healpix.py index b4696ce..9f66c74 100644 --- a/ugali/utils/healpix.py +++ b/ugali/utils/healpix.py @@ -520,7 +520,7 @@ def merge_likelihood_headers(filenames, outfile, **kwargs): data_dict = fileio.load_headers(filenames,ext=ext,keys=keys,mulitproc=8) names = data_dict.dtype.names data_dict.dtype.names = ['PIXEL' if n=='LKDPIX' else n for n in names] - + write_partial_map(outfile, data_dict, nside) return data_dict From d26ad717333c0cad13d0a8ad1771b0ec2bb5da4d Mon Sep 17 00:00:00 2001 From: Alex Drlica-Wagner Date: Tue, 10 Sep 2024 11:11:07 -0500 Subject: [PATCH 12/15] move to healpix module --- ugali/utils/plotting.py | 5 ++--- ugali/utils/skymap.py | 3 ++- 2 files changed, 4 insertions(+), 4 deletions(-) diff --git a/ugali/utils/plotting.py b/ugali/utils/plotting.py index c289bf3..e180cdb 100644 --- a/ugali/utils/plotting.py +++ b/ugali/utils/plotting.py @@ -24,7 +24,6 @@ import ugali.utils.config import ugali.observation.roi import ugali.observation.catalog -import ugali.utils.skymap import ugali.utils.projector import ugali.utils.healpix import ugali.isochrone @@ -153,7 +152,7 @@ def sparseHealpixFiles(title, infiles, field='MAGLIM',**kwargs): Inputs: field """ #map = ugali.utils.skymap.readSparseHealpixMaps(infiles,field) - map = ugali.utils.skymap.read_partial_map(infiles,field) + nside,pixels,hpxmap = ugali.utils.healpix.read_partial_map(infiles,field,fullsky=True) ax = hp.mollview(map=map, title=title, **kwargs) return ax, map @@ -602,7 +601,7 @@ def drawMembership(self, ax=None, radius=None, zidx=0, mc_source_id=1): pix = ang2pix(self.nside, self.lon, self.lat) - likelihood_pix = ugali.utils.skymap.superpixel(pix,self.nside,self.config.params['coords']['nside_likelihood']) + likelihood_pix = ugali.utils.healpix.superpixel(pix,self.nside,self.config.params['coords']['nside_likelihood']) config = self.config scan = ugali.analysis.scan.Scan(self.config,likelihood_pix) likelihood = scan.likelihood diff --git a/ugali/utils/skymap.py b/ugali/utils/skymap.py index cca24e8..f51fd8d 100644 --- a/ugali/utils/skymap.py +++ b/ugali/utils/skymap.py @@ -6,7 +6,8 @@ import healpy as hp import ugali.utils.projector -from ugali.utils.healpix import superpixel,subpixel,ang2pix,pix2ang,query_disc +from ugali.utils.healpix import superpixel,subpixel +from ugali.utils.healpix import ang2pix,pix2ang,query_disc from ugali.utils.healpix import read_partial_map from ugali.utils.logger import logger from ugali.utils.config import Config From e51f4371956a649b8f57226c2245ff1df4a31acd Mon Sep 17 00:00:00 2001 From: Alex Drlica-Wagner Date: Tue, 10 Sep 2024 14:42:18 -0500 Subject: [PATCH 13/15] actions --- .github/workflows/python-package.yml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/.github/workflows/python-package.yml b/.github/workflows/python-package.yml index b671fef..d20dd9e 100644 --- a/.github/workflows/python-package.yml +++ b/.github/workflows/python-package.yml @@ -18,7 +18,7 @@ jobs: strategy: max-parallel: 3 matrix: - python-version: ["3.7", "3.8", "3.9"] + python-version: ["3.8", "3.9", "3.10"] steps: - uses: actions/checkout@v4 - name: Set up Python ${{ matrix.python-version }} From 5be1740064266994652c897270ccfa7e8a2932bf Mon Sep 17 00:00:00 2001 From: Alex Drlica-Wagner Date: Tue, 10 Sep 2024 14:42:38 -0500 Subject: [PATCH 14/15] import subprocess --- ugali/pipeline/run_03.0_likelihood.py | 1 + 1 file changed, 1 insertion(+) diff --git a/ugali/pipeline/run_03.0_likelihood.py b/ugali/pipeline/run_03.0_likelihood.py index 9947c25..5ae0bd5 100755 --- a/ugali/pipeline/run_03.0_likelihood.py +++ b/ugali/pipeline/run_03.0_likelihood.py @@ -98,6 +98,7 @@ def run(self): print("Writing %s..."%outfile) plt.savefig(outfile) # Make the movie + import subprocess cmd = "convert -delay 30 -loop 0 %s/*.png %s.gif"%(outdir,basename) print(cmd) subprocess.check_call(cmd,shell=True) From 38756d6ba74d906cdfd1e9c12df777720d4458e4 Mon Sep 17 00:00:00 2001 From: Alex Drlica-Wagner Date: Tue, 10 Sep 2024 14:48:33 -0500 Subject: [PATCH 15/15] actions --- .github/workflows/python-package.yml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/.github/workflows/python-package.yml b/.github/workflows/python-package.yml index d20dd9e..ab93a18 100644 --- a/.github/workflows/python-package.yml +++ b/.github/workflows/python-package.yml @@ -18,7 +18,7 @@ jobs: strategy: max-parallel: 3 matrix: - python-version: ["3.8", "3.9", "3.10"] + python-version: ["3.8", "3.9", "3.10", "3.11"] steps: - uses: actions/checkout@v4 - name: Set up Python ${{ matrix.python-version }}