<jupyter_start><jupyter_text>Check Serial Times<jupyter_code>hptile = oss.HealpixTiles(nside=256,
                          preComputedMap='/Users/rbiswas/data/LSST/OpSimData/healpixelized_MINION_1016_256.db')
!rm simTiles.hdf
tstart = time.time()
tiles = [0]
for tileID in tiles:
    simTile = snsims.SimulationTile(snsims.GMM_SALT2Params,
                                    NSIDE=NSIDE,
                                    tileID=tileID,
                                    hpOpSim='/Users/rbiswas/data/LSST/OpSimData/healpixelized_MINION_1016_256.db',
                                    rate=snsims.PowerLawRates,
                                    allPointings=opsimOut,
                                    )
    simTile.writeTile(fileName='simTiles.hdf')
tend = time.time()
simTile.minPeakTime
simTile.snParams().t0.hist(histtype='step', bins=20)
plt.axvline(59800., ls='--', color='k', lw=2)
plt.axvline(59800. + 3650., ls='--', color='k', lw=2)
!du -h simTiles.hdf
df =  pd.read_hdf('simTiles.hdf')
df.snid.unique()
print('Writing the {0} tiles took {1} minutes'.format(len(tiles), (tend -tstart) / 60. ))
import line_profiler
%lprun 
import multiprocessing
from joblib import Parallel, delayed
help(snsims.SimulationTile)
simTile = snsims.SimulationTile(snsims.SimpleSALTDist, NSIDE=256, tileID=0, 
                      hpOpSim='/Users/rbiswas/data/LSST/OpSimData/healpixelized_MINION_1016_256.db',
                      rate=snsims.PowerLawRates,
                      allPointings=opsimOut)
fields = np.array(list(len(simTile.lc(i).lightCurve.fieldID.unique()) for i in range(125)))
fields.sort()
fields
np.bincount(fields)
simTile.writeSNParams('simTiles_0.hdf')
simTile.writeSN(0, 'simTile_SN_0.hdf')
def func(i):
    tileID = i
    simTile = snsims.SimulationTile(snsims.SimpleSALTDist,
                                    NSIDE=NSIDE,
                                    tileID=tileID,
                                    hpOpSim='/Users/rbiswas/data/LSST/OpSimData/healpixelized_MINION_1016_256.db',
                                    rate=snsims.PowerLawRates,
                                    allPointings=opsimOut)
    simTile.writeTile(fileName='simTiles_new_{}.hdf'.format(i))<jupyter_output><empty_output><jupyter_text>- WHERE is snid- WHy is fieldID there twice<jupyter_code>func(0)
!ls -ltr  simTiles_0.hdf
df = pd.read_hdf('simTiles_new_0.hdf')
df["SNR"] = df['flux']/df['fluxerr']
t = df.groupby(['snid', 'band']).agg(dict(SNR=max))
t.i
t
t.pivot('band', 'SNR')
snids = df.query('SNR > 5.').snid.unique()
paramsdf = pd.read_hdf('simTiles_new_0_params.hdf')
paramsdf.ix[snids].query('z > 0.5')
import analyzeSN as ans
import sncosmo
sncosmo.plot_lc(simTile.lc(0).snCosmoLC(), color='k')
df.head()
df.
%load_ext line_profiler

%prun func(0)
from math import sqrt
Parallel(n_jobs=2)(delayed(sqrt)(i ** 2) for i in range(10))

!rm simTiles*.hdf
tstart = time.time()
Parallel(n_jobs=1)(delayed(func)(i) for i in range(2))
tend = time.time()
print(tend -  tstart)/60.
simTile = snsims.SimulationTile(snsims.SimpleSALTDist,
                                NSIDE=NSIDE,
                                tileID=2,
                                hpOpSim='/Users/rbiswas/data/LSST/OpSimData/healpixelized_MINION_1016_256.db',
                                rate=snsims.PowerLawRates,
                                allPointings=opsimOut)
simTile.snParamTable.head()
lc = simTile.lc(2097153)

lc.lightCurve
lc.coaddedLC(coaddTimes=1.)
lc.snCosmoLC()
import analyzeSN as ans
simTile.writeTile(fileName='SimTile-2.hdf')
tend = time.time()
print(tend - tstart)
len(simTile.snParamTable)
df = pd.read_hdf('simTiles.hdf')
df.fieldID.unique()
df['SNR'] = df.flux / df.fluxerrs
totalCounts, edges = np.histogram(simTile.snParamTable.ix[df.query('SNR > 0').snid.unique()].z, bins=np.arange(0., 1.4, 0.1))
fiveSigmaCounts, edges = np.histogram(simTile.snParamTable.ix[df.query('SNR > 5').snid.unique()].z, bins=np.arange(0., 1.4, 0.1))
fig, ax = plt.subplots()
ax.plot(0.5 * (edges[: -1] + edges[1:] ), fiveSigmaCounts / totalCounts.astype(np.float), 'o')<jupyter_output><empty_output><jupyter_text>Scratch<jupyter_code>#/ np.float(totalCounts)
df.snid.unique().size * 200 /10000
df.expMJD.max() - df.expMJD.min()
class SimulationTile(snsims.Universe):
    def __init__(self, paramDist, rate, NSIDE, tileID, hpOpSim, allPointings=None, timeRange=None):
        self._randomState = None
        self.Tiling = snsims.HealpixTiles(nside=NSIDE, preComputedMap=hpOpSim)
        self.tileID = tileID
        self.fieldArea = self.Tiling.area(tileID)
        self.zdist = rate(rng=self.randomState, fieldArea=self.fieldArea)
        self.zsamples = self.zdist.zSamples
        self.numSN = len(self.zsamples)
        self.positions = self.Tiling.positions(self.tileID, self.numSN, rng=self.randomState)
        self._snParamTable = None
        self.columns = ('expMJD', 'filter', 'fieldID', 'fiveSigmaDepth')
        self.tilePointings = self.Tiling.pointingSequenceForTile(self.tileID, 
                                                                 allPointings=allPointings, columns=self.columns)
        self._timeRange = timeRange
        self.bandPasses = BandpassDict.loadTotalBandpassesFromFiles()
        
    @property
    def minPeakTime(self):
        if  self._timeRange is None:
            minTime = self.tilePointings.expMJD.min()
        else:
            minTime = self._timeRange[0]
        return minTime
    
    @property
    def maxPeakTime(self):
        if  self._timeRange is None:
            maxTime = self.tilePointings.expMJD.max()
        else:
            maxTime = self._timeRange[1]
        return maxTime
    
    @property
    def snParamTable(self):
        if self._snParamTable is None:
            self.snParams()
        return self._snParamTable
    @property
    def randomState(self):
        if self._randomState is None:
            self._randomState = np.random.RandomState(self.tileID)
        return self._randomState
    def snParams(self):
        zsamples = self.zdist.zSamples
        numSN = len(zsamples)
        positions = self.Tiling.positions(self.tileID, numSN, rng=self.randomState)
        ra = self.positions[0]
        dec = - self.positions[1] + 45.0 
        # Why do we need numSN
        sp = snsims.SimpleSALTDist(numSN=numSN, rng=self.randomState, zSamples=self.zsamples).paramSamples
        sp['ra'] = self.positions[0]
        sp['dec'] = self.positions[1]
        sp['snid'] = np.left_shift(self.tileID, 20) + np.arange(numSN)
        sp.set_index('snid', inplace=True)
        self._snParamTable = sp
        
        if self.minPeakTime is None or self.maxPeakTime is None:
            pass
        else:
            sp['t0'] = self.minPeakTime + (self.maxPeakTime - self.minPeakTime) * sp['t0']
        return sp
    
    @staticmethod
    def getSNCosmoParamDict(odict, SNCosmoModel):
        mydict = dict()
        param_names = SNCosmoModel.param_names
        for param in odict.index.values:
            if param in param_names:
                mydict[param] = odict[param]
        return mydict
                
    def SN(self, snid, timeRange='model'):
        mySNParams = self.snParamTable.ix[snid]
        sn = SNObject(ra=mySNParams.ra, dec=mySNParams.dec)
        #print mySNParams
        sncosmo_params = self.getSNCosmoParamDict(mySNParams, sn)
        #print(sncosmo_params)
        sn.set(**sncosmo_params)
        z = sn.get('z')
        t0 = sn.get('t0')
        lcMinTime = t0 - 20. * (1.0 + z)
        lcMaxTime = t0 + 50. * (1.0 + z )
        df = self.tilePointings.query('expMJD < @lcMaxTime and expMJD > @lcMinTime')
        df['snid'] = snid
        fluxes = []
        fluxerrs = []
        for rows in df.iterrows():
            row = rows[1]
            # print(row['expMJD'], row['filter'], row['fiveSigmaDepth'])
            bp = self.bandPasses[row['filter']]
            flux = sn.catsimBandFlux(bandpassobject=bp, time=row['expMJD'])
            fluxerr = sn.catsimBandFluxError(time=row['expMJD'], bandpassobject=bp, fluxinMaggies=flux,
                                             m5=row['fiveSigmaDepth'])
            fluxes.append(flux)
            fluxerrs.append(fluxerr)
        df['flux'] = fluxes
        df['fluxerrs'] = fluxerrs
        return sn, df
    def writeTile(self, fileName, timeRange='model'):
        count = 0
        for snid in self.snParamTable.index.values:
            self.writeSN(snid, fileName, timeRange=timeRange)
            if count % 50 == 0:
                print('another 50', snid)
            count += 1
    def writeSN(self, snid, fileName, timeRange='model'):
        sn, df = self.SN(snid, timeRange)
        df['filter'] = df['filter'].astype(str)
        with pd.get_store(fileName) as store:
        # with pd.HDFStore(fileName, mode='a') as hdf:
        #    hdf.append(df, format='table')
             store.append('tile_{}'.format(self.tileID), df)
        #df.to_hdf(fileName, key='tile_{}'.format(self.tileID), mode='a', format='t')
hpOpSim = oss.HealPixelizedOpSim.fromOpSimDB(opsimdb, NSIDE=256)
hpTiles = snsims.HealpixTiles(nside=16, healpixelizedOpSim=hpOpSim)
hpTiles.positions(2, 100, rng=np.random.RandomState(2))
mapvalues = np.ones(hp.nside2npix(256)) * hp.UNSEEN
mapvalues[100] = 100
r, t = hp.pix2ang(nside=256, ipix=100, nest=True)
np.degrees(r)
np.degrees(t)
hp.mollview(mapvalues)
hp.projscatter(r, t)
hp.projscatter(theta, phi)
phi, theta = hpTiles.positions(100, 100)
phi = np.radians(phi)
theta = np.radians(theta)
hp
hpTiles.area(tileID=23)/ 41252.96
hpTiles.positions()
import healpy as hp
zdist.zSamples
zdist.numSN()
zdist.zbinEdges
zdist.zSampleSize()
from future.builtins import range
sd = snsims.SimpleSALTDist(rng=np.random.RandomState(43), numSN=3, zSamples=[0.03, 0.2, 0.5])
sd.numSN
sd = snsims.SimpleSALTDist(3, zSamples=[0.2, 0.5, 0.7], rng=np.random.RandomState(2))
sd.paramSamples
hp.nside2npix(256)
np.binary_repr(hp.nside2npix(256))
np.binary_repr(1000)
np.binary_repr(x)
x = hp.nside2npix(16) + np.left_shift(1000, 20)
hp.nside2npix(256)
np.log2(hp.nside2npix(256))
from lsst.sims.catUtils.supernovae import SNObject
NSIDE=256
tileSim = snsims.HealpixTiles(nside=NSIDE, healpixelizedOpSim=hpOpSim)
tileSim.positions(2, 5, rng=np.random.RandomState(2))
SNObject(ra = sn)

from lsst.sims.photUtils import BandpassDict
NSIDE=256
opsimout = oss.OpSimOutput.fromOpSimDB(opsimdb)
!rm simulatedTile-2.hdf
import time
tstart = time.time()
ss = simTile.snParamTable
ss.tail()
simTile.writeTile(fileName='simulatedTile-2.hdf')
# simTile.writeSN(snid=2097152, fileName='simulatedTile-2.hdf')
pd.read_hdf('simulatedTile-2.hdf')
sn, lc = simTile.SN(2097152)
tend = time.time()
tend - tstart
len(simTile.snParamTable)
data = pd.read_hdf('simulatedTile-2.hdf', 'tile_2')
len(data.snid.unique())
len(data)
glc
glc.fluxerrs
fig, ax = plt.subplots()
glc = lc.query('filter == "i"')
print(len(glc))
ax.errorbar(glc.expMJD, glc.flux, yerr=glc.fluxerrs.values, fmt='o')
mySN = SNObject(ra=44.78, dec=0.28)
simTile.snParamTable.ix[2097152]
for param in simTile.snParamTable.ix[2097152]:
    print param
mySN.param_names
simTile.getSNCosmoParamDict(simTile.snParamTable.ix[2097152], mySN)
ht = snsims.HealpixTiles(nside=256, 
                         preComputedMap='/Users/rbiswas/.local/lib/python2.7/site-packages/opsimsummary/example_data/healpixels_micro.db')
subset = ht.pointingSequenceForTile(allPointings=opsimOut, columns=('expMJD', 'filter', 'fiveSigmaDepth'), tileID=2)
subset.head()
row = next(subset.iterrows())
    
for rows in subset.head().iterrows():
    row = rows[1]
    print(row['expMJD'], row['filter'], row['fiveSigmaDepth'])
row[1]['expMJD']
opsimdb
hpopSim = oss.HealPixelizedOpSim
from lsst.sims.photUtils import BandpassDict
opsimOut
ss.head()
from lsst.sims.catUtils.supernovae import SNObject
sn, df = simTile.SN(2097152)
sn.SNstate
simTile.tilePointings.head()
simTile.tilePointings
opsimout.summary.ix[pointings].fieldID.unique()
opsimout.summary.ix[pointings][['fieldRA', 'fieldDec', 'fieldID', 'fiveSigmaDepth', 'expMJD', 'filter']]
pointings = simTile.Tiling.pointingSequenceForTile(2, allPointings=None)
ss.head()
snObj = SNObject(ra=ss.ix[2097152].ra, dec=ss.ix[2097152].dec)
snObj.set(ss.ix[209])
mydict = dict()
for key in ss.ix[2097152].index.values:
    if key in snObj.param_names:
        mydict[key] = ss.ix[2097152][key]
ss.ix[2097152].index.values
mydict
snObj.set(**mydict)
ss.z.hist(histtype='step')
plt.hexbin(ss.ra, ss.dec)
ss.dec.hist(histtype='step')
ss.ra.ix[0]
sp.paramSamples
simTile.numSN
simTile.snParams
simTile.zsamples.max()
plt.hist(simTile.zsamples, bins=20, histtype='step')
mapvals = np.ones(hp.nside2npix(NSIDE))* hp.UNSEEN
mapvals[2] =1
hp.mollview(mapvals, nest=True)
hp.projscatter(np.radians(simTile.positions[1][:50]), np.radians(simTile.positions[0][:50]), **dict(s=0.01))
class TestUniverse(snsims.Universe()):
    pass
zdist = snsims.PowerLawRates(rng=np.random.)
np.pi * 4 * np.degrees(1.)**2
1.0/ 12.0/ 256**2
rng = np.random.RandomState()
hpTiles.nside
p, t = hpTiles.positions(23, 20, rng)
p
t
hp.pix2ang(256, 23, nest=True)
hpTiles._angularSamples(phi_c, theta_c, 0.13, 20, 23, rng)
phi,theta = hpTiles.samplePatchOnSphere(phi_c, theta_c, 0.13, 20, rng, degrees=False)
phi
phi_c
hp.ang2pix(256, theta, phi, nest=True)
import healpy as hp
theta_c, phi_c = np.degrees(hp.pix2ang(nside=256, ipix=23, nest=True))
theta_c
np.sqrt(hp.nside2pixarea(256, degrees=True) / np.pi)
import healpy as hp
hp.nside2pixarea(256, degrees=True)
hpTiles.pointingSequenceForTile(tileID=23, allPointings=None)
import healpy as hp
from lsst.sims.utils import samplePatchOnSphere
samplePatchOnSphere()
def samplePatchOnSphere(phi, theta, delta, size):
    """
    Uniformly distributes samples on a patch on a sphere between phi \pm delta,
    and theta \pm delta on a sphere. Uniform distribution implies that the
    number of points in a patch of sphere is proportional to the area of the
    patch. Here, the coordinate system is the usual
    spherical coordinate system but with the azimuthal angle theta going from
    90 degrees at the North Pole, to -90 degrees at the South Pole, through
    0. at the equator. 
    
    This function is not equipped to handle wrap-around the ranges of theta
    phi and therefore does not work at the poles.
 
    Parameters
    ----------
    phi: float, mandatory, degrees
	center of the spherical patch in ra with range 
    theta: float, mandatory, degrees
    delta: float, mandatory, degrees
    size: int, mandatory
        number of samples
    seed : int, optional, defaults to 1
        random Seed used for generating values
    Returns
    -------
    tuple of (phivals, thetavals) where phivals and thetavals are arrays of 
        size size in degrees.
    """
    u = rng.uniform(size=size)
    v = rng.uniform(size=size)
    phi = np.radians(phi)
    theta = np.radians(theta)
    delta = np.radians(delta)

    phivals = 2. * delta * u + (phi - delta)
    phivals = np.where ( phivals >= 0., phivals, phivals + 2. * np.pi)
    
    # use conventions in spherical coordinates
    theta = np.pi/2.0 - theta
 
    thetamax = theta + delta
    thetamin = theta - delta

    if thetamax > np.pi or thetamin < 0. :
        raise ValueError('Function not implemented to cover wrap around poles')

    # Cumulative Density Function is cos(thetamin) - cos(theta) / cos(thetamin) - cos(thetamax)
    a = np.cos(thetamin) - np.cos(thetamax)
    thetavals = np.arccos(-v * a + np.cos(thetamin))

    # Get back to -pi/2 to pi/2 range of decs
    thetavals = np.pi/2.0 - thetavals 
    return np.degrees(phivals) , np.degrees(thetavals)
class CoordSamples(snsims.PositionSamples, snsims.HealpixTiles):
    def __init__(self, nside, hpOpSim, rng):
        self.nside = nside
        super(self.__class__, self).__init__(nside=nside, healpixelizedOpSim=hpOpSim)
        self._rng = rng
    @property
    def randomState(self):
        if self._rng is None:
            raise ValueError('self._rng should not be None')
        return self._rng
    def positions(self, tileID, numSamples):
        res_phi = 
        ang = hp.pix2ang(nside=self.nside, ipix=tileID, nest=True)
        radius = np.degrees(np.sqrt(hp.nside2pixarea / np.pi))
        phi_c, theta_c = np.degrees(ang[::-1])
        phi, theta = super(self.__class__).samplePatchOnSphere(phi_c, theta_c, delta=radius, size=numSamples)
        
        return ang
ang[::-1]
rng = np.random.RandomState(23)
cs = CoordSamples(nside=256, hpOpSim=hpOpSim, rng=rng)
ang = cs.positions(23, 20)
area = hp.nside2pixarea(256)
side = np.sqrt(area / np.pi) 
res = samplePatchOnSphere(phi=ang[1], theta=ang[0], delta=side, size=100)
hp.ang2pix(256, theta=res[1], phi=res[0], nest=True)
snsims.SimulatedTile(tileID=23)
class SimTile(snsims.SimulatedTile):
    pass
simulatedTile = SimTile(tileID=23)
from future.utils import with_metaclass
import abc
class TestDist(with_metaclass(abc.ABCMeta, snsims.Universe)):
    @property
    def randomState(self):
        if self._rng is None:
            self._rng = np.randomState(self.tileID)
        return self._rng
class ConcreteImplementation(TestDist):
    def __init__(self, tileID, NSIDE=1):
        hpOpSim = oss.HealPixelizedOpSim.fromOpSimDB(NSIDE=NSIDE, opSimDBpath=opsimdb)
        self.hpTile = snsims.HealpixTiles(nside=NSIDE, healpixelizedOpSim=hpOpSim)
        self.tileID = tileID
        self._rng = None
    @property
    def randomState(self):
        if self._rng is None:
            self._rng = np.randomState(self.tileID)
        return self._rng

    @property
    def snParams(self):
        
imp = ConcreteImplementation()
import numpy as np
rng = np.random.RandomState(1)
rng.poisson(lam=np.array([3405.3, 34., 21.2]))
state =rng.get_state()
state
rng.uniform()
zbin_edges = np.array([0., 1.0, 2.0, 3.0])
zbin_edges
x = zbin_edges[:-1]
y = zbin_edges[1:]
y
zip(x , y)
zz = (34, 23, 56)
map(len, list(rng.uniform(low=xx, high=yy, size=zzz) for (xx, yy, zzz) in zip(x,y, zz)))
arr = list(rng.uniform(low=xx, high=yy, size=zzz).tolist() for (xx, yy, zzz) in zip(x,y, zz))
list(_x for lst in arr for _x in lst)
mylist = []
for lst in arr:
    for x in lst:
        mylist.append(x)
len(mylist)
mylist
list(x for x in )
np.shape(arr)
np.shape(arr)
class coordSamples(snsims.HealpixTiles):
    pass
from astropy.cosmology import Planck15
Planck15.comoving_volume(z=1.4) - Planck15.comoving_volume(z=1.3)
Planck15.comoving_volume(z=0.1) - Planck15.comoving_volume(z=0.01)
from lsst.sims.catUtils.supernovae import SNObject
snobj = SNObject(ra=54., dec=-47.)
snobj.set(**dict(x0=1.4e-6, x1-))
t = ('expMJD', 'filter', 'fiveSigmaDepth')
list(t)
import healpy as hp
hp.ud_grade(np.arange(hp.nside2npix(16)), nside_out=256)
8189/60/60
8202/ 60.
189216030.0 /365/60/60 /24
94608030.0 /365/60/60 /24 
189216060.0 / 365/60/60 /24 
283824060.0 / 365/60/60 /24 
94608000.0/ 365/60/60 /24 
l = [3, 4, 5]<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.13307964695312088, 'fraction_numerical': 0.034194809999471484, 'mean_word_length': 3.163256325632563, 'pattern_counts': {'":': 0, '<': 15, '<?xml version=': 0, '>': 20, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 15, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '3', 'dir': 'jupyter-structured-clean-dedup', 'id': '535929', 'license': 'MIT', 'n_tokens_mistral': 7484, 'n_tokens_neox': 6977, 'n_words': 1327, 'path': 'examples/TestingSims.ipynb', 'repo_name': 'rbiswas4/SNsims', 'size': '18921'}