Skip to content

Commit

Permalink
Merge pull request #10 from kadrlica/cleanup
Browse files Browse the repository at this point in the history
Ok, this is probably enough for right now.
  • Loading branch information
kadrlica authored May 7, 2018
2 parents 349d82f + 41e51b5 commit 17d78e2
Show file tree
Hide file tree
Showing 16 changed files with 533 additions and 871 deletions.
15 changes: 7 additions & 8 deletions ugali/analysis/kernel.py
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -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).
Expand All @@ -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):
Expand All @@ -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)
Expand Down Expand Up @@ -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):
Expand All @@ -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)
Expand Down Expand Up @@ -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']:
Expand Down
61 changes: 30 additions & 31 deletions ugali/candidate/associate.py
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand All @@ -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/")
#
Expand All @@ -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)
#
Expand Down Expand Up @@ -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")
Expand All @@ -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):
Expand Down Expand Up @@ -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)
Expand All @@ -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']),'*')
Expand Down Expand Up @@ -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'])
Expand All @@ -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)
Expand All @@ -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)
Expand All @@ -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)
Expand All @@ -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
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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']
Expand All @@ -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)
Expand All @@ -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']
Expand Down
11 changes: 5 additions & 6 deletions ugali/isochrone/model.py
Original file line number Diff line number Diff line change
Expand Up @@ -37,7 +37,6 @@
import glob
from functools import wraps

import numpy
import numpy as np
import scipy.interpolate
import scipy.stats
Expand Down Expand Up @@ -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)
Expand Down Expand Up @@ -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())
Expand All @@ -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):
Expand Down
Loading

0 comments on commit 17d78e2

Please sign in to comment.