In [1]:

    

import numpy as np
import matplotlib.pyplot as plt
import wisps
import wisps.simulations as wispsim
import matplotlib as mpl
import astropy.units as u
from astropy.coordinates import SkyCoord
import theano
import theano.tensor as tt
import pandas as pd
import pymc3 as pm
import seaborn as sns 
%matplotlib inline


    

    




Warning: spectrum object has a flux vector of zero length - maybe empty?

Warning: normalize is attempting to divide by nan; ignoring

Warning: spectrum object has a flux vector of zero length - maybe empty?

Warning: normalize is attempting to divide by nan; ignoring

In [2]:

    

COORDS=SkyCoord([p.coord for p in wisps.OBSERVED_POINTINGS ])


    

In [3]:

    

galc=COORDS.transform_to('galactic')


    

In [4]:

    

LS=galc.l.radian


    

In [5]:

    

BS=galc.b.radian


    

In [6]:

    

Rsun=8300.
Zsun=27.


    

In [7]:

    

#-------------------------------------------
def density_function(r, z, h):
    
    """
    A custom juric density function that only uses numpy arrays for speed
    All units are in pc
    """
    l = 2600. # radial length scale of exponential thin disk 
    zpart=np.exp(-abs(z-Zsun)/h)
    rpart=np.exp(-(r-Rsun)/l)
    return zpart*rpart


    

In [8]:

    

np.min(abs(BS))


    

    
Out[8]:





0.35233193655075173

In [9]:

    

def sample_distances(dmax, nsample=1000):
    """
    sample the galaxy given a scale height
    
    """

    def logprior(l, b):
        return tt.switch(( abs(b) < 0.35),-np.inf, 0)
        
    def loglike(r, z, d, h):
        return np.log10((d**2)*density_function(r, z, h))
    
    def logp(l, b, r, z, d, h):
        return np.log10((d**2)*density_function(r, z, h))

    with pm.Model() as model:
        h=350.
        l=pm.Uniform('l', lower=-2*np.pi, upper=2*np.pi, testval=np.pi/2, observed=LS)
        b=pm.Uniform('b', lower=-2*np.pi, upper=2*np.pi, testval=np.pi/3, observed=BS)
    
        d=pm.Uniform('d', lower=0., upper=dmax, testval=dmax/2,shape=BS.shape)
        
        x=pm.Deterministic('x',  Rsun-d*np.cos(b)*np.cos(l))
        y=pm.Deterministic('y', -d*np.cos(b)*np.sin(l))
        r=pm.Deterministic('r', (x**2+y**2)**0.5 )
        z=pm.Deterministic('z', Zsun+ d * np.sin(b))
        
        like = pm.DensityDist('likelihood', logp, observed={'l':l, 'b':b,
                             'r': r, 'z': z, 'd':d, 'h':h})

        trace = pm.sample(draws=int(nsample), cores=2)
    return trace, model


    

In [10]:

    

tr, model=sample_distances(16000, nsample=100)


    

    




Only 100 samples in chain.
Auto-assigning NUTS sampler...
Initializing NUTS using jitter+adapt_diag...
Multiprocess sampling (2 chains in 2 jobs)
NUTS: [d]
Sampling 2 chains: 100%|██████████| 1200/1200 [00:17<00:00, 69.11draws/s] 

In [11]:

    

fig, ax=plt.subplots(ncols=2, figsize=(8, 4))
h=ax[0].hist(tr.r.flatten(), histtype='step', normed=True, bins='auto', alpha=0.5)
h=ax[1].hist(tr.z.flatten(), histtype='step', normed=True, bins='auto', alpha=0.5)
ax[0].set_xlabel('r (pc)',fontsize=18)
ax[1].set_xlabel('z (pc)',fontsize=18)


plt.tight_layout()


    

In [12]:

    

tr.d.shape


    

    
Out[12]:





(200, 533)

In [13]:

    

gc=SkyCoord(l=LS*u.rad, b=BS*u.rad, distance=tr.d*u.pc, frame='galactic')


    

In [14]:

    

fig = plt.figure(figsize=(8,6))
ax = fig.add_subplot(111, projection="mollweide")
ax.scatter(gc.l.wrap_at(180*u.degree).radian,gc.b.wrap_at(90*u.degree).radian, marker='+', alpha=0.1)
ax.set_xlabel('l (deg)', fontsize=18)
ax.set_ylabel('b (deg)', fontsize=18)
plt.grid()


    

In [15]:

    

kde1=wisps.kernel_density([tr.r.flatten(), tr.z.flatten()])
kde2=wisps.kernel_density([tr.x.flatten(), tr.y.flatten()])


    

In [17]:

    

#dens1=kde1.pdf([tr.r.flatten(), tr.z.flatten()])
#dens2=kde2.pdf([tr.x.flatten(), tr.y.flatten()])


    

In [ ]:

    

tr.r.flatten().shape


    

In [ ]:

    

#cmap= sns.diverging_palette(100, 300, s=80, l=55, n=19, as_cmap=True)


    

In [19]:

    

#fig, ax=plt.subplots(ncols=2, figsize=(12, 4))

#c=ax[0].scatter(tr.x.flatten(), tr.y.flatten(), c=dens2, s=15.,  alpha=0.5)
#ax[0].set_xlabel('x (pc)', fontsize=18)
#ax[0].set_ylabel('y (pc)', fontsize=18)

#c=ax[1].scatter(tr.r.flatten(), tr.z.flatten(),c=dens1, s=15.,  alpha=0.5)
#ax[1].set_xlabel('r (pc)', fontsize=18)
#ax[1].set_ylabel('z (pc)', fontsize=18)


#ax[0].minorticks_on()
#ax[1].minorticks_on()



#plt.tight_layout()


    

In [20]:

    

import itertools


    

In [21]:

    

dmaxses=[50, 400, 1000, 2000, 5000]
lbses=[(0, 0), (0, np.pi/2), (0, -np.pi/2), 
       (-np.pi, 0.0), (np.pi, np.pi/2), (np.pi, -np.pi/2),
      (np.pi, 0.0), (np.pi, np.pi/2), (np.pi, -np.pi/2)]


    

In [22]:

    

c = list(itertools.product(dmaxses, lbses))


    

In [23]:

    

c[0]


    

    
Out[23]:





(50, (0, 0))

In [24]:

    

import corner


    

In [25]:

    

import astropy.units as u


    

In [26]:

    

from scipy.stats import gaussian_kde


    

In [27]:

    

def sample_distances_coord(dmax, lb, nsample=1000):
    """
    sample the galaxy given a scale height
    
    """

    def logprior(l, b):
        return tt.switch(( abs(b) < 0.35),-np.inf, 0)
        
    def loglike(r, z, d, h):
        return np.log10((d**2)*density_function(r, z, h))
    
    def logp(l, b, r, z, d, h):
        return np.log10((d**2)*density_function(r, z, h))

    with pm.Model() as model:
        l, b=lb
        h=350.
        d=pm.Uniform('d', lower=0., upper=dmax, testval=dmax/2)
        
        x=pm.Deterministic('x',  Rsun-d*np.cos(b)*np.cos(l))
        y=pm.Deterministic('y', -d*np.cos(b)*np.sin(l))
        r=pm.Deterministic('r', (x**2+y**2)**0.5 )
        z=pm.Deterministic('z', Zsun+ d * np.sin(b))
        
        like = pm.DensityDist('likelihood', logp, observed={'l':l, 'b':b,
                             'r': r, 'z': z, 'd':d, 'h':h})

        trace = pm.sample(draws=int(nsample), cores=2)
    return trace, model


    

In [28]:

    

for comb in c:
    tr, model=sample_distances_coord(  comb[0], comb[1], nsample=1000)
    fig, ax=plt.subplots(ncols=2, figsize=(8, 4))
    h=ax[0].hist(tr.d, histtype='step', normed=True, bins='auto')
    h=ax[1].hist(tr.z, histtype='step', normed=True, bins='auto')
    ax[0].set_xlabel('d (pc)',fontsize=18)
    ax[1].set_xlabel('z (pc)',fontsize=18)
    ax[0].spines['top'].set_visible(False)
    ax[1].spines['top'].set_visible(False)
    plt.tight_layout()
    fig.suptitle("dmax = {} l= {} b={}".format(comb[0]*u.pc, (comb[1][0]*u.rad).to(u.deg),(comb[1][1]*u.rad).to(u.deg)), 
                 fontsize=18)
    plt.tight_layout()


    

    




Auto-assigning NUTS sampler...
Initializing NUTS using jitter+adapt_diag...
Multiprocess sampling (2 chains in 2 jobs)
NUTS: [d]
Sampling 2 chains: 100%|██████████| 3000/3000 [00:02<00:00, 1437.22draws/s]
The acceptance probability does not match the target. It is 0.8849133318560523, but should be close to 0.8. Try to increase the number of tuning steps.
Auto-assigning NUTS sampler...
Initializing NUTS using jitter+adapt_diag...
Multiprocess sampling (2 chains in 2 jobs)
NUTS: [d]
Sampling 2 chains: 100%|██████████| 3000/3000 [00:02<00:00, 1305.25draws/s]
The acceptance probability does not match the target. It is 0.8805856705290144, but should be close to 0.8. Try to increase the number of tuning steps.
Auto-assigning NUTS sampler...
Initializing NUTS using jitter+adapt_diag...
Multiprocess sampling (2 chains in 2 jobs)
NUTS: [d]
Sampling 2 chains: 100%|██████████| 3000/3000 [00:02<00:00, 1465.34draws/s]
Auto-assigning NUTS sampler...
Initializing NUTS using jitter+adapt_diag...
Multiprocess sampling (2 chains in 2 jobs)
NUTS: [d]
Sampling 2 chains: 100%|██████████| 3000/3000 [00:01<00:00, 1727.53draws/s]
Auto-assigning NUTS sampler...
Initializing NUTS using jitter+adapt_diag...
Multiprocess sampling (2 chains in 2 jobs)
NUTS: [d]
Sampling 2 chains: 100%|██████████| 3000/3000 [00:01<00:00, 1633.06draws/s]
Auto-assigning NUTS sampler...
Initializing NUTS using jitter+adapt_diag...
Multiprocess sampling (2 chains in 2 jobs)
NUTS: [d]
Sampling 2 chains: 100%|██████████| 3000/3000 [00:02<00:00, 1264.76draws/s]
The acceptance probability does not match the target. It is 0.8861868896601115, but should be close to 0.8. Try to increase the number of tuning steps.
Auto-assigning NUTS sampler...
Initializing NUTS using jitter+adapt_diag...
Multiprocess sampling (2 chains in 2 jobs)
NUTS: [d]
Sampling 2 chains: 100%|██████████| 3000/3000 [00:01<00:00, 1826.16draws/s]
Auto-assigning NUTS sampler...
Initializing NUTS using jitter+adapt_diag...
Multiprocess sampling (2 chains in 2 jobs)
NUTS: [d]
Sampling 2 chains: 100%|██████████| 3000/3000 [00:01<00:00, 1543.07draws/s]
The acceptance probability does not match the target. It is 0.8831658849024592, but should be close to 0.8. Try to increase the number of tuning steps.
Auto-assigning NUTS sampler...
Initializing NUTS using jitter+adapt_diag...
Multiprocess sampling (2 chains in 2 jobs)
NUTS: [d]
Sampling 2 chains: 100%|██████████| 3000/3000 [00:01<00:00, 1565.56draws/s]
Auto-assigning NUTS sampler...
Initializing NUTS using jitter+adapt_diag...
Multiprocess sampling (2 chains in 2 jobs)
NUTS: [d]
Sampling 2 chains: 100%|██████████| 3000/3000 [00:01<00:00, 1693.15draws/s]
Auto-assigning NUTS sampler...
Initializing NUTS using jitter+adapt_diag...
Multiprocess sampling (2 chains in 2 jobs)
NUTS: [d]
Sampling 2 chains: 100%|██████████| 3000/3000 [00:01<00:00, 1842.77draws/s]
Auto-assigning NUTS sampler...
Initializing NUTS using jitter+adapt_diag...
Multiprocess sampling (2 chains in 2 jobs)
NUTS: [d]
Sampling 2 chains: 100%|██████████| 3000/3000 [00:01<00:00, 1770.93draws/s]
Auto-assigning NUTS sampler...
Initializing NUTS using jitter+adapt_diag...
Multiprocess sampling (2 chains in 2 jobs)
NUTS: [d]
Sampling 2 chains: 100%|██████████| 3000/3000 [00:01<00:00, 1788.11draws/s]
Auto-assigning NUTS sampler...
Initializing NUTS using jitter+adapt_diag...
Multiprocess sampling (2 chains in 2 jobs)
NUTS: [d]
Sampling 2 chains: 100%|██████████| 3000/3000 [00:01<00:00, 2154.12draws/s]
The acceptance probability does not match the target. It is 0.8801158297180247, but should be close to 0.8. Try to increase the number of tuning steps.
Auto-assigning NUTS sampler...
Initializing NUTS using jitter+adapt_diag...
Multiprocess sampling (2 chains in 2 jobs)
NUTS: [d]
Sampling 2 chains: 100%|██████████| 3000/3000 [00:01<00:00, 2428.35draws/s]
Auto-assigning NUTS sampler...
Initializing NUTS using jitter+adapt_diag...
Multiprocess sampling (2 chains in 2 jobs)
NUTS: [d]
Sampling 2 chains: 100%|██████████| 3000/3000 [00:01<00:00, 2326.70draws/s]
Auto-assigning NUTS sampler...
Initializing NUTS using jitter+adapt_diag...
Multiprocess sampling (2 chains in 2 jobs)
NUTS: [d]
Sampling 2 chains: 100%|██████████| 3000/3000 [00:01<00:00, 2275.01draws/s]
Auto-assigning NUTS sampler...
Initializing NUTS using jitter+adapt_diag...
Multiprocess sampling (2 chains in 2 jobs)
NUTS: [d]
Sampling 2 chains: 100%|██████████| 3000/3000 [00:01<00:00, 2002.50draws/s]
Auto-assigning NUTS sampler...
Initializing NUTS using jitter+adapt_diag...
Multiprocess sampling (2 chains in 2 jobs)
NUTS: [d]
Sampling 2 chains: 100%|██████████| 3000/3000 [00:02<00:00, 1280.19draws/s]
Auto-assigning NUTS sampler...
Initializing NUTS using jitter+adapt_diag...
Multiprocess sampling (2 chains in 2 jobs)
NUTS: [d]
Sampling 2 chains: 100%|██████████| 3000/3000 [00:01<00:00, 1980.84draws/s]
Auto-assigning NUTS sampler...
Initializing NUTS using jitter+adapt_diag...
Multiprocess sampling (2 chains in 2 jobs)
NUTS: [d]
Sampling 2 chains: 100%|██████████| 3000/3000 [00:01<00:00, 1973.23draws/s]
Auto-assigning NUTS sampler...
Initializing NUTS using jitter+adapt_diag...
Multiprocess sampling (2 chains in 2 jobs)
NUTS: [d]
Sampling 2 chains: 100%|██████████| 3000/3000 [00:01<00:00, 2389.14draws/s]
Auto-assigning NUTS sampler...
Initializing NUTS using jitter+adapt_diag...
Multiprocess sampling (2 chains in 2 jobs)
NUTS: [d]
Sampling 2 chains: 100%|██████████| 3000/3000 [00:01<00:00, 1862.76draws/s]
Auto-assigning NUTS sampler...
Initializing NUTS using jitter+adapt_diag...
Multiprocess sampling (2 chains in 2 jobs)
NUTS: [d]
Sampling 2 chains: 100%|██████████| 3000/3000 [00:01<00:00, 2189.52draws/s]
Auto-assigning NUTS sampler...
Initializing NUTS using jitter+adapt_diag...
Multiprocess sampling (2 chains in 2 jobs)
NUTS: [d]
Sampling 2 chains: 100%|██████████| 3000/3000 [00:01<00:00, 1902.75draws/s]
Auto-assigning NUTS sampler...
Initializing NUTS using jitter+adapt_diag...
Multiprocess sampling (2 chains in 2 jobs)
NUTS: [d]
Sampling 2 chains: 100%|██████████| 3000/3000 [00:01<00:00, 2094.20draws/s]
The acceptance probability does not match the target. It is 0.8795014352002515, but should be close to 0.8. Try to increase the number of tuning steps.
Auto-assigning NUTS sampler...
Initializing NUTS using jitter+adapt_diag...
Multiprocess sampling (2 chains in 2 jobs)
NUTS: [d]
Sampling 2 chains: 100%|██████████| 3000/3000 [00:01<00:00, 2022.38draws/s]
Auto-assigning NUTS sampler...
Initializing NUTS using jitter+adapt_diag...
Multiprocess sampling (2 chains in 2 jobs)
NUTS: [d]
Sampling 2 chains: 100%|██████████| 3000/3000 [00:01<00:00, 1739.55draws/s]
Auto-assigning NUTS sampler...
Initializing NUTS using jitter+adapt_diag...
Multiprocess sampling (2 chains in 2 jobs)
NUTS: [d]
Sampling 2 chains: 100%|██████████| 3000/3000 [00:01<00:00, 1690.90draws/s]
Auto-assigning NUTS sampler...
Initializing NUTS using jitter+adapt_diag...
Multiprocess sampling (2 chains in 2 jobs)
NUTS: [d]
Sampling 2 chains: 100%|██████████| 3000/3000 [00:01<00:00, 2152.30draws/s]
Auto-assigning NUTS sampler...
Initializing NUTS using jitter+adapt_diag...
Multiprocess sampling (2 chains in 2 jobs)
NUTS: [d]
Sampling 2 chains: 100%|██████████| 3000/3000 [00:01<00:00, 2048.63draws/s]
Auto-assigning NUTS sampler...
Initializing NUTS using jitter+adapt_diag...
Multiprocess sampling (2 chains in 2 jobs)
NUTS: [d]
Sampling 2 chains: 100%|██████████| 3000/3000 [00:01<00:00, 2101.32draws/s]
Auto-assigning NUTS sampler...
Initializing NUTS using jitter+adapt_diag...
Multiprocess sampling (2 chains in 2 jobs)
NUTS: [d]
Sampling 2 chains: 100%|██████████| 3000/3000 [00:01<00:00, 2177.93draws/s]
Auto-assigning NUTS sampler...
Initializing NUTS using jitter+adapt_diag...
Multiprocess sampling (2 chains in 2 jobs)
NUTS: [d]
Sampling 2 chains: 100%|██████████| 3000/3000 [00:01<00:00, 2176.02draws/s]
Auto-assigning NUTS sampler...
Initializing NUTS using jitter+adapt_diag...
Multiprocess sampling (2 chains in 2 jobs)
NUTS: [d]
Sampling 2 chains: 100%|██████████| 3000/3000 [00:01<00:00, 1831.72draws/s]
Auto-assigning NUTS sampler...
Initializing NUTS using jitter+adapt_diag...
Multiprocess sampling (2 chains in 2 jobs)
NUTS: [d]
Sampling 2 chains: 100%|██████████| 3000/3000 [00:01<00:00, 2147.45draws/s]
Auto-assigning NUTS sampler...
Initializing NUTS using jitter+adapt_diag...
Multiprocess sampling (2 chains in 2 jobs)
NUTS: [d]
Sampling 2 chains: 100%|██████████| 3000/3000 [00:01<00:00, 2120.89draws/s]
Auto-assigning NUTS sampler...
Initializing NUTS using jitter+adapt_diag...
Multiprocess sampling (2 chains in 2 jobs)
NUTS: [d]
Sampling 2 chains: 100%|██████████| 3000/3000 [00:01<00:00, 2160.66draws/s]
Auto-assigning NUTS sampler...
Initializing NUTS using jitter+adapt_diag...
Multiprocess sampling (2 chains in 2 jobs)
NUTS: [d]
Sampling 2 chains: 100%|██████████| 3000/3000 [00:01<00:00, 2162.56draws/s]
Auto-assigning NUTS sampler...
Initializing NUTS using jitter+adapt_diag...
Multiprocess sampling (2 chains in 2 jobs)
NUTS: [d]
Sampling 2 chains: 100%|██████████| 3000/3000 [00:01<00:00, 2242.22draws/s]
Auto-assigning NUTS sampler...
Initializing NUTS using jitter+adapt_diag...
Multiprocess sampling (2 chains in 2 jobs)
NUTS: [d]
Sampling 2 chains: 100%|██████████| 3000/3000 [00:01<00:00, 2225.90draws/s]
Auto-assigning NUTS sampler...
Initializing NUTS using jitter+adapt_diag...
Multiprocess sampling (2 chains in 2 jobs)
NUTS: [d]
Sampling 2 chains: 100%|██████████| 3000/3000 [00:01<00:00, 2249.19draws/s]
Auto-assigning NUTS sampler...
Initializing NUTS using jitter+adapt_diag...
Multiprocess sampling (2 chains in 2 jobs)
NUTS: [d]
Sampling 2 chains: 100%|██████████| 3000/3000 [00:01<00:00, 2178.26draws/s]
Auto-assigning NUTS sampler...
Initializing NUTS using jitter+adapt_diag...
Multiprocess sampling (2 chains in 2 jobs)
NUTS: [d]
Sampling 2 chains: 100%|██████████| 3000/3000 [00:01<00:00, 2267.80draws/s]
Auto-assigning NUTS sampler...
Initializing NUTS using jitter+adapt_diag...
Multiprocess sampling (2 chains in 2 jobs)
NUTS: [d]
Sampling 2 chains: 100%|██████████| 3000/3000 [00:01<00:00, 2149.20draws/s]






    













    













    













    













    













    













    













    













    













    













    













    













    













    













    













    













    













    













    













    













    













    













    













    













    













    













    













    













    













    













    













    













    













    













    













    













    













    













    













    













    













    













    













    













    

In [29]:

    

#cmap=sns.diverging_palette(150, 275, s=80, l=55, n=9, as_cmap=True)


    

In [30]:

    

np.arctan( Zsun / Rsun)*u.radian.to(u.deg)


    

    
Out[30]:





0.18638320362739255

In [ ]: