diff --git a/tests/test_results.py b/tests/test_results.py new file mode 100644 index 0000000..705b505 --- /dev/null +++ b/tests/test_results.py @@ -0,0 +1,30 @@ +#!/usr/bin/env python +""" +Generic python script. +""" +__author__ = "Alex Drlica-Wagner" + +import numpy as np +import ugali.analysis.results + +def test_surface_brightness(): + + # Draco values from McConnachie et al. 2012 + abs_mag = -8.8 + app_mag = 10.6 + distance = 76 #kpc + distance_modulus = 19.40 + a_half_arcmin = 10.0 # arcmin + a_physical_kpc = 0.221 # kpc + ellip = 0.31 + mu0 = 26.1 # mag/arcsec^2 + + # Convert to azimuthally average + r_half_arcmin = a_half_arcmin * np.sqrt(1-ellip) + r_physical_kpc = a_physical_kpc * np.sqrt(1-ellip) + + mu1 = ugali.analysis.results.surfaceBrightness(abs_mag, r_physical_kpc, distance) + mu2 = ugali.analysis.results.surfaceBrightness2(app_mag, r_half_arcmin) + + np.testing.assert_allclose(mu1,mu0,atol=2e-2) + np.testing.assert_allclose(mu2,mu0,atol=2e-2) diff --git a/ugali/analysis/farm.py b/ugali/analysis/farm.py index 25465c5..3073b9b 100755 --- a/ugali/analysis/farm.py +++ b/ugali/analysis/farm.py @@ -54,52 +54,6 @@ def command(self, outfile, configfile, pix): verbose='-v' if self.verbose else '') cmd = '%(script)s %(config)s %(outfile)s --hpx %(nside)i %(pix)i %(verbose)s'%params return cmd - - # ADW: Should probably be in a utility - def footprint(self, nside=None): - """ - UNTESTED. - Should return a boolean array representing the pixels in the footprint. - """ - if nside is None: - nside = self.nside_pixel - elif nside < self.nside_catalog: - raise Exception('Requested nside=%i is less than catalog_nside'%nside) - elif nside > self.nside_pixel: - raise Exception('Requested nside=%i is greater than pixel_nside'%nside) - pix = np.arange(hp.nside2npix(nside), dtype=int) - map = self.inFootprint(pix,nside) - return map - - # ADW: Should probably be in a utility - def inFootprint(self, pixels, nside=None): - """ - Open each valid filename for the set of pixels and determine the set - of subpixels with valid data. - """ - if np.isscalar(pixels): pixels = np.array([pixels]) - if nside is None: nside = self.nside_likelihood - - inside = np.zeros( len(pixels), dtype='bool') - if not self.nside_catalog: - catalog_pix = [0] - else: - catalog_pix = superpixel(pixels,nside,self.nside_catalog) - catalog_pix = np.intersect1d(catalog_pix,self.catalog_pixels) - - for filenames in self.filenames[catalog_pix]: - # ADW: Need to replace with healpix functions... - #logger.debug("Loading %s"%filenames['mask_1']) - #subpix_1,val_1 = ugali.utils.skymap.readSparseHealpixMap(filenames['mask_1'],'MAGLIM',construct_map=False) - _n,subpix_1,val_1 = read_partial_map(filenames['mask_1'],'MAGLIM',fullsky=False) - #logger.debug("Loading %s"%filenames['mask_2']) - #subpix_2,val_2 = ugali.utils.skymap.readSparseHealpixMap(filenames['mask_2'],'MAGLIM',construct_map=False) - _n,subpix_2,val_2 = read_partial_map(filenames['mask_2'],'MAGLIM',fullsky=False) - subpix = np.intersect1d(subpix_1,subpix_2) - superpix = np.unique(superpixel(subpix,self.nside_pixel,nside)) - inside |= np.in1d(pixels, superpix) - - return inside def submit_all(self, coords=None, queue=None, debug=False): """ @@ -114,18 +68,9 @@ def submit_all(self, coords=None, queue=None, debug=False): if coords is None: pixels = np.arange(hp.nside2npix(self.nside_likelihood)) else: - coords = np.asarray(coords) - if coords.ndim == 1: - coords = np.array([coords]) - if coords.shape[1] == 2: - lon,lat = coords.T - radius = np.zeros(len(lon)) - elif coords.shape[1] == 3: - lon,lat,radius = coords.T - else: - raise Exception("Unrecognized coords shape:"+str(coords.shape)) + lon,lat,radius = coords['lon'],coords['lat'],coords['radius'] - #ADW: targets is still in glon,glat + #ADW: coords are always parsed in GAL, so convert to CEL if necessary if self.config['coords']['coordsys'].lower() == 'cel': lon,lat = gal2cel(lon,lat) @@ -158,11 +103,11 @@ def submit(self, pixels, queue=None, debug=False, configfile=None): """ # For backwards compatibility batch = self.config['scan'].get('batch',self.config['batch']) - queue = batch['cluster'] if queue is None else queue + queue = batch.get('default','medium') if queue is None else queue # Need to develop some way to take command line arguments... - self.batch = ugali.utils.batch.batchFactory(queue,**batch['opts']) - self.batch.max_jobs = batch.get('max_jobs',200) + self.batch = ugali.utils.batch.batchFactory(queue,**batch.get(queue,{})) + self.batch.max_jobs = self.config['scan'].get('max_jobs',200) if np.isscalar(pixels): pixels = np.array([pixels]) @@ -178,7 +123,7 @@ def submit(self, pixels, queue=None, debug=False, configfile=None): lon,lat = pix2ang(self.nside_likelihood,pixels) commands = [] - chunk = batch['chunk'] + chunk = self.config['scan'].get('chunk',25) istart = 0 logger.info('=== Submit Likelihood ===') for ii,pix in enumerate(pixels): @@ -189,7 +134,7 @@ def submit(self, pixels, queue=None, debug=False, configfile=None): # Create outfile name outfile = self.config.likefile%(pix,self.config['coords']['coordsys'].lower()) outbase = os.path.basename(outfile) - jobname = batch['jobname'] + jobname = batch.get('jobname','ugali') # Submission command sub = not os.path.exists(outfile) @@ -226,21 +171,48 @@ def submit(self, pixels, queue=None, debug=False, configfile=None): time.sleep(0.5) def write_script(self, filename, commands): + """ Write a batch submission script. + + Parameters + ---------- + filename : filename of batch script + commands : list of commands to execute + + Returns + ------- + None + """ info = 'echo "{0:=^60}";\n' hline = info.format("") newline = 'echo;\n' + shebang = "#!/usr/bin/env bash" + # Limit the memory based on SLAC 4 GB per node (defined in KB) + # Careful, shell arithmetic is weird. + memory_limit = """ +if [ -n "$LSB_CG_MEMLIMIT" ] & [ -n "$LSB_HOSTS" ]; then + mlimit=$(( $(wc -w <<< $LSB_HOSTS) * $LSB_CG_MEMLIMIT/1024 * 9/10 )) + ulimit -v ${mlimit}; ulimit -H -v ${mlimit}; +fi +""" + memory_usage=r"""free -m | awk 'NR==2{printf "Memory Usage: %.2f/%.2fGB (%.2f%%)\n",$3/1024,$2/1024,$3*100/$2}';""" + memory_usage=r"""ps -U $USER --no-headers -o rss | awk '{sum+=$1} END {print "Memory Usage: " int(sum/1024**2) "GB"}'""" istart, iend = commands[0][0], commands[-1][0] script = open(filename,'w') + script.write(shebang) + #script.write(memory_limit) script.write(hline) - script.write(info.format('Submit Pixels %i to %i'%(istart,iend))) + script.write(info.format('Submit Jobs %i to %i'%(istart,iend))) script.write(hline) script.write(newline) script.write('status=0;\n') for i,cmd,lon,lat,sub in commands: - script.write(info.format('Pixel %i: (%.2f, %.2f)'%(i,lon,lat))) - if sub: script.write('%s; [ $? -ne 0 ] && status=1;\n'%cmd) - else: script.write('echo "%s";\n'%self.skip) + script.write(info.format('Job %i: (%.2f, %.2f)'%(i,lon,lat))) + if sub: + script.write(memory_usage+'\n') + script.write('%s; [ $? -ne 0 ] && status=1;\n'%cmd) + else: + script.write('echo "%s";\n'%self.skip) script.write(hline) script.write(newline) script.write('exit $status;\n') diff --git a/ugali/analysis/kernel.py b/ugali/analysis/kernel.py index c9b13d1..02f8f09 100644 --- a/ugali/analysis/kernel.py +++ b/ugali/analysis/kernel.py @@ -40,29 +40,42 @@ def __init__(self, proj='ait', **kwargs): super(Kernel,self).__init__(**kwargs) def __call__(self, lon, lat): + """ + Return the value of the pdf at a give location. + + Parameters + ---------- + lon : longitude (deg) + lat : latitude (deg) + + Returns + ------- + pdf : normalized, truncated pdf + """ return self.pdf(lon, lat) @abstractmethod def _kernel(self, r): - # unnormalized, untruncated kernel + """Unnormalized, untruncated kernel""" pass def _pdf(self, radius): - # unnormalized, truncated kernel + """Unnormalized, truncated kernel""" return np.where(radius<=self.edge, self._kernel(radius), 0.) @abstractmethod def pdf(self, lon, lat): - # normalized, truncated pdf + """Normalized, truncated pdf""" pass @property def norm(self): - # Numerically integrate the pdf + """Normalization from the integated pdf""" return 1./self.integrate() @property def projector(self): + """Projector used to transform to local sky coordinates.""" if self.proj is None or self.proj.lower()=='none': return None else: @@ -133,12 +146,14 @@ class EllipticalKernel(Kernel): _params = odict( list(Kernel._params.items()) + [ + # This is the semi-major axis ('extension', Parameter(0.1, [0.0001,0.5]) ), - ('ellipticity', Parameter(0.0, [0.0, 0.99]) ), # Default 0 for RadialKernel - ('position_angle',Parameter(0.0, [0.0, 180.0]) ), # Default 0 for RadialKernel - # This is the PA *WEST* of North. + # This is e = 1 - b/a (0 for RadialKernel) + ('ellipticity', Parameter(0.0, [0.0, 0.99]) ), + # This is the PA *WEST* of North (0 for RadialKernel) # to get the conventional PA EAST of North take 90-PA # Would it be better to have bounds [-90,90]? + ('position_angle',Parameter(0.0, [0.0, 180.0]) ), # ]) _mapping = odict([ ('e','ellipticity'), @@ -183,19 +198,36 @@ def sample_radius(self, n): """ Sample the radial distribution (deg) from the 2D stellar density. Output is elliptical radius in true projected coordinates. + + Parameters + ---------- + n : number of stars to sample + + Returns + ------- + radius : distance from centroid (deg) """ + size = int(n) edge = self.edge if self.edge<20*self.extension else 20*self.extension - radius = np.linspace(0, edge, 1.e5) + radius = np.linspace(0, edge, int(10**5)) pdf = self._pdf(radius) * np.sin(np.radians(radius)) cdf = np.cumsum(pdf) cdf /= cdf[-1] fn = scipy.interpolate.interp1d(cdf, list(range(0, len(cdf)))) - index = np.floor(fn(np.random.uniform(size=n))).astype(int) + index = np.floor(fn(np.random.uniform(size=size))).astype(int) return radius[index] def sample_lonlat(self, n): """ Sample 2D distribution of points in lon, lat + + Parameters + ---------- + n : number of stars to sample + + Returns + ------- + lon,lat : longitude and latitude (deg) """ # From http://en.wikipedia.org/wiki/Ellipse#General_parametric_form # However, Martin et al. (2009) use PA theta "from North to East" @@ -374,7 +406,7 @@ def _kernel(self, radius): return ((1./np.sqrt(1.+(radius/self.r_c)**2))-(1./np.sqrt(1.+(self.r_t/self.r_c)**2)))**2 def _cache(self, name=None): - if name in ['extension','ellipticity','truncate']: + if name in [None,'extension','ellipticity','truncate']: self._norm = 1./self.integrate() else: return diff --git a/ugali/analysis/model.py b/ugali/analysis/model.py index e47865c..31f426f 100644 --- a/ugali/analysis/model.py +++ b/ugali/analysis/model.py @@ -25,9 +25,9 @@ def asscalar(a): # Do we want to check that the value is numeric? #if isinstance(value, (int, long, float)): return value try: - return np.asscalar(a) + return a.item() except AttributeError as e: - return np.asscalar(np.asarray(a)) + return np.asarray(a).item() class Model(object): diff --git a/ugali/analysis/results.py b/ugali/analysis/results.py index 27e29b5..c98796c 100644 --- a/ugali/analysis/results.py +++ b/ugali/analysis/results.py @@ -23,13 +23,15 @@ from ugali.utils.config import Config from ugali.utils.logger import logger +_alpha = 0.32 + class Results(object): """ Calculate results from a MCMC chain. """ def __init__(self, config, loglike, samples=None): self.config = Config(config) - self.alpha = self.config['results'].get('alpha',0.10) + self.alpha = self.config['results'].get('alpha',_alpha) self.nwalkers = self.config['mcmc'].get('nwalkers',100) self.nburn = self.config['results'].get('nburn',10) self.coordsys = self.config['coords']['coordsys'].lower() @@ -54,9 +56,22 @@ def get_mle(self): return mle - def estimate(self,param,burn=None,clip=10.0,alpha=0.32): - """ Estimate parameter value and uncertainties """ + def estimate(self,param,burn=None,clip=10.0,alpha=_alpha): + """ Estimate parameter value and uncertainties + + Parameters + ---------- + param : parameter of interest + burn : number of burn in samples + clip : sigma clipping + alpha : confidence interval + + Returns + ------- + [mle, [lo,hi]] : value and interval + """ # FIXME: Need to add age and metallicity to composite isochrone params (currently properties) + if alpha is None: alpha = self.alpha if param not in list(self.samples.names) + list(self.source.params) + ['age','metallicity']: msg = 'Unrecognized parameter: %s'%param raise KeyError(msg) @@ -89,18 +104,20 @@ def estimate(self,param,burn=None,clip=10.0,alpha=0.32): ### msg = "Unrecognized parameter: %s"%param ### raise KeyError(msg) - def estimate_params(self,burn=None,clip=10.0,alpha=0.32): + def estimate_params(self,burn=None,clip=10.0,alpha=None): """ Estimate all source parameters """ + if alpha is None: alpha = self.alpha mle = self.get_mle() out = odict() for param in mle.keys(): out[param] = self.estimate(param,burn=burn,clip=clip,alpha=alpha) return out - def estimate_position_angle(self,param='position_angle',burn=None,clip=10.0,alpha=0.32): + def estimate_position_angle(self,param='position_angle',burn=None,clip=10.0,alpha=None): """ Estimate the position angle from the posterior dealing with periodicity. """ + if alpha is None: alpha = self.alpha # Transform so peak in the middle of the distribution pa = self.samples.get(param,burn=burn,clip=clip) peak = ugali.utils.stats.kde_peak(pa) @@ -141,17 +158,31 @@ def get_results(self,**kwargs): # Extra parameters from the MCMC chain logger.debug('Estimating auxiliary parameters...') + logger.debug("alpha = %.2f"%kwargs['alpha']) + results['alpha'] = kwargs['alpha'] try: results['ra'] = self.estimate('ra',**kwargs) results['dec'] = self.estimate('dec',**kwargs) + results['glon'] = self.estimate('glon',**kwargs) + results['glat'] = self.estimate('glat',**kwargs) except KeyError: - logger.warn("Didn't find 'ra' or 'dec'") - ra,dec = gal2cel(results['lon'][0],results['lat'][0]) - results['ra'] = ugali.utils.stats.interval(ra) - results['dec'] = ugali.utils.stats.interval(dec) - - ra,dec = results['ra'][0],results['dec'][0] - glon,glat = lon,lat = results['lon'][0],results['lat'][0] + logger.warn("Didn't find 'ra' or 'dec' in Samples...") + if self.coordsys == 'gal': + results['glon'] = results['lon'] + results['glat'] = results['lat'] + ra,dec = gal2cel(results['lon'][0],results['lat'][0]) + results['ra'] = ugali.utils.stats.interval(ra) + results['dec'] = ugali.utils.stats.interval(dec) + else: + results['ra'] = results['lon'] + results['dec'] = results['lat'] + glon,glat = cel2gal(results['lon'][0],results['lat'][0]) + results['glon'] = ugali.utils.stats.interval(glon) + results['glat'] = ugali.utils.stats.interval(glat) + + lon,lat = results['lon'][0],results['lat'][0] + ra,dec = results['ra'][0],results['dec'][0] + glon,glat = results['glon'][0],results['glat'][0] results.update(gal=[float(glon),float(glat)]) results.update(cel=[float(ra),float(dec)]) @@ -165,14 +196,6 @@ def get_results(self,**kwargs): ts = 2*self.loglike.value(**params) results['ts'] = ugali.utils.stats.interval(ts,np.nan,np.nan) - #lon,lat = estimate['lon'][0],estimate['lat'][0] - # - #results.update(gal=[float(lon),float(lat)]) - #ra,dec = gal2cel(lon,lat) - #results.update(cel=[float(ra),float(dec)]) - #results['ra'] = ugali.utils.stats.interval(ra,np.nan,np.nan) - #results['dec'] = ugali.utils.stats.interval(dec,np.nan,np.nan) - # Celestial position angle # Break ambiguity in direction with '% 180.' pa,pa_err = results['position_angle'] @@ -256,7 +279,7 @@ def get_results(self,**kwargs): # ADW: WARNING this is very fragile. # Also, this is not quite right, should cut on the CMD available space kwargs = dict(richness=rich,mag_bright=16., mag_faint=23., - n_trials=5000,alpha=self.alpha, seed=0) + n_trials=5000,alpha=kwargs['alpha'], seed=0) martin = self.config['results'].get('martin') if martin: logger.info("Calculating Martin magnitude...") @@ -319,22 +342,65 @@ def write(self,filename): out.close() def surfaceBrightness(abs_mag, r_physical, distance): - """ - Compute the average surface brightness [mag arcsec^-2] within the half-light radius + """Compute the average surface brightness [mag/arcsec^2] within the + azimuthally averaged physical half-light radius. + + The azimuthally averaged physical half-light radius is used + because the surface brightness is calculated as a function of + area. The area of an ellipse is: - abs_mag = absolute magnitude [mag] - r_physical = half-light radius [kpc] - distance = [kpc] + A = pi * a*b = pi * a**2 * (1-e) = pi * r**2 - The factor 2 in the c_v equation below account for half the luminosity - within the half-light radius. The 3600.**2 is conversion from deg^2 to arcsec^2 + The factor 2 in the c_v equation below accounts for the fact that + we are taking half the luminosity within the half-light + radius. The 3600.**2 is conversion from deg^2 to arcsec^2 c_v = 2.5 * np.log10(2.) + 2.5 * np.log10(np.pi * 3600.**2) = 19.78 + + TODO: Distance could be optional + + Parameters + ---------- + abs_mag : absolute V-band magnitude [mag] + r_physical : azimuthally averaged physical half-light radius [kpc] + distance : heliocentric distance [kpc] + + Returns + ------- + mu : surface brightness [mag/arcsec^2] """ + c_v = 19.78 # conversion to mag/arcsec^2 r_angle = np.degrees(np.arctan(r_physical / distance)) - c_v = 19.78 # mag/arcsec^2 return abs_mag + dist2mod(distance) + c_v + 2.5 * np.log10(r_angle**2) +def surfaceBrightness2(app_mag, r_half_arcmin): + """Compute the average surface brightness [mag/arcsec^2] within the + azimuthally averaged angular half-light radius. + + The azimuthally averaged half-light radius is used + because the surface brightness is calculated as a function of + area. The area of an ellipse is: + + A = pi * a*b = pi * a**2 * (1-e) = pi * r**2 + + The factor 2.5*log10(2) accounts for the fact that half the + luminosity is within the half-light radius. + + Parameters + ---------- + app_mag : apparent V-band magnitude [mag] + r_half : azimuthally averaged half-light radius [arcmin] + + Returns + ------- + mu : surface brightness [mag/arcsec^2] + + """ + r_half_arcsec = r_half_arcmin * 60 # arcmin to arcsec + mu = app_mag + 2.5*np.log10(2.) + 2.5*np.log10(np.pi * r_half_arcsec**2) + return mu + + def createResults(config,srcfile,section='source',samples=None): """ Create an MCMC instance """ source = ugali.analysis.source.Source() diff --git a/ugali/analysis/scan.py b/ugali/analysis/scan.py index 7272505..93c1ea3 100755 --- a/ugali/analysis/scan.py +++ b/ugali/analysis/scan.py @@ -17,6 +17,7 @@ from ugali.analysis.source import Source from ugali.utils.parabola import Parabola +from ugali.utils.batch import LSF from ugali.utils.config import Config from ugali.utils.logger import logger from ugali.utils.healpix import superpixel, subpixel, pix2ang, ang2pix @@ -88,105 +89,82 @@ def __init__(self, config, loglike): # What it should be... self.config = Config(config) self.loglike = loglike + self.source = self.loglike.source self.roi = self.loglike.roi self.mask = self.loglike.mask logger.info(str(self.loglike)) - self.stellar_mass_conversion = self.loglike.source.stellar_mass() + self.stellar_mass_conversion = self.source.stellar_mass() self.distance_modulus_array = np.asarray(self.config['scan']['distance_modulus_array']) + self.extension_array = np.asarray(self.config['scan'].get('extension_array',[self.source.extension])) - def precompute(self, distance_modulus_array=None): - """ - DEPRECATED: ADW 20170627 - - Precompute color probabilities for background ('u_background') - and signal ('u_color') for each star in catalog. Precompute - observable fraction in each ROI pixel. # Precompute still - operates over the full ROI, not just the likelihood region - Parameters: - ----------- - distance_modulus_array : Array of distance moduli - - Returns: - -------- - None + def search(self, coords=None, distance_modulus=None, extension=None, tolerance=1.e-2): """ - msg = "'%s.precompute': ADW 2017-09-20"%self.__class__.__name__ - DeprecationWarning(msg) - - if distance_modulus_array is not None: - self.distance_modulus_array = distance_modulus_array - else: - self.distance_modulus_array = sel - - # Observable fraction for each pixel - self.u_color_array = [[]] * len(self.distance_modulus_array) - self.observable_fraction_sparse_array = [[]] * len(self.distance_modulus_array) - - logger.info('Looping over distance moduli in precompute ...') - for ii, distance_modulus in enumerate(self.distance_modulus_array): - logger.info(' (%i/%i) Distance Modulus = %.2f ...'%(ii+1, len(self.distance_modulus_array), distance_modulus)) - - self.u_color_array[ii] = False - if self.config['scan']['color_lut_infile'] is not None: - DeprecationWarning("'color_lut' is deprecated") - logger.info(' Precomputing signal color from %s'%(self.config['scan']['color_lut_infile'])) - self.u_color_array[ii] = ugali.analysis.color_lut.readColorLUT(self.config['scan']['color_lut_infile'], - distance_modulus, - self.loglike.catalog.mag_1, - self.loglike.catalog.mag_2, - self.loglike.catalog.mag_err_1, - self.loglike.catalog.mag_err_2) - if not np.any(self.u_color_array[ii]): - logger.info(' Precomputing signal color on the fly...') - self.u_color_array[ii] = self.loglike.calc_signal_color(distance_modulus) - - # Calculate over all pixels in ROI - self.observable_fraction_sparse_array[ii] = self.loglike.calc_observable_fraction(distance_modulus) - - self.u_color_array = np.array(self.u_color_array) - - def search(self, coords=None, distance_modulus=None, tolerance=1.e-2): - """ - Organize a grid search over ROI target pixels and distance moduli in distance_modulus_array - coords: (lon,lat) - distance_modulus: scalar + Organize a grid search over ROI target pixels, distance + moduli, and extensions. If coords, distance_modulus, or + extension is specified, then the nearest value in the + predefined scan grid is used. ***This may be different than + the input value.** + + Parameters + ---------- + coords : (float,float) + coordinate to search (matched to nearest scan value) + distance_modulus : float + distance modulus to search (matched to nearest scan value) + extension : float + extension to search (matched to nearest scan value) + tolerance : float + tolerance on richness maximization + + Returns + ------- + None """ nmoduli = len(self.distance_modulus_array) - npixels = len(self.roi.pixels_target) - self.log_likelihood_sparse_array = np.zeros([nmoduli, npixels]) - self.richness_sparse_array = np.zeros([nmoduli, npixels]) - self.richness_lower_sparse_array = np.zeros([nmoduli, npixels]) - self.richness_upper_sparse_array = np.zeros([nmoduli, npixels]) - self.richness_upper_limit_sparse_array = np.zeros([nmoduli, npixels]) - self.stellar_mass_sparse_array = np.zeros([nmoduli, npixels]) - self.fraction_observable_sparse_array = np.zeros([nmoduli, npixels]) + npixels = len(self.roi.pixels_target) + self.loglike_array = np.zeros([nmoduli,npixels],dtype='f4') + self.richness_array = np.zeros([nmoduli,npixels],dtype='f4') + 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') # Specific pixel/distance_modulus - coord_idx, distance_modulus_idx = None, None + coord_idx, distance_modulus_idx, extension_idx = None, None, None if coords is not None: # Match to nearest grid coordinate index coord_idx = self.roi.indexTarget(coords[0],coords[1]) if distance_modulus is not None: # Match to nearest distance modulus index distance_modulus_idx=np.fabs(self.distance_modulus_array-distance_modulus).argmin() + if extension is not None: + # Match to nearest extension + extension_idx=np.fabs(self.extension_array-extension).argmin() lon, lat = self.roi.pixels_target.lon, self.roi.pixels_target.lat logger.info('Looping over distance moduli in grid search ...') for ii, distance_modulus in enumerate(self.distance_modulus_array): - - # Specific pixel + # Specific distance if distance_modulus_idx is not None: if ii != distance_modulus_idx: continue logger.info(' (%-2i/%i) Distance Modulus=%.1f ...'%(ii+1,nmoduli,distance_modulus)) + # No objects, continue + if len(self.loglike.catalog) == 0: + logger.warn("No catalog objects") + continue + # Set distance_modulus once to save time self.loglike.set_params(distance_modulus=distance_modulus) - + # Loop over pixels for jj in range(0, npixels): # Specific pixel if coord_idx is not None: @@ -194,84 +172,90 @@ def search(self, coords=None, distance_modulus=None, tolerance=1.e-2): # Set kernel location self.loglike.set_params(lon=lon[jj],lat=lat[jj]) - # Doesn't re-sync distance_modulus each time - self.loglike.sync_params() - - args = (jj+1, npixels, self.loglike.source.lon, self.loglike.source.lat) - msg = ' (%-3i/%i) Candidate at (%.2f, %.2f) ... '%(args) - - self.log_likelihood_sparse_array[ii][jj], self.richness_sparse_array[ii][jj], parabola = self.loglike.fit_richness() - self.stellar_mass_sparse_array[ii][jj] = self.stellar_mass_conversion * self.richness_sparse_array[ii][jj] - self.fraction_observable_sparse_array[ii][jj] = self.loglike.f + + loglike = 0 + # Loop over extensions + for kk,ext in enumerate(self.extension_array): + # Specific extension + if extension_idx is not None: + if kk != extension_idx: continue + + # Set extension + self.loglike.set_params(extension=ext) + + # Doesn't re-sync distance_modulus each time + self.loglike.sync_params() + + # Maximize the likelihood with respect to richness + loglike,rich,p = self.loglike.fit_richness() + + if loglike < self.loglike_array[ii][jj]: + # No loglike increase, continue + continue + + self.loglike_array[ii][jj],self.richness_array[ii][jj], parabola = loglike,rich,p + self.stellar_mass_array[ii][jj] = self.stellar_mass_conversion*self.richness_array[ii][jj] + self.fraction_observable_array[ii][jj] = self.loglike.f + self.extension_fit_array[ii][jj] = self.source.extension + + # ADW: Careful, we are leaving the extension at the + # last value in the array, not at the maximum... + + # Debug output + args = (jj+1, npixels, lon[jj], lat[jj], + 2.*self.loglike_array[ii][jj], + self.stellar_mass_array[ii][jj], + self.fraction_observable_array[ii][jj], + self.extension_fit_array[ii][jj] + ) + msg = ' (%-3i/%i) Candidate at (%.2f, %.2f) ... ' + msg += 'TS=%.1f, Mstar=%.2g, ObsFrac=%.2g, Ext=%.2g' + logger.debug(msg%args) + + """ + # This is debugging output if self.config['scan']['full_pdf']: - #n_pdf_points = 100 - #richness_range = parabola.profileUpperLimit(delta=25.) - self.richness_sparse_array[ii][jj] - #richness = np.linspace(max(0., self.richness_sparse_array[ii][jj] - richness_range), - # self.richness_sparse_array[ii][jj] + richness_range, - # n_pdf_points) - #if richness[0] > 0.: - # richness = np.insert(richness, 0, 0.) - # n_pdf_points += 1 - # - #log_likelihood = np.zeros(n_pdf_points) - #for kk in range(0, n_pdf_points): - # log_likelihood[kk] = self.loglike.value(richness=richness[kk]) - #parabola = ugali.utils.parabola.Parabola(richness, 2.*log_likelihood) - #self.richness_lower_sparse_array[ii][jj], self.richness_upper_sparse_array[ii][jj] = parabola.confidenceInterval(0.6827) - self.richness_lower_sparse_array[ii][jj], self.richness_upper_sparse_array[ii][jj] = self.loglike.richness_interval(0.6827) + DeprecationWarning("'full_pdf' is deprecated.") + self.richness_lower_array[ii][jj], self.richness_upper_array[ii][jj] = self.loglike.richness_interval(0.6827) - self.richness_upper_limit_sparse_array[ii][jj] = parabola.bayesianUpperLimit(0.95) + self.richness_ulimit_array[ii][jj] = parabola.bayesianUpperLimit(0.95) args = ( - 2. * self.log_likelihood_sparse_array[ii][jj], - self.stellar_mass_conversion*self.richness_sparse_array[ii][jj], - self.stellar_mass_conversion*self.richness_lower_sparse_array[ii][jj], - self.stellar_mass_conversion*self.richness_upper_sparse_array[ii][jj], - self.stellar_mass_conversion*self.richness_upper_limit_sparse_array[ii][jj] + 2. * self.loglike_array[ii][jj], + self.stellar_mass_conversion*self.richness_array[ii][jj], + self.stellar_mass_conversion*self.richness_lower_array[ii][jj], + self.stellar_mass_conversion*self.richness_upper_array[ii][jj], + self.stellar_mass_conversion*self.richness_ulimit_array[ii][jj] ) - msg += 'TS=%.1f, Stellar Mass=%.1f (%.1f -- %.1f @ 0.68 CL, < %.1f @ 0.95 CL)'%(args) - else: - args = ( - 2. * self.log_likelihood_sparse_array[ii][jj], - self.stellar_mass_conversion * self.richness_sparse_array[ii][jj], - self.fraction_observable_sparse_array[ii][jj] - ) - msg += 'TS=%.1f, Stellar Mass=%.1f, Fraction=%.2g'%(args) - logger.debug(msg) + msg = 'TS=%.1f, Stellar Mass=%.1f (%.1f -- %.1f @ 0.68 CL, < %.1f @ 0.95 CL)'%(args) + logger.debug(msg) + """ - #if coords is not None and distance_modulus is not None: - # results = [self.richness_sparse_array[ii][jj], - # self.log_likelihood_sparse_array[ii][jj], - # self.richness_lower_sparse_array[ii][jj], - # self.richness_upper_sparse_array[ii][jj], - # self.richness_upper_limit_sparse_array[ii][jj], - # richness, log_likelihood, self.loglike.p, self.loglike.f] - # return results - - jj_max = self.log_likelihood_sparse_array[ii].argmax() + jj_max = self.loglike_array[ii].argmax() args = ( jj_max+1, npixels, lon[jj_max], lat[jj_max], - 2. * self.log_likelihood_sparse_array[ii][jj_max], - self.stellar_mass_conversion * self.richness_sparse_array[ii][jj_max] + 2. * self.loglike_array[ii][jj_max], + self.stellar_mass_conversion * self.richness_array[ii][jj_max], + self.extension_fit_array[ii][jj_max] ) - msg = ' (%-3i/%i) Maximum at (%.2f, %.2f) ... TS=%.1f, Stellar Mass=%.1f'%(args) + msg = ' (%-3i/%i) Max at (%.2f, %.2f) : TS=%.1f, Mstar=%.2g, Ext=%.2f'%(args) logger.info(msg) def mle(self): - a = self.log_likelihood_sparse_array + a = self.loglike_array j,k = np.unravel_index(a.argmax(),a.shape) mle = odict() - mle['richness'] = self.richness_sparse_array[j][k] + mle['richness'] = self.richness_array[j][k] mle['lon'] = self.roi.pixels_target.lon[k] mle['lat'] = self.roi.pixels_target.lat[k] mle['distance_modulus'] = self.distance_modulus_array[j] - mle['extension'] = float(self.loglike.source.extension) - mle['ellipticity'] = float(self.loglike.source.ellipticity) - mle['position_angle'] = float(self.loglike.source.position_angle) + mle['extension'] = float(self.source.extension) + mle['ellipticity'] = float(self.source.ellipticity) + mle['position_angle'] = float(self.source.position_angle) # ADW: FIXME! try: - mle['age'] = np.mean(self.loglike.source.age) - mle['metallicity'] = np.mean(self.loglike.source.metallicity) + mle['age'] = np.mean(self.source.age) + mle['metallicity'] = np.mean(self.source.metallicity) except AttributeError: mle['age'] = np.nan mle['metallicity'] = np.nan @@ -288,7 +272,7 @@ def err(self): err.update([(k,np.nan*np.ones(2)) for k in err.keys()]) # Find the maximum likelihood - a = self.log_likelihood_sparse_array + a = self.loglike_array j,k = np.unravel_index(a.argmax(),a.shape) self.loglike.set_params(distance_modulus=self.distance_modulus_array[j], @@ -331,17 +315,17 @@ def write(self, outfile): data['PIXEL']=self.roi.pixels_target # Full data output (too large for survey) if self.config['scan']['full_pdf']: - data['LOG_LIKELIHOOD']=self.log_likelihood_sparse_array.T - data['RICHNESS']=self.richness_sparse_array.T - data['RICHNESS_LOWER']=self.richness_lower_sparse_array.T - data['RICHNESS_UPPER']=self.richness_upper_sparse_array.T - data['RICHNESS_LIMIT']=self.richness_upper_limit_sparse_array.T - #data['STELLAR_MASS']=self.stellar_mass_sparse_array.T - data['FRACTION_OBSERVABLE']=self.fraction_observable_sparse_array.T + data['LOG_LIKELIHOOD']=self.loglike_array.T + data['RICHNESS']=self.richness_array.T + data['RICHNESS_LOWER']=self.richness_lower_array.T + data['RICHNESS_UPPER']=self.richness_upper_array.T + data['RICHNESS_LIMIT']=self.richness_ulimit_array.T + #data['STELLAR_MASS']=self.stellar_mass_array.T + data['FRACTION_OBSERVABLE']=self.fraction_observable_array.T else: - data['LOG_LIKELIHOOD']=self.log_likelihood_sparse_array.T - data['RICHNESS']=self.richness_sparse_array.T - data['FRACTION_OBSERVABLE']=self.fraction_observable_sparse_array.T + data['LOG_LIKELIHOOD']=self.loglike_array.T + data['RICHNESS']=self.richness_array.T + data['FRACTION_OBSERVABLE']=self.fraction_observable_array.T # Convert to 32bit float for k in list(data.keys())[1:]: @@ -386,6 +370,8 @@ def write(self, outfile): parser = ugali.utils.parser.Parser(description=description) parser.add_config() parser.add_argument('outfile',metavar='outfile.fits',help='Output fits file.') + parser.add_argument('-m','--mlimit',nargs='?',default=-1, type=int, + help='Memory limit (KB)') parser.add_debug() parser.add_verbose() parser.add_coords(required=True,radius=False) @@ -395,6 +381,14 @@ def write(self, outfile): raise Exception('Must specify exactly one coordinate.') lon,lat,radius = opts.coords[0] + #if opts.mlimit > 0: + if True: + # Convert from KB to GB + mlimit,_ = LSF.get_memory_limit() + logger.info("Setting memory limit: %.1fGB"%(mlimit/1024.**3)) + soft,hard = LSF.set_memory_limit(mlimit) + logger.info("Memory limit: %.1fGB"%(soft/1024.**3)) + grid = createGridSearch(opts.config,lon,lat) if not opts.debug: result = grid.search() diff --git a/ugali/analysis/search.py b/ugali/analysis/search.py index 1307c8e..85b12a0 100644 --- a/ugali/analysis/search.py +++ b/ugali/analysis/search.py @@ -7,6 +7,7 @@ import copy import subprocess from collections import OrderedDict as odict +import glob import healpy import fitsio @@ -23,6 +24,7 @@ from ugali.utils.healpix import pix2ang, ang2pix from ugali.candidate.associate import SourceCatalog, catalogFactory from ugali.utils.config import Config +from ugali.utils import fileio class CandidateSearch(object): """ @@ -54,16 +56,44 @@ def _load(self): self.loadLikelihood() self.loadROI() - def loadLikelihood(self,filename=None): - if filename is None: filename = self.mergefile - f = fitsio.FITS(self.mergefile) - data = f[1].read() + def loadLikelihood(self,filenames=None): + """Load the merged likelihood healpix map. + + Parameters + ---------- + filename : input filename for sparse healpix map (default: mergefile) + + Returns + ------- + None : sets attributes: `pixels`, `values`, `distances`, `richness` + """ + if filenames is None: + if os.path.exists(self.mergefile): + filenames = self.mergefile + else: + filenames = glob.glob(self.mergefile.split('_%')[0]+'_*.fits') + + filenames = np.atleast_1d(filenames) + data = fileio.load_files(filenames) + self.pixels = data['PIXEL'] if 'PIXEL' in data.dtype.names else data['PIX'] self.values = 2*data['LOG_LIKELIHOOD'] - self.distances = f[2].read()['DISTANCE_MODULUS'] self.richness = data['RICHNESS'] + # Load distances from first file (should all be the same) + self.distances = fileio.load_files(filenames[0],ext=2,columns='DISTANCE_MODULUS') + def loadROI(self,filename=None): + """Load the ROI parameter sparse healpix map. + + Parameters + ---------- + filename : input filename for sparse healpix map (default: roifile) + + Returns + ------- + None : sets attributes: `ninterior`, `nannulus`, `stellar` + """ if filename is None: filename = self.roifile self.ninterior = healpix.read_partial_map(filename,'NINSIDE')[-1] @@ -71,16 +101,29 @@ def loadROI(self,filename=None): self.stellar = healpix.read_partial_map(filename,'STELLAR')[-1] def createLabels2D(self): - """ 2D labeling at zmax """ + """Create a flattened (2D) labeled array at zmax. + + Parameters + ---------- + None + + Returns + ------- + labels, nlables : labeled healpix array. + """ logger.debug(" Creating 2D labels...") - self.zmax = np.argmax(self.values,axis=1) - self.vmax = self.values[np.arange(len(self.pixels),dtype=int),self.zmax] + zmax = np.argmax(self.values,axis=1) + vmax = self.values[np.arange(len(self.pixels),dtype=int),zmax] - kwargs=dict(pixels=self.pixels,values=self.vmax,nside=self.nside, + kwargs=dict(pixels=self.pixels,values=vmax,nside=self.nside, threshold=self.threshold,xsize=self.xsize) labels,nlabels = CandidateSearch.labelHealpix(**kwargs) self.nlabels = nlabels + + # This is wasteful, but returns the same shape array for 2D and 3D labels self.labels = np.repeat(labels,len(self.distances)).reshape(len(labels),len(self.distances)) + #self.labels = labels + return self.labels, self.nlabels def createLabels3D(self): @@ -91,13 +134,20 @@ def createLabels3D(self): return self.labels, self.nlabels def writeLabels(self,filename=None): + """Write the labeled array to a file. + + Parameters + ---------- + filename : output filename (default: labelfile) + + Returns + ------- + None + """ if filename is None: filename = self.labelfile - # ADW: Is it necessary to convert labels? - data_dict = {'PIXEL':self.pixels, - 'LABEL':self.labels.astype(float,copy=False)} + data_dict = {'PIXEL':self.pixels,'LABEL':self.labels} - logger.info("Writing %s..."%filename) healpix.write_partial_map(filename,data_dict,self.nside) fitsio.write(filename, {'DISTANCE_MODULUS':self.distances.astype('f4',copy=False)}, @@ -106,30 +156,56 @@ def writeLabels(self,filename=None): def loadLabels(self,filename=None): + """Load labels from filename. + + Parameters + ---------- + filename : labeled healpix map + + Returns + ------- + labels,nlabels : labels and number of labels + """ if filename is None: filename = self.labelfile - data = fitsio.read(filename) - distances = fitsio.read(filename,ext='DISTANCE_MODULUS')['DISTANCE_MODULUS'] + f = fitsio.FITS(filename) + + data = f['PIX_DATA'].read() if not (self.pixels == data['PIXEL']).all(): raise Exception("Pixels do not match") + self.pixels = data['PIXEL'] # this is just to save memory + + distances = f['DISTANCE_MODULUS'].read(columns='DISTANCE_MODULUS') if not (self.distances == distances).all(): raise Exception("Distance moduli do not match.") + self.distances = distances - self.labels = data['LABEL'].astype(int) + self.labels = data['LABEL'] self.nlabels = self.labels.max() if self.nlabels != (len(np.unique(self.labels)) - 1): raise Exception("Incorrect number of labels found.") + return self.labels, self.nlabels - + def createObjects(self): + """Create list of objects from the labeled array + + Parameters + ---------- + None + + Returns + ------- + objects : array of objects + """ logger.debug(" Creating objects...") - hist,edges,rev = reverseHistogram(self.labels,bins=numpy.arange(self.nlabels+2)) + hist,edges,rev = reverseHistogram(self.labels,bins=np.arange(self.nlabels+2)) self.rev = rev # Make some cut on the minimum size of a labelled object - good, = numpy.where( (hist >= self.minpix) ) + good, = np.where( (hist >= self.minpix) ) # Get rid of zero label (below threshold) - self.good, = numpy.nonzero(good) + self.good, = np.nonzero(good) - kwargs=dict(pixels=self.pixels,values=self.values,nside=self.nside, + kwargs=dict(pixels=self.pixels,values=self.values,nside=self.nside, zvalues=self.distances, rev=self.rev, good=self.good) objects = self.findObjects(**kwargs) self.objects = self.finalizeObjects(objects) @@ -171,13 +247,13 @@ def labelHealpix(pixels, values, nside, threshold=0, xsize=1000): if values.ndim < 2: iterate = [values] else: iterate = values.T for i,value in enumerate(iterate): - logger.debug("Labeling slice %i...") - searchim = numpy.zeros(xx.shape,dtype=bool) + logger.debug("Labeling slice %i..."%i) + searchim = np.zeros(xx.shape,dtype=bool) select = (value > threshold) yidx = ij[0][select]; xidx = ij[1][select] searchim[yidx,xidx] |= True searchims.append( searchim ) - searchims = numpy.array(searchims) + searchims = np.array(searchims) # Full binary structure s = ndimage.generate_binary_structure(searchims.ndim,searchims.ndim) @@ -203,19 +279,19 @@ def findObjects(pixels, values, nside, zvalues, rev, good): Characterize labelled candidates in a multi-dimensional HEALPix map. Parameters: + pixels : Pixel values associated to (sparse) HEALPix array values : (Sparse) HEALPix array of data values nside : HEALPix dimensionality - pixels : Pixel values associated to (sparse) HEALPix array zvalues : Values of the z-dimension (usually distance modulus) rev : Reverse indices for pixels in each "island" good : Array containg labels for each "island" Returns: - objs : numpy.recarray of object characteristics + objs : np.recarray of object characteristics """ ngood = len(good) - objs = numpy.recarray((ngood,), + objs = np.recarray((ngood,), dtype=[('LABEL','i4'), ('NPIX','i4'), ('VAL_MAX','f4'), @@ -269,7 +345,7 @@ def findObjects(pixels, values, nside, zvalues, rev, good): xpix,ypix = proj.sphereToImage(xval,yval) # Projected centroid - x_cent,y_cent,zval_cent = numpy.average([xpix,ypix,zval],axis=1) + x_cent,y_cent,zval_cent = np.average([xpix,ypix,zval],axis=1) xval_cent, yval_cent = proj.imageToSphere(x_cent,y_cent) objs[i]['X_CENT'] = xval_cent objs[i]['Y_CENT'] = yval_cent @@ -277,7 +353,7 @@ def findObjects(pixels, values, nside, zvalues, rev, good): # Projected barycenter weights=[island,island,island] - x_bary,y_bary,zval_bary = numpy.average([xpix,ypix,zval],weights=weights,axis=1) + x_bary,y_bary,zval_bary = np.average([xpix,ypix,zval],weights=weights,axis=1) xval_bary,yval_bary = proj.imageToSphere(x_bary, y_bary) objs[i]['X_BARY'] = xval_bary objs[i]['Y_BARY'] = yval_bary @@ -286,7 +362,7 @@ def findObjects(pixels, values, nside, zvalues, rev, good): return objs def finalizeObjects(self, objects): - objs = numpy.recarray(len(objects), + objs = np.recarray(len(objects), dtype=[('NAME','S24'), ('TS','f4'), ('GLON','f4'), @@ -334,21 +410,6 @@ def finalizeObjects(self, objects): return out - def loadLabels(self,filename=None): - if filename is None: filename = self.labelfile - data = fitsio.read(filename) - distances = fitsio.read(filename,ext='DISTANCE_MODULUS')['DISTANCE_MODULUS'] - if not (self.pixels == data['PIXEL']).all(): - raise Exception("Pixels do not match") - if not (self.distances == distances).all(): - raise Exception("Distance moduli do not match.") - - self.labels = data['LABEL'].astype(int) - self.nlabels = self.labels.max() - if self.nlabels != (len(np.unique(self.labels)) - 1): - raise Exception("Incorrect number of labels found.") - return self.labels, self.nlabels - def loadObjects(self,filename=None): if filename is None: filename = self.objectfile self.objects = fitsio.read(filename) diff --git a/ugali/isochrone/model.py b/ugali/isochrone/model.py index 7858d20..8f592dc 100644 --- a/ugali/isochrone/model.py +++ b/ugali/isochrone/model.py @@ -355,6 +355,8 @@ def absolute_magnitude(self, richness=1, steps=1e4): [astro-ph/0506022]. + TODO: ADW If richness not specified, should use self.richness + Parameters: ----------- richness : isochrone normalization parameter @@ -1159,6 +1161,17 @@ def download(self,age=None,metallicity=None,outdir=None,force=False): """ Check valid parameter range and download isochrones from: http://stev.oapd.inaf.it/cgi-bin/cmd + + Parameters + ---------- + age : age in (Gyr) + metallicity : Z + outdir : output directory (default to current directory) + force : force overwrite of file + + Returns + ------- + outfile : the output isochrone """ try: from urllib.error import URLError diff --git a/ugali/isochrone/parsec.py b/ugali/isochrone/parsec.py index dd62fa7..5157fec 100644 --- a/ugali/isochrone/parsec.py +++ b/ugali/isochrone/parsec.py @@ -92,10 +92,68 @@ #'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') +# Access prior to 3.1 seems to be gone +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) + +# This seems to maintain old ischrone format +defaults_31 = dict(defaults_cmd,cmd_version=3.1) + +# New query and file format for 3.3... +defaults_33 = {'cmd_version': 3.3, + 'track_parsec': 'parsec_CAF09_v1.2S', + 'track_colibri': 'parsec_CAF09_v1.2S_S35', + 'track_postagb': 'no', + 'n_inTPC': 10, + 'eta_reimers': 0.2, + 'kind_interp': 1, + 'kind_postagb': -1, + 'photsys_file': photsys_dict['des'], + 'photsys_version': 'OBC', + 'dust_sourceM': 'dpmod60alox40', + 'dust_sourceC': 'AMCSIC15', + 'kind_mag': 2, + 'kind_dust': 0, + #'extinction_av': 0.0, + 'extinction_coeff': 'constant', + 'extinction_curve': 'cardelli', + 'imf_file': 'tab_imf/imf_chabrier_lognormal.dat', + 'isoc_isagelog': 0, + 'isoc_agelow': 1.0e9, + 'isoc_ageupp': 1.0e10, + 'isoc_dage': 0.0, + 'isoc_lagelow': 6.6, + 'isoc_lageupp': 10.13, + 'isoc_dlage': 0.0, + 'isoc_ismetlog': 0, + 'isoc_zlow': 0.0152, + 'isoc_zupp': 0.03, + 'isoc_dz': 0.0, + 'isoc_metlow': -2, + 'isoc_metupp': 0.3, + 'isoc_dmet': 0.0, + 'output_kind': 0, + 'output_evstage': 1, + #'lf_maginf': -15, + #'lf_magsup': 20, + #'lf_deltamag': 0.5, + #'sim_mtot': 1.0e4, + 'submit_form': 'Submit', + #'.cgifields': 'dust_sourceC', + #'.cgifields': 'track_colibri', + #'.cgifields': 'extinction_curve', + #'.cgifields': 'output_kind', + #'.cgifields': 'photsys_version', + #'.cgifields': 'isoc_isagelog', + #'.cgifields': 'track_parsec', + #'.cgifields': 'extinction_coeff', + #'.cgifields': 'track_postagb', + #'.cgifields': 'output_gzip', + #'.cgifields': 'isoc_ismetlog', + #'.cgifields': 'dust_sourceM', + } class ParsecIsochrone(Isochrone): """ Base class for PARSEC-style isochrones. """ @@ -154,20 +212,29 @@ def query_server(self,outfile,age,metallicity): epsilon = 1e-4 lage = np.log10(age*1e9) - lage_min,lage_max = params['isoc_lage0'],params['isoc_lage1'] + + lage_min = params.get('isoc_lage0',6.602) + lage_max = params.get('isoc_lage1',10.1303) + if not (lage_min-epsilon < lage = threshold] + + return candidates +def run(self): if 'label' in self.opts.run: logger.info("Running 'label'...") - if exists(search.labelfile) and not self.opts.force: - logger.info(" Found %s; skipping..."%search.labelfile) + if not hasattr(self,'search'): + self.search = CandidateSearch(self.config) + if exists(self.search.labelfile) and not self.opts.force: + logger.info(" Found %s; skipping..."%self.search.labelfile) else: - #search.createLabels3D() - search.createLabels2D() - search.writeLabels() + #self.search.createLabels3D() + #self.search.loadLikelhood() + #self.search.loadROI() + self.search.createLabels2D() + self.search.writeLabels() if 'objects' in self.opts.run: logger.info("Running 'objects'...") - if exists(search.objectfile) and not self.opts.force: - logger.info(" Found %s; skipping..."%search.labelfile) + if not hasattr(self,'search'): + self.search = CandidateSearch(self.config) + if exists(self.search.objectfile) and not self.opts.force: + logger.info(" Found %s; skipping..."%self.search.labelfile) else: - search.loadLabels() - search.createObjects() - search.writeObjects() + self.search.loadLabels() + self.search.createObjects() + self.search.writeObjects() if 'associate' in self.opts.run: logger.info("Running 'associate'...") - if exists(search.assocfile) and not self.opts.force: - logger.info(" Found %s; skipping..."%search.assocfile) + if not hasattr(self,'search'): + self.search = CandidateSearch(self.config) + if exists(self.search.assocfile) and not self.opts.force: + logger.info(" Found %s; skipping..."%self.search.assocfile) else: - search.loadObjects() - search.createAssociations() - search.writeAssociations() + self.search.loadObjects() + self.search.createAssociations() + self.search.writeAssociations() if 'candidate' in self.opts.run: logger.info("Running 'candidate'...") - if exists(search.candfile) and not self.opts.force: - logger.info(" Found %s; skipping..."%search.candfile) + if exists(self.search.candfile) and not self.opts.force: + logger.info(" Found %s; skipping..."%self.search.candfile) else: - search.loadAssociations() - search.writeCandidates() + self.search.loadAssociations() + self.search.writeCandidates() if 'plot' in self.opts.run: + self.opts.run.append('www') logger.info("Running 'plot'...") - import fitsio threshold = self.config['search']['cand_threshold'] outdir = mkdir(self.config['output']['plotdir']) - logdir = mkdir(os.path.join(outdir,'log')) + logdir = mkdir(join(outdir,'log')) # Eventually move this into 'plotting' module - candidates = fitsio.read(self.config.candfile,lower=True,trim_strings=True) - candidates = candidates[candidates['ts'] >= threshold] + candidates = load_candidates(self.config.candfile,threshold) for i,c in enumerate(candidates): - msg = "(%i/%i) Plotting %s (%.2f,%.2f)..."%(i,len(candidates),c['name'],c['ra'],c['dec']) + name = c['name'].replace('(','').replace(')','') + msg = "(%i/%i) Plotting %s (%.2f,%.2f)..."%(i,len(candidates),name,c['ra'],c['dec']) logger.info(msg) - params = (self.opts.config,outdir,c['name'],c['ra'], + params = (self.opts.config,outdir,name,c['ra'], c['dec'],0.5,c['modulus']) cmd = 'ugali/scratch/PlotCandidate.py %s %s -n="%s" --cel %f %f --radius %s -m %.2f' cmd = cmd%params - logger.info(cmd) - jobname = c['name'].lower().replace(' ','_') - logfile = os.path.join(logdir,jobname+'.log') + jobname = name.lower().replace(' ','_') + logfile = join(logdir,jobname+'.log') batch = self.config['search'].get('batch',self.config['batch']) - self.batch.submit(cmd,jobname,logfile,**batch['opts']) - time.sleep(3) + out = [join(outdir,jobname+'.png'), + join(outdir,jobname+'_dist.png'), + join(outdir,jobname+'_scat.png')] + if all([exists(o) for o in out]) and not self.opts.force: + logger.info(" Found plots for %s; skipping..."%name) + else: + logger.info(cmd) + self.batch.submit(cmd,jobname,logfile,**batch.get(self.opts.queue,{})) + time.sleep(3) + + if 'www' in self.opts.run: + logger.info("Running 'www'...") + + threshold = self.config['search']['cand_threshold'] + outdir = mkdir(self.config['output']['plotdir']) + + # Eventually move this into 'plotting' module + candidates = load_candidates(self.config.candfile,threshold) + + from ugali.utils.www import create_index_html + filename = os.path.join(outdir,'index.html') + create_index_html(filename,candidates) + Pipeline.run = run pipeline = Pipeline(__doc__,components) pipeline.parse_args() diff --git a/ugali/pipeline/run_06.0_simulate.py b/ugali/pipeline/run_06.0_simulate.py index e27d614..ab7c88d 100755 --- a/ugali/pipeline/run_06.0_simulate.py +++ b/ugali/pipeline/run_06.0_simulate.py @@ -6,6 +6,7 @@ from os.path import join, splitext, exists import time import glob +import shutil import numpy as np import numpy.lib.recfunctions as recfuncs @@ -20,12 +21,17 @@ from ugali.utils.logger import logger from ugali.utils.healpix import pix2ang -components = ['simulate','analyze','merge','plot'] +#components = ['simulate','analyze','merge','plot'] +components = ['analyze','merge'] def run(self): outdir=mkdir(self.config['output']['simdir']) logdir=mkdir(join(outdir,'log')) + # Actually copy config instead of re-writing + shutil.copy(self.config.filename,outdir) + configfile = join(outdir,os.path.basename(self.config.filename)) + if 'simulate' in self.opts.run: logger.info("Running 'simulate'...") @@ -36,7 +42,7 @@ def run(self): logfile=join(logdir,base+'.log') jobname=base script = self.config['simulator']['script'] - cmd='%s %s %s --seed %i'%(script,self.opts.config,outfile,i) + cmd='%s %s %s --seed %i'%(script,configfile,outfile,i) #cmd='%s %s %s'%(script,self.opts.config,outfile) self.batch.submit(cmd,jobname,logfile) time.sleep(0.1) @@ -51,16 +57,22 @@ def run(self): for i,catfile in enumerate(catfiles): basename = os.path.basename(catfile) outfile = join(outdir,basename) - - if exists(outfile) and not self.opts.force: - logger.info("Found %s; skipping..."%outfile) - continue - base = splitext(os.path.basename(outfile))[0] logfile=join(logdir,base+'.log') jobname=base + + if exists(outfile) and not self.opts.force: + msg = "Found %s;"%outfile + if exists(logfile) and len(self.batch.bfail(logfile)): + msg += " failed." + logger.info(msg) + else: + msg += " skipping..." + logger.info(msg) + continue + script = self.config['simulate']['script'] - cmd='%s %s -p %s -c %s -o %s'%(script,self.opts.config,popfile,catfile,outfile) + cmd='%s %s -m 0 --rerun -p %s -c %s -o %s'%(script,configfile,popfile,catfile,outfile) self.batch.max_jobs = batch.get('max_jobs',200) opts = batch.get(self.opts.queue,dict()) self.batch.submit(cmd,jobname,logfile,**opts) @@ -72,21 +84,15 @@ def run(self): if 'merge' in self.opts.run: logger.info("Running 'merge'...") - filenames=join(outdir,self.config['output']['simfile']).split('_%')[0]+'_*' + filenames=join(outdir,self.config['simulate']['catfile']) infiles=sorted(glob.glob(filenames)) + print("Reading %i files..."%len(infiles)) + data = np.concatenate([fitsio.read(f,ext=1) for f in infiles]) + hdr = fitsio.read_header(infiles[0],ext=1) - f = fitsio.read(infiles[0]) - table = np.empty(0,dtype=data.dtype) - for filename in infiles: - logger.debug("Reading %s..."%filename) - d = fitsio.read(filename) - t = d[~np.isnan(d['ts'])] - table = recfuncs.stack_arrays([table,t],usemask=False,asrecarray=True) - - logger.info("Found %i simulations."%len(table)) - outfile = join(outdir,"merged_sims.fits") + outfile = "./merged_sims.fits" logger.info("Writing %s..."%outfile) - fitsio.write(outfile,table,clobber=True) + fitsio.write(outfile,data,header=hdr,clobber=True) if 'plot' in self.opts.run: logger.info("Running 'plot'...") diff --git a/ugali/preprocess/padova.py b/ugali/preprocess/padova.py index eedfb8d..3d4ad03 100755 --- a/ugali/preprocess/padova.py +++ b/ugali/preprocess/padova.py @@ -23,7 +23,7 @@ import numpy as np from ugali.utils.logger import logger from ugali.utils.shell import mkdir -from ugali.analysis.isochrone import Padova as PadovaIsochrone +from ugali.analysis.isochrone import Padova # survey system photsys_dict = odict([ @@ -171,7 +171,7 @@ def download(self,age,metallicity,outdir=None,force=False): class Padova(Download): defaults = copy.deepcopy(defaults_27) - isochrone = PadovaIsochrone + isochrone = Padova abins = np.arange(1.0, 13.5 + 0.1, 0.1) zbins = np.arange(1e-4,1e-3 + 1e-5,1e-5) @@ -208,7 +208,6 @@ def query_server(self,outfile,age,metallicity): msg = "Output filename not found" raise RuntimeError(msg) - #out = '{0}/~lgirardi/tmp/{1}.dat'.format(server, fname[0]) out = '{0}/tmp/{1}.dat'.format(server, fname[0]) cmd = 'wget %s -O %s'%(out,outfile) logger.debug(cmd) @@ -321,8 +320,8 @@ def factory(name, **kwargs): parser = ugali.utils.parser.Parser(description=description) parser.add_verbose() parser.add_force() - parser.add_argument('-a','--age',default=None,type=float) - parser.add_argument('-z','--metallicity',default=None,type=float) + 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='Bressan2012') parser.add_argument('-s','--survey',default='des') parser.add_argument('-o','--outdir',default=None) @@ -332,10 +331,11 @@ def factory(name, **kwargs): if args.verbose: try: from http.client import HTTPConnection - except ImportError: from httplib import HTTPConnection + except ImportError: + from httplib import HTTPConnection HTTPConnection.debuglevel = 1 - if args.outdir is None: + 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) @@ -349,9 +349,9 @@ def factory(name, **kwargs): logger.info("Metallicities:\n %s"%np.unique(grid[1])) def run(args): - try: + try: p.download(*args) - except Exception as e: + except Exception as e: logger.warn(str(e)) logger.error("Download failed.") diff --git a/ugali/scratch/simulation/farm_simulate_population.py b/ugali/scratch/simulation/farm_simulate_population.py old mode 100644 new mode 100755 index 7481f9f..c5a7875 --- a/ugali/scratch/simulation/farm_simulate_population.py +++ b/ugali/scratch/simulation/farm_simulate_population.py @@ -1,9 +1,9 @@ +#!/usr/bin/env python import os import subprocess -import numpy as np tag = 'ps1_v1' # PS1 -tag = 'v7' # DES +tag = 'des_v7' # DES n_chunk = 100 mc_source_id_start_global = 1 size_batch = 1000 @@ -13,12 +13,26 @@ parser = argparse.ArgumentParser(description=__doc__) parser.add_argument('config') parser.add_argument('-q','--queue',default='local') -parser.add_argument('-t','--tag',default=tag) -parser.add_argument('-k','--nchunk',default=n_chunk,type=int) -parser.add_argument('-b','--size-batch',default=size_batch,type=int) -parser.add_argument('-n','--nbatch',default=number_of_batches,type=int) -parser.add_argument('--mc_source_id',default=mc_source_id_start_global,type=int) -parser.add_argument('-s','--section',default='des',choices=['des','ps1']) +parser.add_argument('-t','--tag',default=tag, + help='tag appended to file name') +parser.add_argument('-k','--nchunk',default=n_chunk,type=int, + help='number of satellites per catalog file') +parser.add_argument('-b','--size-batch',default=size_batch,type=int, + help='number of satellites per population file') +parser.add_argument('-n','--nbatch',default=number_of_batches,type=int, + 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'], + help='section of config file') +parser.add_argument('--dwarfs', dest='dwarfs', action='store_true', + help="Simulate from known dwarfs") +parser.add_argument('--sleep',default=None,type=int, + help='sleep between jobs (seconds)') +parser.add_argument('--njobs',default=10,type=int, + help='number of jobs to run concurrently') +parser.add_argument('--host',default=None,metavar='host1[,host2,...]', + help="submit to comma-delimited list of hosts or 'all'") args = parser.parse_args() @@ -26,21 +40,24 @@ # /home/s1/kadrlica/bin/csub # I try to set it here, but you can always do this in shell # export PATH=/home/s1/kadrlica/bin:$PATH - os.environ['PATH'] = '/home/s1/kadrlica/bin:'+os.environ['PATH'] logdir = '%s/log'%args.tag if not os.path.exists(logdir): os.makedirs(logdir) -for index_batch in np.arange(args.nbatch): +for index_batch in range(args.nbatch): seed = mc_source_id_start = args.mc_source_id + (args.size_batch * index_batch) - command = 'simulate_population.py %s -s %s --tag %s --start %i --size %i --chunk %i --seed %i'%(args.config, args.section, args.tag, mc_source_id_start, args.size_batch, args.nchunk, seed) - + dwarfs = '--dwarfs' if args.dwarfs else '' + command = 'simulate_population.py %s -s %s --tag %s --start %i --size %i --chunk %i --seed %i %s'%(args.config, args.section, args.tag, mc_source_id_start, args.size_batch, args.nchunk, seed, dwarfs) + + # Max number of jobs limited by memory logfile = os.path.join(logdir,'%07i.log'%mc_source_id_start) - submit = 'csub -q %s -o %s -n 10 '%(args.queue,logfile) + command + params = dict(queue=args.queue,logfile=logfile, + sleep='-s %d'%args.sleep if args.sleep else '', + host='--host %s'%args.host if args.host else '', + njobs='-n %s'%args.njobs if args.njobs else '', + command=command) + submit = 'csub -q %(queue)s %(njobs)s %(sleep)s %(host)s -o %(logfile)s %(command)s'%params + print(submit) subprocess.call(submit,shell=True) - #print command - #os.system(command) - - diff --git a/ugali/scratch/simulation/simulate_population.py b/ugali/scratch/simulation/simulate_population.py index 0c1af2f..7cc0847 100755 --- a/ugali/scratch/simulation/simulate_population.py +++ b/ugali/scratch/simulation/simulate_population.py @@ -9,7 +9,7 @@ import yaml import numpy as np import scipy.interpolate -import healpy +import healpy as hp import fitsio import astropy.io.fits as pyfits @@ -44,6 +44,17 @@ def getCompleteness(config): ############################################################ def getPhotoError(config): + """Photometric error model based on the delta-mag from the maglim and + the photometric uncertainty estimated from the data + + Parameters + ---------- + config : configuration dictionary + + Returns + ------- + interp : interpolation function (mag_err as a function of delta_mag) + """ #infile = 'photo_error_model.csv' infile = config['photo_error'] d = np.recfromcsv(infile) @@ -87,7 +98,7 @@ def meanFracdet(map_fracdet, lon_population, lat_population, radius_population): lon, lat, and radius are taken to be arrays of the same length """ - nside_fracdet = healpy.npix2nside(len(map_fracdet)) + nside_fracdet = hp.npix2nside(len(map_fracdet)) map_fracdet_zero = np.where(map_fracdet >= 0., map_fracdet, 0.) fracdet_population = np.empty(len(lon_population)) for ii in range(0, len(lon_population)): @@ -101,11 +112,33 @@ def meanFracdet(map_fracdet, lon_population, lat_population, radius_population): ############################################################ def catsimSatellite(config, lon_centroid, lat_centroid, distance, stellar_mass, r_physical, + ellipticity, position_angle, age, metal, m_maglim_1, m_maglim_2, m_ebv, plot=False, title='test'): - """ - Simulate a single satellite. This is currently only valid for band_1 = g and band_2 = r. - r_physical is azimuthally averaged half-light radius, kpc + """ Simulate a single satellite. This is currently only valid for + band_1 = g and band_2 = r. r_physical is azimuthally averaged + half-light radius, kpc + + Parameters + ---------- + config : configuration + lon_centroid : longitude centroid (deg) + lat_centroid : latitude centroid (deg) + distance : distance (kpc) + stellar_mass : stellar mass (Msun) + r_physical : azimuthally averaged physical half-light radius (kpc) + ellipticity : ellipticity [0, 1] + position_angle : position angle (deg) + age : age (Gyr) + metal : metallicity + m_maglim_1 : mask of magnitude limit in band 2 + m_maglim_2 : mask of magnitude limit in band 2 + m_ebv : mask of E(B-V) values + plot : Plot the output [False] + + Returns + ------- + satellite : ordered dictionary of satellite star output """ # Probably don't want to parse every time @@ -114,40 +147,39 @@ def catsimSatellite(config, lon_centroid, lat_centroid, distance, stellar_mass, s = ugali.analysis.source.Source() - # Following McConnachie 2012, ellipticity = 1 - (b/a) , where a is semi-major axis and b is semi-minor axis - - r_h = np.degrees(np.arcsin(r_physical / distance)) # Azimuthally averaged half-light radius - #ellipticity = 0.3 # Semi-arbitrary default for testing purposes - # See http://iopscience.iop.org/article/10.3847/1538-4357/833/2/167/pdf - # Based loosely on https://arxiv.org/abs/0805.2945 - ellipticity = np.random.uniform(0.1, 0.8) - position_angle = np.random.uniform(0., 180.) # Random position angle (deg) - a_h = r_h / np.sqrt(1. - ellipticity) # semi-major axis (deg) - - - # Elliptical kernels take the "extension" as the semi-major axis - ker = ugali.analysis.kernel.EllipticalPlummer(lon=lon_centroid, lat=lat_centroid, ellipticity=ellipticity, position_angle=position_angle) + # Azimuthally averaged projected half-light radius (deg) + r_h = np.degrees(np.arcsin(r_physical / distance)) + # Elliptical half-light radius along semi-major axis (deg) + a_h = r_h / np.sqrt(1. - ellipticity) + # Create the kernel without extension + ker = ugali.analysis.kernel.EllipticalPlummer(lon=lon_centroid, lat=lat_centroid, + ellipticity=ellipticity, position_angle=position_angle) + # Apply a max extension cut flag_too_extended = False - if a_h >= 1.0: - print('Too extended: a_h = %.2f'%(a_h)) - a_h = 1.0 + max_extension = 5.0 # deg + if a_h >= max_extension: + print 'Too extended: a_h = %.2f'%(a_h) + a_h = max_extension flag_too_extended = True - ker.setp('extension', value=a_h, bounds=[0.0,1.0]) + # 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) - - age = np.random.choice([10., 12.0, 13.5]) - metal_z = np.random.choice([0.0001, 0.0002]) - distance_modulus = ugali.utils.projector.distanceToDistanceModulus(distance) - iso = isochrone_factory('Bressan2012', survey=config['survey'], age=age, z=metal_z, distance_modulus=distance_modulus) + + # Create the isochrone + distance_modulus = ugali.utils.projector.dist2mod(distance) + iso = isochrone_factory('Bressan2012', survey=config['survey'], age=age, z=metal, + distance_modulus=distance_modulus) s.set_isochrone(iso) # Simulate takes stellar mass as an argument, NOT richness mag_1, mag_2 = s.isochrone.simulate(stellar_mass) + # Generate the positions of stars lon, lat = s.kernel.sample_lonlat(len(mag_2)) - nside = healpy.npix2nside(len(m_maglim_1)) # Assuming that the two maglim maps have same resolution - pix = ugali.utils.healpix.angToPix(nside, lon, lat) + nside = hp.npix2nside(len(m_maglim_1)) # Assuming that the two maglim maps have same resolution + pix = ugali.utils.healpix.ang2pix(nside, lon, lat) maglim_1 = m_maglim_1[pix] maglim_2 = m_maglim_2[pix] if config['survey'] == 'des': @@ -201,9 +233,10 @@ def catsimSatellite(config, lon_centroid, lat_centroid, distance, stellar_mass, #abs_mag_martin = s.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 to do would be to sum up the actual simulated stars + # The more clever thing is to sum the simulated stars if config['survey'] == 'des': - v = mag_1 - 0.487*(mag_1 - mag_2) - 0.0249 # See https://github.com/DarkEnergySurvey/ugali/blob/master/ugali/isochrone/model.py + # See https://github.com/DarkEnergySurvey/ugali/blob/master/ugali/isochrone/model.py + v = mag_1 - 0.487*(mag_1 - mag_2) - 0.0249 elif config['survey'] == 'ps1': # https://arxiv.org/pdf/1706.06147.pdf # V - g = C_0 + C_1 * (g - r) @@ -218,18 +251,9 @@ def catsimSatellite(config, lon_centroid, lat_centroid, distance, stellar_mass, flux = np.sum(10**(-v/2.5)) abs_mag = -2.5*np.log10(flux) - distance_modulus - - - #print abs_mag, abs_mag_martin + # Realized surface brightness within azimuthally averaged half-light radius + surface_brightness = ugali.analysis.results.surfaceBrightness(abs_mag, r_physical, distance) - #distance = ugali.utils.projector.distanceModulusToDistance(distance_modulus) - #r_h = extension * np.sqrt(1. - ellipticity) # Azimuthally averaged half-light radius - r_physical = distance * np.tan(np.radians(r_h)) # Azimuthally averaged half-light radius, kpc - #print 'distance = %.3f kpc'%(distance) - #print 'r_physical = %.3f kpc'%(r_physical) - surface_brightness = ugali.analysis.results.surfaceBrightness(abs_mag, r_physical, distance) # Average within azimuthally averaged half-light radius - #print 'surface_brightness = %.3f mag arcsec^-2'%(surface_brightness) - if plot: import pylab pylab.ion() @@ -243,25 +267,40 @@ def catsimSatellite(config, lon_centroid, lat_centroid, distance, stellar_mass, pylab.xlabel('g - r') pylab.ylabel('g') pylab.title('Number of stars with g < 23: %i'%(n_sigma_p)) - pylab.savefig('y3_sat_sim_cmd_%s.png'%(title), dpi=150.) + pylab.savefig('y3_sat_sim_cmd_%s.png'%('test'), dpi=150.) print('n_Sigma_p = %i'%(n_sigma_p)) raw_input('WAIT') - #if flag_too_extended: - # # This is a kludge to remove these satellites. fragile!! - # n_g24 = 1.e6 + satellite=odict(lon=lon[cut_detect], lat=lat[cut_detect], + mag_1=mag_1_meas[cut_detect], mag_2=mag_2_meas[cut_detect], + mag_1_error=mag_1_error[cut_detect], mag_2_error=mag_2_error[cut_detect], + mag_extinction_1=mag_extinction_1[cut_detect], + mag_extinction_2=mag_extinction_2[cut_detect], + n_g22=n_g22, n_g24=n_g24, abs_mag=abs_mag, surface_brightness=surface_brightness, + extension=extension, flag_too_extended=flag_too_extended) - return lon[cut_detect], lat[cut_detect], mag_1_meas[cut_detect], mag_2_meas[cut_detect], mag_1_error[cut_detect], mag_2_error[cut_detect], mag_extinction_1[cut_detect], mag_extinction_2[cut_detect], n_g22, n_g24, abs_mag, surface_brightness, ellipticity, position_angle, age, metal_z, flag_too_extended + return satellite ############################################################ #from memory_profiler import profile #@profile -def catsimPopulation(tag, mc_source_id_start=1, n=5000, n_chunk=100, config='simulate_population.yaml'): - """ - n = Number of satellites to simulation - n_chunk = Number of satellites in a file chunk +def catsimPopulation(config, tag, mc_source_id_start=1, n=5000, n_chunk=100, + known_dwarfs=False): + """ Create a population of satellites and then simulate the stellar distributions for each. + + Parameters + ---------- + config : configuration file or dictionary + tag : output name + mc_source_id_start : starting value of source id + n : number of satellites to simulate [5000] + n_chunk : number of satellites written in a file chunk + + Returns + ------- + None """ assert mc_source_id_start >= 1, "Starting mc_source_id must be >= 1" @@ -279,14 +318,19 @@ def catsimPopulation(tag, mc_source_id_start=1, n=5000, n_chunk=100, config='sim infile_maglim_r = config['maglim_r'] infile_density = config['stellar_density'] - range_distance = config.get('range_distance',[5., 500.]) - range_stellar_mass = config.get('range_stellar_mass',[1.e1, 1.e6]) - range_r_physical = config.get('range_r_physical',[1.e-3, 2.0]) - + range_distance = config.get('range_distance',[5., 500.]) + range_stellar_mass = config.get('range_stellar_mass',[1.e1, 1.e6]) + range_r_physical = config.get('range_r_physical',[1.e-3, 2.0]) + range_ellipticity = config.get('range_ellipticity',[0.1, 0.8]) + range_position_angle= config.get('range_position_angle',[0.0, 180.0]) + choice_age = config.get('choice_age',[10., 12.0, 13.5]) + choice_metal = config.get('choice_metal',[0.00010,0.00020]) + dwarf_file = config.get('known_dwarfs',None) + m_density = np.load(infile_density) - nside_density = healpy.npix2nside(len(m_density)) + nside_density = hp.npix2nside(len(m_density)) m_fracdet = read_map(infile_fracdet, nest=False) #.astype(np.float16) - nside_fracdet = healpy.npix2nside(len(m_fracdet)) + nside_fracdet = hp.npix2nside(len(m_fracdet)) m_maglim_g = read_map(infile_maglim_g, nest=False) #.astype(np.float16) m_maglim_r = read_map(infile_maglim_r, nest=False) #.astype(np.float16) @@ -297,22 +341,31 @@ def catsimPopulation(tag, mc_source_id_start=1, n=5000, n_chunk=100, config='sim mask = (m_fracdet > 0.5) - kwargs = dict(range_distance = range_distance, - range_stellar_mass = range_stellar_mass, - range_r_physical = range_r_physical) - print(kwargs) - # r_physical is azimuthally-averaged half-light radius, kpc - simulation_area, lon_population, lat_population, distance_population, stellar_mass_population, r_physical_population = ugali.simulation.population.satellitePopulation(mask, nside_pix, n, **kwargs) + if known_dwarfs: + # Simulate from known dwarfs + if dwarf_file is None: raise Exception("Must provide known_dwarf file") + print("Simulating dwarfs from: %s"%dwarf_file) + area, population = ugali.simulation.population.knownPopulation(dwarf_file, mask, nside_pix, n) + else: + # r_physical is azimuthally-averaged half-light radius, kpc + kwargs = dict(range_distance = range_distance, + range_stellar_mass = range_stellar_mass, + range_r_physical = range_r_physical, + range_ellipticity=[0.1, 0.8], + range_position_angle=[0.0, 180.0], + choice_age=[10., 12.0, 13.5], + choice_metal=[0.00010,0.00020], + plot=False) + area, population = ugali.simulation.population.satellitePopulation(mask, nside_pix, n, **kwargs) + + population['id'] += mc_source_id_start + simulation_area = area + n_g22_population = np.tile(np.nan, n) n_g24_population = np.tile(np.nan, n) abs_mag_population = np.tile(np.nan, n) surface_brightness_population = np.tile(np.nan, n) - ellipticity_population = np.tile(np.nan, n) - position_angle_population = np.tile(np.nan, n) - age_population = np.tile(np.nan, n) - metal_z_population = np.tile(np.nan, n) - mc_source_id_population = np.arange(mc_source_id_start, mc_source_id_start + n) - #cut_difficulty_population = np.tile(False, n) + extension_population = np.tile(np.nan, n) difficulty_population = np.tile(0, n) lon_array = [] @@ -324,44 +377,40 @@ def catsimPopulation(tag, mc_source_id_start=1, n=5000, n_chunk=100, config='sim mag_extinction_1_array = [] mag_extinction_2_array = [] mc_source_id_array = [] - for ii, mc_source_id in enumerate(mc_source_id_population): - print(' Simulating satellite (%i/%i) ... MC_SOURCE_ID = %i'%(ii + 1, n, mc_source_id)) - print(' distance=%.2e, stellar_mass=%.2e, rhalf=%.2e'%(distance_population[ii],stellar_mass_population[ii],r_physical_population[ii])) - lon, lat, mag_1, mag_2, mag_1_error, mag_2_error, mag_extinction_1, mag_extinction_2, n_g22, n_g24, abs_mag, surface_brightness, ellipticity, position_angle, age, metal_z, flag_too_extended = catsimSatellite(config, - lon_population[ii], - lat_population[ii], - distance_population[ii], - stellar_mass_population[ii], - r_physical_population[ii], - m_maglim_g, - m_maglim_r, - m_ebv) - print(' ', len(lon)) + 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]) + satellite = catsimSatellite(config, population[ii]['lon'], population[ii]['lat'], + population[ii]['distance'], population[ii]['stellar_mass'], + population[ii]['r_physical'],population[ii]['ellipticity'], + population[ii]['position_angle'],population[ii]['age'], + population[ii]['metallicity'], + m_maglim_g, m_maglim_r, m_ebv) - n_g22_population[ii] = n_g22 - n_g24_population[ii] = n_g24 - abs_mag_population[ii] = abs_mag - surface_brightness_population[ii] = surface_brightness - ellipticity_population[ii] = ellipticity - position_angle_population[ii] = position_angle - age_population[ii] = age - metal_z_population[ii] = metal_z - - #print "Difficulty masking..." + n_g22_population[ii] = satellite['n_g22'] + n_g24_population[ii] = satellite['n_g24'] + abs_mag_population[ii] = satellite['abs_mag'] + extension_population[ii] = satellite['extension'] + surface_brightness_population[ii] = satellite['surface_brightness'] # These objects are too extended and are not simulated - if (flag_too_extended): + if (satellite['flag_too_extended']): difficulty_population[ii] |= 0b0001 # We assume that these objects would be easily detected and # remove them to reduce data volume - if (surface_brightness_population[ii]<25.)&(n_g22_population[ii]>1e2): + if (surface_brightness_population[ii]<23.5)&(n_g22_population[ii]>1e3): difficulty_population[ii] |= 0b0010 - if (surface_brightness_population[ii]<28.)&(n_g22_population[ii]>1e4): - difficulty_population[ii] |= 0b0100 - if (surface_brightness_population[ii]<30.)&(n_g22_population[ii]>1e5): - difficulty_population[ii] |= 0b1000 + + # ADW 2019-08-31: I don't think these were implmented + #if (surface_brightness_population[ii]<25.)&(n_g22_population[ii]>1e2): + # difficulty_population[ii] |= 0b0010 + #if (surface_brightness_population[ii]<28.)&(n_g22_population[ii]>1e4): + # difficulty_population[ii] |= 0b0100 + #if (surface_brightness_population[ii]<30.)&(n_g22_population[ii]>1e5): + # difficulty_population[ii] |= 0b1000 + # ADW: 2019-08-31: These were Keith's original cuts, which were too aggressive #cut_easy = (surface_brightness_population[ii]<25.)&(n_g22_population[ii]>1.e2) \ # | ((surface_brightness_population[ii] < 30.) & (n_g24_population[ii] > 1.e4)) \ # | ((surface_brightness_population[ii] < 31.) & (n_g24_population[ii] > 1.e5)) @@ -374,19 +423,21 @@ def catsimPopulation(tag, mc_source_id_start=1, n=5000, n_chunk=100, config='sim #if flag_too_extended: # difficulty_population[ii] += 3 # TOO EXTENDED + # Only write satellites that aren't flagged if difficulty_population[ii] == 0: - lon_array.append(lon) - lat_array.append(lat) - mag_1_array.append(mag_1) - mag_2_array.append(mag_2) - mag_1_error_array.append(mag_1_error) - mag_2_error_array.append(mag_2_error) - mag_extinction_1_array.append(mag_extinction_1) - mag_extinction_2_array.append(mag_extinction_2) - mc_source_id_array.append(np.tile(mc_source_id, len(lon))) - - # Concatenate all the arrays - + lon_array.append(satellite['lon']) + lat_array.append(satellite['lat']) + mag_1_array.append(satellite['mag_1']) + mag_2_array.append(satellite['mag_2']) + mag_1_error_array.append(satellite['mag_1_error']) + mag_2_error_array.append(satellite['mag_2_error']) + mag_extinction_1_array.append(satellite['mag_extinction_1']) + 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] + + # Concatenate the arrays print("Concatenating arrays...") lon_array = np.concatenate(lon_array) lat_array = np.concatenate(lat_array) @@ -399,9 +450,8 @@ def catsimPopulation(tag, mc_source_id_start=1, n=5000, n_chunk=100, config='sim mc_source_id_array = np.concatenate(mc_source_id_array) # Now do the masking all at once - print("Fracdet masking...") - pix_array = ugali.utils.healpix.angToPix(nside_fracdet, lon_array, lat_array) + pix_array = ugali.utils.healpix.ang2pix(nside_fracdet, lon_array, lat_array) cut_fracdet = (np.random.uniform(size=len(lon_array)) < m_fracdet[pix_array]) lon_array = lon_array[cut_fracdet] @@ -414,184 +464,167 @@ def catsimPopulation(tag, mc_source_id_start=1, n=5000, n_chunk=100, config='sim mag_extinction_2_array = mag_extinction_2_array[cut_fracdet] mc_source_id_array = mc_source_id_array[cut_fracdet] - # Cut out the entries that are easily detectable - """ - lon_population = lon_population[cut_difficulty_population] - lat_population = lat_population[cut_difficulty_population] - distance_population = distance_population[cut_difficulty_population] - stellar_mass_population = stellar_mass_population[cut_difficulty_population] - r_physical_population = r_physical_population[cut_difficulty_population] - n_g24_population = n_g24_population[cut_difficulty_population] - abs_mag_population = abs_mag_population[cut_difficulty_population] - surface_brightness_population = surface_brightness_population[cut_difficulty_population] - ellipticity_population = ellipticity_population[cut_difficulty_population] - position_angle_population = position_angle_population[cut_difficulty_population] - age_population = age_population[cut_difficulty_population] - metal_z_population = metal_z_population[cut_difficulty_population] - mc_source_id_population = mc_source_id_population[cut_difficulty_population] - """ - # Create bonus columns - print("Creating bonus columns...") - distance_modulus_population = ugali.utils.projector.distanceToDistanceModulus(distance_population) - hpix_32_population = ugali.utils.healpix.angToPix(32, lon_population, lat_population) # Make sure this matches the dataset + distance_modulus_population = ugali.utils.projector.dist2mod(population['distance']) + hpix_32_population = ugali.utils.healpix.ang2pix(32, population['lon'], population['lat']) # Make sure this matches the dataset # Local stellar density - pixarea = healpy.nside2pixarea(nside_density, degrees=True) * 60.**2 # arcmin^2 - density_population = m_density[ugali.utils.healpix.angToPix(nside_density, lon_population, lat_population)] / pixarea # arcmin^-2 + pixarea = hp.nside2pixarea(nside_density, degrees=True) * 60.**2 # arcmin^2 + density_population = m_density[ugali.utils.healpix.ang2pix(nside_density, population['lon'], population['lat'])] / pixarea # arcmin^-2 # Average fracdet within the azimuthally averaged half-light radius #m_fracdet_zero = np.where(m_fracdet >= 0., m_fracdet, 0.) #m_fracdet_zero = m_fracdet - r_half = np.degrees(np.arctan2(r_physical_population, distance_population)) # Azimuthally averaged half-light radius in degrees - fracdet_half_population = meanFracdet(m_fracdet, lon_population, lat_population, r_half) - fracdet_core_population = meanFracdet(m_fracdet, lon_population, lat_population, 0.1) - fracdet_wide_population = meanFracdet(m_fracdet, lon_population, lat_population, 0.5) + + # Azimuthally averaged half-light radius in degrees + r_half = np.degrees(np.arctan2(population['r_physical'], population['distance'])) + fracdet_half_population = meanFracdet(m_fracdet, population['lon'], population['lat'], r_half) + fracdet_core_population = meanFracdet(m_fracdet, population['lon'], population['lat'], 0.1) + fracdet_wide_population = meanFracdet(m_fracdet, population['lon'], population['lat'], 0.5) # Magnitude limits - nside_maglim = healpy.npix2nside(len(m_maglim_g)) - pix_population = ugali.utils.healpix.angToPix(nside_maglim, lon_population, lat_population) + nside_maglim = hp.npix2nside(len(m_maglim_g)) + pix_population = ugali.utils.healpix.ang2pix(nside_maglim, population['lon'], population['lat']) maglim_g_population = m_maglim_g[pix_population] maglim_r_population = m_maglim_r[pix_population] # E(B-V) - nside_ebv = healpy.npix2nside(len(m_ebv)) - pix_population = ugali.utils.healpix.angToPix(nside_ebv, lon_population, lat_population) + nside_ebv = hp.npix2nside(len(m_ebv)) + pix_population = ugali.utils.healpix.ang2pix(nside_ebv, population['lon'], population['lat']) ebv_population = m_ebv[pix_population] # Survey - survey_population = np.tile(config['survey'], len(lon_population)) + survey_population = np.tile(config['survey'], len(population)) # Number of surviving catalog stars - n_catalog_population = np.histogram(mc_source_id_array, bins=np.arange(mc_source_id_population[0] - 0.5, mc_source_id_population[-1] + 0.51))[0] + n_catalog_population = np.histogram(mc_source_id_array, bins=np.arange(population['id'][0] - 0.5, population['id'][-1] + 0.51))[0] # Faked-up coadd_object_ids coadd_object_id_array = [] - for mc_source_id in mc_source_id_population: + for mc_source_id in population['id']: coadd_object_id_array.append((1000000 * mc_source_id) + 1 + np.arange(np.sum(mc_source_id == mc_source_id_array))) - coadd_object_id_array = -1 * np.concatenate(coadd_object_id_array) # Assign negative numbers to distinguish from real objects + # Assign negative numbers to distinguish from real objects + coadd_object_id_array = -1 * np.concatenate(coadd_object_id_array) - # Catalog output file + # Object ID assignment can get messed up if there are duplicate population ids + assert len(coadd_object_id_array) == len(mc_source_id_array) + # Population metadata output file + tbhdu = pyfits.BinTableHDU.from_columns([ + pyfits.Column(name='RA', format='E', array=population['lon'], unit='deg'), + pyfits.Column(name='DEC', format='E', array=population['lat'], unit='deg'), + pyfits.Column(name='DISTANCE', format='E', array=population['distance'], unit='kpc'), + pyfits.Column(name='DISTANCE_MODULUS', format='E', array=distance_modulus_population, unit='kpc'), + pyfits.Column(name='STELLAR_MASS', format='E', array=population['stellar_mass'], unit='Msun'), + pyfits.Column(name='R_PHYSICAL', format='E', array=population['r_physical'], unit='kpc'), + pyfits.Column(name='N_G22', format='J', array=n_g22_population, unit=''), + pyfits.Column(name='N_G24', format='J', array=n_g24_population, unit=''), + pyfits.Column(name='N_CATALOG', format='J', array=n_catalog_population, unit=''), + pyfits.Column(name='DIFFICULTY', format='J', array=difficulty_population, unit=''), + pyfits.Column(name='ABS_MAG', format='E', array=abs_mag_population, unit='mag'), + pyfits.Column(name='SURFACE_BRIGHTNESS', format='E', array=surface_brightness_population, unit='mag arcsec^-2'), + pyfits.Column(name='EXTENSION', format='E', array=extension_population, unit='deg'), + pyfits.Column(name='ELLIPTICITY', format='E', array=population['ellipticity'], unit=''), + pyfits.Column(name='POSITION_ANGLE', format='E', array=population['position_angle'], unit='deg'), + pyfits.Column(name='AGE', format='E', array=population['age'], unit='Gyr'), + pyfits.Column(name='METAL_Z', format='E', array=population['metallicity'], unit=''), + pyfits.Column(name='MC_SOURCE_ID', format='K', array=population['id'], unit=''), + pyfits.Column(name='HPIX_32', format='E', array=hpix_32_population, unit=''), + pyfits.Column(name='DENSITY', format='E', array=density_population, unit='arcmin^-2'), + pyfits.Column(name='FRACDET_HALF', format='E', array=fracdet_half_population, unit=''), + pyfits.Column(name='FRACDET_CORE', format='E', array=fracdet_core_population, unit=''), + pyfits.Column(name='FRACDET_WIDE', format='E', array=fracdet_wide_population, unit=''), + pyfits.Column(name='MAGLIM_G', format='E', array=maglim_g_population, unit='mag'), + pyfits.Column(name='MAGLIM_R', format='E', array=maglim_r_population, unit='mag'), + pyfits.Column(name='EBV', format='E', array=ebv_population, unit='mag'), + pyfits.Column(name='SURVEY', format='A12', array=survey_population, unit=''), + ]) + tbhdu.header.set('AREA', simulation_area, 'Simulation area (deg^2)') + print("Writing population metadata file...") + filename = '%s/sim_population_%s_mc_source_id_%07i-%07i.fits'%(tag, tag, mc_source_id_start, mc_source_id_start + n - 1) + tbhdu.writeto(filename, overwrite=True) + + # Write simulated catalogs + + # Simulated catalog output needs to match the real data + # https://github.com/sidneymau/simple/blob/master/search_algorithm.py + # https://github.com/sidneymau/simple/blob/master/config.yaml + # /home/s1/kadrlica/projects/y3a2/dsphs/v2/skim/ + # e.g., /home/s1/kadrlica/projects/y3a2/dsphs/v2/skim/y3a2_ngmix_cm_11755.fits # for ii in range(0, len(d.formats)): print '\'%s\': [ , \'%s\'],'%(d.names[ii], d.formats[ii]) - - # See: - # https://github.com/sidneymau/simple/blob/master/search_algorithm.py - # https://github.com/sidneymau/simple/blob/master/config.yaml - # /home/s1/kadrlica/projects/y3a2/dsphs/v2/skim/ , e.g., /home/s1/kadrlica/projects/y3a2/dsphs/v2/skim/y3a2_ngmix_cm_11755.fits - #default_array = np.tile(np.nan, len(mc_source_id_array)) # To recognize that those values are synthetic filler default_array = np.tile(-9999., len(mc_source_id_array)) - - """ - # Column name, data, fits format - # Y3A2 pre-Gold - key_map = {'CM_MAG_ERR_G': [mag_1_error_array, 'D'], - 'CM_MAG_ERR_R': [mag_2_error_array, 'D'], - 'CM_MAG_G': [mag_1_array, 'D'], - 'CM_MAG_R': [mag_2_array, 'D'], - 'CM_T': [default_array, 'D'], - 'CM_T_ERR': [default_array, 'D'], - 'COADD_OBJECT_ID': [coadd_object_id_array, 'K'], - 'DEC': [lat_array, 'D'], - 'FLAGS': [default_array, 'K'], - 'PSF_MAG_ERR_G': [mag_1_error_array, 'D'], - 'PSF_MAG_ERR_R': [mag_2_error_array, 'D'], - 'PSF_MAG_G': [mag_1_array, 'D'], - 'PSF_MAG_R': [mag_2_array, 'D'], - 'RA': [lon_array, 'D'], - 'SEXTRACTOR_FLAGS_G': [np.tile(0, len(mc_source_id_array)), 'I'], - 'SEXTRACTOR_FLAGS_R': [np.tile(0, len(mc_source_id_array)), 'I'], - 'WAVG_MAG_PSF_G': [mag_1_array, 'E'], - 'WAVG_MAG_PSF_R': [mag_2_array, 'E'], - 'WAVG_MAGERR_PSF_G': [mag_1_error_array, 'E'], - 'WAVG_MAGERR_PSF_R': [mag_2_error_array, 'E'], - 'WAVG_SPREAD_MODEL_I': [default_array, 'E'], - 'WAVG_SPREADERR_MODEL_I': [default_array, 'E'], - 'EXT_SFD98_G': [default_array, 'E'], - 'EXT_SFD98_R': [default_array, 'E'], - 'CM_MAG_SFD_G': [mag_1_array, 'D'], - 'CM_MAG_SFD_R': [mag_2_array, 'D'], - 'FLAG_FOOTPRINT': [np.tile(1, len(mc_source_id_array)), 'J'], - 'FLAG_FOREGROUND': [np.tile(0, len(mc_source_id_array)), 'J'], - 'EXTENDED_CLASS_MASH': [np.tile(0, len(mc_source_id_array)), 'K'], - 'PSF_MAG_SFD_G': [mag_1_array, 'D'], - 'PSF_MAG_SFD_R': [mag_2_array, 'D'], - 'WAVG_MAG_PSF_SFD_G': [mag_1_array, 'E'], - 'WAVG_MAG_PSF_SFD_R': [mag_2_array, 'E']} - """ - if config['survey'] == 'des': # Y3 Gold v2.0 key_map = odict([ - ('COADD_OBJECT_ID', [coadd_object_id_array, 'K']), - ('RA', [lon_array, 'D']), - ('DEC', [lat_array, 'D']), + ('COADD_OBJECT_ID', [coadd_object_id_array, 'K']), + ('RA', [lon_array, 'D']), + ('DEC', [lat_array, 'D']), ('SOF_PSF_MAG_CORRECTED_G', [mag_1_array, 'D']), ('SOF_PSF_MAG_CORRECTED_R', [mag_2_array, 'D']), - ('SOF_PSF_MAG_ERR_G', [mag_1_error_array, 'D']), - ('SOF_PSF_MAG_ERR_R', [mag_2_error_array, 'D']), - ('A_SED_SFD98_G', [mag_extinction_1_array, 'E']), - ('A_SED_SFD98_R', [mag_extinction_2_array, 'E']), - ('WAVG_MAG_PSF_G', [mag_1_array+mag_extinction_1_array, 'E']), - ('WAVG_MAG_PSF_R', [mag_2_array+mag_extinction_2_array, 'E']), - ('WAVG_MAGERR_PSF_G', [mag_1_error_array, 'E']), - ('WAVG_MAGERR_PSF_R', [mag_2_error_array, 'E']), - ('WAVG_SPREAD_MODEL_I', [default_array, 'E']), - ('WAVG_SPREADERR_MODEL_I', [default_array, 'E']), - ('SOF_CM_T', [default_array, 'D']), - ('SOF_CM_T_ERR', [default_array, 'D']), - ('FLAGS_GOLD', [np.tile(0, len(mc_source_id_array)), 'J']), + ('SOF_PSF_MAG_ERR_G', [mag_1_error_array, 'D']), + ('SOF_PSF_MAG_ERR_R', [mag_2_error_array, 'D']), + ('A_SED_SFD98_G', [mag_extinction_1_array, 'E']), + ('A_SED_SFD98_R', [mag_extinction_2_array, 'E']), + ('WAVG_MAG_PSF_G', [mag_1_array+mag_extinction_1_array, 'E']), + ('WAVG_MAG_PSF_R', [mag_2_array+mag_extinction_2_array, 'E']), + ('WAVG_MAGERR_PSF_G', [mag_1_error_array, 'E']), + ('WAVG_MAGERR_PSF_R', [mag_2_error_array, 'E']), + ('WAVG_SPREAD_MODEL_I', [default_array, 'E']), + ('WAVG_SPREADERR_MODEL_I', [default_array, 'E']), + ('SOF_CM_T', [default_array, 'D']), + ('SOF_CM_T_ERR', [default_array, 'D']), + ('FLAGS_GOLD', [np.tile(0, len(mc_source_id_array)), 'J']), ('EXTENDED_CLASS_MASH_SOF', [np.tile(0, len(mc_source_id_array)), 'I']), ]) elif config['survey'] == 'ps1': # PS1 key_map = odict([ - ('OBJID', [coadd_object_id_array, 'K']), - ('RA', [lon_array, 'D']), - ('DEC', [lat_array, 'D']), - #('UNIQUEPSPSOBID', [coadd_object_id_array, 'K']), - #('OBJINFOFLAG', [default_array, 'E']), - #('QUALITYFLAG', [np.tile(16, len(mc_source_id_array)), 'I']), - #('NSTACKDETECTIONS', [np.tile(99, len(mc_source_id_array)), 'I']), - #('NDETECTIONS', [np.tile(99, len(mc_source_id_array)), 'I']), - #('NG', [default_array, 'E']), - #('NR', [default_array, 'E']), - #('NI', [default_array, 'E']), - ('GFPSFMAG', [mag_1_array+mag_extinction_1_array, 'E']), - ('RFPSFMAG', [mag_2_array+mag_extinction_2_array, 'E']), - #('IFPSFMAG', [np.tile(0., len(mc_source_id_array)), 'E'], # Too pass star selection - ('GFPSFMAGERR', [mag_1_error_array, 'E']), - ('RFPSFMAGERR', [mag_2_error_array, 'E']), - #('IFPSFMAGERR', [default_array, 'E']), - #('GFKRONMAG', [mag_1_array, 'E']), - #('RFKRONMAG', [mag_2_array, 'E']), - #('IFKRONMAG', [np.tile(0., len(mc_source_id_array)), 'E'], # Too pass star selection - #('GFKRONMAGERR', [mag_1_error_array, 'E']), - #('RFKRONMAGERR', [mag_2_error_array, 'E']), - #('IFKRONMAGERR', [default_array, 'E']), - #('GFLAGS', [np.tile(0, len(mc_source_id_array)), 'I']), - #('RFLAGS', [np.tile(0, len(mc_source_id_array)), 'I']), - #('IFLAGS', [np.tile(0, len(mc_source_id_array)), 'I']), - #('GINFOFLAG', [np.tile(0, len(mc_source_id_array)), 'I']), - #('RINFOFLAG', [np.tile(0, len(mc_source_id_array)), 'I']), - #('IINFOFLAG', [np.tile(0, len(mc_source_id_array)), 'I']), - #('GINFOFLAG2', [np.tile(0, len(mc_source_id_array)), 'I']), - #('RINFOFLAG2', [np.tile(0, len(mc_source_id_array)), 'I']), - #('IINFOFLAG2', [np.tile(0, len(mc_source_id_array)), 'I']), - #('GINFOFLAG3', [np.tile(0, len(mc_source_id_array)), 'I']), - #('RINFOFLAG3', [np.tile(0, len(mc_source_id_array)), 'I']), - #('IINFOFLAG3', [np.tile(0, len(mc_source_id_array)), 'I']), - #('PRIMARYDETECTION', [default_array, 'E']), - #('BESTDETECTION', [default_array, 'E']), - #('EBV', [default_array, 'E']), - #('EXTSFD_G', [mag_extinction_1_array 'E']), - #('EXTSFD_R', [mag_extinction_2_array, 'E']), - #('EXTSFD_I', [default_array, 'E']), - ('GFPSFMAG_SFD', [mag_1_array, 'E']), - ('RFPSFMAG_SFD', [mag_2_array, 'E']), - ('EXTENDED_CLASS', [np.tile(0, len(mc_source_id_array)), 'I']), + ('OBJID', [coadd_object_id_array, 'K']), + ('RA', [lon_array, 'D']), + ('DEC', [lat_array, 'D']), + #('UNIQUEPSPSOBID', [coadd_object_id_array, 'K']), + #('OBJINFOFLAG', [default_array, 'E']), + #('QUALITYFLAG', [np.tile(16, len(mc_source_id_array)), 'I']), + #('NSTACKDETECTIONS', [np.tile(99, len(mc_source_id_array)), 'I']), + #('NDETECTIONS', [np.tile(99, len(mc_source_id_array)), 'I']), + #('NG', [default_array, 'E']), + #('NR', [default_array, 'E']), + #('NI', [default_array, 'E']), + ('GFPSFMAG', [mag_1_array+mag_extinction_1_array, 'E']), + ('RFPSFMAG', [mag_2_array+mag_extinction_2_array, 'E']), + #('IFPSFMAG', [np.tile(0., len(mc_source_id_array)), 'E'], # Pass star selection + ('GFPSFMAGERR', [mag_1_error_array, 'E']), + ('RFPSFMAGERR', [mag_2_error_array, 'E']), + #('IFPSFMAGERR', [default_array, 'E']), + #('GFKRONMAG', [mag_1_array, 'E']), + #('RFKRONMAG', [mag_2_array, 'E']), + #('IFKRONMAG', [np.tile(0., len(mc_source_id_array)), 'E'], # Pass star selection + #('GFKRONMAGERR', [mag_1_error_array, 'E']), + #('RFKRONMAGERR', [mag_2_error_array, 'E']), + #('IFKRONMAGERR', [default_array, 'E']), + #('GFLAGS', [np.tile(0, len(mc_source_id_array)), 'I']), + #('RFLAGS', [np.tile(0, len(mc_source_id_array)), 'I']), + #('IFLAGS', [np.tile(0, len(mc_source_id_array)), 'I']), + #('GINFOFLAG', [np.tile(0, len(mc_source_id_array)), 'I']), + #('RINFOFLAG', [np.tile(0, len(mc_source_id_array)), 'I']), + #('IINFOFLAG', [np.tile(0, len(mc_source_id_array)), 'I']), + #('GINFOFLAG2', [np.tile(0, len(mc_source_id_array)), 'I']), + #('RINFOFLAG2', [np.tile(0, len(mc_source_id_array)), 'I']), + #('IINFOFLAG2', [np.tile(0, len(mc_source_id_array)), 'I']), + #('GINFOFLAG3', [np.tile(0, len(mc_source_id_array)), 'I']), + #('RINFOFLAG3', [np.tile(0, len(mc_source_id_array)), 'I']), + #('IINFOFLAG3', [np.tile(0, len(mc_source_id_array)), 'I']), + #('PRIMARYDETECTION', [default_array, 'E']), + #('BESTDETECTION', [default_array, 'E']), + #('EBV', [default_array, 'E']), + #('EXTSFD_G', [mag_extinction_1_array 'E']), + #('EXTSFD_R', [mag_extinction_2_array, 'E']), + #('EXTSFD_I', [default_array, 'E']), + ('GFPSFMAG_SFD', [mag_1_array, 'E']), + ('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 key_map = odict([ @@ -609,64 +642,22 @@ def catsimPopulation(tag, mc_source_id_start=1, n=5000, n_chunk=100, config='sim key_map['MC_SOURCE_ID'] = [mc_source_id_array, 'K'] print("Writing catalog files...") - columns = [] - for key in key_map: - columns.append(pyfits.Column(name=key, format=key_map[key][1], array=key_map[key][0])) - tbhdu = pyfits.BinTableHDU.from_columns(columns) - tbhdu.header.set('AREA', simulation_area, 'Simulation area (deg^2)') - - for mc_source_id_chunk in np.split(np.arange(mc_source_id_start, mc_source_id_start + n), n / n_chunk): - print(' writing MC_SOURCE_ID values from %i to %i'%(mc_source_id_chunk[0], mc_source_id_chunk[-1])) - cut_chunk = np.in1d(mc_source_id_array, mc_source_id_chunk) + for mc_source_id_chunk in np.split(np.arange(mc_source_id_start, mc_source_id_start + n), n//n_chunk): outfile = '%s/sim_catalog_%s_mc_source_id_%07i-%07i.fits'%(tag, tag, mc_source_id_chunk[0], mc_source_id_chunk[-1]) - header = copy.deepcopy(tbhdu.header) - header.set('IDMIN',mc_source_id_chunk[0], 'Minimum MC_SOURCE_ID') - header.set('IDMAX',mc_source_id_chunk[-1], 'Maximum MC_SOURCE_ID') - pyfits.writeto(outfile, tbhdu.data[cut_chunk], header, clobber=True) - - # Population metadata output file - - print("Writing population metadata file...") - tbhdu = pyfits.BinTableHDU.from_columns([ - pyfits.Column(name='RA', format='E', array=lon_population, unit='deg'), - pyfits.Column(name='DEC', format='E', array=lat_population, unit='deg'), - pyfits.Column(name='DISTANCE', format='E', array=distance_population, unit='kpc'), - pyfits.Column(name='DISTANCE_MODULUS', format='E', array=distance_modulus_population, unit='kpc'), - pyfits.Column(name='STELLAR_MASS', format='E', array=stellar_mass_population, unit='m_solar'), - pyfits.Column(name='R_PHYSICAL', format='E', array=r_physical_population, unit='kpc'), - pyfits.Column(name='N_G22', format='J', array=n_g22_population, unit=''), - pyfits.Column(name='N_G24', format='J', array=n_g24_population, unit=''), - pyfits.Column(name='N_CATALOG', format='J', array=n_catalog_population, unit=''), - pyfits.Column(name='DIFFICULTY', format='J', array=difficulty_population, unit=''), - pyfits.Column(name='ABS_MAG', format='E', array=abs_mag_population, unit='mag'), - pyfits.Column(name='SURFACE_BRIGHTNESS', format='E', array=surface_brightness_population, unit='mag arcsec^-2'), - pyfits.Column(name='ELLIPTICITY', format='E', array=ellipticity_population, unit=''), - pyfits.Column(name='POSITION_ANGLE', format='E', array=position_angle_population, unit='deg'), - pyfits.Column(name='AGE', format='E', array=age_population, unit='deg'), - pyfits.Column(name='METAL_Z', format='E', array=metal_z_population, unit=''), - pyfits.Column(name='MC_SOURCE_ID', format='K', array=mc_source_id_population, unit=''), - pyfits.Column(name='HPIX_32', format='E', array=hpix_32_population, unit=''), - pyfits.Column(name='DENSITY', format='E', array=density_population, unit='arcmin^-2'), - pyfits.Column(name='FRACDET_HALF', format='E', array=fracdet_half_population, unit=''), - pyfits.Column(name='FRACDET_CORE', format='E', array=fracdet_core_population, unit=''), - pyfits.Column(name='FRACDET_WIDE', format='E', array=fracdet_wide_population, unit=''), - pyfits.Column(name='MAGLIM_G', format='E', array=maglim_g_population, unit='mag'), - pyfits.Column(name='MAGLIM_R', format='E', array=maglim_r_population, unit='mag'), - pyfits.Column(name='EBV', format='E', array=ebv_population, unit='mag'), - pyfits.Column(name='SURVEY', format='A12', array=survey_population, unit=''), - ]) - tbhdu.header.set('AREA', simulation_area, 'Simulation area (deg^2)') - tbhdu.writeto('%s/sim_population_%s_mc_source_id_%07i-%07i.fits'%(tag, tag, mc_source_id_start, mc_source_id_start + n - 1), clobber=True) - - # 5284.2452461023322 + print(' '+outfile) + sel = np.in1d(mc_source_id_array, mc_source_id_chunk) + columns = [pyfits.Column(name=k, format=v[1], array=v[0][sel]) for k,v in key_map.items()] + tbhdu = pyfits.BinTableHDU.from_columns(columns) + tbhdu.header.set('AREA', simulation_area, 'Simulation area (deg^2)') + tbhdu.header.set('IDMIN',mc_source_id_chunk[0], 'Minimum MC_SOURCE_ID') + tbhdu.header.set('IDMAX',mc_source_id_chunk[-1], 'Maximum MC_SOURCE_ID') + tbhdu.writeto(outfile, overwrite=True) # Mask output file - print("Writing population mask file...") - outfile_mask = '%s/sim_mask_%s_cel_nside_%i.fits'%(tag, tag, healpy.npix2nside(len(mask))) + outfile_mask = '%s/sim_mask_%s_cel_nside_%i.fits'%(tag, tag, hp.npix2nside(len(mask))) if not os.path.exists(outfile_mask): - #healpy.write_map(outfile_mask, mask.astype(int), nest=True, coord='C', overwrite=True) - healpy.write_map(outfile_mask, mask.astype(int), nest=True, coord='C') # overwrite no longer appeasr to be a kwarg? + hp.write_map(outfile_mask, mask.astype(int), nest=True, coord='C', overwrite=True) os.system('gzip -f %s'%(outfile_mask)) ############################################################ @@ -688,9 +679,10 @@ def catsimPopulation(tag, mc_source_id_start=1, n=5000, n_chunk=100, config='sim help="Number of MC_SOURCE_ID's per catalog output file") parser.add_argument('--seed', dest='seed', type=int, default=None, help="Random seed") + parser.add_argument('--dwarfs', dest='dwarfs', action='store_true', + help="Simulate from known dwarfs") args = parser.parse_args() - #print args if args.seed is not None: print("Setting random seed: %i"%args.seed) np.random.seed(args.seed) @@ -699,7 +691,7 @@ def catsimPopulation(tag, mc_source_id_start=1, n=5000, n_chunk=100, config='sim config = yaml.load(open(args.config))[args.section] #catsimPopulation(tag, mc_source_id_start=mc_source_id_start, n=n, n_chunk=n_chunk) - catsimPopulation(args.tag, mc_source_id_start=args.mc_source_id_start, n=args.n, n_chunk=args.n_chunk,config=config) + 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 f7ac4c8..4a74752 100644 --- a/ugali/scratch/simulation/simulate_population.yaml +++ b/ugali/scratch/simulation/simulate_population.yaml @@ -1,8 +1,9 @@ des: survey: des - range_distance: [5.0, 1.0e+3] # kpc + range_distance: [5.0, 1.0e+3] # kpc range_stellar_mass: [1.0e+1, 1.0e+6] # Lsun - range_r_physical: [1.0e-3, 2.0] # pc + range_r_physical: [1.0e-3, 2.0] # kpc + known_dwarfs: /home/s1/kadrlica/projects/mw_substructure/y3a2/sim_population/inputs/recovered_satellites_des.csv stellar_density: /home/s1/kadrlica/projects/mw_substructure/y3a2/sim_population/inputs/des_y3a2_stellar_density_map_g_23_cel_nside_128.npy fracdet: /home/s1/kadrlica/projects/mw_substructure/y3a2/sim_population/inputs/y3a2_griz_o.4096_t.32768_coverfoot_EQU_ring.fits maglim_g: /home/s1/kadrlica/projects/mw_substructure/y3a2/sim_population/inputs/y3a2_gold_1.0_cmv02-001_v1_nside4096_ring_g_depth.fits @@ -13,15 +14,20 @@ des: ps1: survey: ps1 - range_distance: [5.0, 1.0e+3] # kpc + range_distance: [5.0, 1.0e+3] # kpc range_stellar_mass: [1.0e+1, 1.0e+6] # Lsun - range_r_physical: [1.0e-3, 2.0] # pc + range_r_physical: [1.0e-3, 2.0] # pc + known_dwarfs: /home/s1/kadrlica/projects/mw_substructure/y3a2/sim_population/inputs/recovered_satellites_ps1.csv stellar_density: /home/s1/kadrlica/projects/mw_substructure/y3a2/sim_population/inputs/ps1_dr1_stellar_density_map_cel_nside_256.npy fracdet: /home/s1/kadrlica/projects/mw_substructure/y3a2/sim_population/inputs/ps1_dr1_fracdet_n2048.fits.gz maglim_g: /home/s1/bechtol/des60.b/projects/mw_substructure/ps1/sim_population/inputs/ps1_nside32_nest_g_maglim.fits maglim_r: /home/s1/bechtol/des60.b/projects/mw_substructure/ps1/sim_population/inputs/ps1_nside32_nest_r_maglim.fits ebv: /home/s1/bechtol/des60.b/projects/mw_substructure/y3a2/sim_population/inputs/ebv_sfd98_fullres_nside_4096_nest_equatorial.fits.gz - completeness: /home/s1/bechtol/des60.b/projects/mw_substructure/ps1/sim_population/inputs/ps1_stellar_classification_summary_r.csv + #ADW: Ring map (should be faster, but not yet tested) + #ebv: /home/s1/kadrlica/projects/mw_substructure/y3a2/sim_population/inputs/ebv_sfd98_fullres_nside_4096_ring_equatorial.fits + #ADW: Completeness with SNR cut + #completeness: /home/s1/bechtol/des60.b/projects/mw_substructure/ps1/sim_population/inputs/ps1_stellar_classification_summary_r.csv + completeness: /home/s1/kadrlica/projects/mw_substructure/y3a2/sim_population/inputs/ps1_stellar_classification_summary_r_no_snr_cut.csv photo_error: /home/s1/bechtol/des60.b/projects/mw_substructure/ps1/sim_population/inputs/ps1_photo_error_model_r.csv kbechtol_des: diff --git a/ugali/simulation/analyzer.py b/ugali/simulation/analyzer.py index 2825a48..1d0f8d3 100755 --- a/ugali/simulation/analyzer.py +++ b/ugali/simulation/analyzer.py @@ -5,6 +5,9 @@ __author__ = "Alex Drlica-Wagner" import copy import os +import time +import resource, psutil +from collections import OrderedDict as odict import numpy import numpy as np @@ -27,9 +30,33 @@ from ugali.utils import mlab # Analysis flags -FLAG_EBV = 0b001 -FLAG_DENSITY = 0b010 -FLAG_FRACDET = 0b100 +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 + +def update_header_flags(filename): + fits = fitsio.FITS(filename,'rw') + for k,v in FLAGS.items(): + fits[1].write_key(k,v) + +# Results dtypes +# FIXME: KERNEL should be removed for next run +DTYPES=[('TS','>f4'), + ('FIT_KERNEL','S18'),('FIT_EXTENSION','>f4'), + ('FIT_MASS','>f4'),('FIT_MASS_ERR','>f4'), + ('FIT_DISTANCE','>f4'),('FIT_DISTANCE_ERR','>f4'),('FLAG','>i4'), + ('RUNTIME','>f4'),('MEMORY','>i8')] + +KB = 1024**1 +MB = 1024**2 +GB = 1024**3 class Analyzer(object): """ @@ -37,12 +64,74 @@ class Analyzer(object): """ def __init__(self, config, catfile=None, popfile=None): self.config = Config(config) - self.catalog = self.read_catalog(catfile) self.population = self.read_population(popfile) - + self.catalog = self.read_catalog(catfile) + self.mlimit = -1 + + def get_memory_usage(self): + """Get the memory usage of this process. + + Parameters + ---------- + """ + process = psutil.Process() + mem = process.memory_info()[0] + return mem + + ## https://stackoverflow.com/a/7669482/4075339 + ## peak memory usage (kilobytes on Linux) + #usage = resource.getrusage(resource.RUSAGE_SELF).ru_maxrss * 1024 + #return usage + + def get_memory_limit(self): + """Get the hard memory limit from LSF. + + Parameters + ---------- + None + + Returns + ------- + mlimit : memory limit (bytes) + """ + rsrc = resource.RLIMIT_AS + soft, hard = resource.getrlimit(rsrc) + if os.getenv('LSB_CG_MEMLIMIT') and os.getenv('LSB_HOSTS'): + # Memory limit per core + memlimit = int(os.getenv('LSB_CG_MEMLIMIT'), 16) + # Number of cores + ncores = len(os.getenv('LSB_HOSTS').split()) + #soft = ncores * memlimit - 100*1024**2 + soft = ncores * memlimit - GB + return soft + + def set_memory_limit(self, mlimit): + """Set the (soft) memory limit through setrlimit. + + Parameters + ---------- + mlimit : soft memory limit (bytes) + + Returns + ------- + soft, hard : memory limits (bytes) + """ + rsrc = resource.RLIMIT_AS + resource.setrlimit(rsrc, (mlimit, mlimit)) + self.mlimit = mlimit + return resource.getrlimit(rsrc) + + #rsrc = resource.RLIMIT_AS + #soft, hard = resource.getrlimit(rsrc) + #resource.setrlimit(rsrc, (mlimit, hard)) + #self.mlimit, hard = resource.getrlimit(rsrc) + #return (self.mlimit,hard) + + def read_population(self, filename=None): if not filename: filename = os.path.join(self.config['simulate']['dirname'],self.config['simulate']['popfile']) + logger.info("Reading population file: %s"%filename) pop = ugali.utils.fileio.read(filename) pop.dtype.names = list(map(str.upper,pop.dtype.names)) return pop @@ -50,130 +139,275 @@ def read_population(self, filename=None): def read_catalog(self, filename=None): if not filename: filename = os.path.join(self.config['simulate']['dirname'],self.config['simulate']['catfile']) + logger.info("Reading catalog file: %s"%filename) catalog = ugali.observation.catalog.Catalog(self.config,filenames=filename) catalog.data = mlab.rec_append_fields(catalog.data, names=['PIX8','PIX4096'], arrs=np.zeros((2,len(catalog.lon)),dtype='>i8')) return catalog - #from memory_profiler import profile - #@profile - def run(self, population=None, catalog=None, outfile=None, mc_source_id=None): - if not population: population = self.population - if not catalog: catalog = self.catalog - - if mc_source_id is not None: - sel = np.in1d(population['MC_SOURCE_ID'],mc_source_id) - if not sel.sum(): - msg = "Requested MC_SOURCE_ID not found: %i"%mc_source_id - logger.warn(msg) - return - else: - # Select only systems that are in the catalog - sel = np.in1d(population['MC_SOURCE_ID'],catalog.mc_source_id) - population = population[sel] - + def read_results(self, filename): + logger.info("Reading results file: %s"%filename) + results = ugali.utils.fileio.read(filename) + return results + + def create_results(self, population=None): + """ Create the results array of len(population) + + Parameters: + ----------- + population : population array (self.population if None) + + Returns: + -------- + results : array of len(population) + """ + if population is None: population = self.population size = len(population) - dtype=[('KERNEL','S18'),('TS','>f4'),('FIT_KERNEL','S18'), - ('FIT_EXTENSION','>f4'),('FIT_MASS','>f4'),('FIT_MASS_ERR','>f4'), - ('FIT_DISTANCE','>f4'),('FIT_DISTANCE_ERR','>f4'),('FLAG','>i8')] - results = np.array(np.nan*np.ones(size),dtype=dtype) - results = recfuncs.merge_arrays([population,results],flatten=True,asrecarray=False,usemask=False) + results = np.array(np.nan*np.ones(size),dtype=DTYPES) + results = recfuncs.merge_arrays([population,results], + flatten=True, + asrecarray=False,usemask=False) results['TS'] = -np.inf - results['FLAG'] = 0 - self.results = results + results['FLAG'] = FLAG_NOPROC + results['MEMORY'] = -1 + + return results + + def write_results(self, filename,**kwargs): + """ Write results array to a file. + + Parameters: + ----------- + filename : output file name + kwargs : arguments passed to fileio.write + + Returns: + -------- + None + """ + ugali.utils.fileio.write(filename,self.results,**kwargs) + update_header_flags(filename) + + def runall(self,outfile=None,mc_source_id=None,rerun=False): + """Run all sources in population. + + Parameters: + ----------- + outfile : file to write output to + mc_source_id : list of sources to process (None is all sources) + + Returns: + -------- + results : processing results + """ + if mc_source_id is None: + mc_source_id = np.unique(self.catalog.mc_source_id) + + # Select only systems that are in the catalog + sel = np.in1d(self.population['MC_SOURCE_ID'],mc_source_id) + + if not sel.sum(): + msg = "Requested MC_SOURCE_IDs not found in population." + raise ValueError(msg) + if os.path.exists(outfile) and rerun: + # read the results from the existing outfile + self.results = self.read_results(outfile) + else: + # create the results + self.results = self.create_results(population=self.population[sel]) + + if not np.in1d(mc_source_id,self.results['MC_SOURCE_ID']).all(): + msg = "Requested MC_SOURCE_IDs not found in results." + raise ValueError(msg) + if outfile: - ugali.utils.fileio.write(outfile,results,clobber=True) - - for i,d in enumerate(results): - params = dict(list(zip(results.dtype.names,d))) - lon,lat = params['RA'],params['DEC'] - distance_modulus = params['DISTANCE_MODULUS'] - mc_source_id = params['MC_SOURCE_ID'] - extension=np.degrees(np.arctan(params['R_PHYSICAL']/params['DISTANCE'])) - logger.info('\n(%i/%i); (id, lon, lat, mod, ext) = (%i, %.2f, %.2f, %.1f, %.3f)'%(i+1,size,mc_source_id,lon,lat,distance_modulus,extension)) - - if params['EBV'] > 0.2: - msg = "High reddening region E(B-V) > 0.2; skipping..." + logger.info("Writing %s..."%outfile) + self.write_results(outfile,clobber=True) + + for i,r in enumerate(self.results): + # Skip if not in mc_source_id list + if self.results[i]['MC_SOURCE_ID'] not in mc_source_id: + msg = "%i skipped."%self.results[i]['MC_SOURCE_ID'] logger.info(msg) - results[i]['FLAG'] |= FLAG_EBV - results[i]['TS'] = np.nan continue - source = ugali.analysis.loglike.createSource(self.config,section='scan',lon=lon,lat=lat) - logger.info("Reading data catalog...") - obs = ugali.analysis.loglike.createObservation(self.config,lon=lon,lat=lat) - - # Select just the simulated target of interest - logger.info("Merging simulated catalog...") - data = catalog.data[catalog.mc_source_id == mc_source_id].copy() - data = np.array(data[list(obs.catalog.data.dtype.names)],dtype=obs.catalog.data.dtype) - obs.catalog = ugali.observation.catalog.Catalog(self.config, data=np.concatenate([obs.catalog.data,data])) + # Rerun in NOPROC or MEM flagged + if (self.results[i]['FLAG'] & (FLAG_NOPROC | FLAG_MEM)) == 0: + msg = "%i already processed."%self.results[i]['MC_SOURCE_ID'] + logger.info(msg) + continue - # Index of closest distance modulus - loglike = ugali.analysis.loglike.LogLikelihood(self.config,obs,source) - grid = ugali.analysis.scan.GridSearch(self.config,loglike) + start_time = time.time() - #if len(loglike.catalog) > 1.0e5: - if len(loglike.catalog) > 5.0e5: - msg = "High stellar density (N_CATALOG = %i); skipping..."%len(loglike.catalog) + try: + self.runone(i) + except MemoryError as e: + msg = "Memory usage exceeded %.3f GB"%(self.mlimit/GB) logger.warn(msg) - results[i]['FLAG'] |= FLAG_DENSITY - results[i]['TS'] = np.nan - continue - - self.grid = grid - self.loglike = self.grid.loglike - - pix = self.loglike.roi.indexTarget(lon,lat) - - # ADW: Should allow fit_distance to float in order to model search procedure? - # Currently we are evaluating at the true distance. - distance_idx = np.fabs(grid.distance_modulus_array-distance_modulus).argmin() - fit_distance = grid.distance_modulus_array[distance_idx] - fit_extension = [0.02, 0.07, 0.15] # half light radii (deg) - - for ext in fit_extension: - grid.loglike.set_params(extension = ext) - try: - grid.search(coords=(lon,lat),distance_modulus=fit_distance) - results[i]['FLAG'] &= ~FLAG_FRACDET - except ValueError as e: - logger.error(str(e)) - results[i]['FLAG'] |= FLAG_FRACDET - continue - mle = grid.mle() + self.results[i]['FLAG'] |= FLAG_MEM + except Exception as e: + logger.error(str(e)) + self.results[i]['FLAG'] |= FLAG_FIT - ts = 2*grid.log_likelihood_sparse_array[distance_idx][pix] - if ts <= results[i]['TS']: - logger.info("No TS increase; continuing...") - continue - - results[i]['FIT_KERNEL'] = grid.loglike.kernel.name - results[i]['TS'] = ts - results[i]['FIT_MASS'] = grid.stellar_mass_conversion*mle['richness'] - results[i]['FIT_DISTANCE'] = fit_distance #mle['distance_modulus'] - results[i]['FIT_EXTENSION'] = ext - - err = grid.err() - richness_err = (err['richness'][1]-err['richness'][0])/2. - results[i]['FIT_MASS_ERR'] = grid.stellar_mass_conversion*richness_err - - distance_modulus_err = (err['distance_modulus'][1]-err['distance_modulus'][0])/2. - results[i]['FIT_DISTANCE_ERR'] = distance_modulus_err + runtime = time.time() - start_time + self.results[i]['MEMORY'] = self.get_memory_usage() + self.results[i]['RUNTIME'] = runtime logger.info("Fit parameter values:") - for d in dtype: - logger.info('\t%s: %s'%(d[0], results[i][d[0]])) + for d in DTYPES: + logger.info('\t%s: %s'%(d[0], self.results[i][d[0]])) + + logger.info("Memory usage: %.3f GB"%(self.get_memory_usage()/GB)) if (i%self.config['simulate']['save'])==0 and outfile: - ugali.utils.fileio.write(outfile,results) - - if outfile: ugali.utils.fileio.write(outfile,results,clobber=True) + logger.info("Writing %s..."%outfile) + self.write_results(outfile,clobber=True) + + if outfile: + logger.info("Writing %s..."%outfile) + self.write_results(outfile,clobber=True) + + return self.results + + #from memory_profiler import profile + #@profile + def runone(self, i): + """ Run one simulation. + + Parameters: + ----------- + i : index of the simulation to run + + Returns: + -------- + results : result array + """ + results = self.results + results[i]['FLAG'] = FLAG_PROC + + params = dict(list(zip(results[i].dtype.names,results[i]))) + size = len(results) + lon,lat = params['RA'],params['DEC'] + distance_modulus = params['DISTANCE_MODULUS'] + mc_source_id = params['MC_SOURCE_ID'] + extension=np.degrees(np.arctan(params['R_PHYSICAL']/params['DISTANCE'])) + + logger.info('\n(%i/%i); (id, lon, lat, mod, ext) = (%i, %.2f, %.2f, %.1f, %.3f)'%(i+1,size,mc_source_id,lon,lat,distance_modulus,extension)) + + if params['EBV'] > 0.2: + msg = "High reddening region; skipping..." + logger.warn(msg) + results[i]['FLAG'] |= FLAG_EBV + #raise Exception(msg) + #results[i]['TS'] = np.nan + #return + + # This links to the parameters in the data scan + section = 'scan' + # This links to the parameters in the simulate section + #section = 'simulate' + + source=ugali.analysis.loglike.createSource(self.config,section=section,lon=lon,lat=lat) + + logger.info("Reading data catalog...") + obs = ugali.analysis.loglike.createObservation(self.config,lon=lon,lat=lat) + + + # Select just the simulated target of interest + logger.info("Merging simulated catalog...") + data = self.catalog.data[self.catalog.mc_source_id == mc_source_id].copy() + data = np.array(data[list(obs.catalog.data.dtype.names)], + dtype=obs.catalog.data.dtype) + obs.catalog = ugali.observation.catalog.Catalog(self.config, data=np.concatenate([obs.catalog.data,data])) + + loglike = ugali.analysis.loglike.LogLikelihood(self.config,obs,source) + + # Mitigate memory overflow issues by cutting objects with + # too many catalog stars + if len(loglike.catalog) > 5e5: # 1e5 + msg = "Large catalog (N_CATALOG = %i)."%len(loglike.catalog) + logger.warn(msg) + results[i]['FLAG'] |= FLAG_NOBJ + + grid = ugali.analysis.scan.GridSearch(self.config,loglike) + self.grid = grid + self.loglike = self.grid.loglike + + pix = self.loglike.roi.indexTarget(lon,lat) + + # ADW: Should fit_distance be free in order to model search procedure? + if self.config['simulate'].get('fit_distance',False): + fit_distance = None + else: + idx = np.fabs(grid.distance_modulus_array-distance_modulus).argmin() + fit_distance = grid.distance_modulus_array[idx] + + try: + grid.search(coords=(lon,lat),distance_modulus=fit_distance) + results[i]['FLAG'] &= ~FLAG_FIT + except ValueError as e: + logger.error(str(e)) + results[i]['FLAG'] |= FLAG_FIT + + mle = grid.mle() + + distance_idx = np.fabs(grid.distance_modulus_array-mle['distance_modulus']).argmin() + + ts = 2*grid.loglike_array[distance_idx][pix] + + results[i]['FIT_KERNEL'] = grid.loglike.kernel.name + results[i]['TS'] = ts + results[i]['FIT_MASS'] = grid.stellar_mass_conversion*mle['richness'] + results[i]['FIT_DISTANCE'] = mle['distance_modulus'] + results[i]['FIT_EXTENSION'] = grid.loglike.kernel.extension + + err = grid.err() + richness_err = (err['richness'][1]-err['richness'][0])/2. + results[i]['FIT_MASS_ERR'] = grid.stellar_mass_conversion*richness_err + + distance_modulus_err = (err['distance_modulus'][1]-err['distance_modulus'][0])/2. + results[i]['FIT_DISTANCE_ERR'] = distance_modulus_err + + """ + fit_extension = [0.02, 0.07, 0.15] # half light radii (deg) + + # ADW: This won't work since we add extension inside the search loop + for ext in fit_extension: + grid.loglike.set_params(extension = ext) + # Fit failures are often due to fracdet = 0 + try: + grid.search(coords=(lon,lat),distance_modulus=fit_distance) + results[i]['FLAG'] &= ~FLAG_FIT + except ValueError as e: + logger.error(str(e)) + results[i]['FLAG'] |= FLAG_FIT + continue + mle = grid.mle() + ts = 2*grid.loglike_array[distance_idx][pix] + if ts <= results[i]['TS']: + logger.info("No TS increase; continuing...") + continue + + results[i]['FIT_KERNEL'] = grid.loglike.kernel.name + results[i]['TS'] = ts + results[i]['FIT_MASS'] = grid.stellar_mass_conversion*mle['richness'] + results[i]['FIT_DISTANCE'] = fit_distance #mle['distance_modulus'] + results[i]['FIT_EXTENSION'] = ext + + err = grid.err() + richness_err = (err['richness'][1]-err['richness'][0])/2. + results[i]['FIT_MASS_ERR'] = grid.stellar_mass_conversion*richness_err + + distance_modulus_err = (err['distance_modulus'][1]-err['distance_modulus'][0])/2. + results[i]['FIT_DISTANCE_ERR'] = distance_modulus_err + """ return results + run = runall + if __name__ == "__main__": import ugali.utils.parser parser = ugali.utils.parser.Parser(description=__doc__) @@ -186,16 +420,30 @@ def run(self, population=None, catalog=None, outfile=None, mc_source_id=None): help='output results file') parser.add_argument('-i','--mc-source-id',default=None,type=int,action='append', help='specific source id to run') - parser.add_force() - parser.add_debug() + parser.add_argument('-m','--mlimit',default=None,type=int, + help='limit memory usage') + parser.add_argument('-r','--rerun',action='store_true', + help='rerun failed jobs') + + #parser.add_force() + #parser.add_debug() parser.add_verbose() args = parser.parse_args() analyzer = Analyzer(args.config,args.catfile,args.popfile) - + + if args.mlimit is not None: + if args.mlimit == 0: + mlimit = analyzer.get_memory_limit() + else: + mlimit = args.mlimit * GB + soft,hard = analyzer.set_memory_limit(mlimit) + logger.info("Setting memory limit to %.3f GB"%(soft/GB)) + if args.mc_source_id is None: basename = os.path.splitext(args.catfile)[0] imin,imax = list(map(int,basename.rsplit('_',1)[-1].split('-'))) args.mc_source_id = np.arange(imin,imax+1) - analyzer.run(outfile=args.outfile,mc_source_id=args.mc_source_id) + analyzer.run(outfile=args.outfile,mc_source_id=args.mc_source_id, + rerun=args.rerun) diff --git a/ugali/simulation/population.py b/ugali/simulation/population.py index 811ec2d..0bc4bf7 100644 --- a/ugali/simulation/population.py +++ b/ugali/simulation/population.py @@ -1,94 +1,110 @@ """ Tool to generate a population of simulated satellite properties. """ - +from collections import OrderedDict as odict import numpy as np import pylab +import pandas as pd import ugali.utils.config import ugali.utils.projector import ugali.utils.skymap import ugali.analysis.kernel import ugali.observation.catalog - -pylab.ion() +import ugali.isochrone ############################################################ +NAMES = ['id','name','lon','lat','distance','stellar_mass','r_physical', + 'ellipticity','position_angle','age','metallicity'] +DTYPE = [('id',int)] + [('name','S8')] + [(n,float) for n in NAMES[2:]] -def satellitePopulation(mask, nside_pix, n, +def satellitePopulation(mask, nside_pix, size, range_distance=[5., 500.], range_stellar_mass=[1.e1, 1.e6], range_r_physical=[1.e-3, 2.], + range_ellipticity=[0.1,0.8], + range_position_angle=[0,180], + choice_age=[10., 12.0, 13.5], + choice_metal=[0.00010,0.00020], plot=False): """ - Create a population of n randomly placed satellites within a survey mask. - Satellites are distributed uniformly in log(distance) (kpc), uniformly in log(stellar_mass) (M_sol), and uniformly in - physical half-light radius log(r_physical) (kpc). The ranges can be set by the user. + Create a population of randomly placed satellites within a + survey mask. Satellites are distributed uniformly in + log(distance/kpc), uniformly in log(stellar_mass/Msun), and + uniformly in physical half-light radius log(r_physical/kpc). The + ranges can be set by the user. - Returns the simulated area (deg^2) as well as the - lon (deg), lat (deg), distance modulus, stellar mass (M_sol), and half-light radius (deg) for each satellite + Returns the simulated area (deg^2) as well as the lon (deg), lat + (deg), distance modulus, stellar mass (M_sol), and half-light + radius (deg) for each satellite Parameters: ----------- - mask : the survey mask of available area + mask : the survey mask of available area nside_pix : coarse resolution npix for avoiding small gaps in survey - n : number of satellites to simulate - + size : number of satellites to simulate + range_distance : heliocentric distance range (kpc) + range_stellar_mass : stellar mass range (Msun) + range_r_physical : projected physical half-light radius (kpc) + range_ellipticity : elliptictiy [0,1] + range_position_angle : position angle (deg) + choice_age : choices for age + choice_metal : choices for metallicity + Returns: -------- - area, lon, lat, distance, stellar_mass, r_physical + area, population : area of sky covered and population of objects """ - - distance = 10**np.random.uniform(np.log10(range_distance[0]), - np.log10(range_distance[1]), - n) - stellar_mass = 10**np.random.uniform(np.log10(range_stellar_mass[0]), - np.log10(range_stellar_mass[1]), - n) + population = np.recarray(size,dtype=DTYPE) + population.fill(np.nan) + population['name'] = '' + + # Source ID + population['id'] = np.arange(size) + + # Distance (kpc) + population['distance'] = 10**np.random.uniform(np.log10(range_distance[0]), + np.log10(range_distance[1]), + size) + # Stellar mass (Msun) + population['stellar_mass'] = 10**np.random.uniform(np.log10(range_stellar_mass[0]), + np.log10(range_stellar_mass[1]), + size) # Physical half-light radius (kpc) - r_physical = 10**np.random.uniform(np.log10(range_r_physical[0]), - np.log10(range_r_physical[1]), - n) + population['r_physical'] = 10**np.random.uniform(np.log10(range_r_physical[0]), + np.log10(range_r_physical[1]), + size) - # Call positions last because while loop has a variable number of calls to np.random (thus not preserving seed information) - lon, lat, simulation_area = ugali.utils.skymap.randomPositions(mask, nside_pix, n=n) + # Ellipticity [e = 1 - (b/a)] with semi-major axis a and semi-minor axis b + # See http://iopscience.iop.org/article/10.3847/1538-4357/833/2/167/pdf + # Based loosely on https://arxiv.org/abs/0805.2945 + population['ellipticity'] = np.random.uniform(range_ellipticity[0], + range_ellipticity[1], + size) - #half_light_radius = np.degrees(np.arcsin(half_light_radius_physical \ - # / ugali.utils.projector.distanceModulusToDistance(distance_modulus))) + # Random position angle (deg) + population['position_angle'] = np.random.uniform(range_position_angle[0], + range_position_angle[1], + size) - # 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))) + # Age (Gyr) + population['age'] = np.random.choice(choice_age,size) + # Metallicity + population['metallicity'] = np.random.choice(choice_metal,size) - 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)') + # Call positions last because while loop has a variable number of + # calls to np.random (thus not preserving seed information) + lon, lat, area = ugali.utils.skymap.randomPositions(mask, nside_pix, n=size) + population['lon'] = lon + population['lat'] = lat - 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, stellar_mass, r_physical + return area, population ############################################################ +# ADW: 2019-09-01 DEPRECATED def satellitePopulationOrig(config, n, range_distance_modulus=[16.5, 24.], range_stellar_mass=[1.e2, 1.e5], @@ -96,12 +112,26 @@ def satellitePopulationOrig(config, n, 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) (M_sol), and uniformly in - log(r_physical) (kpc). The ranges can be set by the user. + 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 - lon (deg), lat (deg), distance modulus, stellar mass (M_sol), and half-light radius (deg) for each satellite + 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 """ if type(config) == str: @@ -162,3 +192,106 @@ def satellitePopulationOrig(config, n, 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) + + stellar_mass,abs_mag = [],[] + for richness in np.logspace(1,8,25): + stellar_mass += [iso.stellar_mass()*richness] + if stellar_mass[-1] < 1e3: + abs_mag += [iso.absolute_magnitude_martin(richness)[0]] + else: + abs_mag += [iso.absolute_magnitude(richness)] + + return abs_mag,stellar_mass + +def knownPopulation(dwarfs, mask, nside_pix, size): + """ Sample parameters from a known population . + + Parameters + ---------- + dwarfs : known dwarfs to sample at + mask : the survey mask of available area + nside_pix : coarse resolution npix for avoiding small gaps in survey + size : number of satellites to simulate; will be broken into + size//len(dwarfs) per dwarf + + Returns + ------- + area, population + """ + # generated from the interpolation function above. + abs_mag_interp = [5.7574, 4.5429, 3.5881, 2.7379, 1.8594, 0.9984, 0.0245, + -0.851, -1.691, -2.495, -3.343, -4.072, -4.801, -5.530, + -6.259, -6.988, -7.718, -8.447, -9.176, -9.905, -10.63, + -11.36, -12.09, -12.82, -13.55][::-1] + + stellar_mass_interp = [2.363705510, 4.626579555, 9.055797468, 17.72529075, + 34.69445217, 67.90890082, 132.9209289, 260.1716878, + 509.2449149, 996.7663490, 1951.012421, 3818.798128, + 7474.693131, 14630.52917, 28636.94603, 56052.29096, + 109713.4910, 214746.8000, 420332.8841, 822735.1162, + 1610373.818, 3152051.957, 6169642.994, 12076100.01, + 23637055.10][::-1] + + if isinstance(dwarfs,str): + # Note that the random seed will be dependent on the order of dwarfs + dwarfs = pd.read_csv(dwarfs).to_records(index=False) + + nsim = size // len(dwarfs) + nlast = nsim + size % len(dwarfs) + + print('Calculating coarse footprint mask...') + coarse_mask = ugali.utils.skymap.coarseFootprint(mask, nside_pix) + + results = [] + for i,dwarf in enumerate(dwarfs): + print(dwarf['name']) + kwargs = dict() + kwargs['range_distance'] = [dwarf['distance'],dwarf['distance']] + r_physical = dwarf['r12']/1000. + kwargs['range_r_physical'] = [r_physical,r_physical] + stellar_mass = np.interp(dwarf['M_V'],abs_mag_interp,stellar_mass_interp) + print("WARNING: Zeroing stellar mass...") + stellar_mass = 0.01 + kwargs['range_stellar_mass'] = [stellar_mass, stellar_mass] + kwargs['range_ellipticity'] = [dwarf['ellipticity'],dwarf['ellipticity']] + kwargs['range_position_angle']=[0.0, 180.0] + kwargs['choice_age']=[10., 12.0, 13.5] + kwargs['choice_metal']=[0.00010,0.00020] + + num = nsim if i < (len(dwarfs)-1) else nlast + area,pop = satellitePopulation(coarse_mask, nside_pix, num, **kwargs) + pop['name'] = dwarf['abbreviation'] + print("WARNING: Fixing location...") + pop['lon'] = dwarf['ra'] + pop['lat'] = dwarf['dec'] + results += [pop] + + 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 6f06727..e8cb4d7 100644 --- a/ugali/utils/batch.py +++ b/ugali/utils/batch.py @@ -2,18 +2,20 @@ """ Cross-platform batch computing interface. """ - +import os import subprocess, subprocess as sub import getpass from collections import OrderedDict as odict from itertools import chain import copy import time +import resource import numpy as np from ugali.utils.logger import logger +GB=1024**3 CLUSTERS = odict([ ('local' ,['local']), ('lsf' ,['lsf','slac','kipac']), @@ -33,7 +35,7 @@ RUNLIMITS = odict([ #Hard limits (None ,'4:00'), # Default value ('express' ,'0:04'), # 0:04 - ('short' ,'0:30'), # 1:00 + ('short' ,'0:59'), # 1:00 ('medium' ,'1:00'), # 48:00 ('long' ,'4:00'), # 120:00 ('xlong' ,'72:00'), # 72:00 @@ -47,16 +49,20 @@ # General queues now only have 4GB of RAM per CPU # To get more memory for a general job, you need to request more cores: # i.e., '-n 2 -R "span[hosts=1]"' for 8GB of RAM -# https://confluence.slac.stanford.edu/display/SCSPub/Batch+Compute+Best+Practices+and+Other+Info +# https://confluence.slac.stanford.edu/x/moNdCw -# These are options for MPI jobs +# These are options for MPI jobs +# '-R "span[ptile=16]"' indicates the number of processors on each +# host that should be allocated to the job MPIOPTS = odict([ (None ,' -R "span[ptile=4]"'), ('local' ,''), ('short' ,''), ('medium' ,''), + ('long' ,''), ('kipac-ibq',' -R "span[ptile=8]"'), - ('bulletmpi',' -R "span[ptile=16]"'), + #('bulletmpi',' -R "span[ptile=16]"'), + ('bulletmpi',''), ]) def factory(queue,**kwargs): @@ -109,6 +115,40 @@ def njobs(self): jobs = self.jobs() return len(jobs.strip().split('\n'))-1 if jobs else 0 + def bfail(self, path): + """ Check logfile(s) for failure. + + Parameters: + ----------- + path : path to logfile(s) + + Returns: + -------- + files : failed logfiles and message + """ + cmd='grep -r "^Exited" %s'%path + out = self.popen(cmd) + stdout = out.communicate()[0] + # There is an extra newline at the end + return stdout.split('\n')[:-1] + + def bcomp(self, path): + """ Return completed logfile(s). + + Parameters: + ----------- + path : path to logfile(s) + + Returns: + -------- + files : completed logfiles and message + """ + cmd='grep -r "^Successfully completed." %s'%path + out = self.popen(cmd) + stdout = out.communicate()[0] + # There is an extra newline at the end + return stdout.split('\n')[:-1] + def throttle(self,max_jobs=None,sleep=60): if max_jobs is None: max_jobs = self.max_jobs if max_jobs is None: return @@ -147,6 +187,9 @@ def submit(self, command, jobname=None, logfile=None, **opts): self.call(cmd) return cmd + def get_memory_limit(): pass + def set_memory_limit(mlimit): pass + class Local(Batch): def __init__(self,**kwargs): super(Local,self).__init__(**kwargs) @@ -202,11 +245,51 @@ def parse_options(self, **opts): #options.update(OPTIONS[options.get('q')]) # User specified options options.update(opts) - if 'n' in list(options.keys()): - #options['a'] = 'mpirun' - options['R'] += self.mpiopts(options.get('q')) + #if 'n' in list(options.keys()): + # #options['a'] = 'mpirun' + # options['R'] += self.mpiopts(options.get('q')) options.setdefault('W',self.runlimit(options.get('q'))) return ''.join('-%s %s '%(k,v) for k,v in options.items()) + + @classmethod + def set_memory_limit(cls, mlimit): + """Set the (soft) memory limit for setrlimit. + + Parameters: + ----------- + mlimit : soft memory limit (bytes) + + Returns: + -------- + soft, hard : memory limits (bytes) + """ + rsrc = resource.RLIMIT_AS + resource.setrlimit(rsrc, (mlimit, mlimit)) + return resource.getrlimit(rsrc) + + @classmethod + def get_memory_limit(cls): + """Get the hard memory limit from LSF. + + Parameters + ---------- + None + + Returns + ------- + mlimit : memory limit (bytes) + """ + rsrc = resource.RLIMIT_AS + soft, hard = resource.getrlimit(rsrc) + if os.getenv('LSB_CG_MEMLIMIT') and os.getenv('LSB_HOSTS'): + # Memory limit per core + memlimit = int(os.getenv('LSB_CG_MEMLIMIT'), 16) + # Number of cores + ncores = len(os.getenv('LSB_HOSTS').split()) + #soft = ncores * memlimit - 100*1024**2 + soft = int( 0.95 * ncores * memlimit ) + return soft,hard + class Slurm(Batch): _defaults = odict([ diff --git a/ugali/utils/binning.py b/ugali/utils/binning.py index a574175..915783d 100644 --- a/ugali/utils/binning.py +++ b/ugali/utils/binning.py @@ -318,8 +318,10 @@ def reverseHistogram(data,bins=None): """ if bins is None: bins = np.arange(data.max()+2) hist, edges = np.histogram(data, bins=bins) - digi = np.digitize(data.flat,bins=np.unique(data)).argsort() - rev = np.hstack( (len(edges), len(edges) + np.cumsum(hist), digi) ) + #digi = np.digitize(data.flat,bins=np.unique(data)).argsort() + #rev = np.hstack((len(edges), len(edges) + np.cumsum(hist), digi)) + rev = np.hstack((len(edges), len(edges)+np.cumsum(hist), + np.digitize(data.flat,bins=np.unique(data)).argsort())) return hist,edges,rev def binnedMedian(x,y,xbins=None): diff --git a/ugali/utils/config.py b/ugali/utils/config.py index 95be93c..6a0d17e 100644 --- a/ugali/utils/config.py +++ b/ugali/utils/config.py @@ -1,3 +1,4 @@ + """ Class for storing and updating config dictionaries. """ @@ -10,13 +11,15 @@ import numpy as np import healpy as hp +import yaml from ugali.utils.logger import logger import ugali.utils.config # To recognize own type from ugali.utils.mlab import isstring -try: import yaml -except ImportError: logger.warning("YAML not found") +#Yaml is a firm dependency +#try: import yaml +#except ImportError: logger.warning("YAML not found") class Config(dict): """ @@ -132,8 +135,10 @@ def _formatFilepaths(self): """ likedir=self['output']['likedir'] self.likefile = join(likedir,self['output']['likefile']) - self.mergefile = join(likedir,self['output']['mergefile']) - self.roifile = join(likedir,self['output']['roifile']) + + mergedir=self['output'].get('mergedir',likedir) + self.mergefile = join(mergedir,self['output']['mergefile']) + self.roifile = join(mergedir,self['output']['roifile']) searchdir=self['output']['searchdir'] self.labelfile = join(searchdir,self['output']['labelfile']) diff --git a/ugali/utils/fileio.py b/ugali/utils/fileio.py index ba484d4..9f1caa1 100644 --- a/ugali/utils/fileio.py +++ b/ugali/utils/fileio.py @@ -15,6 +15,7 @@ import healpy as hp from ugali.utils.logger import logger +import warnings from ugali.utils.mlab import isstring def read(filename,**kwargs): @@ -84,6 +85,49 @@ def add_column(filename,column,formula,force=False): insert_columns(filename,col,force=force) return True +def load_header(kwargs): + """ Load a FITS header with kwargs. + + Parameters + ---------- + kwargs : dict + keyword arguments passed to fitsio.read_header + + Returns + ------- + data : ndarray + fits catalog data + """ + logger.debug("Loading %(filename)s..."%kwargs) + try: + keys = kwargs.pop('keys',None) + hdr = fitsio.read_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) + raise(e) + +def load_headers(filenames,multiproc=False,**kwargs): + """ Load a set of FITS headers with kwargs. """ + filenames = np.atleast_1d(filenames) + logger.debug("Loading %s files..."%len(filenames)) + + kwargs = [dict(filename=f,**kwargs) for f in filenames] + + if multiproc: + from multiprocessing import Pool + processes = multiproc if multiproc > 0 else None + pool = Pool(processes,maxtasksperchild=1) + out = pool.map(load_header,kwargs) + del pool + else: + out = [load_header(kw) for kw in kwargs] + + logger.debug('Concatenating headers...') + return np.rec.fromrecords([o.values() for o in out],names=out[0].keys()) + + def load_file(kwargs): """ Load a FITS file with kwargs. @@ -92,8 +136,13 @@ def load_file(kwargs): Returns: ndarray : fits catalog data """ - logger.debug("Loading %s..."%kwargs['filename']) - return fitsio.read(**kwargs) + logger.debug("Loading %(filename)s..."%kwargs) + try: + return fitsio.read(**kwargs) + except Exception as e: + logger.error("Failed to load file: %(filename)s"%kwargs) + raise(e) + def load_files(filenames,multiproc=False,**kwargs): """ Load a set of FITS files with kwargs. """ @@ -105,8 +154,9 @@ def load_files(filenames,multiproc=False,**kwargs): if multiproc: from multiprocessing import Pool processes = multiproc if multiproc > 0 else None - p = Pool(processes,maxtasksperchild=1) - out = p.map(load_file,kwargs) + pool = Pool(processes,maxtasksperchild=1) + out = pool.map(load_file,kwargs) + del pool else: out = [load_file(kw) for kw in kwargs] diff --git a/ugali/utils/healpix.py b/ugali/utils/healpix.py index bdd368c..d92b4aa 100644 --- a/ugali/utils/healpix.py +++ b/ugali/utils/healpix.py @@ -11,6 +11,7 @@ import healpy import ugali.utils.projector +from ugali.utils import fileio from ugali.utils.logger import logger ############################################################ @@ -181,7 +182,7 @@ def healpixMap(nside, lon, lat, fill_value=0., nest=False): if max_angsep < 10: # More efficient histograming for small regions of sky m = np.tile(fill_value, hp.nside2npix(nside)) - pix_subset = ugali.utils.healpix.angToDisc(nside, lon_median, lat_median, max_angsep, nest=nest) + pix_subset = ang2disc(nside, lon_median, lat_median, max_angsep, nest=nest) bins = np.arange(np.min(pix_subset), np.max(pix_subset) + 1) m_subset = np.histogram(pix, bins=bins - 0.5)[0].astype(float) m[bins[0:-1]] = m_subset @@ -382,7 +383,6 @@ def write_partial_map(filename, data, nside, coord=None, nest=False, None """ import fitsio - import ugali.utils.fileio # ADW: Do we want to make everything uppercase? @@ -400,8 +400,8 @@ def write_partial_map(filename, data, nside, coord=None, nest=False, if header is not None: for k,v in header.items(): fitshdr.add_record({'name':k,'value':v}) - - logger.info("Writing %s"%filename) + # ADW: Should this be a debug? + logger.info("Writing %s..."%filename) fitsio.write(filename,data,extname='PIX_DATA',header=fitshdr,clobber=True) def read_partial_map(filenames, column, fullsky=True, **kwargs): @@ -422,7 +422,6 @@ def read_partial_map(filenames, column, fullsky=True, **kwargs): (nside,pix,map) : pixel array and healpix map (partial or fullsky) """ import fitsio - import ugali.utils.fileio # Make sure that PIXEL is in columns #kwargs['columns'] = ['PIXEL',column] @@ -430,7 +429,7 @@ def read_partial_map(filenames, column, fullsky=True, **kwargs): filenames = np.atleast_1d(filenames) header = fitsio.read_header(filenames[0],ext=kwargs.get('ext',1)) - data = ugali.utils.fileio.load_files(filenames,**kwargs) + data = fileio.load_files(filenames,**kwargs) pix = data['PIXEL'] value = data[column] @@ -467,13 +466,11 @@ def merge_partial_maps(filenames,outfile,**kwargs): None """ import fitsio - import ugali.utils.fileio - filenames = np.atleast_1d(filenames) header = fitsio.read_header(filenames[0],ext=kwargs.get('ext',1)) nside = header['NSIDE'] - data = ugali.utils.fileio.load_files(filenames,**kwargs) + data = fileio.load_files(filenames,**kwargs) pix = data['PIXEL'] ndupes = len(pix) - len(np.unique(pix)) @@ -482,7 +479,8 @@ def merge_partial_maps(filenames,outfile,**kwargs): raise Exception(msg) extname = 'DISTANCE_MODULUS' - distance = ugali.utils.fileio.load_files(filenames,ext=extname)[extname] + kwargs['ext'] = extname + distance = fileio.load_files(filenames,**kwargs)[extname] unique_distance = np.unique(distance) # Check if distance moduli are the same... if np.any(distance[:len(unique_distance)] != unique_distance): @@ -490,11 +488,43 @@ def merge_partial_maps(filenames,outfile,**kwargs): msg += '\n'+str(distance[:len(unique_distance)]) msg += '\n'+str(unique_distance) raise Exception(msg) + del distance - write_partial_map(outfile,data=data,nside=nside,clobber=True) - fitsio.write(outfile,{extname:unique_distance},extname=extname) + # Writing partial maps + if outfile is not None: + write_partial_map(outfile,data=data,nside=nside,clobber=True) + fitsio.write(outfile,{extname:unique_distance},extname=extname) -def merge_likelihood_headers(filenames, outfile): + return nside,data,unique_distance + +def merge_likelihood_headers(filenames, outfile, **kwargs): + """ + Merge header information from likelihood files. + + Parameters: + ----------- + filenames : input filenames + oufile : the merged file to write + + Returns: + -------- + data : the data being written + """ + import fitsio + filenames = np.atleast_1d(filenames) + + ext='PIX_DATA' + nside = fitsio.read_header(filenames[0],ext=ext)['LKDNSIDE'] + + keys = ['LKDPIX','STELLAR','NINSIDE','NANNULUS'] + data_dict = fileio.load_headers(filenames,ext=ext,keys=keys,mulitproc=8) + names = data_dict.dtype.names + data_dict.dtype.names = ['PIXEL' if n=='LKDPIX' else n for n in names] + + write_partial_map(outfile, data_dict, nside) + return data_dict + +def merge_likelihood_headers2(filenames, outfile, **kwargs): """ Merge header information from likelihood files. @@ -593,21 +623,19 @@ def read_map(filename, nest=False, hdu=None, h=False, verbose=True): if (not healpy.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 not nest is None: + if nest is not None: if nest and ordering.startswith('RING'): idx = healpy.nest2ring(nside,np.arange(m.size,dtype=np.int32)) - if verbose: print('Ordering converted to NEST') m = m[idx] - return 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)) m = m[idx] if verbose: print('Ordering converted to RING') - if h: - return m, header - else: - return m + if h: return m, hdr + else: return m + if __name__ == "__main__": import argparse diff --git a/ugali/utils/parser.py b/ugali/utils/parser.py index 3fec993..de4ac18 100644 --- a/ugali/utils/parser.py +++ b/ugali/utils/parser.py @@ -100,8 +100,9 @@ def _parse_coords(self,opts): gal = pix2ang(*opts.hpx) if gal is not None: - opts.coords = [(gal[0],gal[1],radius)] - opts.names = [vars(opts).get('name','')] + dtype = [('lon',float),('lat',float),('radius',float)] + opts.coords = np.array([(gal[0],gal[1],radius)],dtype=dtype) + opts.names = np.array([vars(opts).get('name','')]) else: opts.coords = None opts.names = None diff --git a/ugali/utils/plotting.py b/ugali/utils/plotting.py index 45d4ee3..ff2953b 100644 --- a/ugali/utils/plotting.py +++ b/ugali/utils/plotting.py @@ -359,6 +359,19 @@ def get_galaxies(self,select=None): self.galaxies = catalog.applyCut(cut) return self.galaxies + def get_likelihood(self,select=None): + nside = self.config.params['coords']['nside_merge'] + pixel = ang2pix(nside, self.lon, self.lat) + pixels = np.append([pixel],hp.get_all_neighbours(nside,pixel)) + + filenames = [] + for p in pixels: + f = self.config.mergefile%p + if os.path.exists(f): filenames.append(f) + + return ugali.utils.healpix.merge_partial_maps(filenames,None) + + def drawSmoothCatalog(self, catalog, label=None, **kwargs): ax = plt.gca() ra,dec = catalog.ra_dec @@ -489,13 +502,12 @@ def drawFracdet(self,ax=None, mask=None): def drawTS(self,ax=None, filename=None, zidx=0): if not ax: ax = plt.gca() - if not filename: - #dirname = self.config.params['output2']['searchdir'] - #basename = self.config.params['output2']['mergefile'] - #filename = os.path.join(dirname,basename) - filename = self.config.mergefile - data = fitsio.read(filename) + if filename: + data = fitsio.read(filename) + else: + data = self.get_likelihood()[1] + pixels = data['PIXEL'] values = 2*data['LOG_LIKELIHOOD'] @@ -507,13 +519,7 @@ def drawTS(self,ax=None, filename=None, zidx=0): ts_map[pixels] = values[:,zidx] - #im = hp.gnomview(ts_map,**self.gnom_kwargs) - #hp.graticule(dpar=1,dmer=1,color='0.5',verbose=False) - #plt.close() - #im = ax.imshow(im,origin='bottom') - im = drawHealpixMap(ts_map,self.lon,self.lat,self.radius,coord=self.coord) - try: ax.cax.colorbar(im) except: plt.colorbar(im) ax.annotate("TS",**self.label_kwargs) @@ -543,14 +549,21 @@ def drawSpatial(self, ax=None): def drawCMD(self, ax=None, radius=None, zidx=None): + """ Draw color magnitude diagram with isochrone + + Parameters: + ax : matplotlib axis + radius : selection radius + zidx : distance modulus index + + Returns: + None + """ if not ax: ax = plt.gca() import ugali.isochrone if zidx is not None: - filename = self.config.mergefile - logger.debug("Opening %s..."%filename) - f = fitsio.FITS(filename) - distance_modulus = f[2].read()['DISTANCE_MODULUS'][zidx] + distance_modulus = self.get_likelihood()[2][zidx] iso = ugali.isochrone.Padova(age=12,z=0.0002,mod=distance_modulus) #drawIsochrone(iso,ls='',marker='.',ms=1,c='k') @@ -581,10 +594,7 @@ def drawMembership(self, ax=None, radius=None, zidx=0, mc_source_id=1): if not ax: ax = plt.gca() import ugali.analysis.scan - filename = self.config.mergefile - logger.debug("Opening %s..."%filename) - f = fitsio.FITS(filenaem) - distance_modulus = f[2].read()['DISTANCE_MODULUS'][zidx] + distance_modulus = self.get_likelihood()[2] for ii, name in enumerate(self.config.params['isochrone']['infiles']): logger.info('%s %s'%(ii, name)) @@ -621,13 +631,10 @@ def drawMembership(self, ax=None, radius=None, zidx=0, mc_source_id=1): except: plt.colorbar(sc) def plotDistance(self): - filename = self.config.mergefile - logger.debug("Opening %s..."%filename) - f = fitsio.FITS(filename) - d = f[1].read() + _,d,distances = self.get_likelihood() + pixels,values = d['PIXEL'],2*d['LOG_LIKELIHOOD'] if values.ndim == 1: values = values.reshape(-1,1) - distances = f[2].read()['DISTANCE_MODULUS'] if distances.ndim == 1: distances = distances.reshape(-1,1) ts_map = hp.UNSEEN * np.ones(hp.nside2npix(self.nside)) @@ -657,8 +664,13 @@ def plotDistance(self): for i,val in enumerate(values.T): ax = axes[i] - im = ax.imshow(images[i],origin='bottom',vmin=vmin,vmax=vmax) - ax.cax.colorbar(im) + + #https://github.com/matplotlib/matplotlib/issues/9720/ + im = ax.imshow(images[i].data,origin='bottom',vmin=vmin,vmax=vmax) + try: + ax.cax.colorbar(im) + except TypeError as e: + print(e) #ax.annotate(r"$\mu = %g$"%distances[i],**self.label_kwargs) ax.annotate(r"$d = %.0f$ kpc"%mod2dist(distances[i]),**self.label_kwargs) diff --git a/ugali/utils/projector.py b/ugali/utils/projector.py index 3b8538e..4e6f838 100644 --- a/ugali/utils/projector.py +++ b/ugali/utils/projector.py @@ -38,7 +38,7 @@ def setReference(self, lon_ref, lat_ref, zenithal=False): cos_phi,sin_phi = np.cos(phi),np.sin(phi) cos_theta,sin_theta = np.cos(theta),np.sin(theta) - self.rotation_matrix = np.matrix([ + self.rotation_matrix = np.array([ [cos_psi * cos_phi - cos_theta * sin_phi * sin_psi, cos_psi * sin_phi + cos_theta * cos_phi * sin_psi, sin_psi * sin_theta], @@ -481,16 +481,30 @@ def deg2sr(solid_angle): ############################################################ def distanceToDistanceModulus(distance): + """ Return distance modulus for a given distance (kpc). + + Parameters + ---------- + distance : distance (kpc) + + Returns + ------- + mod : distance modulus """ - Return distance modulus for a given distance (kpc). - """ - return 5. * (np.log10(np.array(distance) * 1.e3) - 1.) + return 5. * (np.log10(np.array(distance)) + 2.) dist2mod = distanceToDistanceModulus def distanceModulusToDistance(distance_modulus): - """ - Return distance (kpc) for a given distance modulus. + """ Return distance (kpc) for a given distance modulus. + + Parameters + ---------- + distance_modulus : distance modulus + + Returns + ------- + distance : distance (kpc) """ return 10**((0.2 * np.array(distance_modulus)) - 2.) diff --git a/ugali/utils/skymap.py b/ugali/utils/skymap.py index c6b3e2b..5937ff3 100644 --- a/ugali/utils/skymap.py +++ b/ugali/utils/skymap.py @@ -1,7 +1,5 @@ """ -Tools for making maps of the sky with healpix. - -ADW 2018-05-07: Largely deprecated? +Tools for making maps of the sky with healpix. Used by simulations. """ import numpy as np @@ -9,6 +7,7 @@ import ugali.utils.projector from ugali.utils.healpix import superpixel,subpixel,ang2pix,pix2ang,query_disc +from ugali.utils.healpix import read_partial_map from ugali.utils.logger import logger from ugali.utils.config import Config @@ -33,7 +32,23 @@ def inFootprint(config, pixels, nside=None): """ Open each valid filename for the set of pixels and determine the set of subpixels with valid data. + + Parameters + ---------- + config : config + Configuration (file or object) + pixels : array or int + List of pixels to create footprint for + nside : int, optional + Healpix nside + + Returns + ------- + inside : array + Boolean array of whether pixel is in footprint """ + logger.info("Calculating survey footprint...") + config = Config(config) nside_catalog = config['coords']['nside_catalog'] nside_likelihood = config['coords']['nside_likelihood'] @@ -52,17 +67,21 @@ def inFootprint(config, pixels, nside=None): catalog_pix = superpixel(pixels,nside,nside_catalog) catalog_pix = np.intersect1d(catalog_pix,catalog_pixels) - for fnames in filenames[catalog_pix]: - logger.debug("Loading %s"%filenames['mask_1']) - #subpix_1,val_1 = ugali.utils.skymap.readSparseHealpixMap(fnames['mask_1'],'MAGLIM',construct_map=False) - _nside,subpix_1,val_1 = ugali.utils.healpix.read_partial_map(fnames['mask_1'],'MAGLIM',fullsky=False) - logger.debug("Loading %s"%fnames['mask_2']) - #subpix_2,val_2 = ugali.utils.skymap.readSparseHealpixMap(fnames['mask_2'],'MAGLIM',construct_map=False) - _nside,subpix_2,val_2 = ugali.utils.healpix.read_partial_map(fnames['mask_2'],'MAGLIM',fullsky=False) - subpix = np.intersect1d(subpix_1,subpix_2) - superpix = np.unique(superpixel(subpix,nside_pixel,nside)) - inside |= np.in1d(pixels, superpix) - + fnames = filenames[catalog_pix] + + # Load the first mask + logger.debug("Loading %s"%fnames['mask_1']) + _nside,subpix1,val1 = read_partial_map(fnames['mask_1'],'MAGLIM', + fullsky=False,multiproc=8) + # Load the second mask + logger.debug("Loading %s"%fnames['mask_2']) + _nside,subpix2,val2 = read_partial_map(fnames['mask_2'],'MAGLIM', + fullsky=False,multiproc=8) + # Run the subpixels + subpix = np.intersect1d(subpix1,subpix2) + superpix = np.unique(superpixel(subpix,nside_pixel,nside)) + inside |= np.in1d(pixels, superpix) + return inside def footprint(config, nside=None): @@ -70,6 +89,7 @@ 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'] @@ -90,21 +110,20 @@ def allSkyCoordinates(nside): lon,lat = pix2ang(nside, np.arange(0, hp.nside2npix(nside))) return lon, lat - -def randomPositions(input, nside_pix, n=1): +def coarseFootprint(input, nside_pix): """ - Generate n random positions within a full HEALPix mask of booleans, or a set of (lon, lat) coordinates. + Generate a coarse healpix mask of booleans from a finer healpix + mask or a set of (lon, lat) coordinates. Parameters: ----------- input : (1) full HEALPix mask of booleans, or (2) a set of (lon, lat) coordinates for catalog objects that define the occupied pixels. - nside_pix : nside_pix is meant to be at coarser resolution than the input mask or catalog object positions + nside_pix : nside_pix is meant to be at coarser (or equivalent) resolution than the input mask or catalog object positions so that gaps from star holes, bleed trails, cosmic rays, etc. are filled in. Returns: -------- lon,lat,area : Return the longitude and latitude of the random positions (deg) and the total area (deg^2). - """ input = np.array(input) if len(input.shape) == 1: @@ -124,6 +143,27 @@ def randomPositions(input, nside_pix, n=1): # Create mask at the coarser resolution mask = np.tile(False, hp.nside2npix(nside_pix)) mask[pix] = True + + return mask + +def randomPositions(input, nside_pix, n=1): + """ + Generate n random positions within a full HEALPix mask of booleans, or a set of (lon, lat) coordinates. + + Parameters: + ----------- + input : (1) full HEALPix mask of booleans, or (2) a set of (lon, lat) coordinates for catalog objects that define the occupied pixels. + nside_pix : nside_pix is meant to be at coarser resolution than the input mask or catalog object positions + so that gaps from star holes, bleed trails, cosmic rays, etc. are filled in. + n : number of random points + + Returns: + -------- + lon,lat,area : Return the longitude and latitude of the random positions (deg) and the total area (deg^2). + + """ + mask = coarseFootprint(input, nside_pix) + area = mask.sum() * hp.nside2pixarea(nside_pix, degrees=True) # Estimate the number of points that need to be thrown based off # coverage fraction of the HEALPix mask