diff --git a/ugali/analysis/kernel.py b/ugali/analysis/kernel.py index e9fa08a..2ee487c 100644 --- a/ugali/analysis/kernel.py +++ b/ugali/analysis/kernel.py @@ -9,11 +9,10 @@ from collections import OrderedDict as odict import copy -import numpy import numpy as np +import healpy as hp import scipy.integrate import scipy.interpolate -import healpy import ugali.utils.projector from ugali.utils.projector import Projector, angsep @@ -69,7 +68,7 @@ def projector(self): else: return Projector(self.lon, self.lat, self.proj) - def integrate(self, rmin=0, rmax=numpy.inf): + def integrate(self, rmin=0, rmax=np.inf): """ Calculate the 2D integral of the 1D surface brightness profile (i.e, the flux) between rmin and rmax (elliptical radii). @@ -84,7 +83,7 @@ def integrate(self, rmin=0, rmax=numpy.inf): integral : Solid angle integral (deg^2) """ if rmin < 0: raise Exception('rmin must be >= 0') - integrand = lambda r: self._pdf(r) * 2*numpy.pi * r + integrand = lambda r: self._pdf(r) * 2*np.pi * r return scipy.integrate.quad(integrand,rmin,rmax,full_output=True,epsabs=0)[0] class ToyKernel(Kernel): @@ -101,7 +100,7 @@ class ToyKernel(Kernel): ]) def _cache(self, name=None): - pixel_area = healpy.nside2pixarea(self.nside,degrees=True) + pixel_area = hp.nside2pixarea(self.nside,degrees=True) vec = ang2vec(self.lon, self.lat) self.pix = query_disc(self.nside,vec,self.extension) self._norm = 1./(len(self.pix)*pixel_area) @@ -190,7 +189,7 @@ def sample_radius(self, n): cdf = np.cumsum(pdf) cdf /= cdf[-1] fn = scipy.interpolate.interp1d(cdf, list(range(0, len(cdf)))) - index = numpy.floor(fn(numpy.random.uniform(size=n))).astype(int) + index = np.floor(fn(np.random.uniform(size=n))).astype(int) return radius[index] def sample_lonlat(self, n): @@ -206,7 +205,7 @@ def sample_lonlat(self, n): radius = self.sample_radius(n) a = radius; b = self.jacobian * radius - t = 2. * np.pi * numpy.random.rand(n) + t = 2. * np.pi * np.random.rand(n) cost,sint = np.cos(t),np.sin(t) phi = np.pi/2. - np.deg2rad(self.theta) cosphi,sinphi = np.cos(phi),np.sin(phi) @@ -326,7 +325,7 @@ class EllipticalPlummer(EllipticalKernel): ]) def _kernel(self, radius): - return 1./(numpy.pi*self.r_h**2 * (1.+(radius/self.r_h)**2)**2) + return 1./(np.pi*self.r_h**2 * (1.+(radius/self.r_h)**2)**2) def _cache(self, name=None): if name in [None,'extension','ellipticity','truncate']: diff --git a/ugali/candidate/associate.py b/ugali/candidate/associate.py index 346e7a4..8e33ea9 100644 --- a/ugali/candidate/associate.py +++ b/ugali/candidate/associate.py @@ -4,10 +4,9 @@ import inspect from collections import OrderedDict as odict -import numpy import numpy as np -import fitsio from numpy.lib.recfunctions import stack_arrays +import fitsio import ugali.utils.projector from ugali.utils.projector import gal2cel, cel2gal @@ -16,7 +15,7 @@ from ugali.utils.shell import get_ugali_dir from ugali.utils.logger import logger -#class Catalog(numpy.recarray): +#class Catalog(np.recarray): # # DATADIR=os.path.join(os.path.split(os.path.abspath(__file__))[0],"../data/catalogs/") # @@ -27,23 +26,23 @@ # ('dec',float), # ('glon',float), # ('glat',float)] -# self = numpy.recarray(0,dtype=dtype).view(cls) +# self = np.recarray(0,dtype=dtype).view(cls) # self._load(filename) # return self # # def __add__(self, other): -# return numpy.concatenate([self,other]) +# return np.concatenate([self,other]) # # def __getitem__(self, key): # """ # Support indexing, slicing and direct access. # """ # try: -# return numpy.recarray.__getitem__(key) +# return np.recarray.__getitem__(key) # except ValueError, message: # if key in self.name: # idx = (self.name == key) -# return numpy.recarray.__getitem__(idx) +# return np.recarray.__getitem__(idx) # else: # raise ValueError(message) # @@ -75,7 +74,7 @@ def __init__(self, filename=None): ('dec',float), ('glon',float), ('glat',float)] - self.data = numpy.recarray(0,dtype=columns) + self.data = np.recarray(0,dtype=columns) self._load(filename) if np.isnan([self.data['glon'],self.data['glat']]).any(): raise ValueError("Incompatible values") @@ -94,7 +93,7 @@ def __getitem__(self, key): def __add__(self, other): ret = SourceCatalog() - ret.data = numpy.concatenate([self.data,other.data]) + ret.data = np.concatenate([self.data,other.data]) return ret def __len__(self): @@ -126,10 +125,10 @@ def _load(self,filename): filename = os.path.join(self.DATADIR,"J_AJ_144_4/NearbyGalaxies2012.dat") self.filename = filename - raw = numpy.genfromtxt(filename,delimiter=[19,3,3,5,3,3,3],usecols=range(7),dtype=['|S19']+6*[float],skip_header=36) + raw = np.genfromtxt(filename,delimiter=[19,3,3,5,3,3,3],usecols=range(7),dtype=['|S19']+6*[float],skip_header=36) self.data.resize(len(raw)) - self.data['name'] = numpy.char.strip(raw['f0']) + self.data['name'] = np.char.strip(raw['f0']) ra = raw[['f1','f2','f3']].view(float).reshape(len(raw),-1) dec = raw[['f4','f5','f6']].view(float).reshape(len(raw),-1) @@ -152,7 +151,7 @@ def _load(self,filename): filename = os.path.join(self.DATADIR,"J_AJ_144_4/NearbyGalaxies.dat") self.filename = filename - raw = numpy.genfromtxt(filename,delimiter=[19,3,3,5,3,3,3],usecols=list(range(7)),dtype=['|S19']+6*[float],skip_header=36) + raw = np.genfromtxt(filename,delimiter=[19,3,3,5,3,3,3],usecols=list(range(7)),dtype=['|S19']+6*[float],skip_header=36) self.data.resize(len(raw)) self.data['name'] = np.char.lstrip(np.char.strip(raw['f0']),'*') @@ -180,7 +179,7 @@ def _load(self, filename): raw = fitsio.read(filename,lower=True) self.data.resize(len(raw)) - self.data['name'] = numpy.char.mod("RedMaPPer %d",raw['mem_match_id']) + self.data['name'] = np.char.mod("RedMaPPer %d",raw['mem_match_id']) self.data['ra'] = raw['ra'] self.data['dec'] = raw['dec'] glon,glat = cel2gal(raw['ra'],raw['dec']) @@ -203,10 +202,10 @@ def _load(self,filename): self.filename = filename kwargs = dict(delimiter=[12,12,3,3,6,5,3,6,8,8,6],dtype=2*['S12']+7*[float],skip_header=72,skip_footer=363) - raw = numpy.genfromtxt(filename,**kwargs) + raw = np.genfromtxt(filename,**kwargs) self.data.resize(len(raw)) - self.data['name'] = numpy.char.strip(raw['f0']) + self.data['name'] = np.char.strip(raw['f0']) ra = raw[['f2','f3','f4']].view(float).reshape(len(raw),-1) dec = raw[['f5','f6','f7']].view(float).reshape(len(raw),-1) @@ -228,18 +227,18 @@ def _load(self,filename): for basename in ['VII_239A/ngcpos.dat','VII_239A/icpos.dat']: filename = os.path.join(self.DATADIR,basename) raw.append(np.genfromtxt(filename,**kwargs)) - raw = numpy.concatenate(raw) + raw = np.concatenate(raw) else: - raw = numpy.genfromtxt(filename,**kwargs) + raw = np.genfromtxt(filename,**kwargs) self.filename = filename # Some entries are missing... - raw['f4'] = numpy.where(numpy.isnan(raw['f4']),0,raw['f4']) - raw['f7'] = numpy.where(numpy.isnan(raw['f7']),0,raw['f7']) + raw['f4'] = np.where(np.isnan(raw['f4']),0,raw['f4']) + raw['f7'] = np.where(np.isnan(raw['f7']),0,raw['f7']) self.data.resize(len(raw)) - names = numpy.where(raw['f0'] == 'N', 'NGC %04i', 'IC %04i') - self.data['name'] = numpy.char.mod(names,raw['f1']) + names = np.where(raw['f0'] == 'N', 'NGC %04i', 'IC %04i') + self.data['name'] = np.char.mod(names,raw['f1']) ra = raw[['f2','f3','f4']].view(float).reshape(len(raw),-1) dec = raw[['f5','f6','f7']].view(float).reshape(len(raw),-1) @@ -258,16 +257,16 @@ def _load(self,filename): # if filename is None: # filename = os.path.join(self.DATADIR,"NI2013.csv") # -# raw = numpy.genfromtxt(filename,delimiter=',',usecols=[5,6]+range(13,20),dtype=['S1',int]+3*[float]+['S1']+3*[float]) +# raw = np.genfromtxt(filename,delimiter=',',usecols=[5,6]+range(13,20),dtype=['S1',int]+3*[float]+['S1']+3*[float]) # # self.data.resize(len(raw)) -# names = numpy.where(raw['f0'] == 'N', 'NGC %04i', 'IC %04i') -# self.data['name'] = numpy.char.mod(names,raw['f1']) +# names = np.where(raw['f0'] == 'N', 'NGC %04i', 'IC %04i') +# self.data['name'] = np.char.mod(names,raw['f1']) # -# sign = numpy.where(raw['f5'] == '-',-1,1) +# sign = np.where(raw['f5'] == '-',-1,1) # ra = raw[['f2','f3','f4']].view(float).reshape(len(raw),-1) # dec = raw[['f6','f7','f8']].view(float).reshape(len(raw),-1) -# dec[:,0] = numpy.copysign(dec[:,0], sign) +# dec[:,0] = np.copysign(dec[:,0], sign) # # self.data['ra'] = ugali.utils.projector.hms2dec(ra) # self.data['dec'] = ugali.utils.projector.dms2dec(dec) @@ -288,7 +287,7 @@ def _load(self,filename): raw = np.genfromtxt(filename,delimiter=[3,7,2,4,3,2],dtype=['S3']+['S7']+4*[float]) self.data.resize(len(raw)) - self.data['name'] = numpy.char.mod('UGC %s',numpy.char.strip(raw['f1'])) + self.data['name'] = np.char.mod('UGC %s',np.char.strip(raw['f1'])) ra = raw[['f2','f3']].view(float).reshape(len(raw),-1) ra = np.vstack([ra.T,np.zeros(len(raw))]).T @@ -318,7 +317,7 @@ def _load(self,filename): for basename in ['VII_151/table1a.dat','VII_151/table1c.dat']: filename = os.path.join(self.DATADIR,basename) raw.append(np.genfromtxt(filename,**kwargs)) - raw = numpy.concatenate(raw) + raw = np.concatenate(raw) else: raw = np.genfromtxt(filename,**kwargs) self.filename = filename @@ -355,7 +354,7 @@ def _load(self,filename): raw = np.genfromtxt(filename,**kwargs) self.data.resize(len(raw)) - self.data['name'] = numpy.char.strip(raw['f0']) + self.data['name'] = np.char.strip(raw['f0']) self.data['glon'] = raw['f1'] self.data['glat'] = raw['f2'] @@ -378,7 +377,7 @@ def _load(self,filename): raw = np.genfromtxt(filename,**kwargs) self.data.resize(len(raw)) - self.data['name'] = numpy.char.strip(raw['f0']) + self.data['name'] = np.char.strip(raw['f0']) ra = raw[['f1','f2','f3']].view(float).reshape(len(raw),-1) dec = raw[['f4','f5','f6']].view(float).reshape(len(raw),-1) @@ -402,7 +401,7 @@ def _load(self,filename): raw = np.genfromtxt(filename,**kwargs) self.data.resize(len(raw)) - self.data['name'] = numpy.char.strip(raw['f0']) + self.data['name'] = np.char.strip(raw['f0']) self.data['glon'] = raw['f1'] self.data['glat'] = raw['f2'] diff --git a/ugali/isochrone/model.py b/ugali/isochrone/model.py index a993bcf..36e3cd0 100644 --- a/ugali/isochrone/model.py +++ b/ugali/isochrone/model.py @@ -37,7 +37,6 @@ import glob from functools import wraps -import numpy import numpy as np import scipy.interpolate import scipy.stats @@ -163,7 +162,7 @@ def sample(self, mode='data', mass_steps=1000, mass_min=0.1, full_data_range=Fal # ADW: Assume that the isochrones are pre-sorted by mass_init # This avoids some numerical instability from points that have the same # mass_init value (discontinuities in the isochrone). - # ADW: Might consider using numpy.interp for speed + # ADW: Might consider using np.interp for speed mass_act_interpolation = scipy.interpolate.interp1d(mass_init, mass_act,assume_sorted=True) mag_1_interpolation = scipy.interpolate.interp1d(mass_init, mag_1,assume_sorted=True) mag_2_interpolation = scipy.interpolate.interp1d(mass_init, mag_2,assume_sorted=True) @@ -406,12 +405,12 @@ def sumMag(mag_1, mag_2): # Analytic part mass_init, mass_pdf, mass_act, mag_1, mag_2 = self.sample(mass_steps = steps) g,r = (mag_1,mag_2) if self.band_1 == 'g' else (mag_2,mag_1) - #cut = numpy.logical_not((g > mag_bright) & (g < mag_faint) & (r > mag_bright) & (r < mag_faint)) + #cut = np.logical_not((g > mag_bright) & (g < mag_faint) & (r > mag_bright) & (r < mag_faint)) cut = ((g + self.distance_modulus) > mag_faint) if self.band_1 == 'g' else ((r + self.distance_modulus) > mag_faint) mag_unobs = visual(g[cut], r[cut], richness * mass_pdf[cut]) # Stochastic part - abs_mag_obs_array = numpy.zeros(n_trials) + abs_mag_obs_array = np.zeros(n_trials) for ii in range(0, n_trials): if ii%100==0: logger.debug('%i absolute magnitude trials'%ii) g, r = self.simulate(richness * self.stellar_mass()) @@ -421,11 +420,11 @@ def sumMag(mag_1, mag_2): abs_mag_obs_array[ii] = sumMag(mag_obs, mag_unobs) # ADW: This shouldn't be necessary - #abs_mag_obs_array = numpy.sort(abs_mag_obs_array)[::-1] + #abs_mag_obs_array = np.sort(abs_mag_obs_array)[::-1] # ADW: Careful, fainter abs mag is larger (less negative) number q = [100*alpha/2., 50, 100*(1-alpha/2.)] - hi,med,lo = numpy.percentile(abs_mag_obs_array,q) + hi,med,lo = np.percentile(abs_mag_obs_array,q) return ugali.utils.stats.interval(med,lo,hi) def simulate(self, stellar_mass, distance_modulus=None, **kwargs): diff --git a/ugali/observation/mask.py b/ugali/observation/mask.py index c988bef..1c1577d 100644 --- a/ugali/observation/mask.py +++ b/ugali/observation/mask.py @@ -9,11 +9,9 @@ """ import os -import numpy import numpy as np -import scipy.signal -import healpy import healpy as hp +import scipy.signal #import ugali.utils.plotting import ugali.utils.binning @@ -154,14 +152,14 @@ def _pruneMMD(self, minimum_solid_angle): self.solid_angle_mmd = solid_angle_mmd # Compute which magnitudes the clipping correspond to - index_mag_1, index_mag_2 = numpy.nonzero(self.solid_angle_mmd) + index_mag_1, index_mag_2 = np.nonzero(self.solid_angle_mmd) self.mag_1_clip = self.roi.bins_mag[1:][np.max(index_mag_1)] self.mag_2_clip = self.roi.bins_mag[1:][np.max(index_mag_2)] logger.info('Clipping mask 1 at %.2f mag'%(self.mag_1_clip) ) logger.info('Clipping mask 2 at %.2f mag'%(self.mag_2_clip) ) - self.mask_1.mask_roi_sparse = numpy.clip(self.mask_1.mask_roi_sparse, 0., self.mag_1_clip) - self.mask_2.mask_roi_sparse = numpy.clip(self.mask_2.mask_roi_sparse, 0., self.mag_2_clip) + self.mask_1.mask_roi_sparse = np.clip(self.mask_1.mask_roi_sparse, 0., self.mag_1_clip) + self.mask_2.mask_roi_sparse = np.clip(self.mask_2.mask_roi_sparse, 0., self.mag_2_clip) def _solidAngleCMD(self): """ @@ -169,7 +167,7 @@ def _solidAngleCMD(self): function of color and magnitude. """ - self.solid_angle_cmd = numpy.zeros([len(self.roi.centers_mag), + self.solid_angle_cmd = np.zeros([len(self.roi.centers_mag), len(self.roi.centers_color)]) for index_mag in np.arange(len(self.roi.centers_mag)): @@ -196,7 +194,7 @@ def _solidAngleCMD(self): mag_2 = mag + (0.5 * self.roi.delta_mag) # ADW: Is there a problem here? - #self.solid_angle_cmd[index_mag, index_color] = self.roi.area_pixel * numpy.sum((self.mask_1.mask > mag_1) * (self.mask_2.mask > mag_2)) + #self.solid_angle_cmd[index_mag, index_color] = self.roi.area_pixel * np.sum((self.mask_1.mask > mag_1) * (self.mask_2.mask > mag_2)) # ADW: I think we want to keep pixels that are >= mag unmasked_mag_1 = (self.mask_1.mask_annulus_sparse >= mag_1) @@ -225,7 +223,7 @@ def _solidAngleCMD(self): solid_angle_cmd : 2d array """ - self.solid_angle_cmd = numpy.zeros([len(self.roi.centers_mag), + self.solid_angle_cmd = np.zeros([len(self.roi.centers_mag), len(self.roi.centers_color)]) idx_mag,idx_color=np.where(self.solid_angle_cmd == 0) @@ -287,13 +285,13 @@ def _pruneCMD(self, minimum_solid_angle): else: mag_2 = mag mag_1 = color + mag_2 - self.mag_1_clip = numpy.max(mag_1) + (0.5 * self.roi.delta_color) - self.mag_2_clip = numpy.max(mag_2) + (0.5 * self.roi.delta_mag) + self.mag_1_clip = np.max(mag_1) + (0.5 * self.roi.delta_color) + self.mag_2_clip = np.max(mag_2) + (0.5 * self.roi.delta_mag) logger.info('Clipping mask 1 at %.2f mag'%(self.mag_1_clip) ) logger.info('Clipping mask 2 at %.2f mag'%(self.mag_2_clip) ) - self.mask_1.mask_roi_sparse = numpy.clip(self.mask_1.mask_roi_sparse, 0., self.mag_1_clip) - self.mask_2.mask_roi_sparse = numpy.clip(self.mask_2.mask_roi_sparse, 0., self.mag_2_clip) + self.mask_1.mask_roi_sparse = np.clip(self.mask_1.mask_roi_sparse, 0., self.mag_1_clip) + self.mask_2.mask_roi_sparse = np.clip(self.mask_2.mask_roi_sparse, 0., self.mag_2_clip) def completeness(self, delta, method='step'): @@ -363,39 +361,39 @@ def photo_err_2(delta): # Band 1 mag_1_thresh = self.mask_1.mask_roi_sparse[catalog.pixel_roi_index] - catalog.mag_1 - sorting_indices = numpy.argsort(mag_1_thresh) + sorting_indices = np.argsort(mag_1_thresh) mag_1_thresh_sort = mag_1_thresh[sorting_indices] mag_err_1_sort = catalog.mag_err_1[sorting_indices] - # ADW: Can't this be done with numpy.median(axis=?) + # ADW: Can't this be done with np.median(axis=?) mag_1_thresh_medians = [] mag_err_1_medians = [] for i in range(0, int(len(mag_1_thresh) / float(n_per_bin))): - mag_1_thresh_medians.append(numpy.median(mag_1_thresh_sort[n_per_bin * i: n_per_bin * (i + 1)])) - mag_err_1_medians.append(numpy.median(mag_err_1_sort[n_per_bin * i: n_per_bin * (i + 1)])) + mag_1_thresh_medians.append(np.median(mag_1_thresh_sort[n_per_bin * i: n_per_bin * (i + 1)])) + mag_err_1_medians.append(np.median(mag_err_1_sort[n_per_bin * i: n_per_bin * (i + 1)])) if mag_1_thresh_medians[0] > 0.: - mag_1_thresh_medians = numpy.insert(mag_1_thresh_medians, 0, -99.) - mag_err_1_medians = numpy.insert(mag_err_1_medians, 0, mag_err_1_medians[0]) + mag_1_thresh_medians = np.insert(mag_1_thresh_medians, 0, -99.) + mag_err_1_medians = np.insert(mag_err_1_medians, 0, mag_err_1_medians[0]) photo_err_1 = scipy.interpolate.interp1d(mag_1_thresh_medians, mag_err_1_medians, bounds_error=False, fill_value=mag_err_1_medians[-1]) # Band 2 mag_2_thresh = self.mask_2.mask_roi_sparse[catalog.pixel_roi_index] - catalog.mag_2 - sorting_indices = numpy.argsort(mag_2_thresh) + sorting_indices = np.argsort(mag_2_thresh) mag_2_thresh_sort = mag_2_thresh[sorting_indices] mag_err_2_sort = catalog.mag_err_2[sorting_indices] mag_2_thresh_medians = [] mag_err_2_medians = [] for i in range(0, int(len(mag_2_thresh) / float(n_per_bin))): - mag_2_thresh_medians.append(numpy.median(mag_2_thresh_sort[n_per_bin * i: n_per_bin * (i + 1)])) - mag_err_2_medians.append(numpy.median(mag_err_2_sort[n_per_bin * i: n_per_bin * (i + 1)])) + mag_2_thresh_medians.append(np.median(mag_2_thresh_sort[n_per_bin * i: n_per_bin * (i + 1)])) + mag_err_2_medians.append(np.median(mag_err_2_sort[n_per_bin * i: n_per_bin * (i + 1)])) if mag_2_thresh_medians[0] > 0.: - mag_2_thresh_medians = numpy.insert(mag_2_thresh_medians, 0, -99.) - mag_err_2_medians = numpy.insert(mag_err_2_medians, 0, mag_err_2_medians[0]) + mag_2_thresh_medians = np.insert(mag_2_thresh_medians, 0, -99.) + mag_err_2_medians = np.insert(mag_err_2_medians, 0, mag_err_2_medians[0]) photo_err_2 = scipy.interpolate.interp1d(mag_2_thresh_medians, mag_err_2_medians, bounds_error=False, fill_value=mag_err_2_medians[-1]) @@ -478,8 +476,8 @@ def backgroundMMD(self, catalog, method='cloud-in-cells', weights=None): elif method == 'bootstrap': # Not implemented raise ValueError("Bootstrap method not implemented") - mag_1 + (mag_1_err * numpy.random.normal(0, 1., len(mag_1))) - mag_2 + (mag_2_err * numpy.random.normal(0, 1., len(mag_2))) + mag_1 + (mag_1_err * np.random.normal(0, 1., len(mag_1))) + mag_2 + (mag_2_err * np.random.normal(0, 1., len(mag_2))) elif method == 'histogram': # Apply raw histogram @@ -516,12 +514,12 @@ def backgroundMMD(self, catalog, method='cloud-in-cells', weights=None): # ADW: This accounts for leakage to faint magnitudes # But what about the objects that spill out to red colors?? # Maximum obsevable magnitude index for each color (uses the fact that - # numpy.argmin returns first minimum (zero) instance found. + # np.argmin returns first minimum (zero) instance found. # NOTE: More complicated maps may have holes causing problems observable = (self.solid_angle_mmd > self.minimum_solid_angle) index_mag_1 = observable.argmin(axis=0) - 1 - index_mag_2 = numpy.arange(len(self.roi.centers_mag)) + index_mag_2 = np.arange(len(self.roi.centers_mag)) # Add the cumulative leakage back into the last bin of the CMD leakage = (mmd_background * ~observable).sum(axis=0) ### mmd_background[[index_mag_1,index_mag_2]] += leakage @@ -584,8 +582,8 @@ def backgroundCMD(self, catalog, mode='cloud-in-cells', weights=None): mag_1_array = catalog.mag_1 mag_2_array = catalog.mag_2 - catalog.mag_1 + (catalog.mag_1_err * numpy.random.normal(0, 1., len(catalog.mag_1))) - catalog.mag_2 + (catalog.mag_2_err * numpy.random.normal(0, 1., len(catalog.mag_2))) + catalog.mag_1 + (catalog.mag_1_err * np.random.normal(0, 1., len(catalog.mag_1))) + catalog.mag_2 + (catalog.mag_2_err * np.random.normal(0, 1., len(catalog.mag_2))) elif mode == 'histogram': # Apply raw histogram @@ -621,12 +619,12 @@ def backgroundCMD(self, catalog, mode='cloud-in-cells', weights=None): # ADW: This accounts for leakage to faint magnitudes # But what about the objects that spill out to red colors?? # Maximum obsevable magnitude index for each color (uses the fact that - # numpy.argmin returns first minimum (zero) instance found. + # np.argmin returns first minimum (zero) instance found. # NOTE: More complicated maps may have holes causing problems observable = (self.solid_angle_cmd > self.minimum_solid_angle) index_mag = observable.argmin(axis=0) - 1 - index_color = numpy.arange(len(self.roi.centers_color)) + index_color = np.arange(len(self.roi.centers_color)) # Add the cumulative leakage back into the last bin of the CMD leakage = (cmd_background * ~observable).sum(axis=0) cmd_background[[index_mag,index_color]] += leakage @@ -682,7 +680,7 @@ def restrictCatalogToObservableSpaceMMD(self, catalog): catalog.mag_2, catalog.mag_1, self.roi.bins_mag, self.roi.bins_mag) > 0. - sel = numpy.all([sel_roi,sel_mag_1,sel_mag_2,sel_mmd], axis=0) + sel = np.all([sel_roi,sel_mag_1,sel_mag_2,sel_mmd], axis=0) return sel def restrictCatalogToObservableSpaceCMD(self, catalog): @@ -712,9 +710,9 @@ def restrictCatalogToObservableSpaceCMD(self, catalog): ### # Check that the objects fall in the color-magnitude space of the ROI ### # ADW: I think this is degenerate with the cut_cmd - ### sel_mag = numpy.logical_and(catalog.mag > self.roi.bins_mag[0], + ### sel_mag = np.logical_and(catalog.mag > self.roi.bins_mag[0], ### catalog.mag < self.roi.bins_mag[-1]) - ### sel_color = numpy.logical_and(catalog.color > self.roi.bins_color[0], + ### sel_color = np.logical_and(catalog.color > self.roi.bins_color[0], ### catalog.color < self.roi.bins_color[-1]) # and are observable in the ROI-specific mask for both bands @@ -730,7 +728,7 @@ def restrictCatalogToObservableSpaceCMD(self, catalog): catalog.color, catalog.mag, self.roi.bins_color, self.roi.bins_mag) > 0. - sel = numpy.all([sel_roi,sel_mag_1,sel_mag_2,sel_cmd], axis=0) + sel = np.all([sel_roi,sel_mag_1,sel_mag_2,sel_cmd], axis=0) return sel # FIXME: Need to parallelize CMD and MMD formulation @@ -847,9 +845,9 @@ def plot(self): import ugali.utils.plotting - mask = healpy.UNSEEN * numpy.ones(healpy.nside2npix(self.nside)) + mask = hp.UNSEEN * np.ones(hp.nside2npix(self.nside)) mask[self.roi.pixels] = self.mask_roi_sparse - mask[mask == 0.] = healpy.UNSEEN + mask[mask == 0.] = hp.UNSEEN ugali.utils.plotting.zoomedHealpixMap('Completeness Depth', mask, self.roi.lon, self.roi.lat, @@ -867,7 +865,7 @@ def __init__(self, infiles, roi): """ self.roi = roi mask = ugali.utils.skymap.readSparseHealpixMaps(infiles, field='COVERAGE') - self.nside = healpy.npix2nside(len(mask)) + self.nside = hp.npix2nside(len(mask)) # Sparse maps of pixels in various ROI regions self.mask_roi_sparse = mask[self.roi.pixels] @@ -909,7 +907,7 @@ def __init__(self, maglim, roi): """ self.roi = roi mask = maglim*np.ones(hp.nside2npix(self.roi.config['coords']['nside_pixel'])) - self.nside = healpy.npix2nside(len(mask)) + self.nside = hp.npix2nside(len(mask)) # Sparse maps of pixels in various ROI regions self.mask_roi_sparse = mask[self.roi.pixels] @@ -922,13 +920,13 @@ def simpleMask(config): # De-project the bin centers to get magnitude depths - mesh_x, mesh_y = numpy.meshgrid(roi.centers_x, roi.centers_y) - r = numpy.sqrt(mesh_x**2 + mesh_y**2) # Think about x, y conventions here + mesh_x, mesh_y = np.meshgrid(roi.centers_x, roi.centers_y) + r = np.sqrt(mesh_x**2 + mesh_y**2) # Think about x, y conventions here #z = (0. * (r > 1.)) + (21. * (r < 1.)) #z = 21. - r #z = (21. - r) * (mesh_x > 0.) * (mesh_y < 0.) - z = (21. - r) * numpy.logical_or(mesh_x > 0., mesh_y > 0.) + z = (21. - r) * np.logical_or(mesh_x > 0., mesh_y > 0.) return MaskBand(z, roi) @@ -945,7 +943,7 @@ def readMangleFile(infile, lon, lat, index = None): DeprecationWarning(msg) if index is None: - index = numpy.random.randint(0, 1.e10) + index = np.random.randint(0, 1.e10) coordinate_file = 'temp_coordinate_%010i.dat'%(index) maglim_file = 'temp_maglim_%010i.dat'%(index) @@ -981,7 +979,7 @@ def readMangleFile(infile, lon, lat, index = None): logger.error(msg) break - maglim = numpy.array(maglim) + maglim = np.array(maglim) return maglim ############################################################ @@ -1003,14 +1001,14 @@ def scale(mask, mag_scale, outfile=None): """ msg = "'mask.scale': ADW 2018-05-05" DeprecationWarning(msg) - mask_new = healpy.UNSEEN * numpy.ones(len(mask)) + mask_new = hp.UNSEEN * np.ones(len(mask)) mask_new[mask == 0.] = 0. mask_new[mask > 0.] = mask[mask > 0.] + mag_scale if outfile is not None: - pix = numpy.nonzero(mask_new > 0.)[0] + pix = np.nonzero(mask_new > 0.)[0] data_dict = {'MAGLIM': mask_new[pix]} - nside = healpy.npix2nside(len(mask_new)) + nside = hp.npix2nside(len(mask_new)) ugali.utils.skymap.writeSparseHealpixMap(pix, data_dict, nside, outfile) return mask_new diff --git a/ugali/observation/roi.py b/ugali/observation/roi.py index 955bf1c..734db91 100644 --- a/ugali/observation/roi.py +++ b/ugali/observation/roi.py @@ -11,9 +11,8 @@ """ -import numpy import numpy as np -import healpy +import healpy as hp import ugali.utils.binning import ugali.utils.projector @@ -85,7 +84,7 @@ def __init__(self, config, lon, lat): # Pixels in the outer annulus pix = query_disc(self.config['coords']['nside_pixel'], vec, self.config['coords']['roi_radius_annulus']) - pix = numpy.setdiff1d(self.pixels, pix) + pix = np.setdiff1d(self.pixels, pix) self.pixels_annulus = PixelRegion(self.config['coords']['nside_pixel'],pix) # Pixels within target healpix region @@ -95,18 +94,18 @@ def __init__(self, config, lon, lat): # Boolean arrays for selecting given pixels # (Careful, this works because pixels are pre-sorted by query_disc before in1d) - self.pixel_interior_cut = numpy.in1d(self.pixels, self.pixels_interior) + self.pixel_interior_cut = np.in1d(self.pixels, self.pixels_interior) # ADW: Updated for more general ROI shapes #self.pixel_annulus_cut = ~self.pixel_interior_cut - self.pixel_annulus_cut = numpy.in1d(self.pixels, self.pixels_annulus) + self.pixel_annulus_cut = np.in1d(self.pixels, self.pixels_annulus) # # These should be unnecessary now # self.centers_lon, self.centers_lat = self.pixels.lon, self.pixels.lat # self.centers_lon_interior,self.centers_lat_interior = self.pixels_interior.lon,self.pixels_interior.lat # self.centers_lon_target, self.centers_lat_target = self.pixels_target.lon, self.pixels_target.lat - self.area_pixel = healpy.nside2pixarea(self.config.params['coords']['nside_pixel'],degrees=True) # deg^2 + self.area_pixel = hp.nside2pixarea(self.config.params['coords']['nside_pixel'],degrees=True) # deg^2 """ self.centers_x = self._centers(self.bins_x) @@ -125,11 +124,11 @@ def __init__(self, config, lon, lat): # ADW: These are really bin edges, should be careful and consistent # It would be cleaner to separate the CMD from ROI - self.bins_mag = numpy.linspace(self.config.params['mag']['min'], + self.bins_mag = np.linspace(self.config.params['mag']['min'], self.config.params['mag']['max'], self.config.params['mag']['n_bins'] + 1) - self.bins_color = numpy.linspace(self.config.params['color']['min'], + self.bins_color = np.linspace(self.config.params['color']['min'], self.config.params['color']['max'], self.config.params['color']['n_bins'] + 1) @@ -159,8 +158,8 @@ def plot(self, value=None, pixel=None): # DEPRECATED import ugali.utils.plotting - map_roi = numpy.array(healpy.UNSEEN \ - * numpy.ones(healpy.nside2npix(self.config.params['coords']['nside_pixel']))) + map_roi = np.array(hp.UNSEEN \ + * np.ones(hp.nside2npix(self.config.params['coords']['nside_pixel']))) if value is None: #map_roi[self.pixels] = ugali.utils.projector.angsep(self.lon, self.lat, self.centers_lon, self.centers_lat) @@ -223,9 +222,9 @@ def getCatalogPixels(self): nside_pixel = self.config.params['coords']['nside_pixel'] # All possible catalog pixels spanned by the ROI superpix = ugali.utils.skymap.superpixel(self.pixels,nside_pixel,nside_catalog) - superpix = numpy.unique(superpix) + superpix = np.unique(superpix) # Only catalog pixels that exist in catalog files - pixels = numpy.intersect1d(superpix, filenames['pix'].compressed()) + pixels = np.intersect1d(superpix, filenames['pix'].compressed()) return pixels ############################################################ diff --git a/ugali/preprocess/database.py b/ugali/preprocess/database.py index 9c100ee..13404bf 100755 --- a/ugali/preprocess/database.py +++ b/ugali/preprocess/database.py @@ -1,12 +1,11 @@ #!/usr/bin/env python -import numpy -import numpy as np -import sys -import subprocess +import os, sys import re -import os import io +import subprocess + +import numpy as np try: import http.client as httpcl @@ -303,8 +302,8 @@ def footprint(self,nside): import healpy import ugali.utils.projector if nside > 2**9: raise Exception("Overflow error: nside must be <=2**9") - pix = numpy.arange(healpy.nside2npix(nside),dtype='int') - footprint = numpy.zeros(healpy.nside2npix(nside),dtype='bool') + pix = np.arange(healpy.nside2npix(nside),dtype='int') + footprint = np.zeros(healpy.nside2npix(nside),dtype='bool') ra,dec = ugali.utils.projector.pixToAng(nside,pix) table_name = 'Pix%i'%nside self.upload(np.array([pix,ra,dec]).T, ['pix','ra','dec'], name=table_name) diff --git a/ugali/preprocess/maglims.py b/ugali/preprocess/maglims.py index 16275f1..f4b1038 100644 --- a/ugali/preprocess/maglims.py +++ b/ugali/preprocess/maglims.py @@ -5,17 +5,15 @@ import os from os.path import join import shutil - -import fitsio -import numpy -import numpy as np - -import numpy.lib.recfunctions as recfuncs import tempfile import subprocess -import healpy from collections import Counter from collections import OrderedDict as odict + +import fitsio +import numpy as np +import numpy.lib.recfunctions as recfuncs +import healpy as np from scipy.interpolate import interp1d from scipy.optimize import brentq @@ -98,18 +96,18 @@ def calculate(self, infile, field=1, simple=False): data = fitsio.read(infile,columns=[pixel_pix_name]) - #mask_pixels = numpy.arange( healpy.nside2npix(self.nside_mask), dtype='int') - mask_maglims = numpy.zeros(healpy.nside2npix(self.nside_mask)) + #mask_pixels = np.arange( hp.nside2npix(self.nside_mask), dtype='int') + mask_maglims = np.zeros(hp.nside2npix(self.nside_mask)) - out_pixels = numpy.zeros(0,dtype='int') - out_maglims = numpy.zeros(0) + out_pixels = np.zeros(0,dtype='int') + out_maglims = np.zeros(0) # Find the objects in each pixel pixel_pix = data[pixel_pix_name] mask_pix = ugali.utils.skymap.superpixel(pixel_pix,self.nside_pixel,self.nside_mask) count = Counter(mask_pix) pixels = sorted(count.keys()) - pix_digi = numpy.digitize(mask_pix,pixels).argsort() + pix_digi = np.digitize(mask_pix,pixels).argsort() idx = 0 min_num = 500 signal_to_noise = 10. @@ -132,15 +130,15 @@ def calculate(self, infile, field=1, simple=False): # Estimate the magnitude limit as suggested by: # https://deswiki.cosmology.illinois.edu/confluence/display/DO/SVA1+Release+Document # (https://desweb.cosmology.illinois.edu/confluence/display/Operations/SVA1+Doc) - maglim = numpy.median(mag[(magerr>0.9*magerr_lim)&(magerr<1.1*magerr_lim)]) + maglim = np.median(mag[(magerr>0.9*magerr_lim)&(magerr<1.1*magerr_lim)]) # Alternative method to estimate the magnitude limit by fitting median #mag_min, mag_max = mag.min(),mag.max() - #mag_bins = numpy.arange(mag_min,mag_max,0.1) #0.1086? + #mag_bins = np.arange(mag_min,mag_max,0.1) #0.1086? #x,y = ugali.utils.binning.binnedMedian(mag,magerr,mag_bins) - #x,y = x[~numpy.isnan(y)],y[~numpy.isnan(y)] + #x,y = x[~np.isnan(y)],y[~np.isnan(y)] #magerr_med = interp1d(x,y) - #mag0 = numpy.median(x) + #mag0 = np.median(x) #maglim = brentq(lambda a: magerr_med(a)-magerr_lim,x.min(),x.max(),disp=False) # Median from just objects near magerr cut @@ -148,13 +146,13 @@ def calculate(self, infile, field=1, simple=False): logger.debug("%i (n=%i): maglim=%g"%(pix,num,mask_maglims[pix])) subpix = ugali.utils.skymap.subpixel(pix, self.nside_mask, self.nside_pixel) - maglims = numpy.zeros(len(subpix)) + mask_maglims[pix] - out_pixels = numpy.append(out_pixels,subpix) - out_maglims = numpy.append(out_maglims,maglims) + maglims = np.zeros(len(subpix)) + mask_maglims[pix] + out_pixels = np.append(out_pixels,subpix) + out_maglims = np.append(out_maglims,maglims) # Remove empty pixels logger.info("Removing empty pixels") - idx = numpy.nonzero(out_maglims > 0)[0] + idx = np.nonzero(out_maglims > 0)[0] out_pixels = out_pixels[idx] out_maglims = out_maglims[idx] @@ -164,11 +162,11 @@ def calculate(self, infile, field=1, simple=False): glon,glat = pix2ang(self.nside_pixel,out_pixels) ra,dec = gal2cel(glon,glat) footprint = inFootprint(self.footprint,ra,dec) - idx = numpy.nonzero(footprint)[0] + idx = np.nonzero(footprint)[0] out_pixels = out_pixels[idx] out_maglims = out_maglims[idx] - logger.info("MAGLIM = %.3f +/- %.3f"%(numpy.mean(out_maglims),numpy.std(out_maglims))) + logger.info("MAGLIM = %.3f +/- %.3f"%(np.mean(out_maglims),np.std(out_maglims))) return out_pixels,out_maglims def inFootprint(footprint,ra,dec): @@ -188,9 +186,9 @@ def inFootprint(footprint,ra,dec): try: if isinstance(footprint,str) and os.path.exists(footprint): filename = footprint - #footprint = healpy.read_map(filename,verbose=False) + #footprint = hp.read_map(filename,verbose=False) footprint = fitsio.read(filename)['I'].ravel() - nside = healpy.npix2nside(len(footprint)) + nside = hp.npix2nside(len(footprint)) pix = ang2pix(nside,ra,dec) inside = (footprint[pix] > 0) except IOError: @@ -201,7 +199,7 @@ def inFootprint(footprint,ra,dec): def inMangle(polyfile,ra,dec): coords = tempfile.NamedTemporaryFile(suffix='.txt',delete=False) logger.debug("Writing coordinates to %s"%coords.name) - numpy.savetxt(coords, numpy.array( [ra,dec] ).T, fmt='%.6g' ) + np.savetxt(coords, np.array( [ra,dec] ).T, fmt='%.6g' ) coords.close() weights = tempfile.NamedTemporaryFile(suffix='.txt',delete=False) @@ -214,7 +212,7 @@ def inMangle(polyfile,ra,dec): logger.debug(cmd) subprocess.call(cmd,shell=True) - data = numpy.loadtxt(tmp.name,unpack=True,skiprows=1)[-1] + data = np.loadtxt(tmp.name,unpack=True,skiprows=1)[-1] for f in [coords,weights,tmp]: logger.debug("Removing %s"%f.name) os.remove(f.name) @@ -278,23 +276,23 @@ def split(config,dirname='split',force=False): manglefile_1 = join(mangledir,config['mangle']['filename_1']) logger.info("Reading %s..."%manglefile_1) - mangle1 = healpy.read_map(manglefile_1) + mangle1 = hp.read_map(manglefile_1) manglefile_2 = join(mangledir,config['mangle']['filename_2']) logger.info("Reading %s..."%manglefile_2) - mangle2 = healpy.read_map(manglefile_2) + mangle2 = hp.read_map(manglefile_2) # Read the footprint footfile = config['data']['footprint'] logger.info("Reading %s..."%footfile) - footprint = healpy.read_map(footfile) + footprint = hp.read_map(footfile) # Output mask names mask1 = os.path.basename(config['mask']['basename_1']) mask2 = os.path.basename(config['mask']['basename_2']) for band,mangle,base in [(band1,mangle1,mask1),(band2,mangle2,mask2)]: - nside_mangle = healpy.npix2nside(len(mangle)) + nside_mangle = hp.npix2nside(len(mangle)) if nside_mangle != nside_pixel: msg = "Mask nside different from pixel nside" logger.warning(msg) diff --git a/ugali/preprocess/pixelize.py b/ugali/preprocess/pixelize.py index a28e572..237ae58 100644 --- a/ugali/preprocess/pixelize.py +++ b/ugali/preprocess/pixelize.py @@ -7,15 +7,13 @@ import os from os.path import join +import glob +import collections import fitsio -import numpy import numpy as np import numpy.lib.recfunctions as recfuncs -import collections -import healpy import healpy as hp -import glob from matplotlib import mlab #import ugali.utils.binning @@ -138,14 +136,14 @@ def pixelizeDensity(config, nside=None, force=False): def stellarDensity(infile, nside=256, lon_field='RA', lat_field='DEC'): - area = healpy.nside2pixarea(nside,degrees=True) + area = hp.nside2pixarea(nside,degrees=True) logger.debug("Reading %s"%infile) data = fitsio.read(infile,columns=[lon_field,lat_field]) lon,lat = data[lon_field],data[lat_field] pix = ang2pix(nside,lon,lat) counts = collections.Counter(pix) - pixels, number = numpy.array(sorted(counts.items())).T + pixels, number = np.array(sorted(counts.items())).T density = number/area return pixels, density diff --git a/ugali/simulation/simulator.py b/ugali/simulation/simulator.py index 6e5efa0..d5cc366 100755 --- a/ugali/simulation/simulator.py +++ b/ugali/simulation/simulator.py @@ -6,11 +6,10 @@ import copy import os -import numpy import numpy as np import scipy.interpolate import astropy.io.fits as pyfits -import healpy +import healpy as hp import numpy.lib.recfunctions as recfuncs import fitsio @@ -249,39 +248,39 @@ def _photometricErrors(self, n_per_bin=100, plot=False): # Band 1 mag_1_thresh = self.mask.mask_1.mask_roi_sparse[self.catalog.pixel_roi_index] - self.catalog.mag_1 - sorting_indices = numpy.argsort(mag_1_thresh) + sorting_indices = np.argsort(mag_1_thresh) mag_1_thresh_sort = mag_1_thresh[sorting_indices] mag_err_1_sort = self.catalog.mag_err_1[sorting_indices] - # ADW: Can't this be done with numpy.median(axis=?) + # ADW: Can't this be done with np.median(axis=?) mag_1_thresh_medians = [] mag_err_1_medians = [] for i in range(0, int(len(mag_1_thresh) / float(n_per_bin))): - mag_1_thresh_medians.append(numpy.median(mag_1_thresh_sort[n_per_bin * i: n_per_bin * (i + 1)])) - mag_err_1_medians.append(numpy.median(mag_err_1_sort[n_per_bin * i: n_per_bin * (i + 1)])) + mag_1_thresh_medians.append(np.median(mag_1_thresh_sort[n_per_bin * i: n_per_bin * (i + 1)])) + mag_err_1_medians.append(np.median(mag_err_1_sort[n_per_bin * i: n_per_bin * (i + 1)])) if mag_1_thresh_medians[0] > 0.: - mag_1_thresh_medians = numpy.insert(mag_1_thresh_medians, 0, -99.) - mag_err_1_medians = numpy.insert(mag_err_1_medians, 0, mag_err_1_medians[0]) + mag_1_thresh_medians = np.insert(mag_1_thresh_medians, 0, -99.) + mag_err_1_medians = np.insert(mag_err_1_medians, 0, mag_err_1_medians[0]) self.photo_err_1 = scipy.interpolate.interp1d(mag_1_thresh_medians, mag_err_1_medians, bounds_error=False, fill_value=mag_err_1_medians[-1]) # Band 2 mag_2_thresh = self.mask.mask_2.mask_roi_sparse[self.catalog.pixel_roi_index] - self.catalog.mag_2 - sorting_indices = numpy.argsort(mag_2_thresh) + sorting_indices = np.argsort(mag_2_thresh) mag_2_thresh_sort = mag_2_thresh[sorting_indices] mag_err_2_sort = self.catalog.mag_err_2[sorting_indices] mag_2_thresh_medians = [] mag_err_2_medians = [] for i in range(0, int(len(mag_2_thresh) / float(n_per_bin))): - mag_2_thresh_medians.append(numpy.median(mag_2_thresh_sort[n_per_bin * i: n_per_bin * (i + 1)])) - mag_err_2_medians.append(numpy.median(mag_err_2_sort[n_per_bin * i: n_per_bin * (i + 1)])) + mag_2_thresh_medians.append(np.median(mag_2_thresh_sort[n_per_bin * i: n_per_bin * (i + 1)])) + mag_err_2_medians.append(np.median(mag_err_2_sort[n_per_bin * i: n_per_bin * (i + 1)])) if mag_2_thresh_medians[0] > 0.: - mag_2_thresh_medians = numpy.insert(mag_2_thresh_medians, 0, -99.) - mag_err_2_medians = numpy.insert(mag_err_2_medians, 0, mag_err_2_medians[0]) + mag_2_thresh_medians = np.insert(mag_2_thresh_medians, 0, -99.) + mag_err_2_medians = np.insert(mag_err_2_medians, 0, mag_err_2_medians[0]) self.photo_err_2 = scipy.interpolate.interp1d(mag_2_thresh_medians, mag_err_2_medians, bounds_error=False, fill_value=mag_err_2_medians[-1]) @@ -301,7 +300,7 @@ def _setup_subpix(self,nside=2**16): logger.info("Setup subpixels...") self.nside_pixel = self.config['coords']['nside_pixel'] self.nside_subpixel = self.nside_pixel * 2**4 # Could be config parameter - epsilon = np.degrees(healpy.max_pixrad(self.nside_pixel)) # Pad roi radius to cover edge healpix + epsilon = np.degrees(hp.max_pixrad(self.nside_pixel)) # Pad roi radius to cover edge healpix subpix = ugali.utils.healpix.query_disc(self.nside_subpixel,self.roi.vec,self.roi_radius+epsilon) superpix = ugali.utils.healpix.superpixel(subpix,self.nside_subpixel,self.nside_pixel) self.subpix = subpix[np.in1d(superpix,self.roi.pixels)] @@ -376,10 +375,10 @@ def toy_background(self,mc_source_id=2,seed=None): mag_2 = mag_1 - color # There is probably a better way to do this step without creating the full HEALPix map - mask = -1. * numpy.ones(healpy.nside2npix(self.nside_pixel)) + mask = -1. * np.ones(hp.nside2npix(self.nside_pixel)) mask[self.roi.pixels] = self.mask.mask_1.mask_roi_sparse mag_lim_1 = mask[pix] - mask = -1. * numpy.ones(healpy.nside2npix(self.nside_pixel)) + mask = -1. * np.ones(hp.nside2npix(self.nside_pixel)) mask[self.roi.pixels] = self.mask.mask_2.mask_roi_sparse mag_lim_2 = mask[pix] @@ -387,7 +386,7 @@ def toy_background(self,mc_source_id=2,seed=None): #mag_err_2 = 1.0*np.ones(len(pix)) mag_err_1 = self.photo_err_1(mag_lim_1 - mag_1) mag_err_2 = self.photo_err_2(mag_lim_2 - mag_2) - mc_source_id = mc_source_id * numpy.ones(len(mag_1)) + mc_source_id = mc_source_id * np.ones(len(mag_1)) select = (mag_lim_1>mag_1)&(mag_lim_2>mag_2) @@ -429,7 +428,7 @@ def background(self,mc_source_id=2,seed=None): self._setup_cmd() # Randomize the number of stars per bin according to Poisson distribution - nstar_per_bin = numpy.random.poisson(lam=self.bkg_lambda) + nstar_per_bin = np.random.poisson(lam=self.bkg_lambda) nstar = nstar_per_bin.sum() logger.info("Simulating %i background stars..."%nstar) @@ -443,10 +442,10 @@ def background(self,mc_source_id=2,seed=None): # Distribute points within each color-mag bins xx,yy = np.meshgrid(self.bkg_centers_color,self.bkg_centers_mag) - color = numpy.repeat(xx.flatten(),repeats=nstar_per_bin.flatten()) - color += numpy.random.uniform(-delta_color/2.,delta_color/2.,size=nstar) - mag_1 = numpy.repeat(yy.flatten(),repeats=nstar_per_bin.flatten()) - mag_1 += numpy.random.uniform(-delta_mag/2.,delta_mag/2.,size=nstar) + color = np.repeat(xx.flatten(),repeats=nstar_per_bin.flatten()) + color += np.random.uniform(-delta_color/2.,delta_color/2.,size=nstar) + mag_1 = np.repeat(yy.flatten(),repeats=nstar_per_bin.flatten()) + mag_1 += np.random.uniform(-delta_mag/2.,delta_mag/2.,size=nstar) else: # Uniform color-magnitude distribution logger.info("Generating uniform CMD.") @@ -464,20 +463,20 @@ def background(self,mc_source_id=2,seed=None): pix = ang2pix(nside_pixel, lon, lat) # There is probably a better way to do this step without creating the full HEALPix map - mask = -1. * numpy.ones(healpy.nside2npix(nside_pixel)) + mask = -1. * np.ones(hp.nside2npix(nside_pixel)) mask[self.roi.pixels] = self.mask.mask_1.mask_roi_sparse mag_lim_1 = mask[pix] - mask = -1. * numpy.ones(healpy.nside2npix(nside_pixel)) + mask = -1. * np.ones(hp.nside2npix(nside_pixel)) mask[self.roi.pixels] = self.mask.mask_2.mask_roi_sparse mag_lim_2 = mask[pix] mag_err_1 = self.photo_err_1(mag_lim_1 - mag_1) mag_err_2 = self.photo_err_2(mag_lim_2 - mag_2) - mc_source_id = mc_source_id * numpy.ones(len(mag_1)) + mc_source_id = mc_source_id * np.ones(len(mag_1)) # ADW: Should magnitudes be randomized by the erros? - #mag_1 += (numpy.random.normal(size=len(mag_1)) * mag_err_1) - #mag_2 += (numpy.random.normal(size=len(mag_2)) * mag_err_2) + #mag_1 += (np.random.normal(size=len(mag_1)) * mag_err_1) + #mag_2 += (np.random.normal(size=len(mag_2)) * mag_err_2) select = (mag_lim_1>mag_1)&(mag_lim_2>mag_2) @@ -512,10 +511,10 @@ def satellite(self,stellar_mass,distance_modulus,mc_source_id=1,seed=None,**kwar pix = ang2pix(self.config['coords']['nside_pixel'], lon, lat) # There is probably a better way to do this step without creating the full HEALPix map - mask = -1. * numpy.ones(healpy.nside2npix(self.config['coords']['nside_pixel'])) + mask = -1. * np.ones(hp.nside2npix(self.config['coords']['nside_pixel'])) mask[self.roi.pixels] = self.mask.mask_1.mask_roi_sparse mag_lim_1 = mask[pix] - mask = -1. * numpy.ones(healpy.nside2npix(self.config['coords']['nside_pixel'])) + mask = -1. * np.ones(hp.nside2npix(self.config['coords']['nside_pixel'])) mask[self.roi.pixels] = self.mask.mask_2.mask_roi_sparse mag_lim_2 = mask[pix] @@ -523,12 +522,12 @@ def satellite(self,stellar_mass,distance_modulus,mc_source_id=1,seed=None,**kwar mag_err_2 = self.photo_err_2(mag_lim_2 - mag_2) # Randomize magnitudes by their errors - mag_obs_1 = mag_1+numpy.random.normal(size=len(mag_1))*mag_err_1 - mag_obs_2 = mag_2+numpy.random.normal(size=len(mag_2))*mag_err_2 + mag_obs_1 = mag_1+np.random.normal(size=len(mag_1))*mag_err_1 + mag_obs_2 = mag_2+np.random.normal(size=len(mag_2))*mag_err_2 #mag_obs_1 = mag_1 #mag_obs_2 = mag_2 - #select = numpy.logical_and(mag_obs_1 < mag_lim_1, mag_obs_2 < mag_lim_2) + #select = np.logical_and(mag_obs_1 < mag_lim_1, mag_obs_2 < mag_lim_2) select = (mag_lim_1>mag_obs_1)&(mag_lim_2>mag_obs_2) # Make sure objects lie within the original cmd (should also be done later...) @@ -536,7 +535,7 @@ def satellite(self,stellar_mass,distance_modulus,mc_source_id=1,seed=None,**kwar #return mag_1_obs[cut], mag_2_obs[cut], lon[cut], lat[cut] logger.info("Clipping %i simulated satellite stars..."%(~select).sum()) - mc_source_id = mc_source_id * numpy.ones(len(mag_1)) + mc_source_id = mc_source_id * np.ones(len(mag_1)) hdu = ugali.observation.catalog.makeHDU(self.config,mag_obs_1[select],mag_err_1[select], mag_obs_2[select],mag_err_2[select], @@ -564,10 +563,10 @@ def satellite2(self,stellar_mass,distance_modulus,mc_source_id=1,seed=None,**kwa pix = ang2pix(self.config['coords']['nside_pixel'], lon, lat) # There is probably a better way to do this step without creating the full HEALPix map - mask = -1. * numpy.ones(healpy.nside2npix(self.config['coords']['nside_pixel'])) + mask = -1. * np.ones(hp.nside2npix(self.config['coords']['nside_pixel'])) mask[self.roi.pixels] = self.mask.mask_1.mask_roi_sparse mag_lim_1 = mask[pix] - mask = -1. * numpy.ones(healpy.nside2npix(self.config['coords']['nside_pixel'])) + mask = -1. * np.ones(hp.nside2npix(self.config['coords']['nside_pixel'])) mask[self.roi.pixels] = self.mask.mask_2.mask_roi_sparse mag_lim_2 = mask[pix] @@ -589,10 +588,10 @@ def satellite2(self,stellar_mass,distance_modulus,mc_source_id=1,seed=None,**kwa accept = comp > 1 - np.random.uniform(size=len(mag_1)) # Randomize magnitudes by their errors - mag_obs_1 = mag_1 + (numpy.random.normal(size=len(mag_1))*mag_err_1) - mag_obs_2 = mag_2 + (numpy.random.normal(size=len(mag_2))*mag_err_2) + mag_obs_1 = mag_1 + (np.random.normal(size=len(mag_1))*mag_err_1) + mag_obs_2 = mag_2 + (np.random.normal(size=len(mag_2))*mag_err_2) - #select = numpy.logical_and(mag_obs_1 < mag_lim_1, mag_obs_2 < mag_lim_2) + #select = np.logical_and(mag_obs_1 < mag_lim_1, mag_obs_2 < mag_lim_2) select = (mag_lim_1>mag_obs_1)&(mag_lim_2>mag_obs_2)&accept ### # Make sure objects lie within the original cmd (should also be done later...) @@ -601,7 +600,7 @@ def satellite2(self,stellar_mass,distance_modulus,mc_source_id=1,seed=None,**kwa #return mag_1_obs[cut], mag_2_obs[cut], lon[cut], lat[cut] logger.info("Clipping %i simulated satellite stars..."%(~select).sum()) - mc_source_id = mc_source_id * numpy.ones(len(mag_1)) + mc_source_id = mc_source_id * np.ones(len(mag_1)) hdu = ugali.observation.catalog.makeHDU(self.config,mag_obs_1[select],mag_err_1[select], mag_obs_2[select],mag_err_2[select], diff --git a/ugali/utils/bayesian_efficiency.py b/ugali/utils/bayesian_efficiency.py index 8c538dc..b1de4ae 100644 --- a/ugali/utils/bayesian_efficiency.py +++ b/ugali/utils/bayesian_efficiency.py @@ -1,9 +1,13 @@ """ -Documentation. +Utility functions for calculating efficiency uncertainties in a bayesian manner. + +A good reference is probably: +http://home.fnal.gov/~paterno/images/effic.pdf + +ADW: This should be moved into stats.py """ import scipy.special -import numpy import numpy as np ############################################################ @@ -12,8 +16,8 @@ def gammalnStirling(z): """ Uses Stirling's approximation for the log-gamma function suitable for large arguments. """ - return (0.5 * (numpy.log(2. * numpy.pi) - numpy.log(z))) \ - + (z * (numpy.log(z + (1. / ((12. * z) - (1. / (10. * z))))) - 1.)) + return (0.5 * (np.log(2. * np.pi) - np.log(z))) \ + + (z * (np.log(z + (1. / ((12. * z) - (1. / (10. * z))))) - 1.)) ############################################################ @@ -24,12 +28,12 @@ def confidenceInterval(n, k, alpha = 0.68, errorbar=False): try: e = float(k) / float(n) except ZeroDivisionError: - return numpy.nan, [numpy.nan, numpy.nan] + return np.nan, [np.nan, np.nan] bins = 1000001 dx = 1. / bins - efficiency = numpy.linspace(0, 1, bins) + efficiency = np.linspace(0, 1, bins) # MODIFIED FOR LARGE NUMBERS if n + 2 > 1000: @@ -46,25 +50,25 @@ def confidenceInterval(n, k, alpha = 0.68, errorbar=False): c = scipy.special.gammaln(n - k + 1) if k == 0: - p = numpy.concatenate([[numpy.exp(a - b - c)], - numpy.exp(a - b - c + (k * numpy.log(efficiency[1: -1])) + (n - k) * numpy.log(1. - efficiency[1: -1])), + p = np.concatenate([[np.exp(a - b - c)], + np.exp(a - b - c + (k * np.log(efficiency[1: -1])) + (n - k) * np.log(1. - efficiency[1: -1])), [0.]]) elif k == n: - p = numpy.concatenate([[0.], - numpy.exp(a - b - c + (k * numpy.log(efficiency[1: -1])) + (n - k) * numpy.log(1. - efficiency[1: -1])), - [numpy.exp(a - b - c)]]) + p = np.concatenate([[0.], + np.exp(a - b - c + (k * np.log(efficiency[1: -1])) + (n - k) * np.log(1. - efficiency[1: -1])), + [np.exp(a - b - c)]]) else: - p = numpy.concatenate([[0.], - numpy.exp(a - b - c + (k * numpy.log(efficiency[1: -1])) + (n - k) * numpy.log(1. - efficiency[1: -1])), + p = np.concatenate([[0.], + np.exp(a - b - c + (k * np.log(efficiency[1: -1])) + (n - k) * np.log(1. - efficiency[1: -1])), [0.]]) - i = numpy.argsort(p)[::-1] - p_i = numpy.take(p, i) + i = np.argsort(p)[::-1] + p_i = np.take(p, i) - s = i[numpy.cumsum(p_i * dx) < alpha] + s = i[np.cumsum(p_i * dx) < alpha] - low = min(numpy.min(s) * dx, e) - high = max(numpy.max(s) * dx, e) + low = min(np.min(s) * dx, e) + high = max(np.max(s) * dx, e) if not errorbar: return e, [low, high] @@ -82,7 +86,7 @@ def binomialInterval(n, k, alpha = 0.68): Given n tests and k successes, return efficiency and confidence interval. """ e = float(k)/n - delta_e = 1/float(n) * numpy.sqrt(e * (1 - e) * float(n)) * alpha/0.68 + delta_e = 1/float(n) * np.sqrt(e * (1 - e) * float(n)) * alpha/0.68 return e, [e - delta_e, e + delta_e] ############################################################ diff --git a/ugali/utils/binning.py b/ugali/utils/binning.py index 7fb614d..a574175 100644 --- a/ugali/utils/binning.py +++ b/ugali/utils/binning.py @@ -2,7 +2,6 @@ Tools for binning data. """ -import numpy import numpy as np import collections @@ -28,27 +27,27 @@ def take2D(histogram, x, y, bins_x, bins_y): bins_x : the xbin edges, including upper edge (n-dim) bins_y : the ybin edges, including upper edge (m-dim) """ - histogram = numpy.array(histogram) + histogram = np.array(histogram) - if numpy.isscalar(x): + if np.isscalar(x): x = [x] - if numpy.isscalar(y): + if np.isscalar(y): y = [y] bins_x[-1] += 1.e-10 * (bins_x[-1] - bins_x[-2]) # Numerical stability bins_y[-1] += 1.e-10 * (bins_y[-1] - bins_y[-2]) - #return numpy.take(histogram, (histogram.shape[1] * (numpy.digitize(y, bins_y) - 1)) + (numpy.digitize(x, bins_x) - 1)) + #return np.take(histogram, (histogram.shape[1] * (np.digitize(y, bins_y) - 1)) + (np.digitize(x, bins_x) - 1)) - # Return numpy.nan for entries which are outside the binning range on either axis - index = (histogram.shape[1] * (numpy.digitize(y, bins_y) - 1)) + (numpy.digitize(x, bins_x) - 1) - index_clipped = numpy.clip(index, 0, (histogram.shape[0] * histogram.shape[1]) - 1) - val = numpy.take(histogram, index_clipped) + # Return np.nan for entries which are outside the binning range on either axis + index = (histogram.shape[1] * (np.digitize(y, bins_y) - 1)) + (np.digitize(x, bins_x) - 1) + index_clipped = np.clip(index, 0, (histogram.shape[0] * histogram.shape[1]) - 1) + val = np.take(histogram, index_clipped) - outlier_x = numpy.logical_or(x < bins_x[0], x > bins_x[-1]) - outlier_y = numpy.logical_or(y < bins_y[0], y > bins_y[-1]) - outlier = numpy.logical_or(outlier_x, outlier_y) - val[outlier] = numpy.nan + outlier_x = np.logical_or(x < bins_x[0], x > bins_x[-1]) + outlier_y = np.logical_or(y < bins_y[0], y > bins_y[-1]) + outlier = np.logical_or(outlier_x, outlier_y) + val[outlier] = np.nan return val @@ -73,27 +72,27 @@ def cloudInCells(x, y, bins, weights=None): # For consistency, the variable names should be changed in this function, but low priority... - x_bins = numpy.array(bins[0]) + x_bins = np.array(bins[0]) delta_x = x_bins[1] - x_bins[0] # Overflow and underflow bins - x_bins = numpy.insert(x_bins, 0, x_bins[0] - delta_x) - x_bins = numpy.append(x_bins, x_bins[-1] + delta_x) - y_bins = numpy.array(bins[1]) + x_bins = np.insert(x_bins, 0, x_bins[0] - delta_x) + x_bins = np.append(x_bins, x_bins[-1] + delta_x) + y_bins = np.array(bins[1]) delta_y = y_bins[1] - y_bins[0] - y_bins = numpy.insert(y_bins, 0, y_bins[0] - delta_y) - y_bins = numpy.append(y_bins, y_bins[-1] + delta_y) + y_bins = np.insert(y_bins, 0, y_bins[0] - delta_y) + y_bins = np.append(y_bins, y_bins[-1] + delta_y) - x_bound_cut = numpy.logical_and(x >= x_bins[0], x <= x_bins[-1]) - y_bound_cut = numpy.logical_and(y >= y_bins[0], y <= y_bins[-1]) - bound_cut = numpy.logical_and(x_bound_cut, y_bound_cut) + x_bound_cut = np.logical_and(x >= x_bins[0], x <= x_bins[-1]) + y_bound_cut = np.logical_and(y >= y_bins[0], y <= y_bins[-1]) + bound_cut = np.logical_and(x_bound_cut, y_bound_cut) - if not numpy.any(weights): - bound_weights = numpy.ones(len(x))[bound_cut] + if not np.any(weights): + bound_weights = np.ones(len(x))[bound_cut] else: - bound_weights = numpy.array(weights)[bound_cut] + bound_weights = np.array(weights)[bound_cut] - x_vals = numpy.array(x)[bound_cut] - y_vals = numpy.array(y)[bound_cut] + x_vals = np.array(x)[bound_cut] + y_vals = np.array(y)[bound_cut] x_width = x_bins[1] - x_bins[0] y_width = y_bins[1] - y_bins[0] @@ -101,8 +100,8 @@ def cloudInCells(x, y, bins, weights=None): x_centers = x_bins[0: -1] + (0.5 * x_width) y_centers = y_bins[0: -1] + (0.5 * y_width) - dx = x_vals - x_centers[numpy.digitize(x_vals, x_bins) - 1] - dy = y_vals - y_centers[numpy.digitize(y_vals, y_bins) - 1] + dx = x_vals - x_centers[np.digitize(x_vals, x_bins) - 1] + dy = y_vals - y_centers[np.digitize(y_vals, y_bins) - 1] ux = ((dx / x_width) * (dx >= 0)) +\ ((1. + (dx / x_width)) * (dx < 0)) @@ -133,15 +132,15 @@ def cloudInCells(x, y, bins, weights=None): new_y_vals.append(y_vals - (0.5 * y_width)) cell_weights.append(bound_weights * lx * ly) - new_x_vals = numpy.concatenate(new_x_vals) - new_y_vals = numpy.concatenate(new_y_vals) - cell_weights = numpy.concatenate(cell_weights) + new_x_vals = np.concatenate(new_x_vals) + new_y_vals = np.concatenate(new_y_vals) + cell_weights = np.concatenate(cell_weights) - result = numpy.histogram2d(new_x_vals, new_y_vals, + result = np.histogram2d(new_x_vals, new_y_vals, bins = [x_bins, y_bins], weights = cell_weights)[0] - result = numpy.transpose(result[1: result.shape[0] - 1])[1: result.shape[1] - 1] + result = np.transpose(result[1: result.shape[0] - 1])[1: result.shape[1] - 1] return result, x_bins, y_bins @@ -300,14 +299,14 @@ def sum_chunks(data,fx=4,fy=4): def reverseHistogram(data,bins=None): """ - Bins data using numpy.histogram and calculates the + Bins data using np.histogram and calculates the reverse indices for the entries like IDL. Parameters: - data : data to pass to numpy.histogram - bins : bins to pass to numpy.histogram + data : data to pass to np.histogram + bins : bins to pass to np.histogram Returns: - hist : bin content output by numpy.histogram - edges : edges output from numpy.histogram + hist : bin content output by np.histogram + edges : edges output from np.histogram rev : reverse indices of entries in each bin Using Reverse Indices: h,e,rev = histogram(data, bins=bins) @@ -317,10 +316,10 @@ def reverseHistogram(data,bins=None): indices = rev[ rev[i]:rev[i+1] ] # do calculations with data[indices] ... """ - if bins is None: bins = numpy.arange(data.max()+2) - hist, edges = numpy.histogram(data, bins=bins) - digi = numpy.digitize(data.flat,bins=numpy.unique(data)).argsort() - rev = numpy.hstack( (len(edges), len(edges) + numpy.cumsum(hist), digi) ) + 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) ) return hist,edges,rev def binnedMedian(x,y,xbins=None): @@ -329,5 +328,5 @@ def binnedMedian(x,y,xbins=None): med_y = [] for i,n in enumerate(hist): indices = rev[ rev[i]:rev[i+1] ] - med_y.append( numpy.median(y[indices]) ) - return avg_x, numpy.array(med_y) + med_y.append( np.median(y[indices]) ) + return avg_x, np.array(med_y) diff --git a/ugali/utils/parabola.py b/ugali/utils/parabola.py index 8948be3..24a9c65 100644 --- a/ugali/utils/parabola.py +++ b/ugali/utils/parabola.py @@ -1,11 +1,11 @@ """ Class to construct parabolas from 3 points. -ADW: Need to move all of the plotting stuff +ADW: Need to get rid of all of the plotting stuff ADW: Doesn't this all exist in np.poly? """ -import numpy +import numpy as np import scipy.stats import scipy.interpolate @@ -21,11 +21,11 @@ def __init__(self, x, y): """ # Sort the input - argsort = numpy.argsort(x) - self.x = numpy.array(x)[argsort] - self.y = numpy.array(y)[argsort] + argsort = np.argsort(x) + self.x = np.array(x)[argsort] + self.y = np.array(y)[argsort] - index = numpy.argmax(self.y) + index = np.argmax(self.y) if index == 0: index_0 = 0 index_1 = 1 @@ -47,12 +47,12 @@ def __init__(self, x, y): y_2 = self.y[index_2] # Invert matrix - a = numpy.matrix([[x_0**2, x_0, 1.], + a = np.matrix([[x_0**2, x_0, 1.], [x_1**2, x_1, 1.], [x_2**2, x_2, 1.]]) - a_inverse = numpy.linalg.inv(a) - b = numpy.array([y_0, y_1, y_2]) - p = numpy.dot(numpy.array(a_inverse), b) + a_inverse = np.linalg.inv(a) + b = np.array([y_0, y_1, y_2]) + p = np.dot(np.array(a_inverse), b) self.p_2 = p[0] self.p_1 = p[1] @@ -63,7 +63,7 @@ def __init__(self, x, y): self.vertex_y = self.p_0 - (self.p_1**2 / (4. * self.p_2)) def __eq__(self,other): - return numpy.allclose([self.p_0,self.p_1,self.p_2],[other.p_0,other.p_1,other.p_2]) + return np.allclose([self.p_0,self.p_1,self.p_2],[other.p_0,other.p_1,other.p_2]) def __ne__(self,other): return not self.__eq__(other) @@ -88,21 +88,21 @@ def densify(self, factor=10): y = [] for ii in range(0, len(self.x) - 2): p = Parabola(self.x[ii: ii + 3], self.y[ii: ii + 3]) - x.append(numpy.linspace(self.x[ii], self.x[ii + 1], factor)[0: -1]) + x.append(np.linspace(self.x[ii], self.x[ii + 1], factor)[0: -1]) y.append(p(x[-1])) p = Parabola(self.x[len(self.x) - 3:], self.y[len(self.y) - 3:]) - x.append(numpy.linspace(self.x[-2], self.x[-1], factor)[0: -1]) + x.append(np.linspace(self.x[-2], self.x[-1], factor)[0: -1]) y.append(p(x[-1])) x.append([self.x[-1]]) y.append([self.y[-1]]) - #f = scipy.interpolate.interp1d(numpy.concatenate(x), numpy.concatenate(y)) - #x = numpy.linspace(self.x[0], self.x[-1], len(x) * factor) + #f = scipy.interpolate.interp1d(np.concatenate(x), np.concatenate(y)) + #x = np.linspace(self.x[0], self.x[-1], len(x) * factor) #return x, f(x) - return numpy.concatenate(x), numpy.concatenate(y) + return np.concatenate(x), np.concatenate(y) def profileUpperLimit(self, delta = 2.71): """ @@ -118,28 +118,9 @@ def profileUpperLimit(self, delta = 2.71): if b**2 - 4. * a * c < 0.: print('WARNING') print(a, b, c) - - #pylab.figure() - #pylab.scatter(self.x, self.y) - #raw_input('WAIT') return 0. - - - return max((numpy.sqrt(b**2 - 4. * a * c) - b) / (2. * a), (-1. * numpy.sqrt(b**2 - 4. * a * c) - b) / (2. * a)) - - #def bayesianUpperLimit3(self, alpha, steps = 1.e5): - # """ - # Compute one-sided upper limit using Bayesian Method of Helene. - # """ - # # Need a check to see whether limit is reliable - # pdf = scipy.interpolate.interp1d(self.x, numpy.exp(self.y / 2.)) # Convert from 2 * log(likelihood) to likelihood - # x_pdf = numpy.linspace(self.x[0], self.x[-1], steps) - # cdf = numpy.cumsum(pdf(x_pdf)) - # cdf /= cdf[-1] - # cdf_reflect = scipy.interpolate.interp1d(cdf, x_pdf) - # return cdf_reflect(alpha) - # #return self.x[numpy.argmin((cdf - alpha)**2)] + return max((np.sqrt(b**2 - 4. * a * c) - b) / (2. * a), (-1. * np.sqrt(b**2 - 4. * a * c) - b) / (2. * a)) def bayesianUpperLimit(self, alpha, steps=1.e5, plot=False): """ @@ -147,34 +128,23 @@ def bayesianUpperLimit(self, alpha, steps=1.e5, plot=False): Several methods of increasing numerical stability have been implemented. """ x_dense, y_dense = self.densify() - y_dense -= numpy.max(y_dense) # Numeric stability + y_dense -= np.max(y_dense) # Numeric stability f = scipy.interpolate.interp1d(x_dense, y_dense, kind='linear') - x = numpy.linspace(0., numpy.max(x_dense), steps) - pdf = numpy.exp(f(x) / 2.) - cut = (pdf / numpy.max(pdf)) > 1.e-10 + x = np.linspace(0., np.max(x_dense), steps) + pdf = np.exp(f(x) / 2.) + cut = (pdf / np.max(pdf)) > 1.e-10 x = x[cut] pdf = pdf[cut] #pdf /= pdf[0] - #forbidden = numpy.nonzero(pdf < 1.e-10)[0] + #forbidden = np.nonzero(pdf < 1.e-10)[0] #if len(forbidden) > 0: # index = forbidden[0] # Numeric stability # x = x[0: index] # pdf = pdf[0: index] - cdf = numpy.cumsum(pdf) + cdf = np.cumsum(pdf) cdf /= cdf[-1] cdf_reflect = scipy.interpolate.interp1d(cdf, x) - #if plot: - # pylab.figure() - # pylab.plot(x, f(x)) - # pylab.scatter(self.x, self.y, c='red') - # - # pylab.figure() - # pylab.plot(x, pdf) - # - # pylab.figure() - # pylab.plot(cdf, x) - return cdf_reflect(alpha) def bayesianUpperLimit2(self, alpha, steps=1.e5, plot=False): @@ -186,85 +156,44 @@ def bayesianUpperLimit2(self, alpha, steps=1.e5, plot=False): f = scipy.interpolate.interp1d(self.x[cut], self.y[cut], kind='cubic') except: f = scipy.interpolate.interp1d(self.x[cut], self.y[cut], kind='linear') - x = numpy.linspace(0., numpy.max(self.x[cut]), steps) - y = numpy.exp(f(x) / 2.) - #forbidden = numpy.nonzero((y / numpy.exp(self.vertex_y / 2.)) < 1.e-10)[0] - forbidden = numpy.nonzero((y / self.vertex_y) < 1.e-10)[0] + x = np.linspace(0., np.max(self.x[cut]), steps) + y = np.exp(f(x) / 2.) + #forbidden = np.nonzero((y / np.exp(self.vertex_y / 2.)) < 1.e-10)[0] + forbidden = np.nonzero((y / self.vertex_y) < 1.e-10)[0] if len(forbidden) > 0: index = forbidden[0] # Numeric stability x = x[0: index] y = y[0: index] - cdf = numpy.cumsum(y) + cdf = np.cumsum(y) cdf /= cdf[-1] cdf_reflect = scipy.interpolate.interp1d(cdf, x) - #if plot: - # pylab.figure() - # pylab.scatter(self.x, self.y) - # - # pylab.figure() - # pylab.plot(x, f(x)) - # - # pylab.figure() - # pylab.plot(x, y) - # - # pylab.figure() - # pylab.plot(cdf, x) - return cdf_reflect(alpha) - """ - if numpy.isnan(result): - import pylab - - for ii in range(0, len(self.x)): - print('%.3f %.3f'%(self.x[ii], self.y[ii])) - - pylab.figure() - pylab.scatter(self.x, self.y) - pylab.figure() - pylab.scatter(cdf, x) - raw_input('WAIT') - - return result - """ def confidenceInterval(self, alpha=0.6827, steps=1.e5, plot=False): """ Compute two-sided confidence interval by taking x-values corresponding to the largest PDF-values first. """ x_dense, y_dense = self.densify() - y_dense -= numpy.max(y_dense) # Numeric stability + y_dense -= np.max(y_dense) # Numeric stability f = scipy.interpolate.interp1d(x_dense, y_dense, kind='linear') - x = numpy.linspace(0., numpy.max(x_dense), steps) + x = np.linspace(0., np.max(x_dense), steps) # ADW: Why does this start at 0, which often outside the input range? # Wouldn't starting at xmin be better: - #x = numpy.linspace(numpy.min(x_dense), numpy.max(x_dense), steps) - pdf = numpy.exp(f(x) / 2.) - cut = (pdf / numpy.max(pdf)) > 1.e-10 + #x = np.linspace(np.min(x_dense), np.max(x_dense), steps) + pdf = np.exp(f(x) / 2.) + cut = (pdf / np.max(pdf)) > 1.e-10 x = x[cut] pdf = pdf[cut] - sorted_pdf_indices = numpy.argsort(pdf)[::-1] # Indices of PDF in descending value - cdf = numpy.cumsum(pdf[sorted_pdf_indices]) + sorted_pdf_indices = np.argsort(pdf)[::-1] # Indices of PDF in descending value + cdf = np.cumsum(pdf[sorted_pdf_indices]) cdf /= cdf[-1] - sorted_pdf_index_max = numpy.argmin((cdf - alpha)**2) + sorted_pdf_index_max = np.argmin((cdf - alpha)**2) x_select = x[sorted_pdf_indices[0: sorted_pdf_index_max]] - - #if plot: - # cdf = numpy.cumsum(pdf) - # cdf /= cdf[-1] - # print( cdf[numpy.max(sorted_pdf_indices[0: sorted_pdf_index_max])] \ - # - cdf[numpy.min(sorted_pdf_indices[0: sorted_pdf_index_max])] ) - # - # pylab.figure() - # pylab.plot(x, f(x)) - # pylab.scatter(self.x, self.y, c='red') - # - # pylab.figure() - # pylab.plot(x, pdf) - return numpy.min(x_select), numpy.max(x_select) + return np.min(x_select), np.max(x_select) ############################################################ @@ -277,8 +206,8 @@ def upperLimitsDeltaTS(confidence_level, one_sided=True, degrees_of_freedom=1): ts_min = 0 # TS = Test Statistic ts_max = 5 ts_steps = 1000 - x = numpy.linspace(ts_min, ts_max, ts_steps) + x = np.linspace(ts_min, ts_max, ts_steps) y = (0.5 * scipy.stats.chi2.sf(x, degrees_of_freedom) - (1. - confidence_level))**2 - return x[numpy.argmin(y)] + return x[np.argmin(y)] ############################################################ diff --git a/ugali/utils/plotting.py b/ugali/utils/plotting.py index 6196b1f..55aa7a6 100644 --- a/ugali/utils/plotting.py +++ b/ugali/utils/plotting.py @@ -10,11 +10,8 @@ except KeyError: matplotlib.use('Agg') import yaml -import numpy import numpy as np -import pylab import pylab as plt -import healpy import healpy as hp import fitsio import scipy.ndimage as nd @@ -66,11 +63,11 @@ def histogram(title, title_x, title_y, """ Plot a basic histogram. """ - pylab.figure() - pylab.hist(x, bins_x) - pylab.xlabel(title_x) - pylab.ylabel(title_y) - pylab.title(title) + plt.figure() + plt.hist(x, bins_x) + plt.xlabel(title_x) + plt.ylabel(title_y) + plt.title(title) ############################################################ @@ -81,27 +78,27 @@ def twoDimensionalHistogram(title, title_x, title_y, """ Create a two-dimension histogram plot or binned map. - If using the outputs of numpy.histogram2d, remember to transpose the histogram. + If using the outputs of np.histogram2d, remember to transpose the histogram. INPUTS """ - pylab.figure() + plt.figure() - mesh_x, mesh_y = numpy.meshgrid(bins_x, bins_y) + mesh_x, mesh_y = np.meshgrid(bins_x, bins_y) if vmin != None and vmin == vmax: - pylab.pcolor(mesh_x, mesh_y, z) + plt.pcolor(mesh_x, mesh_y, z) else: - pylab.pcolor(mesh_x, mesh_y, z, vmin=vmin, vmax=vmax) - pylab.xlabel(title_x) - pylab.ylabel(title_y) - pylab.title(title) - pylab.colorbar() + plt.pcolor(mesh_x, mesh_y, z, vmin=vmin, vmax=vmax) + plt.xlabel(title_x) + plt.ylabel(title_y) + plt.title(title) + plt.colorbar() if lim_x: - pylab.xlim(lim_x[0], lim_x[1]) + plt.xlim(lim_x[0], lim_x[1]) if lim_y: - pylab.ylim(lim_y[0], lim_y[1]) + plt.ylim(lim_y[0], lim_y[1]) ############################################################ @@ -114,20 +111,20 @@ def twoDimensionalScatter(title, title_x, title_y, INPUTS """ - pylab.figure() + plt.figure() - pylab.scatter(x, y, c=color, s=size, alpha=alpha, edgecolors='none') + plt.scatter(x, y, c=color, s=size, alpha=alpha, edgecolors='none') - pylab.xlabel(title_x) - pylab.ylabel(title_y) - pylab.title(title) + plt.xlabel(title_x) + plt.ylabel(title_y) + plt.title(title) if type(color) is not str: - pylab.colorbar() + plt.colorbar() if lim_x: - pylab.xlim(lim_x[0], lim_x[1]) + plt.xlim(lim_x[0], lim_x[1]) if lim_y: - pylab.ylim(lim_y[0], lim_y[1]) + plt.ylim(lim_y[0], lim_y[1]) ############################################################ @@ -137,7 +134,7 @@ def zoomedHealpixMap(title, map, lon, lat, radius, Inputs: lon (deg), lat (deg), radius (deg) """ reso = 60. * 2. * radius / xsize # Deg to arcmin - healpy.gnomview(map=map, rot=[lon, lat, 0], title=title, xsize=xsize, reso=reso, degree=False, **kwargs) + hp.gnomview(map=map, rot=[lon, lat, 0], title=title, xsize=xsize, reso=reso, degree=False, **kwargs) ############################################################ @@ -146,7 +143,7 @@ def projScatter(lon, lat, **kwargs): Create a scatter plot on HEALPix projected axes. Inputs: lon (deg), lat (deg) """ - healpy.projscatter(lon, lat, lonlat=True, **kwargs) + hp.projscatter(lon, lat, lonlat=True, **kwargs) ############################################################ @@ -157,7 +154,7 @@ def sparseHealpixFiles(title, infiles, field='MAGLIM',**kwargs): """ #map = ugali.utils.skymap.readSparseHealpixMaps(infiles,field) map = ugali.utils.skymap.read_partial_map(infiles,field) - ax = healpy.mollview(map=map, title=title, **kwargs) + ax = hp.mollview(map=map, title=title, **kwargs) return ax, map ############################################################ @@ -187,7 +184,7 @@ def drawHealpixMap(hpxmap, lon, lat, size=1.0, xsize=501, coord='GC', **kwargs): pix = ang2pix(get_nside(hpxmap),llon,llat) values = hpxmap[pix].reshape(xx.shape) - zz = np.ma.array(values,mask=(values==healpy.UNSEEN),fill_value=np.nan) + zz = np.ma.array(values,mask=(values==hp.UNSEEN),fill_value=np.nan) return drawProjImage(xx,yy,zz,coord=coord,**kwargs) @@ -236,7 +233,7 @@ def getSDSSImage(ra,dec,radius=1.0,xsize=800,opt='GML',**kwargs): tmp = tempfile.NamedTemporaryFile(suffix='.jpeg') cmd='wget --progress=dot:mega -O %s "%s"'%(tmp.name,url+query) subprocess.call(cmd,shell=True) - im = pylab.imread(tmp.name) + im = plt.imread(tmp.name) tmp.close() return im @@ -294,7 +291,7 @@ def getDSSImage(ra,dec,radius=1.0,xsize=800,**kwargs): tmp = tempfile.NamedTemporaryFile(suffix='.gif') cmd='wget --progress=dot:mega -O %s "%s"'%(tmp.name,url+query) subprocess.call(cmd,shell=True) - im = pylab.imread(tmp.name) + im = plt.imread(tmp.name) tmp.close() if service == 'stsci' and xsize: im = scipy.misc.imresize(im,size=(xsize,xsize)) @@ -369,8 +366,8 @@ def drawSmoothCatalog(self, catalog, label=None, **kwargs): delta_x = self.radius/100. smoothing = 2*delta_x - bins = numpy.arange(-self.radius, self.radius + 1.e-10, delta_x) - h, xbins, ybins = numpy.histogram2d(x, y, bins=[bins, bins]) + bins = np.arange(-self.radius, self.radius + 1.e-10, delta_x) + h, xbins, ybins = np.histogram2d(x, y, bins=[bins, bins]) blur = nd.filters.gaussian_filter(h.T, smoothing / delta_x) defaults = dict(cmap='gray_r',rasterized=True) @@ -381,12 +378,12 @@ def drawSmoothCatalog(self, catalog, label=None, **kwargs): if label: plt.text(0.05, 0.95, label, fontsize=10, ha='left', va='top', - color='k', transform=pylab.gca().transAxes, + color='k', transform=plt.gca().transAxes, bbox=dict(facecolor='white', alpha=1., edgecolor='none')) def drawROI(self, ax=None, value=None, pixel=None): if not ax: ax = plt.gca() - roi_map = np.array(healpy.UNSEEN*np.ones(healpy.nside2npix(self.nside))) + roi_map = np.array(hp.UNSEEN*np.ones(hp.nside2npix(self.nside))) if value is None: roi_map[self.roi.pixels] = 1 @@ -398,7 +395,7 @@ def drawROI(self, ax=None, value=None, pixel=None): roi_map[pixel] = value else: logger.warning('Unable to parse input') - #im = healpy.gnomview(roi_map,**self.gnom_kwargs) + #im = hp.gnomview(roi_map,**self.gnom_kwargs) im = drawHealpixMap(roi_map,self.lon,self.lat,self.radius,coord=self.coord) return im @@ -442,19 +439,19 @@ def drawStellarDensity(self,ax=None,nside=None): #catalog=ugali.observation.catalog.Catalog(self.config,roi=self.roi) pix = ang2pix(nside, catalog.lon, catalog.lat) counts = collections.Counter(pix) - pixels, number = numpy.array(sorted(counts.items())).T - star_map = hp.UNSEEN * numpy.ones(hp.nside2npix(nside)) + pixels, number = np.array(sorted(counts.items())).T + star_map = hp.UNSEEN * np.ones(hp.nside2npix(nside)) star_map[pixels] = number - star_map[star_map == 0] = healpy.UNSEEN + star_map[star_map == 0] = hp.UNSEEN - #im = healpy.gnomview(star_map,**self.gnom_kwargs) - #healpy.graticule(dpar=1,dmer=1,color='0.5',verbose=False) - #pylab.close() + #im = hp.gnomview(star_map,**self.gnom_kwargs) + #hp.graticule(dpar=1,dmer=1,color='0.5',verbose=False) + #plt.close() im = drawHealpixMap(star_map,self.lon,self.lat,self.radius,coord=self.coord) #im = ax.imshow(im,origin='bottom') try: ax.cax.colorbar(im) - except: pylab.colorbar(im,ax=ax) + except: plt.colorbar(im,ax=ax) ax.annotate("Stars",**self.label_kwargs) return im @@ -477,7 +474,7 @@ def drawMask(self,ax=None, mask=None, mtype='maglim'): im = drawHealpixMap(mask_map,self.lon,self.lat,self.radius,coord=self.coord) try: cbar = ax.cax.colorbar(im) - except: cbar = pylab.colorbar(im) + except: cbar = plt.colorbar(im) cbar.ax.set_xticklabels(cbar.ax.get_xticklabels(),rotation=90) ax.annotate(mtype,**self.label_kwargs) return im @@ -503,22 +500,22 @@ def drawTS(self,ax=None, filename=None, zidx=0): values = 2*data['LOG_LIKELIHOOD'] if values.ndim == 1: values = values.reshape(-1,1) - ts_map = healpy.UNSEEN * numpy.ones(healpy.nside2npix(self.nside)) + ts_map = hp.UNSEEN * np.ones(hp.nside2npix(self.nside)) # Sum through all distance_moduli #ts_map[pixels] = values.sum(axis=1) # Just at maximum slice from object ts_map[pixels] = values[:,zidx] - #im = healpy.gnomview(ts_map,**self.gnom_kwargs) - #healpy.graticule(dpar=1,dmer=1,color='0.5',verbose=False) - #pylab.close() + #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: pylab.colorbar(im) + except: plt.colorbar(im) ax.annotate("TS",**self.label_kwargs) return im @@ -526,7 +523,7 @@ def drawCatalog(self, ax=None): if not ax: ax = plt.gca() # Stellar Catalog self._create_catalog() - healpy.projscatter(self.catalog.lon,self.catalog.lat,c='k',marker='.',lonlat=True,coord=self.gnom_kwargs['coord']) + hp.projscatter(self.catalog.lon,self.catalog.lat,c='k',marker='.',lonlat=True,coord=self.gnom_kwargs['coord']) ax.annotate("Stars",**self.label_kwargs) def drawSpatial(self, ax=None): @@ -621,7 +618,7 @@ def drawMembership(self, ax=None, radius=None, zidx=0, mc_source_id=1): ax.set_xlabel('Color (mag)') ax.set_ylabel('Magnitude (mag)') try: ax.cax.colorbar(sc) - except: pylab.colorbar(sc) + except: plt.colorbar(sc) def plotDistance(self): filename = self.config.mergefile @@ -632,27 +629,27 @@ def plotDistance(self): 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 = healpy.UNSEEN * numpy.ones(healpy.nside2npix(self.nside)) + ts_map = hp.UNSEEN * np.ones(hp.nside2npix(self.nside)) ndim = len(distances) - nrows = int(numpy.sqrt(ndim)) + nrows = int(np.sqrt(ndim)) ncols = ndim // nrows + (ndim%nrows > 0) # Create the healpy images, but close the figures images = [] for i,val in enumerate(values.T): ts_map[pixels] = val - im = healpy.gnomview(ts_map,**self.gnom_kwargs) - pylab.close() + im = hp.gnomview(ts_map,**self.gnom_kwargs) + plt.close() images.append(im) - data = numpy.array(images); mask = (data == healpy.UNSEEN) - images = numpy.ma.array(data=data,mask=mask) - vmin = numpy.ma.min(images) - vmax = numpy.ma.max(images) + data = np.array(images); mask = (data == hp.UNSEEN) + images = np.ma.array(data=data,mask=mask) + vmin = np.ma.min(images) + vmax = np.ma.max(images) # Create the image grid - fig = pylab.figure() + fig = plt.figure() axes = AxesGrid(fig, 111, nrows_ncols = (nrows, ncols), axes_pad=0,label_mode='1', cbar_mode='single',cbar_pad=0,cbar_size='5%', @@ -675,7 +672,7 @@ def plotDistance(self): def plot3(self): - fig = pylab.figure(figsize=(8,4)) + fig = plt.figure(figsize=(8,4)) axes = AxesGrid(fig, 111,nrows_ncols = (1, 3),axes_pad=0.1, cbar_mode='each',cbar_pad=0,cbar_size='5%', cbar_location='top',share_all=True) @@ -691,7 +688,7 @@ def plot3(self): def plot4(self): - fig = pylab.figure(figsize=(10,8)) + fig = plt.figure(figsize=(10,8)) axes = AxesGrid(fig, 111,nrows_ncols=(2, 2), axes_pad=0.35, cbar_mode='each',cbar_pad=0,cbar_size='5%', share_all=True,aspect=True, @@ -825,14 +822,14 @@ def drawHessDiagram(self,catalog=None): color = catalog.color[cut_annulus] mag = catalog.mag[cut_annulus] - h, xbins, ybins = numpy.histogram2d(color, mag, bins=[cbins,mbins]) + h, xbins, ybins = np.histogram2d(color, mag, bins=[cbins,mbins]) blur = nd.filters.gaussian_filter(h.T, 2) kwargs = dict(extent=[xbins.min(),xbins.max(),ybins.min(),ybins.max()], cmap='gray_r', aspect='auto', origin='lower', rasterized=True, interpolation='none') ax.imshow(blur, **kwargs) - pylab.scatter(catalog.color[cut_inner], catalog.mag[cut_inner], + plt.scatter(catalog.color[cut_inner], catalog.mag[cut_inner], c='red', s=7, edgecolor='none')# label=r'$r < %.2f$ deg'%(r_peak)) ugali.utils.plotting.drawIsochrone(self.isochrone, c='b', zorder=10) ax.set_xlim(-0.5, 1.) @@ -840,12 +837,12 @@ def drawHessDiagram(self,catalog=None): plt.xlabel(r'$g - r$') plt.ylabel(r'$g$') plt.xticks([-0.5, 0., 0.5, 1.]) - plt.yticks(numpy.arange(mmax - 1., mmin - 1., -1.)) + plt.yticks(np.arange(mmax - 1., mmin - 1., -1.)) radius_string = (r'${\rm r}<%.1f$ arcmin'%( 60 * r_peak)) - pylab.text(0.05, 0.95, radius_string, + plt.text(0.05, 0.95, radius_string, fontsize=10, ha='left', va='top', color='red', - transform=pylab.gca().transAxes, + transform=plt.gca().transAxes, bbox=dict(facecolor='white', alpha=1., edgecolor='none')) @@ -863,7 +860,7 @@ def drawMembersSpatial(self,data): sel = (x_prob > xmin)&(x_prob < xmax) & (y_prob > ymin)&(y_prob < ymax) sel_prob = data['PROB'][sel] > 5.e-2 - index_sort = numpy.argsort(data['PROB'][sel][sel_prob]) + index_sort = np.argsort(data['PROB'][sel][sel_prob]) plt.scatter(x_prob[sel][~sel_prob], y_prob[sel][~sel_prob], marker='o', s=2, c='0.75', edgecolor='none') @@ -884,7 +881,7 @@ def drawMembersSpatial(self,data): divider = make_axes_locatable(ax) ax_cb = divider.new_horizontal(size="7%", pad=0.1) plt.gcf().add_axes(ax_cb) - pylab.colorbar(sc, cax=ax_cb, orientation='vertical', ticks=[0, 0.2, 0.4, 0.6, 0.8, 1.0], label='Membership Probability') + plt.colorbar(sc, cax=ax_cb, orientation='vertical', ticks=[0, 0.2, 0.4, 0.6, 0.8, 1.0], label='Membership Probability') ax_cb.yaxis.tick_right() def drawMembersCMD(self,data): @@ -907,39 +904,39 @@ def drawMembersCMD(self,data): sel = (x_prob > xmin)&(x_prob < xmax) & (y_prob > ymin)&(y_prob < ymax) sel_prob = data['PROB'][sel] > 5.e-2 - index_sort = numpy.argsort(data['PROB'][sel][sel_prob]) + index_sort = np.argsort(data['PROB'][sel][sel_prob]) plt.scatter(data['COLOR'][sel][~sel_prob], mag_1[sel][~sel_prob], marker='o',s=2,c='0.75',edgecolor='none') - sc = pylab.scatter(data['COLOR'][sel][sel_prob][index_sort], mag_1[sel][sel_prob][index_sort], + sc = plt.scatter(data['COLOR'][sel][sel_prob][index_sort], mag_1[sel][sel_prob][index_sort], c=data['PROB'][sel][sel_prob][index_sort], marker='o', s=10, edgecolor='none', cmap='jet', vmin=0., vmax=1) - pylab.xlim(cmin, cmax) - pylab.ylim(mmax, mmin) - pylab.xlabel(r'$g - r$') - pylab.ylabel(r'$g$') + plt.xlim(cmin, cmax) + plt.ylim(mmax, mmin) + plt.xlabel(r'$g - r$') + plt.ylabel(r'$g$') #axes[1].yaxis.set_major_locator(MaxNLocator(prune='lower')) - pylab.xticks([-0.5, 0., 0.5, 1.]) - pylab.yticks(numpy.arange(mmax - 1., mmin - 1., -1.)) + plt.xticks([-0.5, 0., 0.5, 1.]) + plt.yticks(np.arange(mmax - 1., mmin - 1., -1.)) ugali.utils.plotting.drawIsochrone(self.isochrone, c='k', zorder=10) - pylab.text(0.05, 0.95, r'$\Sigma p_{i} = %i$'%(data['PROB'].sum()), - fontsize=10, horizontalalignment='left', verticalalignment='top', color='k', transform=pylab.gca().transAxes, + plt.text(0.05, 0.95, r'$\Sigma p_{i} = %i$'%(data['PROB'].sum()), + fontsize=10, horizontalalignment='left', verticalalignment='top', color='k', transform=plt.gca().transAxes, bbox=dict(facecolor='white', alpha=1., edgecolor='none')) - divider = make_axes_locatable(pylab.gca()) + divider = make_axes_locatable(plt.gca()) ax_cb = divider.new_horizontal(size="7%", pad=0.1) plt.gcf().add_axes(ax_cb) - pylab.colorbar(sc, cax=ax_cb, orientation='vertical', ticks=[0, 0.2, 0.4, 0.6, 0.8, 1.0], label='Membership Probability') + plt.colorbar(sc, cax=ax_cb, orientation='vertical', ticks=[0, 0.2, 0.4, 0.6, 0.8, 1.0], label='Membership Probability') ax_cb.yaxis.tick_right() def drawDensityProfile(self, catalog=None): rmax = 24. # arcmin - bins = numpy.arange(0, rmax + 1.e-10, 2.) + bins = np.arange(0, rmax + 1.e-10, 2.) centers = 0.5 * (bins[1:] + bins[0:-1]) - area = numpy.pi * (bins[1:]**2 - bins[0:-1]**2) + area = np.pi * (bins[1:]**2 - bins[0:-1]**2) r_peak = self.kernel.extension @@ -949,8 +946,8 @@ def drawDensityProfile(self, catalog=None): angsep_arcmin = angsep * 60 # arcmin cut_iso = self.isochrone_selection(stars) - h = numpy.histogram(angsep_arcmin[(angsep_arcmin < rmax) & cut_iso], bins=bins)[0] - h_out = numpy.histogram(angsep_arcmin[(angsep_arcmin < rmax) & (~cut_iso)], bins=bins)[0] + h = np.histogram(angsep_arcmin[(angsep_arcmin < rmax) & cut_iso], bins=bins)[0] + h_out = np.histogram(angsep_arcmin[(angsep_arcmin < rmax) & (~cut_iso)], bins=bins)[0] gals = self.get_galaxies() if len(gals): @@ -963,25 +960,25 @@ def drawDensityProfile(self, catalog=None): h_gal_out = np.histogram(angsep_gal_arcmin[(angsep_gal_arcmin < rmax) & (~cut_iso_gal)], bins=bins)[0] plt.plot(centers, h/area, c='red', label='Filtered Stars') - plt.errorbar(centers, h/area, yerr=(numpy.sqrt(h) / area), ecolor='red', c='red') + plt.errorbar(centers, h/area, yerr=(np.sqrt(h) / area), ecolor='red', c='red') plt.scatter(centers, h/area, edgecolor='none', c='red', zorder=22) plt.plot(centers, h_out/area, c='gray', label='Unfiltered Stars') - plt.errorbar(centers, h_out/area, yerr=(numpy.sqrt(h_out) / area), ecolor='gray', c='gray') + plt.errorbar(centers, h_out/area, yerr=(np.sqrt(h_out) / area), ecolor='gray', c='gray') plt.scatter(centers, h_out/area, edgecolor='none', c='gray', zorder=21) if len(gals): plt.plot(centers, h_gal/area, c='black', label='Filtered Galaxies') - plt.errorbar(centers, h_gal/area, yerr=(numpy.sqrt(h_gal) / area), ecolor='black', c='black') + plt.errorbar(centers, h_gal/area, yerr=(np.sqrt(h_gal) / area), ecolor='black', c='black') plt.scatter(centers, h_gal/area, edgecolor='none', c='black', zorder=20) plt.xlabel('Angular Separation (arcmin)') plt.ylabel(r'Density (arcmin$^{-2}$)') plt.xlim(0., rmax) - ymax = pylab.ylim()[1] - #pylab.ylim(0, ymax) - pylab.ylim(0, 4) - pylab.legend(loc='upper right', frameon=False, fontsize=10) + ymax = plt.ylim()[1] + #plt.ylim(0, ymax) + plt.ylim(0, 4) + plt.legend(loc='upper right', frameon=False, fontsize=10) def plot6(self, filename, title=None): @@ -1001,11 +998,11 @@ def plot6(self, filename, title=None): self.drawSmoothStars() logger.debug("Drawing density profile...") - pylab.subplot(2, 3, 2) + plt.subplot(2, 3, 2) self.drawDensityProfile() logger.debug("Drawing spatial distribution of members...") - pylab.subplot(2, 3, 3) + plt.subplot(2, 3, 3) self.drawMembersSpatial(filename) logger.debug("Drawing smooth galaxies...") @@ -1017,7 +1014,7 @@ def plot6(self, filename, title=None): self.drawHessDiagram() logger.debug("Drawing CMD of members...") - pylab.subplot(2, 3, 6) + plt.subplot(2, 3, 6) self.drawMembersCMD(filename) def plot_candidates(candidates, config, ts_min=50, outdir='./'): @@ -1127,7 +1124,7 @@ def drawKernelHist(ax, kernel): ax.set_xlim(0,len(x)) #try: ax.cax.colorbar(im) - #except: pylab.colorbar(im) + #except: plt.colorbar(im) #a0 = np.array([0.,0.]) #a1 =kernel.a*np.array([np.sin(np.deg2rad(theta)),-np.cos(np.deg2rad(theta))]) @@ -1441,9 +1438,9 @@ def drawSkymapCatalog(ax,lon,lat,**kwargs): def plotSkymap(skymap, proj='mol', **kwargs): kwargs.setdefault('xsize',1000) if proj.upper() == 'MOL': - im = healpy.mollview(skymap,**kwargs) + im = hp.mollview(skymap,**kwargs) elif proj.upper() == 'CAR': - im = healpy.cartview(skymap,**kwargs) + im = hp.cartview(skymap,**kwargs) return im def plotTriangle(srcfile,samples,burn=0,**kwargs): @@ -1503,8 +1500,8 @@ def makePath(x_path, y_path, epsilon=1.e-10): """ Create closed path. """ - x_path_closed = numpy.concatenate([x_path, x_path[::-1]]) - y_path_closed = numpy.concatenate([y_path, epsilon + y_path[::-1]]) + x_path_closed = np.concatenate([x_path, x_path[::-1]]) + y_path_closed = np.concatenate([y_path, epsilon + y_path[::-1]]) path = matplotlib.path.Path(list(zip(x_path_closed, y_path_closed))) return path @@ -1526,12 +1523,12 @@ def cutIsochronePath(g, r, g_err, r_err, isochrone, radius=0.1, return_all=False return np.array([],dtype=bool) try: - if numpy.all(isochrone.stage == 'Main'): + if np.all(isochrone.stage == 'Main'): # Dotter case index_transition = len(isochrone.stage) else: # Other cases - index_transition = numpy.nonzero(isochrone.stage > 3)[0][0] + 1 + index_transition = np.nonzero(isochrone.stage > 3)[0][0] + 1 except AttributeError: index_transition = 1 @@ -1542,34 +1539,34 @@ def cutIsochronePath(g, r, g_err, r_err, isochrone, radius=0.1, return_all=False # Cut one way... f_isochrone = scipy.interpolate.interp1d(mag_2_rgb, mag_1_rgb - mag_2_rgb, bounds_error=False, fill_value = 999.) - color_diff = numpy.fabs((g - r) - f_isochrone(r)) - cut_2 = (color_diff < numpy.sqrt(0.1**2 + r_err**2 + g_err**2)) + color_diff = np.fabs((g - r) - f_isochrone(r)) + cut_2 = (color_diff < np.sqrt(0.1**2 + r_err**2 + g_err**2)) # ...and now the other f_isochrone = scipy.interpolate.interp1d(mag_1_rgb, mag_1_rgb - mag_2_rgb, bounds_error=False, fill_value = 999.) - color_diff = numpy.fabs((g - r) - f_isochrone(g)) - cut_1 = (color_diff < numpy.sqrt(0.1**2 + r_err**2 + g_err**2)) + color_diff = np.fabs((g - r) - f_isochrone(g)) + cut_1 = (color_diff < np.sqrt(0.1**2 + r_err**2 + g_err**2)) - cut = numpy.logical_or(cut_1, cut_2) + cut = np.logical_or(cut_1, cut_2) # If using Padova isochrone, also include horizontal branch - if not numpy.all(isochrone.stage == 'Main'): - index_transition = numpy.nonzero(isochrone.stage > 3)[0][0] + 1 + if not np.all(isochrone.stage == 'Main'): + index_transition = np.nonzero(isochrone.stage > 3)[0][0] + 1 mag_1_hb = isochrone.mag_1[index_transition:] + isochrone.distance_modulus mag_2_hb = isochrone.mag_2[index_transition:] + isochrone.distance_modulus path_hb = makePath(mag_1_hb, mag_2_hb) cut_hb = path_hb.contains_points(list(zip(g, r)), radius=0.1) logger.debug('Applying HB selection') - logger.debug(numpy.sum(cut)) - cut = numpy.logical_or(cut, cut_hb) - logger.debug(numpy.sum(cut)) + logger.debug(np.sum(cut)) + cut = np.logical_or(cut, cut_hb) + logger.debug(np.sum(cut)) - mag_bins = numpy.arange(16., 24.1, 0.1) + mag_bins = np.arange(16., 24.1, 0.1) mag_centers = 0.5 * (mag_bins[1:] + mag_bins[0:-1]) - magerr = numpy.tile(0., len(mag_centers)) + magerr = np.tile(0., len(mag_centers)) for ii in range(0, len(mag_bins) - 1): cut_mag_bin = (g > mag_bins[ii]) & (g < mag_bins[ii + 1]) - magerr[ii] = numpy.median(numpy.sqrt(0.1**2 + r_err[cut_mag_bin]**2 + g_err[cut_mag_bin]**2)) + magerr[ii] = np.median(np.sqrt(0.1**2 + r_err[cut_mag_bin]**2 + g_err[cut_mag_bin]**2)) if return_all: return cut, mag_centers[f_isochrone(mag_centers) < 100], (f_isochrone(mag_centers) + magerr)[f_isochrone(mag_centers) < 100], (f_isochrone(mag_centers) - magerr)[f_isochrone(mag_centers) < 100] diff --git a/ugali/utils/projector.py b/ugali/utils/projector.py index f85b233..e5099a1 100644 --- a/ugali/utils/projector.py +++ b/ugali/utils/projector.py @@ -5,7 +5,6 @@ http://adsabs.harvard.edu/abs/2002A%26A...395.1077C """ -import numpy import numpy as np from ugali.utils.logger import logger @@ -25,20 +24,20 @@ def __init__(self, lon_ref, lat_ref, zenithal=False): def setReference(self, lon_ref, lat_ref, zenithal=False): if zenithal: - phi = (numpy.pi / 2.) + numpy.radians(lon_ref) - theta = (numpy.pi / 2.) - numpy.radians(lat_ref) + phi = (np.pi / 2.) + np.radians(lon_ref) + theta = (np.pi / 2.) - np.radians(lat_ref) psi = 0. if not zenithal: - phi = (-numpy.pi / 2.) + numpy.radians(lon_ref) - theta = numpy.radians(lat_ref) - psi = numpy.radians(90.) # psi = 90 corresponds to (0, 0) psi = -90 corresponds to (180, 0) + phi = (-np.pi / 2.) + np.radians(lon_ref) + theta = np.radians(lat_ref) + psi = np.radians(90.) # psi = 90 corresponds to (0, 0) psi = -90 corresponds to (180, 0) - cos_psi,sin_psi = numpy.cos(psi),numpy.sin(psi) - cos_phi,sin_phi = numpy.cos(phi),numpy.sin(phi) - cos_theta,sin_theta = numpy.cos(theta),numpy.sin(theta) + cos_psi,sin_psi = np.cos(psi),np.sin(psi) + cos_phi,sin_phi = np.cos(phi),np.sin(phi) + cos_theta,sin_theta = np.cos(theta),np.sin(theta) - self.rotation_matrix = numpy.matrix([ + self.rotation_matrix = np.matrix([ [cos_psi * cos_phi - cos_theta * sin_phi * sin_psi, cos_psi * sin_phi + cos_theta * cos_phi * sin_psi, sin_psi * sin_theta], @@ -50,30 +49,30 @@ def setReference(self, lon_ref, lat_ref, zenithal=False): cos_theta] ]) - self.inverted_rotation_matrix = numpy.linalg.inv(self.rotation_matrix) + self.inverted_rotation_matrix = np.linalg.inv(self.rotation_matrix) def cartesian(self,lon,lat): - lon = numpy.radians(lon) - lat = numpy.radians(lat) + lon = np.radians(lon) + lat = np.radians(lat) - x = numpy.cos(lat) * numpy.cos(lon) - y = numpy.cos(lat) * numpy.sin(lon) - z = numpy.sin(lat) - return numpy.array([x,y,z]) + x = np.cos(lat) * np.cos(lon) + y = np.cos(lat) * np.sin(lon) + z = np.sin(lat) + return np.array([x,y,z]) def rotate(self, lon, lat, invert=False): vec = self.cartesian(lon,lat) if invert: - vec_prime = numpy.dot(numpy.array(self.inverted_rotation_matrix), vec) + vec_prime = np.dot(np.array(self.inverted_rotation_matrix), vec) else: - vec_prime = numpy.dot(numpy.array(self.rotation_matrix), vec) + vec_prime = np.dot(np.array(self.rotation_matrix), vec) - lon_prime = numpy.arctan2(vec_prime[1], vec_prime[0]) - lat_prime = numpy.arcsin(vec_prime[2]) + lon_prime = np.arctan2(vec_prime[1], vec_prime[0]) + lat_prime = np.arcsin(vec_prime[2]) - return (numpy.degrees(lon_prime) % 360.), numpy.degrees(lat_prime) + return (np.degrees(lon_prime) % 360.), np.degrees(lat_prime) ############################################################ @@ -149,15 +148,15 @@ def aitoffSphereToImage(lon, lat): Hammer-Aitoff projection (deg). """ lon = lon - 360.*(lon>180) - lon = numpy.radians(lon) - lat = numpy.radians(lat) + lon = np.radians(lon) + lat = np.radians(lat) half_lon = lon/2. - cos_lat = numpy.cos(lat) + cos_lat = np.cos(lat) - gamma = (180. / numpy.pi) * numpy.sqrt(2. / (1. + (cos_lat * numpy.cos(half_lon)))) - x = 2. * gamma * cos_lat * numpy.sin(half_lon) - y = gamma * numpy.sin(lat) + gamma = (180. / np.pi) * np.sqrt(2. / (1. + (cos_lat * np.cos(half_lon)))) + x = 2. * gamma * cos_lat * np.sin(half_lon) + y = gamma * np.sin(lat) return x, y def aitoffImageToSphere(x, y): @@ -165,12 +164,12 @@ def aitoffImageToSphere(x, y): Inverse Hammer-Aitoff projection (deg). """ x = x - 360.*(x>180) - x = numpy.asarray(numpy.radians(x)) - y = numpy.asarray(numpy.radians(y)) - z = numpy.sqrt(1. - (x / 4.)**2 - (y / 2.)**2) # rad - lon = 2. * numpy.arctan2((2. * z**2) - 1, (z / 2.) * x) - lat = numpy.arcsin( y * z) - return ((180. - numpy.degrees(lon)) % 360.), numpy.degrees(lat) + x = np.asarray(np.radians(x)) + y = np.asarray(np.radians(y)) + z = np.sqrt(1. - (x / 4.)**2 - (y / 2.)**2) # rad + lon = 2. * np.arctan2((2. * z**2) - 1, (z / 2.) * x) + lat = np.arcsin( y * z) + return ((180. - np.degrees(lon)) % 360.), np.degrees(lat) ############################################################ @@ -180,11 +179,11 @@ def gnomonicSphereToImage(lon, lat): """ # Convert angle to [-180, 180] interval lon = lon - 360.*(lon>180) - lon = numpy.radians(lon) - lat = numpy.radians(lat) - r_theta = (180. / numpy.pi) / numpy.tan(lat) - x = r_theta * numpy.cos(lon) - y = r_theta * numpy.sin(lon) + lon = np.radians(lon) + lat = np.radians(lat) + r_theta = (180. / np.pi) / np.tan(lat) + x = r_theta * np.cos(lon) + y = r_theta * np.sin(lon) return x, y def gnomonicImageToSphere(x, y): @@ -193,11 +192,11 @@ def gnomonicImageToSphere(x, y): """ # Convert angle to [-180, 180] interval x = x - 360.*(x>180) - x = numpy.asarray(x) - y = numpy.asarray(y) - lon = numpy.degrees(numpy.arctan2(y, x)) - r_theta = numpy.sqrt(x**2 + y**2) - lat = numpy.degrees(numpy.arctan(180. / (numpy.pi * r_theta))) + x = np.asarray(x) + y = np.asarray(y) + lon = np.degrees(np.arctan2(y, x)) + r_theta = np.sqrt(x**2 + y**2) + lat = np.degrees(np.arctan(180. / (np.pi * r_theta))) return lon, lat ############################################################ @@ -208,11 +207,11 @@ def angsep2(lon_1, lat_1, lon_2, lat_2): """ import healpy - v10, v11, v12 = healpy.ang2vec(numpy.radians(90. - lat_1), numpy.radians(lon_1)).transpose() - v20, v21, v22 = healpy.ang2vec(numpy.radians(90. - lat_2), numpy.radians(lon_2)).transpose() + v10, v11, v12 = healpy.ang2vec(np.radians(90. - lat_1), np.radians(lon_1)).transpose() + v20, v21, v22 = healpy.ang2vec(np.radians(90. - lat_2), np.radians(lon_2)).transpose() val = (v10 * v20) + (v11 * v21) + (v12 * v22) - val = numpy.clip(val, -1., 1.) - return numpy.degrees(numpy.arccos(val)) + val = np.clip(val, -1., 1.) + return np.degrees(np.arccos(val)) def angsep(lon1,lat1,lon2,lat2): """ @@ -228,21 +227,21 @@ def angsep(lon1,lat1,lon2,lat2): [1] http://en.wikipedia.org/wiki/Great-circle_distance """ - lon1,lat1 = numpy.radians([lon1,lat1]) - lon2,lat2 = numpy.radians([lon2,lat2]) + lon1,lat1 = np.radians([lon1,lat1]) + lon2,lat2 = np.radians([lon2,lat2]) - sdlon = numpy.sin(lon2 - lon1) - cdlon = numpy.cos(lon2 - lon1) - slat1 = numpy.sin(lat1) - slat2 = numpy.sin(lat2) - clat1 = numpy.cos(lat1) - clat2 = numpy.cos(lat2) + sdlon = np.sin(lon2 - lon1) + cdlon = np.cos(lon2 - lon1) + slat1 = np.sin(lat1) + slat2 = np.sin(lat2) + clat1 = np.cos(lat1) + clat2 = np.cos(lat2) num1 = clat2 * sdlon num2 = clat1 * slat2 - slat1 * clat2 * cdlon denominator = slat1 * slat2 + clat1 * clat2 * cdlon - return numpy.degrees(numpy.arctan2(numpy.hypot(num1,num2), denominator)) + return np.degrees(np.arctan2(np.hypot(num1,num2), denominator)) ############################################################ @@ -251,26 +250,26 @@ def galToCel(ll, bb): """ Converts Galactic (deg) to Celestial J2000 (deg) coordinates """ - bb = numpy.radians(bb) - sin_bb = numpy.sin(bb) - cos_bb = numpy.cos(bb) + bb = np.radians(bb) + sin_bb = np.sin(bb) + cos_bb = np.cos(bb) - ll = numpy.radians(ll) - ra_gp = numpy.radians(192.85948) - de_gp = numpy.radians(27.12825) - lcp = numpy.radians(122.932) + ll = np.radians(ll) + ra_gp = np.radians(192.85948) + de_gp = np.radians(27.12825) + lcp = np.radians(122.932) - sin_lcp_ll = numpy.sin(lcp - ll) - cos_lcp_ll = numpy.cos(lcp - ll) + sin_lcp_ll = np.sin(lcp - ll) + cos_lcp_ll = np.cos(lcp - ll) - sin_d = (numpy.sin(de_gp) * sin_bb) \ - + (numpy.cos(de_gp) * cos_bb * cos_lcp_ll) - ramragp = numpy.arctan2(cos_bb * sin_lcp_ll, - (numpy.cos(de_gp) * sin_bb) \ - - (numpy.sin(de_gp) * cos_bb * cos_lcp_ll)) - dec = numpy.arcsin(sin_d) - ra = (ramragp + ra_gp + (2. * numpy.pi)) % (2. * numpy.pi) - return numpy.degrees(ra), numpy.degrees(dec) + sin_d = (np.sin(de_gp) * sin_bb) \ + + (np.cos(de_gp) * cos_bb * cos_lcp_ll) + ramragp = np.arctan2(cos_bb * sin_lcp_ll, + (np.cos(de_gp) * sin_bb) \ + - (np.sin(de_gp) * cos_bb * cos_lcp_ll)) + dec = np.arcsin(sin_d) + ra = (ramragp + ra_gp + (2. * np.pi)) % (2. * np.pi) + return np.degrees(ra), np.degrees(dec) gal2cel = galToCel @@ -278,26 +277,26 @@ def celToGal(ra, dec): """ Converts Celestial J2000 (deg) to Calactic (deg) coordinates """ - dec = numpy.radians(dec) - sin_dec = numpy.sin(dec) - cos_dec = numpy.cos(dec) + dec = np.radians(dec) + sin_dec = np.sin(dec) + cos_dec = np.cos(dec) - ra = numpy.radians(ra) - ra_gp = numpy.radians(192.85948) - de_gp = numpy.radians(27.12825) + ra = np.radians(ra) + ra_gp = np.radians(192.85948) + de_gp = np.radians(27.12825) - sin_ra_gp = numpy.sin(ra - ra_gp) - cos_ra_gp = numpy.cos(ra - ra_gp) + sin_ra_gp = np.sin(ra - ra_gp) + cos_ra_gp = np.cos(ra - ra_gp) - lcp = numpy.radians(122.932) - sin_b = (numpy.sin(de_gp) * sin_dec) \ - + (numpy.cos(de_gp) * cos_dec * cos_ra_gp) - lcpml = numpy.arctan2(cos_dec * sin_ra_gp, - (numpy.cos(de_gp) * sin_dec) \ - - (numpy.sin(de_gp) * cos_dec * cos_ra_gp)) - bb = numpy.arcsin(sin_b) - ll = (lcp - lcpml + (2. * numpy.pi)) % (2. * numpy.pi) - return numpy.degrees(ll), numpy.degrees(bb) + lcp = np.radians(122.932) + sin_b = (np.sin(de_gp) * sin_dec) \ + + (np.cos(de_gp) * cos_dec * cos_ra_gp) + lcpml = np.arctan2(cos_dec * sin_ra_gp, + (np.cos(de_gp) * sin_dec) \ + - (np.sin(de_gp) * cos_dec * cos_ra_gp)) + bb = np.arcsin(sin_b) + ll = (lcp - lcpml + (2. * np.pi)) % (2. * np.pi) + return np.degrees(ll), np.degrees(bb) cel2gal = celToGal @@ -398,7 +397,7 @@ def dec2dms(dec): if isinstance(dec,str): dec = float(dec) - sign = numpy.copysign(1.0,dec) + sign = np.copysign(1.0,dec) fdeg = np.abs(dec) deg = int(fdeg) @@ -424,7 +423,7 @@ def hms2dec(hms): SECOND = 3600. if isinstance(hms,str): - hour,minute,second = numpy.array(re.split('[hms]',hms))[:3].astype(float) + hour,minute,second = np.array(re.split('[hms]',hms))[:3].astype(float) else: hour,minute,second = hms.T @@ -446,12 +445,12 @@ def dms2dec(dms): # http://docs.scipy.org/doc/numpy-1.7.0/reference/c-api.coremath.html#NPY_NZERO if isinstance(dms,str): - degree,minute,second = numpy.array(re.split('[dms]',hms))[:3].astype(float) + degree,minute,second = np.array(re.split('[dms]',hms))[:3].astype(float) else: degree,minute,second = dms.T - sign = numpy.copysign(1.0,degree) - decimal = numpy.abs(degree) + minute * 1./MINUTE + second * 1./SECOND + sign = np.copysign(1.0,degree) + decimal = np.abs(degree) + minute * 1./MINUTE + second * 1./SECOND decimal *= sign return decimal @@ -467,7 +466,7 @@ def distanceToDistanceModulus(distance): """ Return distance modulus for a given distance (kpc). """ - return 5. * (numpy.log10(numpy.array(distance) * 1.e3) - 1.) + return 5. * (np.log10(np.array(distance) * 1.e3) - 1.) dist2mod = distanceToDistanceModulus @@ -475,7 +474,7 @@ def distanceModulusToDistance(distance_modulus): """ Return distance (kpc) for a given distance modulus. """ - return 10**((0.2 * numpy.array(distance_modulus)) - 2.) + return 10**((0.2 * np.array(distance_modulus)) - 2.) mod2dist = distanceModulusToDistance @@ -566,10 +565,10 @@ def match(lon1, lat1, lon2, lat2, tol=None, nnearest=1): """ from scipy.spatial import cKDTree - lon1 = numpy.asarray(lon1) - lat1 = numpy.asarray(lat1) - lon2 = numpy.asarray(lon2) - lat2 = numpy.asarray(lat2) + lon1 = np.asarray(lon1) + lat1 = np.asarray(lat1) + lon2 = np.asarray(lon2) + lat2 = np.asarray(lat2) if lon1.shape != lat1.shape: raise ValueError('lon1 and lat1 do not match!') @@ -579,15 +578,15 @@ def match(lon1, lat1, lon2, lat2, tol=None, nnearest=1): rotator = SphericalRotator(0,0) - # This is equivalent, but faster than just doing numpy.array([x1, y1, z1]).T + # This is equivalent, but faster than just doing np.array([x1, y1, z1]).T x1, y1, z1 = rotator.cartesian(lon1.ravel(),lat1.ravel()) - coords1 = numpy.empty((x1.size, 3)) + coords1 = np.empty((x1.size, 3)) coords1[:, 0] = x1 coords1[:, 1] = y1 coords1[:, 2] = z1 x2, y2, z2 = rotator.cartesian(lon2.ravel(),lat2.ravel()) - coords2 = numpy.empty((x2.size, 3)) + coords2 = np.empty((x2.size, 3)) coords2[:, 0] = x2 coords2[:, 1] = y2 coords2[:, 2] = z2 @@ -602,7 +601,7 @@ def match(lon1, lat1, lon2, lat2, tol=None, nnearest=1): ds = angsep(lon1, lat1, lon2[idxs2], lat2[idxs2]) - idxs1 = numpy.arange(lon1.size) + idxs1 = np.arange(lon1.size) if tol is not None: msk = ds < tol diff --git a/ugali/utils/skymap.py b/ugali/utils/skymap.py index 7868e13..9e142b0 100644 --- a/ugali/utils/skymap.py +++ b/ugali/utils/skymap.py @@ -1,13 +1,11 @@ """ Tools for making maps of the sky with healpix. -""" -import sys -import re -import gc +ADW 2018-05-07: Largely deprecated? +""" -import numpy -import healpy +import numpy as np +import healpy as hp import ugali.utils.projector from ugali.utils.healpix import superpixel,subpixel,ang2pix,pix2ang,query_disc @@ -21,7 +19,7 @@ def surveyPixel(lon, lat, nside_pix, nside_subpix = None): Return the set of HEALPix pixels that cover the given coordinates at resolution nside. Optionally return the set of subpixels within those pixels at resolution nside_subpix """ - pix = numpy.unique(ang2pix(nside_pix, lon, lat)) + pix = np.unique(ang2pix(nside_pix, lon, lat)) if nside_subpix is None: return pix else: @@ -29,7 +27,7 @@ def surveyPixel(lon, lat, nside_pix, nside_subpix = None): for ii in range(0, len(pix)): subpix = subpixel(pix[ii], nside_pix, nside_subpix) subpix_array.append(subpix) - return pix, numpy.array(subpix_array) + return pix, np.array(subpix_array) def inFootprint(config, pixels, nside=None): """ @@ -41,18 +39,18 @@ def inFootprint(config, pixels, nside=None): nside_likelihood = config['coords']['nside_likelihood'] nside_pixel = config['coords']['nside_pixel'] - if numpy.isscalar(pixels): pixels = numpy.array([pixels]) + if np.isscalar(pixels): pixels = np.array([pixels]) if nside is None: nside = nside_likelihood filenames = config.getFilenames() catalog_pixels = filenames['pix'].compressed() - inside = numpy.zeros(len(pixels), dtype=bool) + inside = np.zeros(len(pixels), dtype=bool) if not nside_catalog: catalog_pix = [0] else: catalog_pix = superpixel(pixels,nside,nside_catalog) - catalog_pix = numpy.intersect1d(catalog_pix,catalog_pixels) + catalog_pix = np.intersect1d(catalog_pix,catalog_pixels) for fnames in filenames[catalog_pix]: logger.debug("Loading %s"%filenames['mask_1']) @@ -61,9 +59,9 @@ def inFootprint(config, pixels, nside=None): 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 = numpy.intersect1d(subpix_1,subpix_2) - superpix = numpy.unique(superpixel(subpix,nside_pixel,nside)) - inside |= numpy.in1d(pixels, superpix) + subpix = np.intersect1d(subpix_1,subpix_2) + superpix = np.unique(superpixel(subpix,nside_pixel,nside)) + inside |= np.in1d(pixels, superpix) return inside @@ -79,7 +77,7 @@ def footprint(config, nside=None): raise Exception('Requested nside=%i is greater than catalog_nside'%nside) elif nside > config['coords']['nside_pixel']: raise Exception('Requested nside=%i is less than pixel_nside'%nside) - pix = numpy.arange(healpy.nside2npix(nside), dtype=int) + pix = np.arange(hp.nside2npix(nside), dtype=int) return inFootprint(config,pix) @@ -89,261 +87,9 @@ def allSkyCoordinates(nside): """ Generate a set of coordinates at the centers of pixels of resolutions nside across the full sky. """ - lon,lat = pix2ang(nside, np.arange(0, healpy.nside2npix(nside))) + lon,lat = pix2ang(nside, np.arange(0, hp.nside2npix(nside))) return lon, lat -############################################################# -# -#def writeSparseHealpixMap(pix, data_dict, nside, outfile, -# distance_modulus_array = None, -# coordsys = 'NULL', ordering = 'NULL', -# header_dict = None): -# """ -# Sparse HEALPix maps are used to efficiently store maps of the sky by only -# writing out the pixels that contain data. -# -# Three-dimensional data can be saved by supplying a distance modulus array -# which is stored in a separate extension. -# -# coordsys [gal, cel] -# ordering [ring, nest] -# """ -# import astropy.io.fits as pyfits -# -# hdul = pyfits.HDUList() -# -# # Pixel data extension -# columns_array = [pyfits.Column(name = 'PIX', -# format = 'K', -# array = pix)] -# -# for key in data_dict.keys(): -# if data_dict[key].shape[0] != len(pix): -# logger.warning('First dimension of column %s (%i) does not match number of pixels (%i).'%(key, -# data_dict[key].shape[0], -# len(pix))) -# -# if len(data_dict[key].shape) == 1: -# columns_array.append(pyfits.Column(name = key, -# format = 'E', -# array = data_dict[key])) -# elif len(data_dict[key].shape) == 2: -# columns_array.append(pyfits.Column(name = key, -# format = '%iE'%(data_dict[key].shape[1]), -# array = data_dict[key])) -# else: -# logger.warning('Unexpected number of data dimensions for column %s.'%(key)) -# -# hdu_pix_data = pyfits.new_table(columns_array) -# hdu_pix_data.header.update([('NSIDE', nside)]) -# hdu_pix_data.header.update([('COORDSYS', coordsys.upper())]) -# hdu_pix_data.header.update([('ORDERING', ordering.upper())]) -# hdu_pix_data.header.update(header_dict) -# hdu_pix_data.name = 'PIX_DATA' -# hdul.append(hdu_pix_data) -# -# # Distance modulus extension -# if distance_modulus_array is not None: -# hdu_distance_modulus = pyfits.new_table([pyfits.Column(name = 'DISTANCE_MODULUS', -# format = 'E', -# array = distance_modulus_array)]) -# hdu_distance_modulus.name = 'DISTANCE_MODULUS' -# hdul.append(hdu_distance_modulus) -# -# hdul.writeto(outfile, clobber = True) -# -############################################################# -# -#def readSparseHealpixMap(infile, field, extension='PIX_DATA', default_value=healpy.UNSEEN, construct_map=True): -# """ -# Open a sparse HEALPix map fits file. -# Convert the contents into a HEALPix map or simply return the contents. -# Flexibility to handle -# """ -# import astropy.io.fits as pyfits -# -# reader = pyfits.open(infile,memmap=False) -# nside = reader[extension].header['NSIDE'] -# -# # Trying to fix avoid a memory leak -# try: -# pix = numpy.array(reader[extension].data.field('PIX'),copy=True) -# except: -# pix = numpy.array(reader[extension].data.field('PIXEL'),copy=True) -# value = numpy.array(reader[extension].data.field(field),copy=True) -# reader.close() -# -# if construct_map: -# if len(value.shape) == 1: -# map = default_value * numpy.ones(healpy.nside2npix(nside)) -# map[pix] = value -# else: -# map = default_value * numpy.ones([value.shape[1], healpy.nside2npix(nside)]) -# for ii in range(0, value.shape[1]): -# map[ii][pix] = numpy.take(value, [ii], axis=1) -# ret = map -# else: -# if len(value.shape) == 1: -# ret = (pix,value) -# else: -# ret = (pix,value.transpose()) -# return ret -# -############################################################# -# -#def readSparseHealpixMaps(infiles, field, extension='PIX_DATA', default_value=healpy.UNSEEN, construct_map=True): -# """ -# Read multiple sparse healpix maps and output the results -# identically to a single file read. -# """ -# import astropy.io.fits as pyfits -# -# if isinstance(infiles,str): infiles = [infiles] -# -# pix_array = [] -# value_array = [] -# -# # Create a map based on the first file in the list -# map = readSparseHealpixMap(infiles[0], field, extension=extension, default_value=healpy.UNSEEN, construct_map=True) -# -# for ii in range(0, len(infiles)): -# logger.debug('(%i/%i) %s'%(ii+1, len(infiles), infiles[ii])) -# pix_array_current, value_array_current = readSparseHealpixMap(infiles[ii], field, -# extension=extension, -# construct_map=False) -# pix_array.append(pix_array_current) -# value_array.append(value_array_current) -# map[pix_array[ii]] = value_array[ii] -# -# # Check to see whether there are any conflicts -# pix_master = numpy.concatenate(pix_array) -# value_master = numpy.concatenate(value_array) -# -# n_conflicting_pixels = len(pix_master) - len(numpy.unique(pix_master)) -# if n_conflicting_pixels != 0: -# logger.warning('%i conflicting pixels during merge.'%(n_conflicting_pixels)) -# -# if construct_map: -# return map -# else: -# if n_conflicting_pixels == 0: -# pix_master = numpy.sort(pix_master) -# return pix_master, map[pix_master] -# else: -# pix_valid = numpy.nonzero(map != default_value)[0] -# return pix_valid, map[pix_valid] -# -############################################################# -# -#def mergeSparseHealpixMaps(infiles, outfile=None, -# pix_data_extension='PIX_DATA', -# pix_field='PIX', -# distance_modulus_extension='DISTANCE_MODULUS', -# distance_modulus_field='DISTANCE_MODULUS', -# default_value=healpy.UNSEEN): -# """ -# Use the first infile to determine the basic contents to expect for the other files. -# """ -# import astropy.io.fits as pyfits -# -# # Setup -# if isinstance(infiles,str): infiles = [infiles] -# -# distance_modulus_array = None -# pix_array = [] -# data_dict = {} -# -# reader = pyfits.open(infiles[0]) -# nside = reader[pix_data_extension].header['NSIDE'] -# -# for ii in range(0, len(reader)): -# if reader[ii].name == distance_modulus_extension: -# distance_modulus_array = reader[distance_modulus_extension].data.field(distance_modulus_field) -# -# for key in reader[pix_data_extension].data.names: -# if key == pix_field: -# continue -# data_dict[key] = [] -# #if distance_modulus_array is None: -# # data_dict[key] = default_value * numpy.ones(healpy.nside2npix(nside)) -# #else: -# # data_dict[key] = default_value * numpy.ones([len(distance_modulus_array), -# # healpy.nside2npix(nside)]) -# reader.close() -# -# # Now loop over the infiles -# -# for ii in range(0, len(infiles)): -# logger.debug('(%i/%i) %s'%(ii+1, len(infiles), infiles[ii])) -# -# reader = pyfits.open(infiles[ii]) -# distance_modulus_array_current = numpy.array(reader[distance_modulus_extension].data.field(distance_modulus_field),copy=True) -# if not numpy.array_equal(distance_modulus_array_current,distance_modulus_array): -# logger.warning("Distance moduli do not match; skipping...") -# continue -# reader.close() -# -# pix_array_current = readSparseHealpixMap(infiles[ii], pix_field, -# extension=pix_data_extension, construct_map=False)[0] -# pix_array.append(pix_array_current) -# -# for key in data_dict.keys(): -# value_array_current = readSparseHealpixMap(infiles[ii], key, -# extension=pix_data_extension, construct_map=False)[1] -# data_dict[key].append(value_array_current) -# #if distance_modulus_array is None: -# # data_dict[key][pix_array_current] = value -# #else: -# # for jj in range(0, len(distance_modulus_array)): -# # data_dict[key][jj] = value[jj] -# -# gc.collect() -# -# pix_master = numpy.concatenate(pix_array) -# n_conflicting_pixels = len(pix_master) - len(numpy.unique(pix_master)) -# if n_conflicting_pixels != 0: -# logger.warning('%i conflicting pixels during merge.'%(n_conflicting_pixels)) -# -# for key in data_dict.keys(): -# if distance_modulus_array is not None: -# data_dict[key] = numpy.concatenate(data_dict[key], axis=1).transpose() -# else: -# data_dict[key] = numpy.concatenate(data_dict[key]) -# -# if outfile is not None: -# writeSparseHealpixMap(pix_master, data_dict, nside, outfile, -# distance_modulus_array=distance_modulus_array, -# coordsys='NULL', ordering='NULL') -# else: -# return data_dict -# -# -############################################################# -# -#def mergeLikelihoodFiles(infiles, lkhdfile, roifile): -# import astropy.io.fits as pyfits -# -# mergeSparseHealpixMaps(infiles,lkhdfile) -# -# ext='PIX_DATA' -# keys=['STELLAR','NINSIDE','NANNULUS'] -# nside = pyfits.open(infiles[0])[ext].header['LKDNSIDE'] -# -# pix_array = [] -# data_dict = dict([(k,[]) for k in keys]) -# for ii in range(0, len(infiles)): -# logger.debug('(%i/%i) %s'%(ii+1, len(infiles), infiles[ii])) -# reader = pyfits.open(infiles[ii]) -# pix_array.append(reader[ext].header['LKDPIX']) -# for key in data_dict.keys(): -# data_dict[key].append(reader[ext].header[key]) -# -# pix_array = numpy.array(pix_array) -# for key in data_dict.keys(): -# data_dict[key] = numpy.array(data_dict[key]) -# writeSparseHealpixMap(pix_array, data_dict, nside, roifile) -# -############################################################# def randomPositions(input, nside_pix, n=1): """ @@ -357,12 +103,12 @@ def randomPositions(input, nside_pix, n=1): so that gaps from star holes, bleed trails, cosmic rays, etc. are filled in. Return the longitude and latitude of the random positions (deg) and the total area (deg^2). """ - input = numpy.array(input) + input = np.array(input) if len(input.shape) == 1: - if healpy.npix2nside(len(input)) < nside_pix: + if hp.npix2nside(len(input)) < nside_pix: logger.warning('Expected coarser resolution nside_pix in skymap.randomPositions') - subpix = numpy.nonzero(input)[0] # All the valid pixels in the mask at the NSIDE for the input mask - lon, lat = pix2ang(healpy.npix2nside(len(input)), subpix) + subpix = np.nonzero(input)[0] # All the valid pixels in the mask at the NSIDE for the input mask + lon, lat = pix2ang(hp.npix2nside(len(input)), subpix) elif len(input.shape) == 2: lon, lat = input[0], input[1] # All catalog object positions else: @@ -370,28 +116,28 @@ def randomPositions(input, nside_pix, n=1): pix = surveyPixel(lon, lat, nside_pix) # Area with which the random points are thrown - area = len(pix) * healpy.nside2pixarea(nside_pix, degrees=True) + area = len(pix) * hp.nside2pixarea(nside_pix, degrees=True) # Create mask at the coarser resolution - mask = numpy.tile(False, healpy.nside2npix(nside_pix)) + mask = np.tile(False, hp.nside2npix(nside_pix)) mask[pix] = True # Estimate the number of points that need to be thrown based off # coverage fraction of the HEALPix mask - coverage_fraction = float(numpy.sum(mask)) / len(mask) + coverage_fraction = float(np.sum(mask)) / len(mask) n_throw = int(n / coverage_fraction) lon, lat = [], [] count = 0 while len(lon) < n: - lon_throw = numpy.random.uniform(0., 360., n_throw) - lat_throw = numpy.degrees(numpy.arcsin(numpy.random.uniform(-1., 1., n_throw))) + lon_throw = np.random.uniform(0., 360., n_throw) + lat_throw = np.degrees(np.arcsin(np.random.uniform(-1., 1., n_throw))) pix_throw = ugali.utils.healpix.angToPix(nside_pix, lon_throw, lat_throw) cut = mask[pix_throw].astype(bool) - lon = numpy.append(lon, lon_throw[cut]) - lat = numpy.append(lat, lat_throw[cut]) + lon = np.append(lon, lon_throw[cut]) + lat = np.append(lat, lat_throw[cut]) count += 1 if count > 10: @@ -409,25 +155,25 @@ def randomPositionsMask(mask, nside_pix, n): """ npix = len(mask) - nside = healpy.npix2nside(npix) + nside = hp.npix2nside(npix) # Estimate the number of points that need to be thrown based off # coverage fraction of the HEALPix mask - coverage_fraction = float(numpy.sum(mask)) / len(mask) + coverage_fraction = float(np.sum(mask)) / len(mask) n_throw = int(n / coverage_fraction) lon, lat = [], [] latch = True count = 0 while len(lon) < n: - lon_throw = numpy.random.uniform(0., 360., n_throw) - lat_throw = numpy.degrees(numpy.arcsin(numpy.random.uniform(-1., 1., n_throw))) + lon_throw = np.random.uniform(0., 360., n_throw) + lat_throw = np.degrees(np.arcsin(np.random.uniform(-1., 1., n_throw))) pix = ugali.utils.healpix.angToPix(nside, lon_throw, lat_throw) cut = mask[pix].astype(bool) - lon = numpy.append(lon, lon_throw[cut]) - lat = numpy.append(lat, lat_throw[cut]) + lon = np.append(lon, lon_throw[cut]) + lat = np.append(lat, lat_throw[cut]) count += 1 if count > 10: diff --git a/ugali/utils/stats.py b/ugali/utils/stats.py index 22c86f6..857a0f9 100644 --- a/ugali/utils/stats.py +++ b/ugali/utils/stats.py @@ -5,13 +5,13 @@ import copy -import numpy import numpy as np import numpy.lib.recfunctions as recfuncs from matplotlib import mlab import scipy.special import scipy.stats +# These should probably live in this file from ugali.utils.bayesian_efficiency import bayesianInterval, binomialInterval _alpha = 0.32 @@ -45,7 +45,7 @@ def median_interval(data, alpha=_alpha): Median including bayesian credible interval. """ q = [100*alpha/2., 50, 100*(1-alpha/2.)] - lo,med,hi = numpy.percentile(data,q) + lo,med,hi = np.percentile(data,q) return interval(med,lo,hi) def peak(data, bins=_nbins):