From f2c88f3d303c2784cbcb475d2676b3fb967d1ddd Mon Sep 17 00:00:00 2001 From: Alex Drlica-Wagner Date: Tue, 25 Jun 2024 17:42:40 -0500 Subject: [PATCH 01/19] documentation --- ugali/analysis/imf.py | 2 +- ugali/isochrone/dartmouth.py | 4 ++-- ugali/utils/stats.py | 19 +++++++++++++++++++ 3 files changed, 22 insertions(+), 3 deletions(-) diff --git a/ugali/analysis/imf.py b/ugali/analysis/imf.py index 07b4304..ba5491b 100644 --- a/ugali/analysis/imf.py +++ b/ugali/analysis/imf.py @@ -27,7 +27,7 @@ def integrate(self, mass_min, mass_max, log_mode=True, weight=False, steps=10000 ----------- mass_min: minimum mass bound for integration (solar masses) mass_max: maximum mass bound for integration (solar masses) - log_mode[True]: use logarithmic steps in stellar mass as oppose to linear + log_mode[True]: use logarithmic steps in stellar mass as opposed to linear weight[False]: weight the integral by stellar mass steps: number of numerical integration steps diff --git a/ugali/isochrone/dartmouth.py b/ugali/isochrone/dartmouth.py index 3e0a2d1..1e8fc36 100644 --- a/ugali/isochrone/dartmouth.py +++ b/ugali/isochrone/dartmouth.py @@ -103,8 +103,8 @@ class Dotter2008(Isochrone): """ - KCB: currently inheriting from PadovaIsochrone because there are - several useful functions where we would basically be copying code. + Dartmouth isochrones from Dotter et al. 2008 (https://arxiv.org/abs/0804.4473) + http://stellar.dartmouth.edu """ _dirname = os.path.join(get_iso_dir(),'{survey}','dotter2008') #_zsolar = 0.0163 diff --git a/ugali/utils/stats.py b/ugali/utils/stats.py index 7e41e90..8ec502f 100644 --- a/ugali/utils/stats.py +++ b/ugali/utils/stats.py @@ -28,6 +28,16 @@ def mad_clip(data,mad=None,mad_lower=None,mad_upper=None): def interval(best,lo=np.nan,hi=np.nan): """ Pythonized interval for easy output to yaml + + Parameters + ---------- + best : best-fit estimate of the parameter + lo : lower value + hi : higher value + + Returns + ------- + [best, [lo, hi]] : list of values """ return [float(best),[float(lo),float(hi)]] @@ -43,6 +53,15 @@ def mean_interval(data, alpha=_alpha): def median_interval(data, alpha=_alpha): """ Median with bayesian credible interval from percentiles. + + Parameters + ---------- + data : posterior samples + alpha : 1 - confidence interval + + Returns + ------- + [med,[lo, hi]] : median, lower, and upper percentiles """ q = [100*alpha/2., 50, 100*(1-alpha/2.)] lo,med,hi = np.percentile(data,q) From 6c96f650284988b844e71eefacb9299298c6c52b Mon Sep 17 00:00:00 2001 From: Alex Drlica-Wagner Date: Tue, 25 Jun 2024 17:50:56 -0500 Subject: [PATCH 02/19] updates to doc and download --- ugali/isochrone/model.py | 49 ++++++++++++++++++++++++++++------------ 1 file changed, 34 insertions(+), 15 deletions(-) diff --git a/ugali/isochrone/model.py b/ugali/isochrone/model.py index 8f592dc..2ac1af9 100644 --- a/ugali/isochrone/model.py +++ b/ugali/isochrone/model.py @@ -96,9 +96,9 @@ class IsochroneModel(Model): """ Abstract base class for dealing with isochrone models. """ _params = odict([ - ('distance_modulus', Parameter(15.0, [10.0, 30.0]) ), + ('distance_modulus', Parameter(15.0, [10.0, 30.0]) ), # m-M ('age', Parameter(10.0, [0.1, 15.0]) ), # Gyr - ('metallicity', Parameter(0.0002, [0.0,0.02]) ), + ('metallicity', Parameter(0.0002, [0.0,0.02]) ), # Z_ini ]) _mapping = odict([ ('mod','distance_modulus'), @@ -407,7 +407,7 @@ def absolute_magnitude_martin(self, richness=1, steps=1e4, n_trials=1000, mag_br Returns: -------- - med,lo,hi : Total absolute magnitude interval + med,[lo,hi] : Total absolute magnitude interval """ # ADW: This function is not quite right. It should restrict # the catalog to the obsevable space using the mask in each @@ -420,6 +420,10 @@ def absolute_magnitude_martin(self, richness=1, steps=1e4, n_trials=1000, mag_br params.update(band_1='g',band_2='r',survey='sdss') iso = self.__class__(**params) + # ADW: 2022-09-18: There are now direct transformation + # equations from DES to Johnson V-band... + # Appendix B.5: https://arxiv.org/abs/2101.05765 + # Analytic part (below detection threshold) # g, r are absolute magnitudes mass_init, mass_pdf, mass_act, sdss_g, sdss_r = iso.sample(mass_steps = steps) @@ -462,10 +466,15 @@ def simulate(self, stellar_mass, distance_modulus=None, **kwargs): """ if distance_modulus is None: distance_modulus = self.distance_modulus # Total number of stars in system - n = int(round(stellar_mass / self.stellar_mass())) - f_1 = scipy.interpolate.interp1d(self.mass_init, self.mag_1) - f_2 = scipy.interpolate.interp1d(self.mass_init, self.mag_2) - mass_init_sample = self.imf.sample(n, np.min(self.mass_init), np.max(self.mass_init), **kwargs) + richness = stellar_mass / self.stellar_mass() + n = int(round(richness)) + #ADW 2022-08-28: Should this really be a Poisson sample of the richness? + #n = scipy.stats.poisson(richness).rvs() + + mass_init,mag_1,mag_2 = self.mass_init, self.mag_1, self.mag_2 + f_1 = scipy.interpolate.interp1d(mass_init, mag_1) + f_2 = scipy.interpolate.interp1d(mass_init, mag_2) + mass_init_sample = self.imf.sample(n, np.min(mass_init), np.max(mass_init), **kwargs) mag_1_sample, mag_2_sample = f_1(mass_init_sample), f_2(mass_init_sample) return mag_1_sample + distance_modulus, mag_2_sample + distance_modulus @@ -929,8 +938,11 @@ def pdf(self, mag_1, mag_2, mag_err_1, mag_err_2, def raw_separation(self,mag_1,mag_2,steps=10000): - """ - Calculate the separation in magnitude-magnitude space between points and isochrone. Uses a dense sampling of the isochrone and calculates the metric distance from any isochrone sample point. + """ + Calculate the separation in magnitude-magnitude space between + points and isochrone. Uses a dense sampling of the isochrone + and calculates the metric distance from any isochrone sample + point. Parameters: ----------- @@ -1096,11 +1108,14 @@ def create_grid(self,abins=None,zbins=None): if abins is None and zbins is None: filenames = glob.glob(self.get_dirname()+'/%s_*.dat'%(self._prefix)) data = np.array([self.filename2params(f) for f in filenames]) - if not len(data): + if len(data): + arange = np.unique(data[:,0]) + zrange = np.unique(data[:,1]) + else: msg = "No isochrone files found in: %s"%self.get_dirname() - raise Exception(msg) - arange = np.unique(data[:,0]) - zrange = np.unique(data[:,1]) + logger.warn(msg) + arange = np.array([self.age]) + zrange = np.array([self.metallicity]) elif abins is not None and zbins is not None: # Age in units of Gyr arange = np.linspace(abins[0],abins[1],abins[2]+1) @@ -1132,7 +1147,11 @@ def _cache(self,name=None): filename = self.get_filename() if filename != self.filename: self.filename = filename - self._parse(self.filename) + if not os.path.exists(self.filename): + msg = "Filename does not exist: %s"%self.filename + logger.warn(msg) + else: + self._parse(self.filename) def _parse(self,filename): raise Exception("Must be implemented by subclass.") @@ -1165,7 +1184,7 @@ def download(self,age=None,metallicity=None,outdir=None,force=False): Parameters ---------- age : age in (Gyr) - metallicity : Z + metallicity : initial metallicity, Z outdir : output directory (default to current directory) force : force overwrite of file From 15a7a938efc13d7e33834628b32055ce11798c2a Mon Sep 17 00:00:00 2001 From: Alex Drlica-Wagner Date: Tue, 25 Jun 2024 17:53:37 -0500 Subject: [PATCH 03/19] updates to parsec isochrone motivated by LSST DP0 --- ugali/isochrone/parsec.py | 156 ++++++++++++++++++++++++++------------ 1 file changed, 109 insertions(+), 47 deletions(-) diff --git a/ugali/isochrone/parsec.py b/ugali/isochrone/parsec.py index 5157fec..7890ff4 100644 --- a/ugali/isochrone/parsec.py +++ b/ugali/isochrone/parsec.py @@ -34,6 +34,7 @@ ('ps1' ,'tab_mag_odfnew/tab_mag_panstarrs1.dat'), ('acs_wfc' ,'tab_mag_odfnew/tab_mag_acs_wfc.dat'), ('lsst', 'tab_mag_odfnew/tab_mag_lsst.dat'), + ('lsst_dp0', 'tab_mag_odfnew/tab_mag_lsstDP0.dat'), ]) photname_dict = odict([ @@ -42,6 +43,7 @@ ('ps1' ,'Pan-STARRS1'), ('acs_wfc','HST/ACS'), ('lsst', 'LSST'), + ('lsst_dp0', 'LSST'), ]) # Commented options may need to be restored for older version/isochrones. @@ -155,6 +157,8 @@ #'.cgifields': 'dust_sourceM', } +defaults_36 = dict(defaults_33,cmd_version=3.6) + class ParsecIsochrone(Isochrone): """ Base class for PARSEC-style isochrones. """ @@ -242,7 +246,7 @@ def query_server(self,outfile,age,metallicity): urlopen(url,timeout=2) q = urlencode(params).encode('utf-8') - logger.debug(url+'?'+q) + logger.debug("%s?%s"%(url,q)) c = str(urlopen(url, q).read()) aa = re.compile('output\d+') fname = aa.findall(c) @@ -256,49 +260,93 @@ def query_server(self,outfile,age,metallicity): cmd = 'wget --progress dot:binary %s -O %s'%(out,outfile) logger.debug(cmd) stdout = subprocess.check_output(cmd,shell=True,stderr=subprocess.STDOUT) - logger.debug(stdout) + logger.debug(str(stdout)) return outfile @classmethod - def verify(cls, filename, survey, age, metallicity): - age = age*1e9 - nlines=15 + def parse_header(cls, filename, nlines=15): + header = dict( + photname = None, + columns = None + ) + with open(filename,'r') as f: lines = [f.readline() for i in range(nlines)] - if len(lines) < nlines: - msg = "Incorrect file size" - raise Exception(msg) - for i,l in enumerate(lines): - if l.startswith('# Photometric system:'): break - else: - msg = "Incorrect file header" - raise Exception(msg) + if len(lines) < nlines: + msg = "Incorrect file size" + raise Exception(msg) - try: - s = lines[i].split()[3] - assert photname_dict[survey] == s - except: - msg = "Incorrect survey:\n"+lines[i] - raise Exception(msg) + for i,l in enumerate(lines): + if l.startswith('# Photometric system:'): + try: header['photname'] = lines[i].split()[3] + except: header['photname'] = None + if not l.startswith('# '): break - try: - z = lines[-1].split()[0] - assert np.allclose(metallicity,float(z),atol=1e-5) - except: - msg = "Metallicity does not match:\n"+lines[-1] - raise Exception(msg) + header['columns'] = lines[i-1].split()[1:] - try: - a = lines[-1].split()[1] - # Need to deal with age or log-age - assert (np.allclose(age,float(a),atol=1e-2) or - np.allclose(np.log10(age),float(a),atol=1e-2)) - except: - msg = "Age does not match:\n"+lines[-1] + for k,v in header.items(): + if v is None: + msg = "File header missing: '%s'"%k raise Exception(msg) + return header, lines + + @classmethod + def verify(cls, filename, survey, age, metallicity): + """Verify that the isochrone file matches the isochrone + parameters. Used mostly for verifying the integrity of + a download. + + Parameters + ---------- + filename : the downloaded filename + survey : the survey (photometric system) of the download + age : the requested age of the system + metallicity : the requested metallicity + + Returns + ------- + None + """ + age = age*1e9 + nlines=15 + header, lines = cls.parse_header(filename,nlines=nlines) + + try: + assert photname_dict[survey] == header['photname'] + except: + msg = "Incorrect survey:\n"+header['photname'] + raise Exception(msg) + + try: + try: zidx = header['columns'].index('Zini') + except ValueError: zidx = 0 + z = float(lines[-1].split()[zidx]) + assert np.allclose(metallicity,z,atol=1e-5) + except: + msg = "Metallicity does not match:\n"+lines[-1] + raise Exception(msg) + + try: + # Need to deal with age or log-age + names = ['log(age/yr)', 'logAge', 'Age'] + for name in names: + try: + aidx = header['columns'].index(name) + break + except ValueError: + aidx = 1 + if aidx < 0: aidx = 1 + + a = lines[-1].split()[aidx] + assert (np.allclose(age,float(a),atol=1e-2) or + np.allclose(np.log10(age),float(a),atol=1e-2)) + except: + msg = "Age does not match:\n"+lines[-1] + raise Exception(msg) + class Bressan2012(ParsecIsochrone): _dirname = os.path.join(get_iso_dir(),'{survey}','bressan2012') @@ -403,7 +451,7 @@ class Marigo2017(ParsecIsochrone): ('hb_spread',0.1,'Intrinisic spread added to horizontal branch'), ) - download_defaults = copy.deepcopy(defaults_31) + download_defaults = copy.deepcopy(defaults_36) download_defaults['isoc_kind'] = 'parsec_CAF09_v1.2S_NOV13' columns = dict( @@ -442,24 +490,38 @@ class Marigo2017(ParsecIsochrone): (27,('y',float)), (28,('w',float)), ]), - lsst = odict([ - (2, ('mass_init',float)), - (3, ('mass_act',float)), - (4, ('log_lum',float)), - (7, ('stage',int)), - (23,('u',float)), - (24,('g',float)), - (25,('r',float)), - (26,('i',float)), - (27,('z',float)), - (28,('Y',float)), + lsst_dp0 = odict([ + (3, ('mass_init',float)), + (5, ('mass_act',float)), + (6, ('log_lum',float)), + (9, ('stage',int)), + (25,('u',float)), + (26,('g',float)), + (27,('r',float)), + (28,('i',float)), + (29,('z',float)), + (30,('Y',float)), ]), ) - + columns['lsst'] = copy.deepcopy(columns['lsst_dp0']) + + def _find_column_numbers(self): + """ Map from the isochrone column names to the column numbers. """ + header,lines = self.parse_header(self.filename) + columns = header['columns'] + def _parse(self,filename): - """Reads an isochrone file in the Padova (Marigo et al. 2017) + """Reads an isochrone file in the Marigo et al. 2017 format. Creates arrays with the initial stellar mass and corresponding magnitudes for each step along the isochrone. + + Parameters: + ----------- + filename : name of isochrone file to parse + + Returns: + -------- + None """ try: columns = self.columns[self.survey.lower()] @@ -471,7 +533,7 @@ def _parse(self,filename): self.data = np.genfromtxt(filename,**kwargs) # cut out anomalous point: # https://github.com/DarkEnergySurvey/ugali/issues/29 - self.data = self.data[self.data['stage'] != 9] + self.data = self.data[~np.in1d(self.data['stage'], [9])] self.mass_init = self.data['mass_init'] self.mass_act = self.data['mass_act'] From 7931786df204e0ddfab2d8c8365c367879bfd610 Mon Sep 17 00:00:00 2001 From: Alex Drlica-Wagner Date: Tue, 25 Jun 2024 18:02:55 -0500 Subject: [PATCH 04/19] deprecated old preprocessing download --- ugali/preprocess/dotter.py | 3 +++ ugali/preprocess/padova.py | 5 +++++ 2 files changed, 8 insertions(+) diff --git a/ugali/preprocess/dotter.py b/ugali/preprocess/dotter.py index edb8829..99afefa 100755 --- a/ugali/preprocess/dotter.py +++ b/ugali/preprocess/dotter.py @@ -1,5 +1,8 @@ #!/usr/bin/env python """ +DEPRECATED: ADW 2024-06-25 +New download script in ugali/scratch/download_isochrones.py + Script to automate the process of generating an isochrone library using the Dotter isochrones. diff --git a/ugali/preprocess/padova.py b/ugali/preprocess/padova.py index 3d4ad03..823c6af 100755 --- a/ugali/preprocess/padova.py +++ b/ugali/preprocess/padova.py @@ -1,5 +1,8 @@ #!/usr/bin/env python """ +DEPRECATED: ADW 2024-06-25 +New download script in ugali/scratch/download_isochrones.py + Download Padova isochrones from: http://stev.oapd.inaf.it/cgi-bin/cmd @@ -351,9 +354,11 @@ def factory(name, **kwargs): def run(args): try: p.download(*args) + return True except Exception as e: logger.warn(str(e)) logger.error("Download failed.") + return False arguments = [(a,z,args.outdir,args.force) for a,z in zip(*grid)] if args.njobs > 1: From c11b326c65cef51f4a0d1bb24538ddd4b48f38d1 Mon Sep 17 00:00:00 2001 From: Alex Drlica-Wagner Date: Tue, 25 Jun 2024 18:03:27 -0500 Subject: [PATCH 05/19] removing old preprocessing download --- ugali/preprocess/dotter.py | 355 ----------------------------------- ugali/preprocess/padova.py | 369 ------------------------------------- 2 files changed, 724 deletions(-) delete mode 100755 ugali/preprocess/dotter.py delete mode 100755 ugali/preprocess/padova.py diff --git a/ugali/preprocess/dotter.py b/ugali/preprocess/dotter.py deleted file mode 100755 index 99afefa..0000000 --- a/ugali/preprocess/dotter.py +++ /dev/null @@ -1,355 +0,0 @@ -#!/usr/bin/env python -""" -DEPRECATED: ADW 2024-06-25 -New download script in ugali/scratch/download_isochrones.py - -Script to automate the process of generating an isochrone library -using the Dotter isochrones. - -Download isochrones from: -http://stellar.dartmouth.edu/models/isolf_new.html - -New MESA isochrones come from: -http://waps.cfa.harvard.edu/MIST/interp_isos.html -""" - -import os -import re -try: - from urllib.parse import urlencode - from urllib.request import urlopen -except ImportError: - from urllib import urlencode - from urllib2 import urlopen - -import requests -import sys -import copy -import tempfile -import subprocess -import shutil -from multiprocessing import Pool -from collections import OrderedDict as odict - -import numpy -import numpy as np - -import ugali.utils.logger -import ugali.utils.shell -from ugali.utils.logger import logger -from ugali.utils.shell import mkdir -import ugali.analysis.isochrone as iso -from ugali.preprocess.padova import Download - -""" -See Vargas et al. 2013 for the distribution of alpha elements in -dSphs: http://adsabs.harvard.edu/abs/2013ApJ...767..134V - -Josh Simon remarks: For stars at [Fe/H] > -2, [a/Fe] tends to be -around zero. [Note, though, that this paper does not attempt to do -any membership classification, it just accepts the lists from -Simon & Geha 2007. I suspect now that we were not sufficiently -conservative on selecting members in those early days, and so some -of the relatively metal-rich stars may in fact be foreground Milky -Way stars.] More metal-poor stars tend to average more like -[a/Fe] = 0.4-0.5. Fig. 5 of Frebel et al. (2014) shows similar -plots for individual elements from high-resolution spectra. Given -these data, plus the empirical fact that the mean metallicities of -the ultra-faint dwarfs are almost universally [Fe/H] < -2, I guess -I would say [a/Fe] = 0.3 is probably the best compromise. - -From ADW: Other isochrone sets impose [a/Fe] = 0. For an accurate -comparison between isochrones, I suggest that we stick to [a/Fe] = 0 -for the Dotter2008 isochrones as well. - -ADW: The Dotter 2008 isochrones are interpolated from a relative -sparse grid of metallicities [-2.5, -2.0, -1.5, -1.0, -0.5, 0.0]. This -can lead to rather large differences between the input and output Z -> -[Fe/H] conversions. We were initially taking the output Zeff value, -but we now use the input Z value for internal consistency. -""" - - -########################################################### -# Dartmouth Isochrones -# http://stellar.dartmouth.edu/models/isolf_new.php - -dict_clr = {'des' : 14, - 'sdss': 11, - 'ps1' : 12, - } - -dict_hel = {'Y=0.245+1.5*Z' : 1, - 'Y=0.33' : 2, - 'Y=0.40' : 3} - -dict_afe = {'-0.2' : 1, - '0 (scaled-solar)': 2, - '+0.2' : 3, - '+0.4' : 4, - '+0.6' : 5, - '+0.8' : 6} - -dartmouth_defaults = { - 'int':'1', # interpolation: cubic=1, linear=2 (ADW: cubic more accurate) - 'out':'1', # outpur: iso=1, iso+LF=2 - 'age':10, # age [Gyr] - 'feh':-2.0, # feh [-2.5 to 0.0] - 'hel': dict_hel['Y=0.245+1.5*Z'], # initial helium abundance - 'afe': dict_afe['0 (scaled-solar)'], # alpha enhancement - 'clr': dict_clr['des'], # photometric system - 'flt':'', - 'bin':'', - 'imf':1, - 'pls':'', - 'lnm':'', - 'lns':'', - } - -########################################################### -# MESA Isochrones -# http://waps.cfa.harvard.edu/MIST/iso_form.php - -# survey system -dict_output = odict([ - ('des' ,'DECam'), - ('sdss','SDSSugriz'), - ('ps1' ,'PanSTARRS'), - ('lsst','LSST'), -]) - -mesa_defaults = { - 'version':'1.0', - 'v_div_vcrit':'vvcrit0.4', - 'age_scale':'linear', - 'age_type':'single', - 'age_value':10e9, # yr if scale='linear'; log10(yr) if scale='log10' - 'age_range_low':'', - 'age_range_high':'', - 'age_range_delta':'', - 'age_list':'', - 'FeH_value':-3.0, - 'theory_output':'basic', - 'output_option':'photometry', - 'output':'DECam', - 'Av_value':0, -} - -mesa_defaults_10 = dict(mesa_defaults,version='1.0') - -class Dotter2008(Download): - """ Dartmouth isochrones from Dotter et al. 2008: - http://stellar.dartmouth.edu/models/isolf_new.html - """ - defaults = copy.deepcopy(dartmouth_defaults) - isochrone=iso.Dotter2008 - - abins = np.arange(1., 13.5 + 0.1, 0.1) - zbins = np.arange(7e-5,1e-3 + 1e-5,1e-5) - - def query_server(self, outfile, age, metallicity): - z = metallicity - feh = self.isochrone.z2feh(z) - - params = dict(self.defaults) - params['age']=age - params['feh']='%.6f'%feh - params['clr']=dict_clr[self.survey] - - url = 'http://stellar.dartmouth.edu/models/isolf_new.php' - query = url + '?' + urlencode(params) - logger.debug(query) - response = urlopen(query) - page_source = response.read() - try: - file_id = int(page_source.split('tmp/tmp')[-1].split('.iso')[0]) - except Exception as e: - logger.debug(str(e)) - msg = 'Output filename not found' - raise RuntimeError(msg) - - infile = 'http://stellar.dartmouth.edu/models/tmp/tmp%s.iso'%(file_id) - command = 'wget -q %s -O %s'%(infile, outfile) - subprocess.call(command,shell=True) - - ## ADW: Old code to rename the output file based on Zeff ([a/Fe] corrected) - #tmpfile = tempfile.NamedTemporaryFile().name - #tmp = open(tmpfile,'r') - #lines = [tmp.readline() for i in range(4)] - #z_eff = float(lines[3].split()[4]) - #basename = self.params2filename(age,z_eff) - - #logger.info("Writing %s..."%outfile) - #mkdir(outdir) - #shutil.move(tmpfile,outfile) - - @classmethod - def verify(cls, filename, survey, age, metallicity): - nlines=8 - with open(filename,'r') as f: - lines = [f.readline() for i in range(nlines)] - - if len(lines) < nlines: - msg = "Incorrect file size" - raise Exception(msg) - - try: - z = lines[3].split()[4] - assert np.allclose(metallicity,float(z),atol=1e-3) - except: - msg = "Metallicity does not match:\n"+lines[3] - raise Exception(msg) - - try: - s = lines[5].split()[2] - assert dict_output[survey][:4] in s - except: - msg = "Incorrect survey:\n"+lines[5] - raise Exception(msg) - - try: - a = lines[7].split('=')[1].strip().split()[0] - assert np.allclose(age,float(a),atol=1e-5) - except: - msg = "Age does not match:\n"+lines[7] - raise Exception(msg) - -class Dotter2016(Download): - """ MESA isochrones from Dotter 2016: - http://waps.cfa.harvard.edu/MIST/iso_form.php - """ - defaults = copy.deepcopy(mesa_defaults_10) - isochrone = iso.Dotter2016 - - abins = np.arange(1., 13.5 + 0.1, 0.1) - zbins = np.arange(1e-5,1e-3 + 1e-5,1e-5) - - def query_server(self, outfile, age, metallicity): - z = metallicity - feh = self.isochrone.z2feh(z) - - params = dict(self.defaults) - params['output'] = dict_output[self.survey] - params['FeH_value'] = feh - params['age_value'] = age * 1e9 - if params['age_scale'] == 'log10': - params['age_value'] = np.log10(params['age_value']) - - server = 'http://waps.cfa.harvard.edu/MIST' - url = server + '/iso_form.php' - logger.debug("Accessing %s..."%url) - response = requests.post(url,data=params) - - try: - fname = os.path.basename(response.text.split('"')[1]) - except Exception as e: - logger.debug(str(e)) - msg = 'Output filename not found' - raise RuntimeError(msg) - - tmpdir = os.path.dirname(tempfile.NamedTemporaryFile().name) - tmpfile = os.path.join(tmpdir,fname) - - out = '{0}/tmp/{1}'.format(server, fname) - cmd = 'wget %s -P %s'%(out,tmpdir) - logger.debug(cmd) - stdout = subprocess.check_output(cmd,shell=True, - stderr=subprocess.STDOUT) - logger.debug(stdout) - - cmd = 'unzip %s -d %s'%(tmpfile,tmpdir) - logger.debug(cmd) - stdout = subprocess.check_output(cmd,shell=True, - stderr=subprocess.STDOUT) - logger.debug(stdout) - - logger.debug("Creating %s..."%outfile) - shutil.move(tmpfile.replace('.zip','.cmd'),outfile) - os.remove(tmpfile) - - return outfile - - @classmethod - def verify(cls, filename, survey, age, metallicity): - age = age*1e9 - nlines=14 - with open(filename,'r') as f: - lines = [f.readline() for i in range(nlines)] - if len(lines) < nlines: - msg = "Incorrect file size" - raise Exception(msg) - - try: - s = lines[2].split()[-2] - assert dict_output[survey][:4] in s - except: - msg = "Incorrect survey:\n"+lines[2] - raise Exception(msg) - - try: - z = lines[5].split()[2] - assert np.allclose(metallicity,float(z),atol=1e-3) - except: - msg = "Metallicity does not match:\n"+lines[5] - raise Exception(msg) - - try: - a = lines[13].split()[1] - assert np.allclose(age,float(a),atol=1e-5) - except: - msg = "Age does not match:\n"+lines[13] - raise Exception(msg) - -def factory(name, **kwargs): - from ugali.utils.factory import factory - return factory(name, module=__name__, **kwargs) - -if __name__ == "__main__": - import ugali.utils.parser - description = "Download Dotter isochrones" - parser = ugali.utils.parser.Parser(description=description) - parser.add_verbose() - parser.add_force() - parser.add_argument('-a','--age',default=None,type=float,action='append') - parser.add_argument('-z','--metallicity',default=None,type=float,action='append') - parser.add_argument('-k','--kind',default='Dotter2008') - parser.add_argument('-s','--survey',default='des') - parser.add_argument('-o','--outdir',default=None) - parser.add_argument('-n','--njobs',default=1,type=int) - args = parser.parse_args() - - if args.verbose: - try: - from http.client import HTTPConnection - except ImportError: - from httplib import HTTPConnection - HTTPConnection.debuglevel = 1 - - if args.outdir is None: - args.outdir = os.path.join(args.survey.lower(),args.kind.lower()) - logger.info("Writing to output directory: %s"%args.outdir) - - p = factory(args.kind,survey=args.survey) - - abins = args.age if args.age else p.abins - zbins = args.metallicity if args.metallicity else p.zbins - grid = [g.flatten() for g in np.meshgrid(abins,zbins)] - logger.info("Ages:\n %s"%np.unique(grid[0])) - logger.info("Metallicities:\n %s"%np.unique(grid[1])) - - def run(args): - try: - p.download(*args) - except Exception as e: - # power through any exceptions... - logger.warn(str(e)) - - arguments = [(a,z,args.outdir,args.force) for a,z in zip(*grid)] - if args.njobs > 1: - msg = "Multiprocessing does not work for %s download."%args.kind - raise Exception(msg) - #elif args.njobs > 1: - # pool = Pool(processes=args.njobs, maxtasksperchild=100) - # results = pool.map(run,arguments) - else: - results = list(map(run,arguments)) diff --git a/ugali/preprocess/padova.py b/ugali/preprocess/padova.py deleted file mode 100755 index 823c6af..0000000 --- a/ugali/preprocess/padova.py +++ /dev/null @@ -1,369 +0,0 @@ -#!/usr/bin/env python -""" -DEPRECATED: ADW 2024-06-25 -New download script in ugali/scratch/download_isochrones.py - -Download Padova isochrones from: -http://stev.oapd.inaf.it/cgi-bin/cmd - -Adapted from ezpadova by Morgan Fouesneau: -https://github.com/mfouesneau/ezpadova - -""" -import os -try: - from urllib.parse import urlencode - from urllib.request import urlopen -except ImportError: - from urllib import urlencode - from urllib2 import urlopen -import re -import subprocess -from multiprocessing import Pool -from collections import OrderedDict as odict -import copy - -import numpy as np -from ugali.utils.logger import logger -from ugali.utils.shell import mkdir -from ugali.analysis.isochrone import Padova - -# survey system -photsys_dict = odict([ - ('des' ,'tab_mag_odfnew/tab_mag_decam.dat'), - ('sdss','tab_mag_odfnew/tab_mag_sloan.dat'), - ('ps1' ,'tab_mag_odfnew/tab_mag_panstarrs1.dat'), - ('lsst','tab_mag_odfnew/tab_mag_lsst.dat'), -]) - -photname_dict = odict([ - ('des' ,'DECAM'), - ('sdss','SDSS'), - ('ps1' ,'Pan-STARRS1'), - ('lsst','LSST'), -]) - -# Commented options may need to be restored for older version/isochrones. -# The parameters were tracked down by: -# Chrome -> View -> Developer -> Developer Tools -# Network -> Headers -> Request Payload - -defaults_cmd= {#'binary_frac': 0.3, - #'binary_kind': 1, - #'binary_mrinf': 0.7, - #'binary_mrsup': 1, - 'cmd_version': 2.7, - 'dust_source': 'nodust', - 'dust_sourceC': 'nodustC', - 'dust_sourceM': 'nodustM', - 'eta_reimers': 0.2, - #'extinction_av': 0, - #'icm_lim': 4, - 'imf_file': 'tab_imf/imf_chabrier_lognormal.dat', - 'isoc_age': 1e9, - 'isoc_age0': 12.7e9, - 'isoc_dlage': 0.05, - 'isoc_dz': 0.0001, - 'isoc_kind': 'parsec_CAF09_v1.2S', - 'isoc_lage0': 6.602, #Minimum allowed age - 'isoc_lage1': 10.1303, #Maximum allowed age - 'isoc_val': 0, - 'isoc_z0': 0.0001, #Minimum allowed metallicity - 'isoc_z1': 0.03, #Maximum allowed metallicity - 'isoc_zeta': 0.0002, - 'isoc_zeta0': 0.0002, - 'kind_cspecmag': 'aringer09', - 'kind_dust': 0, - 'kind_interp': 1, - 'kind_mag': 2, - 'kind_postagb': -1, - 'kind_pulsecycle': 0, - #'kind_tpagb': 0, - #'lf_deltamag': 0.2, - #'lf_maginf': 20, - #'lf_magsup': -20, - #'mag_lim': 26, - #'mag_res': 0.1, - 'output_evstage': 1, - 'output_gzip': 0, - 'output_kind': 0, - 'photsys_file': photsys_dict['des'], - #'photsys_version': 'yang', - 'submit_form': 'Submit'} - -defaults_27 = dict(defaults_cmd,cmd_version='2.7') -defaults_28 = dict(defaults_cmd,cmd_version='2.8') -defaults_29 = dict(defaults_cmd,cmd_version='2.9') -defaults_30 = dict(defaults_cmd,cmd_version='3.0') - - -class Download(object): - - isochrone = None - - def __init__(self,survey='des',**kwargs): - self.survey=survey.lower() - - def create_grid(self,abins,zbins): - arange = np.linspace(abins[0],abins[1],abins[2]+1) - zrange = np.logspace(np.log10(zbins[0]),np.log10(zbins[1]),zbins[2]+1) - aa,zz = np.meshgrid(arange,zrange) - return aa.flatten(),zz.flatten() - - def print_info(self,age,metallicity): - params = dict(age=age,z=metallicity) - params['name'] = self.__class__.__name__ - params['survey'] = self.survey - params['feh'] = self.isochrone.z2feh(metallicity) - msg = 'Downloading: %(name)s (survey=%(survey)s, age=%(age).1fGyr, Z=%(z).5f, Fe/H=%(feh).3f)'%params - logger.info(msg) - return msg - - def query_server(self,outfile,age,metallicity): - msg = "'query_server' not implemented by base class." - logger.error(msg) - raise RuntimeError(msg) - - @classmethod - def verify(cls,filename,survey,age,metallicity): - msg = "'verify' not implemented by base class." - logger.error(msg) - raise RuntimeError(msg) - - def download(self,age,metallicity,outdir=None,force=False): - """ - Check valid parameter range and download isochrones from: - http://stev.oapd.inaf.it/cgi-bin/cmd - """ - if outdir is None: outdir = './' - basename = self.isochrone.params2filename(age,metallicity) - outfile = os.path.join(outdir,basename) - - if os.path.exists(outfile) and not force: - try: - self.verify(outfile,self.survey,age,metallicity) - logger.info("Found %s; skipping..."%(outfile)) - return - except Exception as e: - msg = "Overwriting corrupted %s..."%(outfile) - logger.warn(msg) - #os.remove(outfile) - - mkdir(outdir) - - self.print_info(age,metallicity) - - try: - self.query_server(outfile,age,metallicity) - except Exception as e: - logger.debug(str(e)) - raise RuntimeError('Bad server response') - - if not os.path.exists(outfile): - raise RuntimeError('Download failed') - - try: - self.verify(outfile,self.survey,age,metallicity) - except Exception as e: - msg = "Output file is corrupted." - logger.error(msg) - #os.remove(outfile) - raise(e) - - return outfile - -class Padova(Download): - defaults = copy.deepcopy(defaults_27) - isochrone = Padova - - abins = np.arange(1.0, 13.5 + 0.1, 0.1) - zbins = np.arange(1e-4,1e-3 + 1e-5,1e-5) - - def query_server(self,outfile,age,metallicity): - epsilon = 1e-4 - lage = np.log10(age*1e9) - lage_min,lage_max = self.defaults['isoc_lage0'],self.defaults['isoc_lage1'] - if not (lage_min-epsilon < lage 1: - pool = Pool(processes=args.njobs, maxtasksperchild=100) - results = pool.map(run,arguments) - else: - results = list(map(run,arguments)) - From 0b84f9fb3e0893d3483e71098f45fe4c46bd66a9 Mon Sep 17 00:00:00 2001 From: Alex Drlica-Wagner Date: Tue, 25 Jun 2024 18:04:02 -0500 Subject: [PATCH 06/19] new isochrone download script --- ugali/scratch/download_isochrones.py | 75 ++++++++++++++++++++++++++++ 1 file changed, 75 insertions(+) create mode 100755 ugali/scratch/download_isochrones.py diff --git a/ugali/scratch/download_isochrones.py b/ugali/scratch/download_isochrones.py new file mode 100755 index 0000000..2c0f573 --- /dev/null +++ b/ugali/scratch/download_isochrones.py @@ -0,0 +1,75 @@ +#!/usr/bin/env python +""" +Script for downloading isochrone grids. +""" +__author__ = "Alex Drlica-Wagner" +import os +import re +import subprocess +from multiprocessing import Pool +from collections import OrderedDict as odict +import copy + +import numpy as np + +from ugali.utils.logger import logger +from ugali.utils.shell import mkdir +from ugali.isochrone import factory as isochrone_factory + +if __name__ == "__main__": + import ugali.utils.parser + description = "Download isochrones" + parser = ugali.utils.parser.Parser(description=description) + parser.add_verbose() + parser.add_force() + parser.add_argument('-a','--age',default=None,type=float,action='append') + parser.add_argument('-z','--metallicity',default=None,type=float,action='append') + parser.add_argument('-k','--kind',default='Marigo2017') + parser.add_argument('-s','--survey',default='lsst_dp0') + parser.add_argument('-o','--outdir',default=None) + parser.add_argument('-n','--njobs',default=1,type=int) + args = parser.parse_args() + + if args.verbose: + from httplib import HTTPConnection + HTTPConnection.debuglevel = 1 + + if args.outdir is None: + args.outdir = os.path.join(args.survey.lower(),args.kind.lower()) + logger.info("Writing to output directory: %s"%args.outdir) + + iso = isochrone_factory(args.kind,survey=args.survey) + + # Defaults + abins = [args.age] if args.age else iso.abins + zbins = [args.metallicity] if args.metallicity else iso.zbins + grid = [g.flatten() for g in np.meshgrid(abins,zbins)] + logger.info("Ages (Gyr):\n %s"%np.unique(grid[0])) + logger.info("Metallicities (Z):\n %s"%np.unique(grid[1])) + + def run(args): + try: + iso.download(*args) + return True + except Exception as e: + logger.warn(str(e)) + logger.error("Download failed.") + return False + + arglist = [(a,z,args.outdir,args.force) for a,z in zip(*grid)] + logger.info("Running %s downloads..."%(len(arglist))) + + if args.njobs > 1 and args.kind.startswith('Dotter'): + msg = "Multiprocessing does not work for %s download."%args.kind + raise Exception(msg) + elif args.njobs > 1: + pool = Pool(processes=args.njobs, maxtasksperchild=100) + results = pool.map(run,arglist) + else: + results = list(map(run,arglist)) + + results = np.array(results) + print("Number of attempted jobs: %s"%len(results)) + print("Number of succesful jobs: %s"%np.sum(results)) + print("Number of failed jobs: %s"%np.sum(~results)) + From 2e420ae33df145f5b1705871c480dde5ab619a5d Mon Sep 17 00:00:00 2001 From: Alex Drlica-Wagner Date: Tue, 25 Jun 2024 18:12:53 -0500 Subject: [PATCH 07/19] updates for LSST DP0 simulations --- .../simulation/farm_simulate_population.py | 7 +- .../scratch/simulation/simulate_population.py | 96 +++++++++++++------ .../simulation/simulate_population.yaml | 37 +++++-- 3 files changed, 99 insertions(+), 41 deletions(-) diff --git a/ugali/scratch/simulation/farm_simulate_population.py b/ugali/scratch/simulation/farm_simulate_population.py index c5a7875..124ba0c 100755 --- a/ugali/scratch/simulation/farm_simulate_population.py +++ b/ugali/scratch/simulation/farm_simulate_population.py @@ -2,8 +2,12 @@ import os import subprocess +# ADW: It'd be better to track this in the config tag = 'ps1_v1' # PS1 tag = 'des_v7' # DES +tag = 'dc2_v1' # LSST DC2 +tag = 'lsst_dc2_v2' # LSST DP0 (DC2) +tag = 'lsst_dp0_v1' # LSST DP0 (DC2) n_chunk = 100 mc_source_id_start_global = 1 size_batch = 1000 @@ -23,7 +27,8 @@ help='number of batches to submit') parser.add_argument('--mc_source_id',default=mc_source_id_start_global,type=int, help='unique identifier') -parser.add_argument('-s','--section',default='des',choices=['des','ps1'], +parser.add_argument('-s','--section',default='des', + choices=['des','ps1','lsst_dc2','lsst_dp0'], help='section of config file') parser.add_argument('--dwarfs', dest='dwarfs', action='store_true', help="Simulate from known dwarfs") diff --git a/ugali/scratch/simulation/simulate_population.py b/ugali/scratch/simulation/simulate_population.py index 7cc0847..d88286a 100755 --- a/ugali/scratch/simulation/simulate_population.py +++ b/ugali/scratch/simulation/simulate_population.py @@ -34,8 +34,15 @@ def getCompleteness(config): x = d['mag_r'] y = d['eff_star'] - x = np.insert(x, 0, 16.) - y = np.insert(y, 0, y[0]) + # Extend to saturation + if min(x) > 16: + x = np.insert(x, 0, 16.) + y = np.insert(y, 0, y[0]) + + # Make efficiency go to zero at faint end + if y[-1] > 0: + x = np.insert(x, -1, x[-1]+1) + y = np.insert(y, -1, 0) f = scipy.interpolate.interp1d(x, y, bounds_error=False, fill_value=0.) @@ -59,11 +66,13 @@ def getPhotoError(config): infile = config['photo_error'] d = np.recfromcsv(infile) - x = d['mag'] + x = d['delta_mag'] y = d['log_mag_err'] - x = np.insert(x, 0, -10.) - y = np.insert(y, 0, y[0]) + # Extend on the bright end + if min(x) > -10.0: + x = np.insert(x, 0, -10.) + y = np.insert(y, 0, y[0]) f = scipy.interpolate.interp1d(x, y, bounds_error=False, fill_value=1.) @@ -145,7 +154,7 @@ def catsimSatellite(config, lon_centroid, lat_centroid, distance, stellar_mass, completeness = getCompleteness(config) log_photo_error = getPhotoError(config) - s = ugali.analysis.source.Source() + src = ugali.analysis.source.Source() # Azimuthally averaged projected half-light radius (deg) r_h = np.degrees(np.arcsin(r_physical / distance)) @@ -159,24 +168,25 @@ def catsimSatellite(config, lon_centroid, lat_centroid, distance, stellar_mass, flag_too_extended = False max_extension = 5.0 # deg if a_h >= max_extension: - print 'Too extended: a_h = %.2f'%(a_h) + print('Too extended: a_h = %.2f'%(a_h)) a_h = max_extension flag_too_extended = True # Elliptical kernels take the "extension" as the semi-major axis extension = a_h # Elliptical half-light radius ker.setp('extension', value=a_h, bounds=[0.0,max_extension]) - s.set_kernel(ker) + src.set_kernel(ker) # Create the isochrone distance_modulus = ugali.utils.projector.dist2mod(distance) - iso = isochrone_factory('Bressan2012', survey=config['survey'], age=age, z=metal, + iso = isochrone_factory(config.get('isochrone','Bressan2012'), + survey=config['survey'], age=age, z=metal, distance_modulus=distance_modulus) - s.set_isochrone(iso) + src.set_isochrone(iso) # Simulate takes stellar mass as an argument, NOT richness - mag_1, mag_2 = s.isochrone.simulate(stellar_mass) + mag_1, mag_2 = src.isochrone.simulate(stellar_mass) # Generate the positions of stars - lon, lat = s.kernel.sample_lonlat(len(mag_2)) + lon, lat = src.kernel.sample_lonlat(len(mag_2)) nside = hp.npix2nside(len(m_maglim_1)) # Assuming that the two maglim maps have same resolution pix = ugali.utils.healpix.ang2pix(nside, lon, lat) @@ -191,12 +201,14 @@ def catsimSatellite(config, lon_centroid, lat_centroid, distance, stellar_mass, # http://iopscience.iop.org/article/10.1088/0004-637X/737/2/103/pdf mag_extinction_1 = 3.172 * m_ebv[pix] mag_extinction_2 = 2.271 * m_ebv[pix] - elif config['survey'] == 'lsst': + elif config['survey'] == 'lsst_dp0': # From Table 6 in Schlafly 2011 with Rv = 3.1 # http://iopscience.iop.org/article/10.1088/0004-637X/737/2/103/pdf mag_extinction_1 = 3.237 * m_ebv[pix] mag_extinction_2 = 2.273 * m_ebv[pix] - + else: + msg = "Unrecognized survey: %s"%config['survey'] + raise Exception(msg) # Photometric uncertainties are larger in the presence of interstellar dust reddening mag_1_error = 0.01 + 10**(log_photo_error((mag_1 + mag_extinction_1) - maglim_1)) @@ -222,15 +234,19 @@ def catsimSatellite(config, lon_centroid, lat_centroid, distance, stellar_mass, cut_detect = (np.random.uniform(size=len(mag_2)) < completeness(mag_2 + mag_extinction_2 + (23.46 - np.clip(maglim_2, 20., 26.)))) elif config['survey'] == 'ps1': cut_detect = (np.random.uniform(size=len(mag_2)) < completeness(mag_2 + mag_extinction_2)) - elif config['survey'] == 'lsst': + elif config['survey'] == 'lsst_dp0': cut_detect = (np.random.uniform(size=len(mag_2)) < completeness(mag_2 + mag_extinction_2 + (25.0 - np.clip(maglim_2, 20., 26.)))) # Using the psuedo mag depth of 25 for now + else: + msg = "Unrecognized survey: %s"%config['survey'] + raise Exception(msg) + n_g22 = np.sum(cut_detect & (mag_1 < 22.)) n_g24 = np.sum(cut_detect & (mag_1 < 24.)) print(' n_sim = %i, n_detect = %i, n_g24 = %i, n_g22 = %i'%(len(mag_1),np.sum(cut_detect),n_g24,n_g22)) - richness = stellar_mass / s.isochrone.stellarMass() - #abs_mag = s.isochrone.absolute_magnitude() - #abs_mag_martin = s.isochrone.absolute_magnitude_martin(richness=richness, n_trials=10)[0] # 100 trials seems to be sufficient for rough estimate + richness = stellar_mass / src.isochrone.stellarMass() + #abs_mag = src.isochrone.absolute_magnitude() + #abs_mag_martin = src.isochrone.absolute_magnitude_martin(richness=richness, n_trials=10)[0] # 100 trials seems to be sufficient for rough estimate #print 'abs_mag_martin = %.2f mag'%(abs_mag_martin) # The more clever thing is to sum the simulated stars @@ -243,11 +259,15 @@ def catsimSatellite(config, lon_centroid, lat_centroid, distance, stellar_mass, C_0 = -0.017 C_1 = -0.508 v = mag_1 + C_0 + C_1 * (mag_1 - mag_2) - elif config['survey'] == 'lsst': + elif config['survey'] == 'lsst_dp0': # Numbers are just placeholders for now, need to figure out exact ones C_0 = -0.02 C_1 = -0.50 v = mag_1 + C_0 + C_1 * (mag_1 - mag_2) + else: + msg = "Unrecognized survey: %s"%config['survey'] + raise Exception(msg) + flux = np.sum(10**(-v/2.5)) abs_mag = -2.5*np.log10(flux) - distance_modulus @@ -310,7 +330,7 @@ def catsimPopulation(config, tag, mc_source_id_start=1, n=5000, n_chunk=100, if not os.path.exists(tag): os.makedirs(tag) if isinstance(config,str): config = yaml.load(open(config)) - assert config['survey'] in ['des', 'ps1', 'lsst'] + assert config['survey'] in ['des', 'ps1', 'lsst_dp0'] infile_ebv = config['ebv'] infile_fracdet = config['fracdet'] @@ -357,7 +377,12 @@ def catsimPopulation(config, tag, mc_source_id_start=1, n=5000, n_chunk=100, choice_metal=[0.00010,0.00020], plot=False) area, population = ugali.simulation.population.satellitePopulation(mask, nside_pix, n, **kwargs) - + # Hack to duplicate RA,DEC + print("Hacking lon,lat...") + for idx in np.array_split(np.arange(len(population)), 10): + population['lon'][idx] = population['lon'][idx[0]] + population['lat'][idx] = population['lat'][idx[0]] + population['id'] += mc_source_id_start simulation_area = area @@ -378,8 +403,8 @@ def catsimPopulation(config, tag, mc_source_id_start=1, n=5000, n_chunk=100, mag_extinction_2_array = [] mc_source_id_array = [] for ii, mc_source_id in enumerate(population['id']): - print 'Simulating satellite (%i/%i) ... mc_source_id = %i'%(ii + 1, n, mc_source_id) - print ' distance=%(distance).2e, stellar_mass=%(stellar_mass).2e, r_physical=%(r_physical).2e'%(population[ii]) + print('Simulating satellite (%i/%i) ... mc_source_id = %i'%(ii + 1, n, mc_source_id)) + print(' distance=%(distance).2e, stellar_mass=%(stellar_mass).2e, r_physical=%(r_physical).2e'%(population[ii])) satellite = catsimSatellite(config, population[ii]['lon'], population[ii]['lat'], population[ii]['distance'], population[ii]['stellar_mass'], population[ii]['r_physical'],population[ii]['ellipticity'], @@ -399,8 +424,12 @@ def catsimPopulation(config, tag, mc_source_id_start=1, n=5000, n_chunk=100, # We assume that these objects would be easily detected and # remove them to reduce data volume - if (surface_brightness_population[ii]<23.5)&(n_g22_population[ii]>1e3): - difficulty_population[ii] |= 0b0010 + if config['survey'] == 'lsst_dp0': + if (surface_brightness_population[ii]<23.5)&(n_g24_population[ii]>1e2): + difficulty_population[ii] |= 0b0010 + else: + if (surface_brightness_population[ii]<23.5)&(n_g22_population[ii]>1e3): + difficulty_population[ii] |= 0b0010 # ADW 2019-08-31: I don't think these were implmented #if (surface_brightness_population[ii]<25.)&(n_g22_population[ii]>1e2): @@ -435,7 +464,7 @@ def catsimPopulation(config, tag, mc_source_id_start=1, n=5000, n_chunk=100, mag_extinction_2_array.append(satellite['mag_extinction_2']) mc_source_id_array.append(np.tile(mc_source_id, len(satellite['lon']))) else: - print ' difficulty=%i; satellite not simulated...'%difficulty_population[ii] + print(' difficulty=%i; satellite not simulated...'%difficulty_population[ii]) # Concatenate the arrays print("Concatenating arrays...") @@ -626,11 +655,11 @@ def catsimPopulation(config, tag, mc_source_id_start=1, n=5000, n_chunk=100, ('RFPSFMAG_SFD', [mag_2_array, 'E']), ('EXTENDED_CLASS', [np.tile(0, len(mc_source_id_array)), 'I']), ]) - elif config['survey'] == 'lsst': # Keys make to match those in the GCRCatalog native_quantities + elif config['survey'] == 'lsst_dp0': # Keys make to match those in the GCRCatalog native_quantities key_map = odict([ ('objectId', [coadd_object_id_array, 'K']), - ('coord_ra', [lon_array, 'D']), - ('coord_dec', [lat_array, 'D']), + ('ra', [lon_array, 'D']), + ('dec', [lat_array, 'D']), ('mag_g', [mag_1_array+mag_extinction_1_array, 'E']), ('mag_r', [mag_2_array+mag_extinction_2_array, 'E']), ('magerr_g', [mag_1_error_array, 'D']), @@ -639,6 +668,10 @@ def catsimPopulation(config, tag, mc_source_id_start=1, n=5000, n_chunk=100, ('mag_corrected_r', [mag_2_array, 'D']), ('extended_class', [np.tile(0, len(mc_source_id_array)), 'I']), ]) + else: + msg = "Unrecognized survey: %s"%config['survey'] + raise Exception(msg) + key_map['MC_SOURCE_ID'] = [mc_source_id_array, 'K'] print("Writing catalog files...") @@ -667,7 +700,8 @@ def catsimPopulation(config, tag, mc_source_id_start=1, n=5000, n_chunk=100, parser = argparse.ArgumentParser(description='Simulate at Milky Way satellite population.') parser.add_argument('config', help='Configuration file') - parser.add_argument('-s','--section',required=True,choices=['des','ps1','lsst'], + parser.add_argument('-s','--section',required=True, + choices=['des','ps1','lsst_dc2','lsst_dp0'], help='Config section for simulation parameters') parser.add_argument('--tag',required=True, help='Descriptive tag for the simulation run') @@ -688,7 +722,7 @@ def catsimPopulation(config, tag, mc_source_id_start=1, n=5000, n_chunk=100, np.random.seed(args.seed) # Load the config and select the survey section - config = yaml.load(open(args.config))[args.section] + config = yaml.safe_load(open(args.config))[args.section] #catsimPopulation(tag, mc_source_id_start=mc_source_id_start, n=n, n_chunk=n_chunk) catsimPopulation(config, args.tag, mc_source_id_start=args.mc_source_id_start, n=args.n, n_chunk=args.n_chunk, known_dwarfs=args.dwarfs) diff --git a/ugali/scratch/simulation/simulate_population.yaml b/ugali/scratch/simulation/simulate_population.yaml index 4a74752..34d1693 100644 --- a/ugali/scratch/simulation/simulate_population.yaml +++ b/ugali/scratch/simulation/simulate_population.yaml @@ -50,12 +50,31 @@ kbechtol_ps1: completeness: /Users/keithbechtol/Documents/DES/external_data/pan-starrs/sim/population/ps1_stellar_classification_summary_r.csv photo_error: /Users/keithbechtol/Documents/DES/external_data/pan-starrs/sim/population/ps1_photo_error_model_r.csv -lsst: - survey: lsst - stellar_density: /global/homes/m/mcnanna/sim_inputs/dc2_pseudo_stellar_density_map_cel_nside_256.npy - fracdet: /global/homes/m/mcnanna/sim_inputs/dc2_pseudo_fracdet_n2048.fits.gz - maglim_g: /global/homes/m/mcnanna/sim_inputs/dc2_pseudo_depth_n32_g_maglim.fits - maglim_r: /global/homes/m/mcnanna/sim_inputs/dc2_pseudo_depth_n32_r_maglim.fits - ebv: /global/homes/m/mcnanna/sim_inputs/ebv_sfd98_fullres_nside_4096_nest_equatorial.fits.gz - completeness: /global/homes/m/mcnanna/sim_inputs/stellar_efficiency.csv - photo_error: /global/homes/m/mcnanna/sim_inputs/photoerror_r.csv +lsst: &lsst + survey: lsst_dp0 + isochrone: 'Marigo2017' + stellar_density: /home/s1/kadrlica/projects/mw_substructure/dc2/sim_population/inputs/mcnanna/dc2_pseudo_stellar_density_map_cel_nside_256.npy + fracdet: /home/s1/kadrlica/projects/mw_substructure/dc2/sim_population/inputs/supreme_dc2_dr6d_v3_gr_fracdet.fits.gz + maglim_g: /home/s1/kadrlica/projects/mw_substructure/dc2/sim_population/inputs/supreme_dc2_dr6d_v3_g_maglim_psf_wmean.fits.gz + maglim_r: /home/s1/kadrlica/projects/mw_substructure/dc2/sim_population/inputs/supreme_dc2_dr6d_v3_r_maglim_psf_wmean.fits.gz + ebv: /home/s1/kadrlica/projects/mw_substructure/y3a2/sim_population/inputs/ebv_sfd98_fullres_nside_4096_ring_equatorial.fits + #completeness: /home/s1/kadrlica/projects/mw_substructure/dc2/sim_population/inputs/dc2_satellite_census/data/stellar_efficiency.csv + completeness: /home/s1/kadrlica/projects/mw_substructure/dc2/sim_population/inputs/dc2_satellite_census/data/stellar_efficiency_truth.csv + photo_error: /home/s1/kadrlica/projects/mw_substructure/dc2/sim_population/inputs/dc2_satellite_census/data/photoerror_r.csv + #completeness: /home/s1/kadrlica/projects/mw_substructure/dc2/sim_population/inputs/mcnanna/stellar_efficiency.csv + #photo_error: /home/s1/kadrlica/projects/mw_substructure/dc2/sim_population/inputs/mcnanna/photoerror_r.csv + +#Mapping... +lsst_dc2: *lsst +lsst_dp0: *lsst + +lsst_dc2_old: + survey: lsst_dp0 + isochrone: 'Bressan2012' + stellar_density: /home/s1/kadrlica/projects/mw_substructure/dc2/sim_population/inputs/mcnanna/dc2_pseudo_stellar_density_map_cel_nside_256.npy + fracdet: /home/s1/kadrlica/projects/mw_substructure/dc2/sim_population/inputs/mcnanna/dc2_pseudo_fracdet_n2048.fits.gz + maglim_g: /home/s1/kadrlica/projects/mw_substructure/dc2/sim_population/inputs/mcnanna/dc2_pseudo_depth_n32_g_maglim.fits + maglim_r: /home/s1/kadrlica/projects/mw_substructure/dc2/sim_population/inputs/mcnanna/dc2_pseudo_depth_n32_r_maglim.fits + ebv: /home/s1/kadrlica/projects/mw_substructure/y3a2/sim_population/inputs/ebv_sfd98_fullres_nside_4096_ring_equatorial.fits + completeness: /home/s1/kadrlica/projects/mw_substructure/dc2/sim_population/inputs/mcnanna/stellar_efficiency.csv + photo_error: /home/s1/kadrlica/projects/mw_substructure/dc2/sim_population/inputs/mcnanna/photoerror_r.csv From 7f3a516eb8b89dac67c5ae18f285e6e90337e7f3 Mon Sep 17 00:00:00 2001 From: Alex Drlica-Wagner Date: Tue, 25 Jun 2024 18:22:10 -0500 Subject: [PATCH 08/19] clean up some coordinate conversion functions --- ugali/utils/projector.py | 125 ++++++++++++++++++--------------------- 1 file changed, 58 insertions(+), 67 deletions(-) diff --git a/ugali/utils/projector.py b/ugali/utils/projector.py index 4e6f838..7de60e9 100644 --- a/ugali/utils/projector.py +++ b/ugali/utils/projector.py @@ -258,37 +258,56 @@ def angsep(lon1,lat1,lon2,lat2): ############################################################ -# ADW: Reduce numpy array operations for speed -def galToCel(ll, bb): - """ - Converts Galactic (deg) to Celestial J2000 (deg) coordinates +def gal2cel(lon, lat): + """Convert from Galactic coordinates (deg) to coordinates in the + Celestial Equatorial J2000 (deg) frame. + + Parameters: + ----------- + lon : Galactic longitude (deg) + lat : Galactic latitude (deg) + + Returns + ------- + ra,dec : Right ascension and declination (deg,deg) + """ - bb = np.radians(bb) - sin_bb = np.sin(bb) - cos_bb = np.cos(bb) + lat = np.radians(lat) + sin_lat = np.sin(lat) + cos_lat = np.cos(lat) - ll = np.radians(ll) + lon = np.radians(lon) ra_gp = np.radians(192.85948) de_gp = np.radians(27.12825) lcp = np.radians(122.932) - sin_lcp_ll = np.sin(lcp - ll) - cos_lcp_ll = np.cos(lcp - ll) + sin_lcp_lon = np.sin(lcp - lon) + cos_lcp_lon = np.cos(lcp - lon) - sin_d = (np.sin(de_gp) * sin_bb) \ - + (np.cos(de_gp) * cos_bb * cos_lcp_ll) - ramragp = np.arctan2(cos_bb * sin_lcp_ll, - (np.cos(de_gp) * sin_bb) \ - - (np.sin(de_gp) * cos_bb * cos_lcp_ll)) + sin_d = (np.sin(de_gp) * sin_lat) \ + + (np.cos(de_gp) * cos_lat * cos_lcp_lon) + ramragp = np.arctan2(cos_lat * sin_lcp_lon, + (np.cos(de_gp) * sin_lat) \ + - (np.sin(de_gp) * cos_lat * cos_lcp_lon)) dec = np.arcsin(sin_d) ra = (ramragp + ra_gp + (2. * np.pi)) % (2. * np.pi) return np.degrees(ra), np.degrees(dec) -gal2cel = galToCel +# DEPRECATED +galToCel = gal2cel -def celToGal(ra, dec): - """ - Converts Celestial J2000 (deg) to Calactic (deg) coordinates +def cel2gal(ra, dec): + """Convert coordinates in the Celestial Equatorial J2000 (deg) frame + to Galactic (deg) coordinates. + + Parameters: + ----------- + ra : Right ascension (deg) + dec : Declination (deg) + + Returns: + -------- + lon, lat : Galactic longitude and latitude (deg, deg) """ dec = np.radians(dec) sin_dec = np.sin(dec) @@ -307,11 +326,12 @@ def celToGal(ra, dec): lcpml = np.arctan2(cos_dec * sin_ra_gp, (np.cos(de_gp) * sin_dec) \ - (np.sin(de_gp) * cos_dec * cos_ra_gp)) - bb = np.arcsin(sin_b) - ll = (lcp - lcpml + (2. * np.pi)) % (2. * np.pi) - return np.degrees(ll), np.degrees(bb) + lat = np.arcsin(sin_b) + lon = (lcp - lcpml + (2. * np.pi)) % (2. * np.pi) + return np.degrees(lon), np.degrees(lat) -cel2gal = celToGal +# DEPRECATED +celToGal = cel2gal def estimate_angle(angle, origin, new_frame, offset=1e-7): """ @@ -339,51 +359,17 @@ def cel2gal_angle(ra,dec,angle,offset=1e-7): origin = SkyCoord(ra,dec,unit=u.deg,frame='fk5') return estimate_angle(angle,origin,'galactic',offset) -### ADW: Probably works, remember 90-pa for kernel convention... -### def gal2cel_angle(glon,glat,angle): -### """ -### WARNING: Not vectorized -### """ -### gal_angle = np.radians(angle) -### galproj = Projector(glon,glat) -### x,y = np.sin(gal_angle),np.cos(gal_angle) -### alon,alat = galproj.imageToSphere(0.1*x,0.1*y) -### -### ra,dec = gal2cel(glon,glat) -### ara,adec = gal2cel(alon,alat) -### celproj = Projector(ra,dec) -### cel_x,cel_y = celproj.sphereToImage(ara,adec) -### cel_angle = np.degrees(np.arctan2(cel_x,cel_y)) -### return cel_angle + 360*(cel_angle<0) - - -### ADW: WARNING DOESN'T WORK YET -### def cel2gal_angle(ra,dec,angle): -### """ -### WARNING: Not vectorized -### """ -### cel_angle = np.radians(angle) -### celproj = Projector(ra,dec) -### x,y = np.sin(cel_angle),np.cos(cel_angle) -### angle_ra,angle_dec = celproj.imageToSphere(1e-2*x,1e-2*y) -### -### glon,glat = gal2cel(ra,dec) -### angle_glon,angle_glat = cel2gal(angle_ra,angle_dec) -### galproj = Projector(glon,glat) -### gal_x,gal_y = galproj.sphereToImage(angle_glon,angle_glat) -### gal_angle = np.degrees(np.arctan2(gal_x,gal_y)) -### return gal_angle + 360*(gal_angle<0) - ############################################################ def dec2hms(dec): - """ - ADW: This should be replaced by astropy... + """Convert from decimal degrees to hours-minutes-seconds. + ADW: This should be replaced by astropy... from astropy.coordinates import Angle hms = Angle(dec*u.deg).hms return (hms.h,hms.m,hms.s) + """ DEGREE = 360. HOUR = 24. @@ -401,12 +387,13 @@ def dec2hms(dec): return (hour, minute, second) def dec2dms(dec): - """ - ADW: This should be replaced by astropy + """Convert from decimal degrees to degrees-minutes-seconds. + ADW: This should be replaced by astropy from astropy.coordinates import Angle dms = Angle(dec*u.deg).dms return (dms.d,dms.m,dms.s) + """ DEGREE = 360. HOUR = 24. @@ -429,11 +416,11 @@ def dec2dms(dec): return (deg, minute, second) def hms2dec(hms): - """ - Convert longitude from hours,minutes,seconds in string or 3-array + """Convert longitude from hours,minutes,seconds in string or 3-array format to decimal degrees. ADW: This really should be replaced by astropy + """ DEGREE = 360. HOUR = 24. @@ -449,9 +436,11 @@ def hms2dec(hms): return decimal def dms2dec(dms): - """ - Convert latitude from degrees,minutes,seconds in string or 3-array - format to decimal degrees. + """Convert latitude from degrees,minutes,seconds in string or 3-array + format to decimal degrees. + + ADW: This really should be replaced by astropy + """ DEGREE = 360. HOUR = 24. @@ -473,9 +462,11 @@ def dms2dec(dms): return decimal def sr2deg(solid_angle): + """ Convert solid angle from steradians to square deg.""" return np.degrees(np.degrees(solid_angle)) def deg2sr(solid_angle): + """ Convert solid angle from square deg to steradians""" return np.radians(np.radians(solid_angle)) ############################################################ From 09d3ae498b59e7f394130fdcd131438acf4abc24 Mon Sep 17 00:00:00 2001 From: Alex Drlica-Wagner Date: Tue, 25 Jun 2024 18:24:31 -0500 Subject: [PATCH 09/19] import fitsio to read_map --- ugali/utils/healpix.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/ugali/utils/healpix.py b/ugali/utils/healpix.py index d92b4aa..a03621d 100644 --- a/ugali/utils/healpix.py +++ b/ugali/utils/healpix.py @@ -565,7 +565,7 @@ def merge_likelihood_headers2(filenames, outfile, **kwargs): return data_dict -def read_map(filename, nest=False, hdu=None, h=False, verbose=True): +def read_map(filename, nest=False, hdu=None, h=False, verbose=False): """Read a healpix map from a fits file. Partial-sky files, if properly identified, are expanded to full size and filled with UNSEEN. Uses fitsio to mirror much (but not all) of the functionality of healpy.read_map @@ -590,7 +590,7 @@ def read_map(filename, nest=False, hdu=None, h=False, verbose=True): m [, header] : array, optionally with header appended The map read from the file, and the header if *h* is True. """ - + import fitsio data,hdr = fitsio.read(filename,header=True,ext=hdu) nside = int(hdr.get('NSIDE')) From cdd871b34f7d7d0546af61628a1adf3b07a46afe Mon Sep 17 00:00:00 2001 From: Alex Drlica-Wagner Date: Wed, 14 Jul 2021 22:21:20 -0500 Subject: [PATCH 10/19] Update README.md [skip ci] --- README.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/README.md b/README.md index 2b1adc5..1966fda 100644 --- a/README.md +++ b/README.md @@ -50,7 +50,7 @@ The `ugali` source code is distributed with several auxiliary libraries for isoc ``` cd $UGALIDIR -wget https://github.com/kadrlica/ugali/releases/download/v1.7.0rc0/ugali-des-bressan2012.tar.gz +wget https://github.com/DarkEnergySurvey/ugali/releases/download/v1.7.0/ugali-des-bressan2012.tar.gz tar -xzf ugali-des-bressan2012.tar.gz ``` From 20a18d827ed707d1a3a6540c1ddf5293f93172a2 Mon Sep 17 00:00:00 2001 From: Alex Drlica-Wagner Date: Wed, 14 Jul 2021 22:13:20 -0500 Subject: [PATCH 11/19] linting --- .github/workflows/python-package.yml | 16 ++-- ugali/analysis/prior.py | 5 ++ ugali/analysis/results.py | 3 +- ugali/candidate/associate.py | 4 +- ugali/isochrone/dartmouth.py | 2 +- ugali/isochrone/empirical.py | 3 + ugali/isochrone/mesa.py | 2 +- ugali/isochrone/padova.py | 17 +++- ugali/isochrone/parsec.py | 6 +- ugali/observation/mask.py | 11 +-- ugali/observation/roi.py | 1 + ugali/preprocess/database.py | 8 +- ugali/scratch/PlotAllSkyHealpix.py | 4 +- ugali/scratch/simulation/assemble.py | 5 +- .../simulation/merge_population_files.py | 10 ++- .../scratch/simulation/simulate_population.py | 2 +- ugali/scratch/simulation/simulation_match.py | 20 ++--- .../simulation/survey_selection_function.py | 21 ++--- .../simulation/validate_sim_population.py | 8 +- ugali/simulation/analyzer.py | 17 ++-- ugali/simulation/population.py | 32 ++------ ugali/utils/batch.py | 2 +- ugali/utils/fileio.py | 79 +++++-------------- ugali/utils/healpix.py | 17 ++-- ugali/utils/idl.py | 8 +- ugali/utils/logger.py | 1 + ugali/utils/plotting.py | 11 --- ugali/utils/projector.py | 10 +-- ugali/utils/stats.py | 1 + 29 files changed, 136 insertions(+), 190 deletions(-) diff --git a/.github/workflows/python-package.yml b/.github/workflows/python-package.yml index fa862a1..c8b745b 100644 --- a/.github/workflows/python-package.yml +++ b/.github/workflows/python-package.yml @@ -7,9 +7,9 @@ name: build on: push: - branches: [ master, actions ] pull_request: - branches: [ master ] + schedule: + - cron: '0 0 1 * *' jobs: build: @@ -38,19 +38,17 @@ jobs: source activate env export UGALIDIR="$HOME/.ugali" python setup.py -q install --isochrones --catalogs - - name: Install test data - run: | - wget https://github.com/DarkEnergySurvey/ugali/releases/download/v1.7.0/ugali-test-data.tar.gz -O ugali-test-data.tar.gz - tar -xzf ugali-test-data.tar.gz - name: Lint with flake8 - if: ${{ false }} + if: ${{ true }} run: | source activate env conda install -q flake8 -c conda-forge # stop the build if there are Python syntax errors or undefined names flake8 . --count --select=E9,F63,F7,F82 --show-source --statistics - # exit-zero treats all errors as warnings. The GitHub editor is 127 chars wide - flake8 . --count --exit-zero --max-complexity=10 --max-line-length=127 --statistics + - name: Install test data + run: | + wget https://github.com/DarkEnergySurvey/ugali/releases/download/v1.7.0/ugali-test-data.tar.gz -O ugali-test-data.tar.gz + tar -xzf ugali-test-data.tar.gz - name: Test with nose run: | source activate env diff --git a/ugali/analysis/prior.py b/ugali/analysis/prior.py index 58ac621..204b93f 100644 --- a/ugali/analysis/prior.py +++ b/ugali/analysis/prior.py @@ -1,4 +1,9 @@ #!/usr/bin/env python +""" +Incomplete module for dealing with priors... +""" + +import scipy.stats class Prior(object): def __call__(self, value): diff --git a/ugali/analysis/results.py b/ugali/analysis/results.py index c98796c..99877d6 100644 --- a/ugali/analysis/results.py +++ b/ugali/analysis/results.py @@ -18,7 +18,8 @@ import ugali.utils.stats from ugali.utils.stats import Samples -from ugali.utils.projector import dist2mod,mod2dist,gal2cel,gal2cel_angle +from ugali.utils.projector import dist2mod,mod2dist +from ugali.utils.projector import cel2gal,gal2cel,gal2cel_angle from ugali.utils.projector import ang2const, ang2iau from ugali.utils.config import Config from ugali.utils.logger import logger diff --git a/ugali/candidate/associate.py b/ugali/candidate/associate.py index 4ca91ce..1667ee8 100644 --- a/ugali/candidate/associate.py +++ b/ugali/candidate/associate.py @@ -54,7 +54,7 @@ # glon, glat = ugali.utils.projector.celToGal(lon,lat) # else: # glon,glat = lon, lat -# return ugali.utils.projector.match(glon,glat,self.data['glon'],self.data['glat'],tol) +# return ugali.utils.projector.match(glon,glat,self['glon'],self['glat'],tol) @@ -452,7 +452,7 @@ def catalogFactory(name, **kwargs): catalogs = odict(inspect.getmembers(sys.modules[__name__], fn)) if name not in list(catalogs.keys()): - msg = "%s not found in catalogs:\n %s"%(name,list(kernels.keys())) + msg = "%s not found in catalogs:\n %s"%(name,list(catalogs.keys())) logger.error(msg) msg = "Unrecognized catalog: %s"%name raise Exception(msg) diff --git a/ugali/isochrone/dartmouth.py b/ugali/isochrone/dartmouth.py index 1e8fc36..944cbf9 100644 --- a/ugali/isochrone/dartmouth.py +++ b/ugali/isochrone/dartmouth.py @@ -165,7 +165,7 @@ def _parse(self,filename): try: columns = self.columns[self.survey.lower()] except KeyError as e: - logger.warning('Unrecognized survey: %s'%(survey)) + logger.warning('Unrecognized survey: %s'%(self.survey)) raise(e) kwargs = dict(comments='#',usecols=list(columns.keys()),dtype=list(columns.values())) diff --git a/ugali/isochrone/empirical.py b/ugali/isochrone/empirical.py index 06b316f..b309a13 100644 --- a/ugali/isochrone/empirical.py +++ b/ugali/isochrone/empirical.py @@ -3,7 +3,10 @@ Generic python script. """ import os +from collections import OrderedDict as odict +from ugali.analysis.model import Model, Parameter +from ugali.utils.shell import mkdir, get_ugali_dir, get_iso_dir from ugali.isochrone.parsec import PadovaIsochrone class EmpiricalPadova(PadovaIsochrone): diff --git a/ugali/isochrone/mesa.py b/ugali/isochrone/mesa.py index 70df82b..28ca103 100644 --- a/ugali/isochrone/mesa.py +++ b/ugali/isochrone/mesa.py @@ -134,7 +134,7 @@ def _parse(self,filename): try: columns = self.columns[self.survey.lower()] except KeyError as e: - logger.warning('Unrecognized survey: %s'%(survey)) + logger.warning('Unrecognized survey: %s'%(self.survey)) raise(e) kwargs = dict(comments='#',usecols=list(columns.keys()),dtype=list(columns.values())) diff --git a/ugali/isochrone/padova.py b/ugali/isochrone/padova.py index 9c919ad..b6236c3 100644 --- a/ugali/isochrone/padova.py +++ b/ugali/isochrone/padova.py @@ -1,7 +1,20 @@ #!/usr/bin/env python """ Older (pre-PARSEC) Padova isochrones. + +Not verified to work... """ +import os.path +from collections import OrderedDict as odict + +import numpy as np + +from ugali.analysis.model import Model, Parameter +from ugali.utils.shell import mkdir, get_ugali_dir, get_iso_dir +from ugali.isochrone.model import Isochrone +from ugali.isochrone.parsec import Marigo2017 as Padova +from ugali.isochrone.parsec import defaults_27 +from ugali.utils.logger import logger class Girardi2002(Padova): defaults = dict(defaults_27) @@ -19,7 +32,7 @@ class Girardi2010b(Padova): defaults = dict(defaults_27) defaults['isoc_kind'] = 'gi10b' -class Girardi2002(PadovaIsochrone): +class Girardi2002(Padova): #_dirname = '/u/ki/kadrlica/des/isochrones/v5/' _dirname = os.path.join(get_iso_dir(),'{survey}','girardi2002') # For use with Marigo et al. (2008) and earlier use Anders & Grevesse 1989 @@ -57,7 +70,7 @@ def _parse(self,filename): try: columns = self.columns[self.survey.lower()] except KeyError as e: - logger.warning('did not recognize survey %s'%(survey)) + logger.warning('did not recognize survey %s'%(self.survey)) raise(e) kwargs = dict(delimiter='\t',usecols=list(columns.keys()),dtype=list(columns.values())) diff --git a/ugali/isochrone/parsec.py b/ugali/isochrone/parsec.py index 7890ff4..af99f24 100644 --- a/ugali/isochrone/parsec.py +++ b/ugali/isochrone/parsec.py @@ -418,7 +418,7 @@ def _parse(self,filename): try: columns = self.columns[self.survey.lower()] except KeyError as e: - logger.warning('Unrecognized survey: %s'%(survey)) + logger.warning('Unrecognized survey: %s'%(self.survey)) raise(e) # delimiter='\t' is used to be compatible with OldPadova... @@ -509,7 +509,7 @@ def _find_column_numbers(self): """ Map from the isochrone column names to the column numbers. """ header,lines = self.parse_header(self.filename) columns = header['columns'] - + def _parse(self,filename): """Reads an isochrone file in the Marigo et al. 2017 format. Creates arrays with the initial stellar mass and @@ -526,7 +526,7 @@ def _parse(self,filename): try: columns = self.columns[self.survey.lower()] except KeyError as e: - logger.warning('Unrecognized survey: %s'%(survey)) + logger.warning('Unrecognized survey: %s'%(self.survey)) raise(e) kwargs = dict(usecols=list(columns.keys()),dtype=list(columns.values())) diff --git a/ugali/observation/mask.py b/ugali/observation/mask.py index 5ad6242..3c41de2 100644 --- a/ugali/observation/mask.py +++ b/ugali/observation/mask.py @@ -474,10 +474,8 @@ def backgroundMMD(self, catalog, method='cloud-in-cells', weights=None): [self.roi.bins_mag,self.roi.bins_mag], weights=number_density)[0] elif method == 'bootstrap': - # Not implemented + # Not implemented; see catalog.bootstrap raise ValueError("Bootstrap method not implemented") - mag_1 + (mag_1_err * np.random.normal(0, 1., len(mag_1))) - mag_2 + (mag_2_err * np.random.normal(0, 1., len(mag_2))) elif method == 'histogram': # Apply raw histogram @@ -577,13 +575,8 @@ def backgroundCMD(self, catalog, mode='cloud-in-cells', weights=None): [self.roi.bins_color,self.roi.bins_mag], weights=number_density)[0] elif mode == 'bootstrap': - # Not implemented + # Not implemented; see catalog.bootstrap raise ValueError("Bootstrap mode not implemented") - mag_1_array = catalog.mag_1 - mag_2_array = catalog.mag_2 - - catalog.mag_1 + (catalog.mag_1_err * np.random.normal(0, 1., len(catalog.mag_1))) - catalog.mag_2 + (catalog.mag_2_err * np.random.normal(0, 1., len(catalog.mag_2))) elif mode == 'histogram': # Apply raw histogram diff --git a/ugali/observation/roi.py b/ugali/observation/roi.py index adadb25..75a5127 100644 --- a/ugali/observation/roi.py +++ b/ugali/observation/roi.py @@ -18,6 +18,7 @@ import ugali.utils.projector import ugali.utils.skymap +from ugali.utils.logger import logger from ugali.utils.config import Config from ugali.utils.healpix import query_disc, ang2pix, pix2ang, ang2vec diff --git a/ugali/preprocess/database.py b/ugali/preprocess/database.py index 13404bf..1d6ceec 100755 --- a/ugali/preprocess/database.py +++ b/ugali/preprocess/database.py @@ -299,15 +299,15 @@ def footprint(self,nside): """ Download the survey footprint for HEALpix pixels. """ - import healpy + import healpy as hp import ugali.utils.projector if nside > 2**9: raise Exception("Overflow error: nside must be <=2**9") - pix = np.arange(healpy.nside2npix(nside),dtype='int') - footprint = np.zeros(healpy.nside2npix(nside),dtype='bool') + pix = np.arange(hp.nside2npix(nside),dtype='int') + footprint = np.zeros(hp.nside2npix(nside),dtype='bool') ra,dec = ugali.utils.projector.pixToAng(nside,pix) table_name = 'Pix%i'%nside self.upload(np.array([pix,ra,dec]).T, ['pix','ra','dec'], name=table_name) - radius = healpy.nside2resol(nside_superpix,arcmin=True) + radius = hp.nside2resol(nside,arcmin=True) query=""" SELECT t.pix, dbo.fInFootprintEq(t.ra, t.dec, %g) diff --git a/ugali/scratch/PlotAllSkyHealpix.py b/ugali/scratch/PlotAllSkyHealpix.py index 35a871d..fae5ad6 100644 --- a/ugali/scratch/PlotAllSkyHealpix.py +++ b/ugali/scratch/PlotAllSkyHealpix.py @@ -1,9 +1,11 @@ #!/usr/bin/env python import healpy import pylab as plt +import numpy + import ugali.utils.skymap from ugali.utils.projector import celToGal -import numpy +from ugali.utils.logger import logger default_kwargs = dict( xytext=(5,5),textcoords='offset points', ha="left",va="center", diff --git a/ugali/scratch/simulation/assemble.py b/ugali/scratch/simulation/assemble.py index f989fe3..7192b89 100644 --- a/ugali/scratch/simulation/assemble.py +++ b/ugali/scratch/simulation/assemble.py @@ -1,3 +1,5 @@ +DeprecationWarning("'assemble.py' should be removed") + import sys import os import shutil @@ -9,7 +11,8 @@ config = yaml.load(open(config_file)) if os.path.exists(outdir): - raw_input('Are you sure you want to continue? [Press ENTER to continue]') + print("Output directory exists: %s"%outdir) + input('Are you sure you want to continue? [Press ENTER to continue]') if not os.path.exists(outdir): os.mkdir(outdir) diff --git a/ugali/scratch/simulation/merge_population_files.py b/ugali/scratch/simulation/merge_population_files.py index 25e67c8..1d5553f 100644 --- a/ugali/scratch/simulation/merge_population_files.py +++ b/ugali/scratch/simulation/merge_population_files.py @@ -1,9 +1,11 @@ +DeprecationWarning("'merge_population_files.py' should be removed") + import sys import glob import numpy as np import astropy.io.fits as pyfits -print sys.argv +print(sys.argv) infiles = sorted(glob.glob(sys.argv[1])) outfile = sys.argv[2] @@ -11,7 +13,7 @@ data_array = [] header_array = [] for infile in infiles: - print infile + print(infile) reader = pyfits.open(infile) data_array.append(reader[1].data) header_array.append(reader[1].header) @@ -20,11 +22,11 @@ data_array = np.concatenate(data_array) -print '\nWill write output to %s\n'%(outfile) +print('\nWill write output to %s\n'%(outfile)) tbhdu = pyfits.BinTableHDU(data_array) tbhdu.header = header_array[0] -raw_input('Continue?') +input('Continue?') tbhdu.writeto(outfile) diff --git a/ugali/scratch/simulation/simulate_population.py b/ugali/scratch/simulation/simulate_population.py index d88286a..a64623b 100755 --- a/ugali/scratch/simulation/simulate_population.py +++ b/ugali/scratch/simulation/simulate_population.py @@ -290,7 +290,7 @@ def catsimSatellite(config, lon_centroid, lat_centroid, distance, stellar_mass, pylab.savefig('y3_sat_sim_cmd_%s.png'%('test'), dpi=150.) print('n_Sigma_p = %i'%(n_sigma_p)) - raw_input('WAIT') + input('WAIT') satellite=odict(lon=lon[cut_detect], lat=lat[cut_detect], mag_1=mag_1_meas[cut_detect], mag_2=mag_2_meas[cut_detect], diff --git a/ugali/scratch/simulation/simulation_match.py b/ugali/scratch/simulation/simulation_match.py index 866f6c9..7cc8476 100644 --- a/ugali/scratch/simulation/simulation_match.py +++ b/ugali/scratch/simulation/simulation_match.py @@ -1,3 +1,5 @@ +DeprecationWarning("'simulation_match.py' should be removed") + import numpy as np import astropy.io.fits as pyfits import pylab @@ -40,7 +42,7 @@ def wrap(x): data_sim = reader_sim[1].data reader_sim.close() -print len(data_sim) +print(len(data_sim)) """ match_search, match_sim, angsep = ugali.utils.projector.match(data_search['ra'], data_search['dec'], data_sim['RA'], data_sim['DEC'], tol=1.) @@ -162,13 +164,13 @@ def wrap(x): if save: pylab.savefig('sim_match_distance_luminosity.pdf') -print '%20s%20s%20s%20s%20s%20s'%('mc_source_id', 'ra', 'dec', 'distance_modulus', 'fracdet', 'density') +print('%20s%20s%20s%20s%20s%20s'%('mc_source_id', 'ra', 'dec', 'distance_modulus', 'fracdet', 'density')) for index in np.nonzero(cut_why_not)[0]: - print '%20i%20.3f%20.3f%20.3f%20.3f%20.3f'%(data_sim['mc_source_id'][index], + print('%20i%20.3f%20.3f%20.3f%20.3f%20.3f'%(data_sim['mc_source_id'][index], data_sim['ra'][index], data_sim['dec'][index], data_sim['distance_modulus'][index], data_sim['fracdet'][index], - data_sim['density'][index]) + data_sim['density'][index])) ############################################################ @@ -190,7 +192,7 @@ def wrap(x): classifier_file = 'trained_classifier.txt' if fit: - print 'Training the machine learning classifier. This may take a while ...' + print('Training the machine learning classifier. This may take a while ...') t_start = time.time() classifier = sklearn.gaussian_process.GaussianProcessClassifier(1.0 * sklearn.gaussian_process.kernels.RBF(0.5)) #classifier = sklearn.neighbors.KNeighborsClassifier(3, weights='uniform') @@ -198,18 +200,18 @@ def wrap(x): #classifier = sklearn.svm.SVC(gamma=2, C=1) classifier.fit(x[cut_train], y[cut_train]) t_end = time.time() - print ' ... training took %.2f seconds'%(t_end - t_start) + print(' ... training took %.2f seconds'%(t_end - t_start)) # Save the trained classifier classifier_data = pickle.dumps(classifier) writer = open(classifier_file, 'w') writer.write(classifier_data) writer.close() - print 'Saving machine learning classifier to %s ...'%(classifier_file) + print('Saving machine learning classifier to %s ...'%(classifier_file)) os.system('gzip %s'%(classifier_file)) else: - print 'Loading machine learning classifier from %s ...'%(classifier_file) + print('Loading machine learning classifier from %s ...'%(classifier_file)) if os.path.exists(classifier_file + '.gz') and not os.path.exists(classifier_file): - print ' Unzipping...' + print(' Unzipping...') os.system('gunzip -k %s.gz'%(classifier_file)) reader = open(classifier_file) classifier_data = ''.join(reader.readlines()) diff --git a/ugali/scratch/simulation/survey_selection_function.py b/ugali/scratch/simulation/survey_selection_function.py index 9719da5..eff06da 100644 --- a/ugali/scratch/simulation/survey_selection_function.py +++ b/ugali/scratch/simulation/survey_selection_function.py @@ -1,5 +1,5 @@ """ - +Create a survey selection function """ import time @@ -9,6 +9,7 @@ import numpy as np import healpy as hp import astropy.io.fits as pyfits +import matplotlib import matplotlib.pyplot as plt import pylab @@ -150,7 +151,7 @@ def __init__(self, config_file): def loadFracdet(self): if self.m_fracdet is None: - print 'Loading fracdet map from %s ...'%(self.config['infile']['fracdet']) + print('Loading fracdet map from %s ...'%(self.config['infile']['fracdet'])) self.m_fracdet = hp.read_map(self.config['infile']['fracdet'], nest=False) def loadPopulationMetadata(self): @@ -165,7 +166,7 @@ def loadSimResults(self): def loadRealResults(self): if self.data_real is None: - print 'Loading real data search results from %s ...'%(self.config[self.algorithm]['real_results']) + print('Loading real data search results from %s ...'%(self.config[self.algorithm]['real_results'])) reader = pyfits.open(self.config[self.algorithm]['real_results']) self.data_real = reader[1].data reader.close() @@ -206,7 +207,7 @@ def trainClassifier(self): # Train random forest classifier if True: - print 'Training the machine learning classifier. This may take a while ...' + print('Training the machine learning classifier. This may take a while ...') t_start = time.time() parameters = {'n_estimators':(500,1000)}#, 'criterion':["gini","entropy"], "min_samples_leaf": [1,2,4]} rf = RandomForestClassifier(oob_score=True) @@ -223,14 +224,14 @@ def trainClassifier(self): print(self.classifier.best_estimator_) print(self.classifier.best_params_) t_end = time.time() - print ' ... training took %.2f seconds'%(t_end - t_start) + print(' ... training took %.2f seconds'%(t_end - t_start)) # Save the trained classifier classifier_data = pickle.dumps(self.classifier) writer = open(self.config[self.algorithm]['classifier'], 'w') writer.write(classifier_data) writer.close() - print 'Saving machine learning classifier to %s ...'%(self.config[self.algorithm]['classifier']) + print('Saving machine learning classifier to %s ...'%(self.config[self.algorithm]['classifier'])) else: self.loadClassifier() @@ -310,7 +311,8 @@ def validateClassifier(self, cut_detect, cut_train, cut_geometry, y_pred): 'actual': 'black', 'hsc': 'black'} - title = 'N_train = %i ; N_test = %i'%(len(cut_train),len(cut_test)) import matplotlib + title = 'N_train = %i ; N_detect = %i'%(len(cut_train),len(cut_detect)) + cmap = matplotlib.colors.ListedColormap(['Gold', 'Orange', 'DarkOrange', 'OrangeRed', 'Red']) pylab.figure() @@ -373,9 +375,9 @@ def validateClassifier(self, cut_detect, cut_train, cut_geometry, y_pred): # pylab.savefig('sim_match_scatter_prediction.pdf')#, dpi=dpi) def loadClassifier(self): - print 'Loading machine learning classifier from %s ...'%(self.config[self.algorithm]['classifier']) + print('Loading machine learning classifier from %s ...'%(self.config[self.algorithm]['classifier'])) if os.path.exists(self.config[self.algorithm]['classifier'] + '.gz') and not os.path.exists(self.config[self.algorithm]['classifier']): - print ' Unzipping...' + print(' Unzipping...') os.system('gunzip -k %s.gz'%(self.config[self.algorithm]['classifier'])) reader = open(self.config[self.algorithm]['classifier']) classifier_data = ''.join(reader.readlines()) @@ -449,7 +451,6 @@ def predict(self, lon, lat, **kwargs): return pred, flags_geometry def validatePredict(self, pred, flags_geometry, lon, lat, r_physical, abs_mag, distance): - import matplotlib cmap = matplotlib.colors.ListedColormap(['Gold', 'Orange', 'DarkOrange', 'OrangeRed', 'Red']) pylab.figure() diff --git a/ugali/scratch/simulation/validate_sim_population.py b/ugali/scratch/simulation/validate_sim_population.py index 207f384..038ef96 100644 --- a/ugali/scratch/simulation/validate_sim_population.py +++ b/ugali/scratch/simulation/validate_sim_population.py @@ -48,11 +48,11 @@ def getCatalog(catalog_dir): catalog_infiles = sorted(glob.glob(catalog_dir + '/*catalog*.fits')) data_array = [] for catalog_infile in catalog_infiles: - print ' Reading %s ...'%(catalog_infile) + print(' Reading %s ...'%(catalog_infile)) reader = pyfits.open(catalog_infile) data_array.append(reader[1].data) reader.close() - print ' Merging ...' + print(' Merging ...') return np.concatenate(data_array) ########## @@ -64,7 +64,7 @@ def getCatalog(catalog_dir): infile_population = 'v4/sim_population_v4.fits' reader_population = pyfits.open(infile_population) data_population = reader_population[1].data -print len(data_population) +print(len(data_population)) data_population = data_population #[0:500] reader_population.close() @@ -207,7 +207,7 @@ def getCatalog(catalog_dir): """ ########## """ -print "Machine learning" +print("Machine learning") save = False diff --git a/ugali/simulation/analyzer.py b/ugali/simulation/analyzer.py index 1d0f8d3..a18b4d0 100755 --- a/ugali/simulation/analyzer.py +++ b/ugali/simulation/analyzer.py @@ -30,16 +30,13 @@ from ugali.utils import mlab # Analysis flags -FLAGS = odict([ - ('FLAG_PROC' , 0 ), # Simulation was processed - ('FLAG_NOPROC' , 1 ), # No processing - ('FLAG_NOBJ' , 2 ), # Too many catalog objects - ('FLAG_FIT' , 4 ), # Fit failure - ('FLAG_EBV' , 8 ), # EBV value too large - ('FLAG_MEM' , 16), # Memory error - ]) -for k,v in FLAGS.items(): - globals()[k] = v +FLAGS = odict([]) +FLAGS['FLAG_PROC' ] = FLAG_PROC = 0 # Simulation was processed +FLAGS['FLAG_NOPROC'] = FLAG_NOPROC = 1 # No processing +FLAGS['FLAG_NOBJ' ] = FLAG_NOBJ = 2 # Too many catalog objects +FLAGS['FLAG_FIT' ] = FLAG_FIT = 4 # Fit failure +FLAGS['FLAG_EBV' ] = FLAG_EBV = 8 # EBV value too large +FLAGS['FLAG_MEM' ] = FLAG_MEM = 16 # Memory error def update_header_flags(filename): fits = fitsio.FITS(filename,'rw') diff --git a/ugali/simulation/population.py b/ugali/simulation/population.py index 0bc4bf7..1765dcb 100644 --- a/ugali/simulation/population.py +++ b/ugali/simulation/population.py @@ -104,7 +104,6 @@ def satellitePopulation(mask, nside_pix, size, ############################################################ -# ADW: 2019-09-01 DEPRECATED def satellitePopulationOrig(config, n, range_distance_modulus=[16.5, 24.], range_stellar_mass=[1.e2, 1.e5], @@ -133,7 +132,9 @@ def satellitePopulationOrig(config, n, lon (deg), lat (deg), distance modulus, stellar mass (Msun), and half-light radius (kpc) for each satellite """ - + msg = "'satellitePopulationOrig': ADW 2019-09-01" + DeprecationWarning(msg) + if type(config) == str: config = ugali.utils.config.Config(config) @@ -155,8 +156,8 @@ def satellitePopulationOrig(config, n, np.log10(range_stellar_mass[1]), n) - half_light_radius_physical = 10**np.random.uniform(np.log10(range_half_light_radius_physical[0]), - np.log10(range_half_light_radius_physical[0]), + half_light_radius_physical = 10**np.random.uniform(np.log10(range_r_physical[0]), + np.log10(range_r_physical[1]), n) # kpc half_light_radius = np.degrees(np.arcsin(half_light_radius_physical \ @@ -272,26 +273,3 @@ def knownPopulation(dwarfs, mask, nside_pix, size): population = np.hstack(results) population['id'] = np.arange(size) return area, population - -def plot_population(population): - # ADW: DEPRECATED: 2019-09-01 - pylab.figure() - #pylab.scatter(lon, lat, c=distance_modulus, s=500 * half_light_radius) - #pylab.colorbar() - pylab.scatter(lon, lat, edgecolors='none') - xmin, xmax = pylab.xlim() # Reverse azimuthal axis - pylab.xlim([xmax, xmin]) - pylab.title('Random Positions in Survey Footprint') - pylab.xlabel('Longitude (deg)') - pylab.ylabel('Latitude (deg)') - - pylab.figure() - pylab.scatter(stellar_mass, distance,c=r_physical, - s=500 * r_physical, edgecolors='none') - pylab.xscale('log') - pylab.yscale('log') - pylab.xlim([0.5 * range_stellar_mass[0], 2. * range_stellar_mass[1]]) - pylab.colorbar() - pylab.title('Half-light Radius (arcmin)') - pylab.xlabel('Stellar Mass (arcmin)') - pylab.ylabel('Distance (kpc)') diff --git a/ugali/utils/batch.py b/ugali/utils/batch.py index e8cb4d7..371c240 100644 --- a/ugali/utils/batch.py +++ b/ugali/utils/batch.py @@ -312,7 +312,7 @@ def parse_options(self, **opts): class Condor(Batch): """ Not implemented yet... """ - def __init__(self): + def __init__(self, **kwargs): super(Condor,self).__init__(**kwargs) logger.warning('Condor cluster is untested') diff --git a/ugali/utils/fileio.py b/ugali/utils/fileio.py index 9f1caa1..04708fa 100644 --- a/ugali/utils/fileio.py +++ b/ugali/utils/fileio.py @@ -66,7 +66,23 @@ def write(filename,data,**kwargs): msg = "Unrecognized file type: %s"%filename raise ValueError(msg) - + +def check_formula(formula): + """ Check that a formula is valid. """ + if not (('data[' in formula) and (']' in formula)): + msg = 'Invalid formula:\n %s'%formula + raise ValueError(msg) + +def parse_formula(formula): + """ Return the columns used in a formula. """ + check_formula(formula) + + columns = [] + for x in formula.split('data[')[1:]: + col = x.split(']')[0].replace('"','').replace("'",'') + columns += [col] + + return columns def add_column(filename,column,formula,force=False): """ Add a column to a FITS file. @@ -246,65 +262,6 @@ def insert_columns(filename,data,ext=1,force=False,colnum=None): # Dealing with FITS files def write_fits(filename,data,header=None,force=False): if os.path.exists(filename) and not force: - found(filename) + logger.found(filename) return fitsio.write(filename,data,header=header,clobber=force) - -# Writing membership files -def write_membership(loglike,filename): - """ - Write a catalog file of the likelihood region including - membership properties. - - Parameters: - ----------- - loglike : input loglikelihood object - filename : output filename - - Returns: - -------- - None - """ - - ra,dec = gal2cel(loglike.catalog.lon,loglike.catalog.lat) - - name_objid = loglike.config['catalog']['objid_field'] - name_mag_1 = loglike.config['catalog']['mag_1_field'] - name_mag_2 = loglike.config['catalog']['mag_2_field'] - name_mag_err_1 = loglike.config['catalog']['mag_err_1_field'] - name_mag_err_2 = loglike.config['catalog']['mag_err_2_field'] - - # Angular and isochrone separations - sep = angsep(loglike.source.lon,loglike.source.lat, - loglike.catalog.lon,loglike.catalog.lat) - isosep = loglike.isochrone.separation(loglike.catalog.mag_1,loglike.catalog.mag_2) - - data = odict() - data[name_objid] = loglike.catalog.objid - data['GLON'] = loglike.catalog.lon - data['GLAT'] = loglike.catalog.lat - data['RA'] = ra - data['DEC'] = dec - data[name_mag_1] = loglike.catalog.mag_1 - data[name_mag_err_1] = loglike.catalog.mag_err_1 - data[name_mag_2] = loglike.catalog.mag_2 - data[name_mag_err_2] = loglike.catalog.mag_err_2 - data['COLOR'] = loglike.catalog.color - data['ANGSEP'] = sep - data['ISOSEP'] = isosep - data['PROB'] = loglike.p - - # HIERARCH allows header keywords longer than 8 characters - header = [] - for param,value in loglike.source.params.items(): - card = dict(name='HIERARCH %s'%param.upper(), - value=value.value, - comment=param) - header.append(card) - card = dict(name='HIERARCH %s'%'TS',value=loglike.ts(), - comment='test statistic') - header.append(card) - card = dict(name='HIERARCH %s'%'TIMESTAMP',value=time.asctime(), - comment='creation time') - header.append(card) - fitsio.write(filename,data,header=header,clobber=True) diff --git a/ugali/utils/healpix.py b/ugali/utils/healpix.py index a03621d..b4696ce 100644 --- a/ugali/utils/healpix.py +++ b/ugali/utils/healpix.py @@ -8,11 +8,12 @@ import numpy import numpy as np import healpy as hp -import healpy +import fitsio import ugali.utils.projector from ugali.utils import fileio from ugali.utils.logger import logger +import ugali.utils.fileio ############################################################ @@ -382,7 +383,6 @@ def write_partial_map(filename, data, nside, coord=None, nest=False, -------- None """ - import fitsio # ADW: Do we want to make everything uppercase? @@ -596,9 +596,9 @@ def read_map(filename, nest=False, hdu=None, h=False, verbose=False): nside = int(hdr.get('NSIDE')) if verbose: print('NSIDE = {0:d}'.format(nside)) - if not healpy.isnsideok(nside): + if not hp.isnsideok(nside): raise ValueError('Wrong nside parameter.') - sz=healpy.nside2npix(nside) + sz=hp.nside2npix(nside) ordering = hdr.get('ORDERING','UNDEF').strip() if verbose: print('ORDERING = {0:s} in fits file'.format(ordering)) @@ -615,28 +615,27 @@ def read_map(filename, nest=False, hdu=None, h=False, verbose=False): # Could be done more efficiently (but complicated) by reordering first if hdr['INDXSCHM'] == 'EXPLICIT': - m = healpy.UNSEEN*np.ones(sz,dtype=data[fields[1]].dtype) + m = hp.UNSEEN*np.ones(sz,dtype=data[fields[1]].dtype) m[data[fields[0]]] = data[fields[1]] else: m = data[fields[0]].ravel() - if (not healpy.isnpixok(m.size) or (sz>0 and sz != m.size)) and verbose: + if (not hp.isnpixok(m.size) or (sz>0 and sz != m.size)) and verbose: print('nside={0:d}, sz={1:d}, m.size={2:d}'.format(nside,sz,m.size)) raise ValueError('Wrong nside parameter.') if nest is not None: if nest and ordering.startswith('RING'): - idx = healpy.nest2ring(nside,np.arange(m.size,dtype=np.int32)) + idx = hp.nest2ring(nside,np.arange(m.size,dtype=np.int32)) m = m[idx] if verbose: print('Ordering converted to NEST') elif (not nest) and ordering.startswith('NESTED'): - idx = healpy.ring2nest(nside,np.arange(m.size,dtype=np.int32)) + idx = hp.ring2nest(nside,np.arange(m.size,dtype=np.int32)) m = m[idx] if verbose: print('Ordering converted to RING') if h: return m, hdr else: return m - if __name__ == "__main__": import argparse description = __doc__ diff --git a/ugali/utils/idl.py b/ugali/utils/idl.py index d6efdcf..14bced8 100644 --- a/ugali/utils/idl.py +++ b/ugali/utils/idl.py @@ -97,8 +97,8 @@ def jprecess(ra, dec, mu_radec=None, parallax=None, rad_vel=None, epoch=None): ra = array(ra) dec = array(dec) else: - ra = array([ra0]) - dec = array([dec0]) + ra = array([ra]) + dec = array([dec]) n = ra.size if rad_vel is None: @@ -111,7 +111,7 @@ def jprecess(ra, dec, mu_radec=None, parallax=None, rad_vel=None, epoch=None): if (mu_radec is not None): if (array(mu_radec).size != 2 * n): - raise Exception('ERROR - MU_RADEC keyword (proper motion) be dimensioned (2,' + strtrim(n, 2) + ')') + raise Exception('ERROR - MU_RADEC keyword (proper motion) be dimensioned (2,%s)'%n) mu_radec = mu_radec * 1. if parallax is None: @@ -315,7 +315,7 @@ def bprecess(ra0, dec0, mu_radec=None, parallax=None, rad_vel=None, epoch=None): if (mu_radec is not None): if (array(mu_radec).size != 2 * n): - raise Exception('ERROR - MU_RADEC keyword (proper motion) be dimensioned (2,' + strtrim(n, 2) + ')') + raise Exception('ERROR - MU_RADEC keyword (proper motion) be dimensioned (2,%s)'%n) mu_radec = mu_radec * 1. if parallax is None: diff --git a/ugali/utils/logger.py b/ugali/utils/logger.py index 8a2d454..8863abc 100644 --- a/ugali/utils/logger.py +++ b/ugali/utils/logger.py @@ -3,6 +3,7 @@ Interface to python logging. For more info see: http://docs.python.org/2/howto/logging.html """ +import os.path import logging class SpecialFormatter(logging.Formatter): diff --git a/ugali/utils/plotting.py b/ugali/utils/plotting.py index ff2953b..c289bf3 100644 --- a/ugali/utils/plotting.py +++ b/ugali/utils/plotting.py @@ -585,8 +585,6 @@ def drawCMD(self, ax=None, radius=None, zidx=None): ax.set_xlabel('Color (mag)') ax.set_ylabel('Magnitude (mag)') - try: ax.cax.colorbar(im) - except: pass ax.annotate("Stars",**self.label_kwargs) @@ -1417,15 +1415,6 @@ def plot_chain(chain,burn=None,clip=None): ################################################### - -def plotSkymapCatalog(lon,lat,**kwargs): - """ - Plot a catalog of coordinates on a full-sky map. - """ - fig = plt.figure() - ax = plt.subplot(111,projection=projection) - drawSkymapCatalog(ax,lon,lat,**kwargs) - def drawSkymapCatalog(ax,lon,lat,**kwargs): mapping = { 'ait':'aitoff', diff --git a/ugali/utils/projector.py b/ugali/utils/projector.py index 7de60e9..9de825a 100644 --- a/ugali/utils/projector.py +++ b/ugali/utils/projector.py @@ -4,7 +4,7 @@ Based on Calabretta & Greisen 2002, A&A, 357, 1077-1122 http://adsabs.harvard.edu/abs/2002A%26A...395.1077C """ - +import re import numpy as np from ugali.utils.logger import logger @@ -31,8 +31,9 @@ def setReference(self, lon_ref, lat_ref, zenithal=False): if not zenithal: phi = (-np.pi / 2.) + np.radians(lon_ref) theta = np.radians(lat_ref) - psi = np.radians(90.) # psi = 90 corresponds to (0, 0) psi = -90 corresponds to (180, 0) - + # psi = 90 corresponds to (0, 0) + # psi = -90 corresponds to (180, 0) + psi = np.radians(90.) cos_psi,sin_psi = np.cos(psi),np.sin(psi) cos_phi,sin_phi = np.cos(phi),np.sin(phi) @@ -440,7 +441,6 @@ def dms2dec(dms): format to decimal degrees. ADW: This really should be replaced by astropy - """ DEGREE = 360. HOUR = 24. @@ -452,7 +452,7 @@ def dms2dec(dms): # http://docs.scipy.org/doc/numpy-1.7.0/reference/c-api.coremath.html#NPY_NZERO if isstring(dms): - degree,minute,second = np.array(re.split('[dms]',hms))[:3].astype(float) + degree,minute,second = np.array(re.split('[dms]',dms))[:3].astype(float) else: degree,minute,second = dms.T diff --git a/ugali/utils/stats.py b/ugali/utils/stats.py index 8ec502f..505530a 100644 --- a/ugali/utils/stats.py +++ b/ugali/utils/stats.py @@ -4,6 +4,7 @@ """ import copy +from collections import OrderedDict as odict import numpy as np import numpy.lib.recfunctions as recfuncs From 9d0692df890fa99444721348a4384e02d34572b6 Mon Sep 17 00:00:00 2001 From: Alex Drlica-Wagner Date: Thu, 15 Jul 2021 12:50:26 -0500 Subject: [PATCH 12/19] Remove deprecated code --- ugali/analysis/scan.py | 11 +- ugali/isochrone/model.py | 22 ---- ugali/observation/mask.py | 185 ++++----------------------------- ugali/observation/roi.py | 14 +-- ugali/preprocess/database.py | 3 +- ugali/simulation/population.py | 90 ---------------- ugali/utils/config.py | 4 +- ugali/utils/skymap.py | 17 --- 8 files changed, 29 insertions(+), 317 deletions(-) diff --git a/ugali/analysis/scan.py b/ugali/analysis/scan.py index 93c1ea3..e3dca85 100755 --- a/ugali/analysis/scan.py +++ b/ugali/analysis/scan.py @@ -130,10 +130,12 @@ def search(self, coords=None, distance_modulus=None, extension=None, tolerance=1 self.stellar_mass_array = np.zeros([nmoduli,npixels],dtype='f4') self.fraction_observable_array = np.zeros([nmoduli,npixels],dtype='f4') self.extension_fit_array = np.zeros([nmoduli,npixels],dtype='f4') - # DEPRECATED: ADW 2019-04-27 - self.richness_lower_array = np.zeros([nmoduli,npixels],dtype='f4') - self.richness_upper_array = np.zeros([nmoduli,npixels],dtype='f4') - self.richness_ulimit_array = np.zeros([nmoduli,npixels],dtype='f4') + if self.config['scan']['full_pdf']: + # DEPRECATED: ADW 2019-04-27 + DeprecationWarning("'full_pdf' is deprecated.") + self.richness_lower_array = np.zeros([nmoduli,npixels],dtype='f4') + self.richness_upper_array = np.zeros([nmoduli,npixels],dtype='f4') + self.richness_ulimit_array = np.zeros([nmoduli,npixels],dtype='f4') # Specific pixel/distance_modulus coord_idx, distance_modulus_idx, extension_idx = None, None, None @@ -315,6 +317,7 @@ def write(self, outfile): data['PIXEL']=self.roi.pixels_target # Full data output (too large for survey) if self.config['scan']['full_pdf']: + DeprecationWarning("'full_pdf' is deprecated.") data['LOG_LIKELIHOOD']=self.loglike_array.T data['RICHNESS']=self.richness_array.T data['RICHNESS_LOWER']=self.richness_lower_array.T diff --git a/ugali/isochrone/model.py b/ugali/isochrone/model.py index 2ac1af9..d673934 100644 --- a/ugali/isochrone/model.py +++ b/ugali/isochrone/model.py @@ -322,27 +322,6 @@ def stellar_luminosity(self, steps=10000): return np.sum(luminosity * d_log_mass * self.imf.pdf(mass, log_mode=True)) - def stellar_luminosity2(self, steps=10000): - """ - DEPRECATED: ADW 2017-09-20 - - Compute the stellar luminosity (L_Sol; average per star). - Uses "sample" to generate mass sample and pdf. The range of - integration only covers the input isochrone data (no - extrapolation used), but this seems like a sub-percent effect - if the isochrone goes to 0.15 Msun for the old and metal-poor - stellar populations of interest. - - Note that the stellar luminosity is very sensitive to the - post-AGB population. - """ - msg = "'%s.stellar_luminosity2': ADW 2017-09-20"%self.__class__.__name__ - DeprecationWarning(msg) - mass_init, mass_pdf, mass_act, mag_1, mag_2 = self.sample(mass_steps=steps) - luminosity_interpolation = scipy.interpolate.interp1d(self.mass_init, self.luminosity,fill_value=0,bounds_error=False) - luminosity = luminosity_interpolation(mass_init) - return np.sum(luminosity * mass_pdf) - # ADW: For temporary backward compatibility stellarMass = stellar_mass stellarLuminosity = stellar_luminosity @@ -354,7 +333,6 @@ def absolute_magnitude(self, richness=1, steps=1e4): g,r -> V transform equations from Jester 2005 [astro-ph/0506022]. - TODO: ADW If richness not specified, should use self.richness Parameters: diff --git a/ugali/observation/mask.py b/ugali/observation/mask.py index 3c41de2..c59146e 100644 --- a/ugali/observation/mask.py +++ b/ugali/observation/mask.py @@ -13,7 +13,6 @@ import healpy as hp import scipy.signal -#import ugali.utils.plotting import ugali.utils.binning import ugali.utils.skymap @@ -403,51 +402,19 @@ def photo_err_2(delta): return self.photo_err_1, self.photo_err_2 - def plotSolidAngleCMD(self): - """ - Solid angle within the mask as a function of color and magnitude. - """ - msg = "'%s.plotSolidAngleCMD': ADW 2018-05-05"%self.__class__.__name__ - DeprecationWarning(msg) - - import ugali.utils.plotting - ugali.utils.plotting.twoDimensionalHistogram('mask', 'color', 'magnitude', - self.solid_angle_cmd, - self.roi.bins_color, - self.roi.bins_mag, - lim_x = [self.roi.bins_color[0], - self.roi.bins_color[-1]], - lim_y = [self.roi.bins_mag[-1], - self.roi.bins_mag[0]]) - - def plotSolidAngleMMD(self): - """ - Solid angle within the mask as a function of color and magnitude. - """ - msg = "'%s.plotSolidAngleMMD': ADW 2018-05-05"%self.__class__.__name__ - DeprecationWarning(msg) - - import ugali.utils.plotting - - ugali.utils.plotting.twoDimensionalHistogram('mask', 'mag_2', 'mag_1', - self.solid_angle_mmd, - self.roi.bins_mag, - self.roi.bins_mag, - lim_x = [self.roi.bins_mag[0], - self.roi.bins_mag[-1]], - lim_y = [self.roi.bins_mag[-1], - self.roi.bins_mag[0]]) - - - def backgroundMMD(self, catalog, method='cloud-in-cells', weights=None): """ Generate an empirical background model in magnitude-magnitude space. - INPUTS: - catalog: Catalog object - OUTPUTS: - background + Parameters + ---------- + catalog: catalog object + method: method for estimated MMD + weights: weights assigned to each catalog object + + Returns + ------- + background: estimate of background magnitude-magnitude distribution """ # Select objects in annulus @@ -545,10 +512,15 @@ def backgroundCMD(self, catalog, mode='cloud-in-cells', weights=None): """ Generate an empirical background model in color-magnitude space. - INPUTS: - catalog: Catalog object - OUTPUTS: - background + Parameters + ---------- + catalog: catalog object + method: method for estimated MMD + weights: weights assigned to each catalog object + + Returns + ------- + background: estimate of background color-magnitude distribution """ # Select objects in annulus @@ -840,24 +812,6 @@ def depth(self, lon, lat): """ pass - def plot(self): - """ - Plot the magnitude depth. - """ - msg = "'%s.plot': ADW 2018-05-05"%self.__class__.__name__ - DeprecationWarning(msg) - - import ugali.utils.plotting - - mask = hp.UNSEEN * np.ones(hp.nside2npix(self.nside)) - mask[self.roi.pixels] = self.mask_roi_sparse - mask[mask == 0.] = hp.UNSEEN - ugali.utils.plotting.zoomedHealpixMap('Completeness Depth', - mask, - self.roi.lon, self.roi.lat, - self.roi.config.params['coords']['roi_radius']) - - class CoverageBand(object): """ Map of coverage fraction for a single observing band. @@ -916,106 +870,3 @@ def __init__(self, maglim, roi): self.mask_roi_sparse = mask[self.roi.pixels] ############################################################ -def simpleMask(config): - - #params = ugali.utils.(config, kwargs) - - roi = ugali.observation.roi.ROI(config) - - # De-project the bin centers to get magnitude depths - - mesh_x, mesh_y = np.meshgrid(roi.centers_x, roi.centers_y) - r = np.sqrt(mesh_x**2 + mesh_y**2) # Think about x, y conventions here - - #z = (0. * (r > 1.)) + (21. * (r < 1.)) - #z = 21. - r - #z = (21. - r) * (mesh_x > 0.) * (mesh_y < 0.) - z = (21. - r) * np.logical_or(mesh_x > 0., mesh_y > 0.) - - return MaskBand(z, roi) - -############################################################ - -def readMangleFile(infile, lon, lat, index = None): - """ - DEPRECATED: 2018-05-04 - Mangle must be set up on your system. - The index argument is a temporary file naming convention to avoid file conflicts. - Coordinates must be given in the native coordinate system of the Mangle file. - """ - msg = "'mask.readMangleFile': ADW 2018-05-05" - DeprecationWarning(msg) - - if index is None: - index = np.random.randint(0, 1.e10) - - coordinate_file = 'temp_coordinate_%010i.dat'%(index) - maglim_file = 'temp_maglim_%010i.dat'%(index) - - writer = open(coordinate_file, 'w') - for ii in range(0, len(lon)): - writer.write('%12.5f%12.5f\n'%(lon[ii], lat[ii])) - writer.close() - - os.system('polyid -W %s %s %s || exit'%(infile, - coordinate_file, - maglim_file)) - - reader = open(maglim_file) - lines = reader.readlines() - reader.close() - - os.remove(maglim_file) - os.remove(coordinate_file) - - maglim = [] - for ii in range(1, len(lines)): - if len(lines[ii].split()) == 3: - maglim.append(float(lines[ii].split()[2])) - elif len(lines[ii].split()) == 2: - maglim.append(0.) # Coordinates outside of the MANGLE ploygon - elif len(lines[ii].split()) > 3: - msg = 'Coordinate inside multiple polygons, masking that coordinate.' - logger.warning(msg) - maglim.append(0.) - else: - msg = 'Cannot parse maglim file, unexpected number of columns.' - logger.error(msg) - break - - maglim = np.array(maglim) - return maglim - -############################################################ - -def allSkyMask(infile, nside): - msg = "'mask.allSkyMask': ADW 2018-05-05" - DeprecationWarning(msg) - lon, lat = ugali.utils.skymap.allSkyCoordinates(nside) - maglim = readMangleFile(infile, lon, lat, index = None) - return maglim - -############################################################ - -def scale(mask, mag_scale, outfile=None): - """ - Scale the completeness depth of a mask such that mag_new = mag + mag_scale. - Input is a full HEALPix map. - Optionally write out the scaled mask as an sparse HEALPix map. - """ - msg = "'mask.scale': ADW 2018-05-05" - DeprecationWarning(msg) - mask_new = hp.UNSEEN * np.ones(len(mask)) - mask_new[mask == 0.] = 0. - mask_new[mask > 0.] = mask[mask > 0.] + mag_scale - - if outfile is not None: - pix = np.nonzero(mask_new > 0.)[0] - data_dict = {'MAGLIM': mask_new[pix]} - nside = hp.npix2nside(len(mask_new)) - ugali.utils.skymap.writeSparseHealpixMap(pix, data_dict, nside, outfile) - - return mask_new - -############################################################ - diff --git a/ugali/observation/roi.py b/ugali/observation/roi.py index 75a5127..f79c49f 100644 --- a/ugali/observation/roi.py +++ b/ugali/observation/roi.py @@ -122,22 +122,12 @@ def __init__(self, config, lon, lat): self.delta_mag = self.bins_mag[1] - self.bins_mag[0] self.delta_color = self.bins_color[1] - self.bins_color[0] - ## Axis labels (DEPRECATED) - #self.label_x = 'x (deg)' - #self.label_y = 'y (deg)' - # - #if self.config.params['catalog']['band_1_detection']: - # self.label_mag = '%s (mag)'%(self.config.params['catalog']['mag_1_field']) - #else: - # self.label_mag = '%s (mag)'%(self.config.params['catalog']['mag_2_field']) - #self.label_color = '%s - %s (mag)'%(self.config.params['catalog']['mag_1_field'], - # self.config.params['catalog']['mag_2_field']) - def plot(self, value=None, pixel=None): """ Plot the ROI """ - # DEPRECATED + # DEPRECATED: ADW 2021-07-15 + DeprecationWarning("'roi.plot' is deprecated and will be removed.") import ugali.utils.plotting map_roi = np.array(hp.UNSEEN \ diff --git a/ugali/preprocess/database.py b/ugali/preprocess/database.py index 1d6ceec..7894aaf 100755 --- a/ugali/preprocess/database.py +++ b/ugali/preprocess/database.py @@ -327,8 +327,7 @@ def __init__(self,release='SVA1_GOLD'): def _setup_desdbi(self): # Function here to setup trivialAccess client... # This should work but it doesn't - import warnings - warnings.warn("desdbi is deprecated", DeprecationWarning) + DeprecationWarning("'desdbi' is deprecated") import despydb.desdbi def generate_query(self, ra_min,ra_max,dec_min,dec_max,filename,db): diff --git a/ugali/simulation/population.py b/ugali/simulation/population.py index 1765dcb..1e8a9b4 100644 --- a/ugali/simulation/population.py +++ b/ugali/simulation/population.py @@ -104,96 +104,6 @@ def satellitePopulation(mask, nside_pix, size, ############################################################ -def satellitePopulationOrig(config, n, - range_distance_modulus=[16.5, 24.], - range_stellar_mass=[1.e2, 1.e5], - range_r_physical=[5.e-3, 1.], - mode='mask', - plot=False): - """ - Create a population of `n` randomly placed satellites within a - survey mask or catalog specified in the config file. Satellites - are distributed uniformly in distance modulus, uniformly in - log(stellar_mass / Msun), and uniformly in - log(r_physical/kpc). The ranges can be set by the user. - - Returns the simulated area (deg^2) as well as the - Parameters - ---------- - config : configuration file or object - n : number of satellites to simulate - range_distance_modulus : range of distance modulus to sample - range_stellar_mass : range of stellar mass to sample (Msun) - range_r_physical : range of azimuthally averaged half light radius (kpc) - mode : how to sample the area['mask','catalog'] - - Returns - ------- - lon (deg), lat (deg), distance modulus, stellar mass (Msun), and - half-light radius (kpc) for each satellite - """ - msg = "'satellitePopulationOrig': ADW 2019-09-01" - DeprecationWarning(msg) - - if type(config) == str: - config = ugali.utils.config.Config(config) - - if mode == 'mask': - mask_1 = ugali.utils.skymap.readSparseHealpixMap(config.params['mask']['infile_1'], 'MAGLIM') - mask_2 = ugali.utils.skymap.readSparseHealpixMap(config.params['mask']['infile_2'], 'MAGLIM') - input = (mask_1 > 0.) * (mask_2 > 0.) - elif mode == 'catalog': - catalog = ugali.observation.catalog.Catalog(config) - input = np.array([catalog.lon, catalog.lat]) - - lon, lat, simulation_area = ugali.utils.skymap.randomPositions(input, - config.params['coords']['nside_likelihood_segmentation'], - n=n) - distance_modulus = np.random.uniform(range_distance_modulus[0], - range_distance_modulus[1], - n) - stellar_mass = 10**np.random.uniform(np.log10(range_stellar_mass[0]), - np.log10(range_stellar_mass[1]), - n) - - half_light_radius_physical = 10**np.random.uniform(np.log10(range_r_physical[0]), - np.log10(range_r_physical[1]), - n) # kpc - - half_light_radius = np.degrees(np.arcsin(half_light_radius_physical \ - / ugali.utils.projector.distanceModulusToDistance(distance_modulus))) - - # One choice of theory prior - #half_light_radius_physical = ugali.analysis.kernel.halfLightRadius(stellar_mass) # kpc - #half_light_radius = np.degrees(np.arcsin(half_light_radius_physical \ - # / ugali.utils.projector.distanceModulusToDistance(distance_modulus))) - - if plot: - pylab.figure() - #pylab.scatter(lon, lat, c=distance_modulus, s=500 * half_light_radius) - #pylab.colorbar() - pylab.scatter(lon, lat, edgecolors='none') - xmin, xmax = pylab.xlim() # Reverse azimuthal axis - pylab.xlim([xmax, xmin]) - pylab.title('Random Positions in Survey Footprint') - pylab.xlabel('Longitude (deg)') - pylab.ylabel('Latitude (deg)') - - pylab.figure() - pylab.scatter(stellar_mass, ugali.utils.projector.distanceModulusToDistance(distance_modulus), - c=(60. * half_light_radius), s=500 * half_light_radius, edgecolors='none') - pylab.xscale('log') - pylab.yscale('log') - pylab.xlim([0.5 * range_stellar_mass[0], 2. * range_stellar_mass[1]]) - pylab.colorbar() - pylab.title('Half-light Radius (arcmin)') - pylab.xlabel('Stellar Mass (arcmin)') - pylab.ylabel('Distance (kpc)') - - return simulation_area, lon, lat, distance_modulus, stellar_mass, half_light_radius - -############################################################ - def interpolate_absolute_magnitude(): iso = ugali.isochrone.factory('Bressan2012',age=12,z=0.00010) diff --git a/ugali/utils/config.py b/ugali/utils/config.py index 6a0d17e..535dac2 100644 --- a/ugali/utils/config.py +++ b/ugali/utils/config.py @@ -200,9 +200,6 @@ def _createFilenames(self,pixels=None): """ nside_catalog = self['coords']['nside_catalog'] - # Deprecated: ADW 2018-06-17 - #if nside_catalog is None: - # pixels = [None] if pixels is not None: pixels = [pixels] if np.isscalar(pixels) else pixels else: @@ -225,6 +222,7 @@ def _createFilenames(self,pixels=None): if pix is None: # DEPRECTATED: ADW 2018-06-17 # This is not really being used anymore + raise ValueError('pix cannot be None') catalog = os.path.join(catalog_dir,catalog_base) mask_1 = os.path.join(mask_dir,mask_base_1) mask_2 = os.path.join(mask_dir,mask_base_2) diff --git a/ugali/utils/skymap.py b/ugali/utils/skymap.py index 5937ff3..cca24e8 100644 --- a/ugali/utils/skymap.py +++ b/ugali/utils/skymap.py @@ -84,23 +84,6 @@ def inFootprint(config, pixels, nside=None): return inside -def footprint(config, nside=None): - """ - UNTESTED. - Should return a boolean array representing the pixels in the footprint. - """ - DeprecationWarning("This function is deprecated.") - config = Config(config) - if nside is None: - nside = config['coords']['nside_pixel'] - elif nside < config['coords']['nside_catalog']: - raise Exception('Requested nside=%i is greater than catalog_nside'%nside) - elif nside > config['coords']['nside_pixel']: - raise Exception('Requested nside=%i is less than pixel_nside'%nside) - pix = np.arange(hp.nside2npix(nside), dtype=int) - return inFootprint(config,pix) - - ############################################################ def allSkyCoordinates(nside): From 8565ad605d1e49aabdb8d87a310a6baca3bb924d Mon Sep 17 00:00:00 2001 From: Alex Drlica-Wagner Date: Fri, 16 Jul 2021 09:27:20 -0500 Subject: [PATCH 13/19] fix typo in mesa isochrone parser (closes #69) --- ugali/isochrone/mesa.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/ugali/isochrone/mesa.py b/ugali/isochrone/mesa.py index 28ca103..0b05fed 100644 --- a/ugali/isochrone/mesa.py +++ b/ugali/isochrone/mesa.py @@ -79,7 +79,7 @@ class Dotter2016(Isochrone): des = odict([ (2, ('mass_init',float)), (3, ('mass_act',float)), - (8, ('log_lum',float)), + (6, ('log_lum',float)), (9, ('u',float)), (10,('g',float)), (11,('r',float)), From 25db1c018b208af48d7c8b6a8088f4574a379084 Mon Sep 17 00:00:00 2001 From: Alex Drlica-Wagner Date: Thu, 15 Jul 2021 23:18:46 -0500 Subject: [PATCH 14/19] adding test data install to setup.py --- .github/workflows/python-package.yml | 7 ++- setup.py | 45 ++++++++++---- tests/config.yaml | 15 ++--- ugali/utils/config.py | 90 +++------------------------- 4 files changed, 53 insertions(+), 104 deletions(-) diff --git a/.github/workflows/python-package.yml b/.github/workflows/python-package.yml index c8b745b..e68b720 100644 --- a/.github/workflows/python-package.yml +++ b/.github/workflows/python-package.yml @@ -18,6 +18,8 @@ jobs: max-parallel: 3 matrix: python-version: [2.7, 3.8] + env: + UGALIDIR: "$HOME/.ugali" steps: - uses: actions/checkout@v2 - name: Set up Python ${{ matrix.python-version }} @@ -36,8 +38,8 @@ jobs: - name: Install package run: | source activate env - export UGALIDIR="$HOME/.ugali" - python setup.py -q install --isochrones --catalogs + #export UGALIDIR="$HOME/.ugali" + python setup.py -q install --isochrones --catalogs --tests - name: Lint with flake8 if: ${{ true }} run: | @@ -46,6 +48,7 @@ jobs: # stop the build if there are Python syntax errors or undefined names flake8 . --count --select=E9,F63,F7,F82 --show-source --statistics - name: Install test data + if: ${{ false }} run: | wget https://github.com/DarkEnergySurvey/ugali/releases/download/v1.7.0/ugali-test-data.tar.gz -O ugali-test-data.tar.gz tar -xzf ugali-test-data.tar.gz diff --git a/setup.py b/setup.py index 201c725..b789caa 100644 --- a/setup.py +++ b/setup.py @@ -24,25 +24,27 @@ def find_packages(): HERE = os.path.abspath(os.path.dirname(__file__)) URL = 'https://github.com/DarkEnergySurvey/ugali' DESC = "Ultra-faint galaxy likelihood toolkit." -LONG_DESC = "See %s"%URL +LONG_DESC = "%s\n%s"%(DESC,URL) CLASSIFIERS = """\ -Development Status :: 3 - Alpha +Development Status :: 4 - Beta Intended Audience :: Science/Research Intended Audience :: Developers License :: OSI Approved :: MIT License Natural Language :: English Operating System :: MacOS :: MacOS X Operating System :: POSIX :: Linux -Programming Language :: Python :: 2.7 +Programming Language :: Python :: 2 +Programming Language :: Python :: 3 Topic :: Scientific/Engineering Topic :: Scientific/Engineering :: Astronomy Topic :: Scientific/Engineering :: Physics """ -RELEASE = URL+'/releases/download/v1.7.0' +RELEASE_URL = URL+'/releases/download/v1.7.0' UGALIDIR = os.getenv("UGALIDIR","$HOME/.ugali") ISOSIZE = "~1MB" CATSIZE = "~20MB" +TSTSIZE = "~1MB" # Could find file size dynamically, but it's a bit slow... # int(urllib.urlopen(ISOCHRONES).info().getheaders("Content-Length")[0])/1024**2 SURVEYS = ['des','ps1','sdss','lsst'] @@ -77,7 +79,7 @@ class TarballCommand(distutils.cmd.Command,object): 'force installation (overwrite any existing files)') ] boolean_options = ['force'] - release = RELEASE + release_url = RELEASE_URL _tarball = None _dirname = None @@ -114,7 +116,7 @@ def install_tarball(self, tarball): os.makedirs(self.ugali_dir) os.chdir(self.ugali_dir) - url = os.path.join(self.release,tarball) + url = os.path.join(self.release_url,tarball) print("downloading %s..."%url) if urlopen(url).getcode() >= 400: @@ -158,6 +160,12 @@ class CatalogCommand(TarballCommand): _tarball = 'ugali-catalogs.tar.gz' _dirname = 'catalogs' +class TestsCommand(TarballCommand): + """ Command for downloading catalog files """ + description = "install test data" + _tarball = 'ugali-test-data.tar.gz' + _dirname = 'testdata' + class IsochroneCommand(TarballCommand): """ Command for downloading isochrone files """ description = "install isochrone files" @@ -223,6 +231,7 @@ class install(_install): user_options = _install.user_options + [ ('isochrones',None,"install isochrone files (%s)"%ISOSIZE), ('catalogs',None,"install catalog files (%s)"%CATSIZE), + ('tests',None,"install test data (%s)"%TSTSIZE), ('ugali-dir=',None,"install file directory [default: %s]"%UGALIDIR), ] boolean_options = _install.boolean_options + ['isochrones','catalogs'] @@ -232,6 +241,7 @@ def initialize_options(self): self.ugali_dir = os.path.expandvars(UGALIDIR) self.isochrones = False self.catalogs = False + self.tests = False def run(self): # run superclass install @@ -246,7 +256,10 @@ def run(self): if self.catalogs: self.install_catalogs() - + + if self.tests: + self.install_tests() + def install_isochrones(self): """ Call to isochrone install command: @@ -266,10 +279,21 @@ def install_catalogs(self): cmd_obj.force = self.force if self.ugali_dir: cmd_obj.ugali_dir = self.ugali_dir self.run_command('catalogs') - + + def install_tests(self): + """ + Call to catalog install command: + http://stackoverflow.com/a/24353921/4075339 + """ + cmd_obj = self.distribution.get_command_obj('tests') + cmd_obj.force = self.force + if self.ugali_dir: cmd_obj.ugali_dir = self.ugali_dir + self.run_command('tests') + CMDCLASS = versioneer.get_cmdclass() CMDCLASS['isochrones'] = IsochroneCommand CMDCLASS['catalogs'] = CatalogCommand +CMDCLASS['tests'] = TestsCommand CMDCLASS['install'] = install setup( @@ -278,9 +302,11 @@ def install_catalogs(self): cmdclass=CMDCLASS, url=URL, author='Keith Bechtol & Alex Drlica-Wagner', - author_email='bechtol@kicp.uchicago.edu, kadrlica@fnal.gov', + author_email='bechtol@wisc.edu, kadrlica@fnal.gov', scripts = [], install_requires=[ + 'astropy', + 'matplotlib', 'numpy >= 1.9.0', 'scipy >= 0.14.0', 'healpy >= 1.6.0', @@ -288,7 +314,6 @@ def install_catalogs(self): 'emcee >= 2.1.0', 'corner >= 1.0.0', 'pyyaml >= 3.10', - # Add astropy, fitsio, matplotlib, ... ], packages=find_packages(), description=DESC, diff --git a/tests/config.yaml b/tests/config.yaml index b2278e3..13a5564 100644 --- a/tests/config.yaml +++ b/tests/config.yaml @@ -27,7 +27,7 @@ coords: catalog: # Color always defined as mag_1 - mag_2 - dirname: ./healpix # directory of catalog files + dirname: ${UGALIDIR}/healpix # directory of catalog files basename: "catalog_hpx%04i.fits" # catalog file basename format objid_field : COADD_OBJECT_ID # unique object identifier lon_field : RA # name of longitude field @@ -45,24 +45,19 @@ catalog: selection : "(np.abs(self.data['WAVG_SPREAD_MODEL_I']) < 0.003)" data: - dirname: ./raw + dirname: ${UGALIDIR}/raw script : ./ugali/preprocess/database.py survey : des release: y3a2 - density: ./density/density_hpx%04i.fits - footprint: ./maps/footprint_nside4096_equ.fits.gz + density: ${UGALIDIR}/density/density_hpx%04i.fits + footprint: ${UGALIDIR}/maps/footprint_nside4096_equ.fits.gz mask: - dirname : ./mask + dirname : ${UGALIDIR}/mask basename_1 : "maglim_g_hpx%04i.fits" basename_2 : "maglim_r_hpx%04i.fits" minimum_solid_angle: 0.1 # deg^2 -mangle: - dirname : ./maps - filename_1 : 'y3a2_g_o.4096_t.32768_maglim_EQU.fits.gz' - filename_2 : 'y3a2_r_o.4096_t.32768_maglim_EQU.fits.gz' - #ADW: Depricated in favor of 'binning' color: min : &cmin -0.5 diff --git a/ugali/utils/config.py b/ugali/utils/config.py index 535dac2..78d62de 100644 --- a/ugali/utils/config.py +++ b/ugali/utils/config.py @@ -172,86 +172,6 @@ def write(self, filename): raise Exception('Unrecognized config format: %s'%ext) writer.close() - def _createFilenames(self,pixels=None): - """ - Create a masked records array of all filenames for the given set of - pixels and store the existence of those files in the mask values. - - Examples: - f = getFilenames([1,2,3]) - # All possible catalog files - f['catalog'].data - # All existing catalog files - f['catalog'][~f.mask['catalog']] - # or - f['catalog'].compressed() - # All missing mask_1 files - f['mask_1'][f.mask['mask_1']] - # Pixels where all files exist - f['pix'][~f.mask['pix']] - - Parameters: - ----------- - pixels : If pixels is None, grab all pixels of 'nside_catalog'. - - Returns: - -------- - recarray : pixels and mask value - """ - nside_catalog = self['coords']['nside_catalog'] - - if pixels is not None: - pixels = [pixels] if np.isscalar(pixels) else pixels - else: - pixels = np.arange(hp.nside2npix(nside_catalog)) - - npix = len(pixels) - - catalog_dir = self['catalog']['dirname'] - catalog_base = self['catalog']['basename'] - - mask_dir = self['mask']['dirname'] - mask_base_1 = self['mask']['basename_1'] - mask_base_2 = self['mask']['basename_2'] - - data = np.ma.empty(npix,dtype=[('pix',int), ('catalog',object), - ('mask_1',object), ('mask_2',object)]) - mask = np.ma.empty(npix,dtype=[('pix',bool), ('catalog',bool), - ('mask_1',bool), ('mask_2',bool)]) - for ii,pix in enumerate(pixels): - if pix is None: - # DEPRECTATED: ADW 2018-06-17 - # This is not really being used anymore - raise ValueError('pix cannot be None') - catalog = os.path.join(catalog_dir,catalog_base) - mask_1 = os.path.join(mask_dir,mask_base_1) - mask_2 = os.path.join(mask_dir,mask_base_2) - else: - catalog = os.path.join(catalog_dir,catalog_base%pix) - mask_1 = os.path.join(mask_dir,mask_base_1%pix) - mask_2 = os.path.join(mask_dir,mask_base_2%pix) - data[ii]['pix'] = pix if pix is not None else -1 - data[ii]['catalog'] = catalog - data[ii]['mask_1'] = mask_1 - data[ii]['mask_2'] = mask_2 - - mask[ii]['catalog'] = not os.path.exists(catalog) - mask[ii]['mask_1'] = not os.path.exists(mask_1) - mask[ii]['mask_2'] = not os.path.exists(mask_2) - - for name in ['catalog','mask_1','mask_2']: - if np.all(mask[name]): logger.warn("All '%s' files masked"%name) - - # mask 'pix' if all files not present - mask['pix'] = mask['catalog'] | mask['mask_1'] | mask['mask_2'] - - if np.all(mask['pix']): logger.warn("All pixels masked") - - #return np.ma.mrecords.MaskedArray(data, mask, fill_value=[-1,None,None,None]) - #return np.ma.mrecords.MaskedArray(data, mask, fill_value=[-1,'','','']) - return np.ma.MaskedArray(data, mask, fill_value=[-1,'','','']) - - def _createFilenames(self): """ Create a masked records array of all filenames for the given set of @@ -269,11 +189,17 @@ def _createFilenames(self): npix = hp.nside2npix(nside_catalog) pixels = np.arange(npix) - catalog_dir = self['catalog']['dirname'] + catalog_dir = os.path.expandvars(self['catalog']['dirname']) + if not os.path.isdir(catalog_dir): + msg = "Directory does not exist: %s"%catalog_dir + raise IOError(msg) catalog_base = self['catalog']['basename'] catalog_path = os.path.join(catalog_dir,catalog_base) - mask_dir = self['mask']['dirname'] + mask_dir = os.path.expandvars(self['mask']['dirname']) + if not os.path.isdir(mask_dir): + msg = "Directory does not exist: %s"%mask_dir + raise IOError(msg) mask_base_1 = self['mask']['basename_1'] mask_base_2 = self['mask']['basename_2'] mask_path_1 = os.path.join(mask_dir,mask_base_1) From a8c69c8dd8249c1783d4555220a178092de51abb Mon Sep 17 00:00:00 2001 From: Alex Drlica-Wagner Date: Thu, 15 Jul 2021 23:26:12 -0500 Subject: [PATCH 15/19] github actions --- .github/workflows/python-package.yml | 7 ++++--- 1 file changed, 4 insertions(+), 3 deletions(-) diff --git a/.github/workflows/python-package.yml b/.github/workflows/python-package.yml index e68b720..7b67fda 100644 --- a/.github/workflows/python-package.yml +++ b/.github/workflows/python-package.yml @@ -18,8 +18,6 @@ jobs: max-parallel: 3 matrix: python-version: [2.7, 3.8] - env: - UGALIDIR: "$HOME/.ugali" steps: - uses: actions/checkout@v2 - name: Set up Python ${{ matrix.python-version }} @@ -34,10 +32,13 @@ jobs: conda info - name: Create conda environment run: | - conda create -q -n env python=${{ matrix.python-version }} numpy scipy matplotlib astropy healpy pyyaml emcee fitsio corner -c conda-forge -c kadrlica + conda create -y -q -n env python=${{ matrix.python-version }} numpy scipy matplotlib astropy healpy pyyaml emcee fitsio corner -c conda-forge -c kadrlica + # Add UGALIDIR to environment + conda env config vars set UGALIDIR=$HOME/.ugali -n env - name: Install package run: | source activate env + conda env config vars list #export UGALIDIR="$HOME/.ugali" python setup.py -q install --isochrones --catalogs --tests - name: Lint with flake8 From 451a2cd524d2d94c51bf55437513c6e421bc7e56 Mon Sep 17 00:00:00 2001 From: Alex Drlica-Wagner Date: Fri, 16 Jul 2021 08:56:02 -0500 Subject: [PATCH 16/19] tests --- .github/workflows/python-package.yml | 6 ------ setup.py | 4 ++-- tests/config.yaml | 8 ++++---- tests/test_config.py | 12 ++++++------ tests/test_observation.py | 7 ++++++- ugali/observation/mask.py | 2 +- ugali/utils/fileio.py | 2 +- 7 files changed, 20 insertions(+), 21 deletions(-) diff --git a/.github/workflows/python-package.yml b/.github/workflows/python-package.yml index 7b67fda..90d992b 100644 --- a/.github/workflows/python-package.yml +++ b/.github/workflows/python-package.yml @@ -39,7 +39,6 @@ jobs: run: | source activate env conda env config vars list - #export UGALIDIR="$HOME/.ugali" python setup.py -q install --isochrones --catalogs --tests - name: Lint with flake8 if: ${{ true }} @@ -48,11 +47,6 @@ jobs: conda install -q flake8 -c conda-forge # stop the build if there are Python syntax errors or undefined names flake8 . --count --select=E9,F63,F7,F82 --show-source --statistics - - name: Install test data - if: ${{ false }} - run: | - wget https://github.com/DarkEnergySurvey/ugali/releases/download/v1.7.0/ugali-test-data.tar.gz -O ugali-test-data.tar.gz - tar -xzf ugali-test-data.tar.gz - name: Test with nose run: | source activate env diff --git a/setup.py b/setup.py index b789caa..217cf1f 100644 --- a/setup.py +++ b/setup.py @@ -40,7 +40,7 @@ def find_packages(): Topic :: Scientific/Engineering :: Physics """ -RELEASE_URL = URL+'/releases/download/v1.7.0' +RELEASE_URL = URL+'/releases/download/v1.8.0' UGALIDIR = os.getenv("UGALIDIR","$HOME/.ugali") ISOSIZE = "~1MB" CATSIZE = "~20MB" @@ -64,7 +64,7 @@ def read(self, size): def progress_bar(count, block_size, total_size): block = 100*block_size/float(total_size) progress = count*block - if progress % 1 < 1.01*block: + if progress % 5 < 1.01*block: msg = '\r[{:51}] ({:d}%)'.format(int(progress//2)*'='+'>',int(progress)) sys.stdout.write(msg) sys.stdout.flush() diff --git a/tests/config.yaml b/tests/config.yaml index 13a5564..e8116e0 100644 --- a/tests/config.yaml +++ b/tests/config.yaml @@ -27,7 +27,7 @@ coords: catalog: # Color always defined as mag_1 - mag_2 - dirname: ${UGALIDIR}/healpix # directory of catalog files + dirname: ${UGALIDIR}/testdata/healpix # directory of catalog files basename: "catalog_hpx%04i.fits" # catalog file basename format objid_field : COADD_OBJECT_ID # unique object identifier lon_field : RA # name of longitude field @@ -49,11 +49,11 @@ data: script : ./ugali/preprocess/database.py survey : des release: y3a2 - density: ${UGALIDIR}/density/density_hpx%04i.fits - footprint: ${UGALIDIR}/maps/footprint_nside4096_equ.fits.gz + density: ./density/density_hpx%04i.fits + footprint: ./maps/footprint_nside4096_equ.fits.gz mask: - dirname : ${UGALIDIR}/mask + dirname : ${UGALIDIR}/testdata/mask basename_1 : "maglim_g_hpx%04i.fits" basename_2 : "maglim_r_hpx%04i.fits" minimum_solid_angle: 0.1 # deg^2 diff --git a/tests/test_config.py b/tests/test_config.py index 172b993..305ea28 100644 --- a/tests/test_config.py +++ b/tests/test_config.py @@ -15,12 +15,12 @@ def test_config(): np.testing.assert_equal(len(config.filenames),768) np.testing.assert_equal(config.filenames['pix'].compressed()[0],687) - np.testing.assert_equal(config.filenames['catalog'].compressed()[0], - './healpix/catalog_hpx0687.fits') - np.testing.assert_equal(config.filenames['mask_1'].compressed()[0], - './mask/maglim_g_hpx0687.fits') - np.testing.assert_equal(config.filenames['mask_2'].compressed()[0], - './mask/maglim_r_hpx0687.fits') + catfile = os.path.basename(config.filenames['catalog'].compressed()[0]) + np.testing.assert_equal(catfile,'catalog_hpx0687.fits') + maskfile = os.path.basename(config.filenames['mask_1'].compressed()[0]) + np.testing.assert_equal(maskfile,'maglim_g_hpx0687.fits') + maskfile = os.path.basename(config.filenames['mask_2'].compressed()[0]) + np.testing.assert_equal(maskfile,'maglim_r_hpx0687.fits') np.testing.assert_equal(config.likefile,'./scan/scan_%08i_%s.fits') np.testing.assert_equal(config.mergefile,'./scan/merged_scan.fits') diff --git a/tests/test_observation.py b/tests/test_observation.py index e97887b..c9e136a 100644 --- a/tests/test_observation.py +++ b/tests/test_observation.py @@ -44,8 +44,13 @@ def test_catalog(): """ Test ugali.observation.catalog """ import ugali.observation.roi import ugali.observation.catalog + import ugali.utils.config - filename='healpix/catalog_hpx0687.fits' + # Get catalog filename + config = ugali.utils.config.Config(CONFIG) + filename = config.filenames['catalog'].compressed()[0] + + # Create catalog from filename catalog = ugali.observation.catalog.Catalog(CONFIG,filenames=filename) roi = ugali.observation.roi.ROI(CONFIG, LON, LAT) diff --git a/ugali/observation/mask.py b/ugali/observation/mask.py index c59146e..611f5ff 100644 --- a/ugali/observation/mask.py +++ b/ugali/observation/mask.py @@ -592,7 +592,7 @@ def backgroundCMD(self, catalog, mode='cloud-in-cells', weights=None): index_color = np.arange(len(self.roi.centers_color)) # Add the cumulative leakage back into the last bin of the CMD leakage = (cmd_background * ~observable).sum(axis=0) - cmd_background[[index_mag,index_color]] += leakage + cmd_background[(index_mag,index_color)] += leakage # Zero out all non-observable bins cmd_background *= observable diff --git a/ugali/utils/fileio.py b/ugali/utils/fileio.py index 04708fa..b62bd6b 100644 --- a/ugali/utils/fileio.py +++ b/ugali/utils/fileio.py @@ -121,7 +121,7 @@ def load_header(kwargs): if not keys: keys = hdr.keys() return {k:hdr[k] for k in keys} except Exception as e: - logger.error("Failed to load file: %(filename)s"%kwargs) + logger.error("Failed to load header from file: %(filename)s"%kwargs) raise(e) def load_headers(filenames,multiproc=False,**kwargs): From 841fc9d58b8086b0014ca66f4c74bd764ed41bbb Mon Sep 17 00:00:00 2001 From: Alex Drlica-Wagner Date: Fri, 16 Jul 2021 13:21:06 -0500 Subject: [PATCH 17/19] README [skip ci] --- README.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/README.md b/README.md index 1966fda..740f056 100644 --- a/README.md +++ b/README.md @@ -50,7 +50,7 @@ The `ugali` source code is distributed with several auxiliary libraries for isoc ``` cd $UGALIDIR -wget https://github.com/DarkEnergySurvey/ugali/releases/download/v1.7.0/ugali-des-bressan2012.tar.gz +wget https://github.com/DarkEnergySurvey/ugali/releases/download/v1.8.0/ugali-des-bressan2012.tar.gz tar -xzf ugali-des-bressan2012.tar.gz ``` From 53a94fd05ec78ac07ba8112c6ff751b854d58cf0 Mon Sep 17 00:00:00 2001 From: Alex Drlica-Wagner Date: Fri, 16 Jul 2021 10:56:17 -0500 Subject: [PATCH 18/19] move deprecated color_lut to scratch --- ugali/{analysis => scratch}/color_lut.py | 0 1 file changed, 0 insertions(+), 0 deletions(-) rename ugali/{analysis => scratch}/color_lut.py (100%) diff --git a/ugali/analysis/color_lut.py b/ugali/scratch/color_lut.py similarity index 100% rename from ugali/analysis/color_lut.py rename to ugali/scratch/color_lut.py From 25777e42b459e689beab40a1dcda7293a4916c78 Mon Sep 17 00:00:00 2001 From: Alex Drlica-Wagner Date: Tue, 25 Jun 2024 18:42:45 -0500 Subject: [PATCH 19/19] add workflow_dispatch trigger --- .github/workflows/python-package.yml | 1 + 1 file changed, 1 insertion(+) diff --git a/.github/workflows/python-package.yml b/.github/workflows/python-package.yml index 90d992b..acef1c4 100644 --- a/.github/workflows/python-package.yml +++ b/.github/workflows/python-package.yml @@ -8,6 +8,7 @@ name: build on: push: pull_request: + workflow_dispatch: schedule: - cron: '0 0 1 * *'