diff --git a/.github/workflows/python-package.yml b/.github/workflows/python-package.yml index 43affee..b671fef 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 * *' 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 7cf713b..944cbf9 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/isochrone/model.py b/ugali/isochrone/model.py index 535eb3f..d673934 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'), @@ -385,7 +385,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 @@ -398,6 +398,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) @@ -440,10 +444,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 @@ -907,8 +916,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: ----------- @@ -1074,11 +1086,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) @@ -1110,7 +1125,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.") @@ -1143,7 +1162,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 diff --git a/ugali/isochrone/parsec.py b/ugali/isochrone/parsec.py index 5f784e3..af99f24 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'] diff --git a/ugali/preprocess/dotter.py b/ugali/preprocess/dotter.py deleted file mode 100755 index edb8829..0000000 --- a/ugali/preprocess/dotter.py +++ /dev/null @@ -1,352 +0,0 @@ -#!/usr/bin/env python -""" -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 fd9fa87..0000000 --- a/ugali/preprocess/padova.py +++ /dev/null @@ -1,364 +0,0 @@ -#!/usr/bin/env python -""" -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)) - 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)) + 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 fd6a364..a64623b 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)) @@ -165,18 +174,19 @@ def catsimSatellite(config, lon_centroid, lat_centroid, distance, stellar_mass, # 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 @@ -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): @@ -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 diff --git a/ugali/utils/healpix.py b/ugali/utils/healpix.py index e4e721b..b4696ce 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 + import fitsio data,hdr = fitsio.read(filename,header=True,ext=hdu) nside = int(hdr.get('NSIDE')) diff --git a/ugali/utils/projector.py b/ugali/utils/projector.py index 75e9b7b..9de825a 100644 --- a/ugali/utils/projector.py +++ b/ugali/utils/projector.py @@ -259,37 +259,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) @@ -308,11 +327,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): """ @@ -340,51 +360,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. @@ -402,12 +388,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. @@ -430,11 +417,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. @@ -475,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)) ############################################################ diff --git a/ugali/utils/stats.py b/ugali/utils/stats.py index aec8a15..505530a 100644 --- a/ugali/utils/stats.py +++ b/ugali/utils/stats.py @@ -29,6 +29,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)]] @@ -44,6 +54,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)