In [1]:

    

import pymc3 as pm
import numpy as np
import matplotlib.pyplot as plt
import corner
import pandas as pd
import wisps
import wisps.simulations as wispsim
from theano.compile.ops import as_op
import theano.tensor as tt
import scipy.integrate as integrate
import astropy.units as u
%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]:

    

lf=pd.read_pickle(wisps.OUTPUT_FILES+'/lf.pkl')


    

In [3]:

    

domega=4.1*(u.arcmin**2).to(u.radian**2)


    

In [4]:

    

Rsun=83000.
Zsun=27.

def density_function(r, z, h=300.):
    
    """
    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

def custom_volume(l,b,dmin, dmax, h):
    nsamp=1000
    ds = np.linspace(dmin,dmax,nsamp)
    rd=np.sqrt( (ds * np.cos( b ) )**2 + Rsun * (Rsun - 2 * ds * np.cos( b ) * np.cos( l ) ) )
    zd=Zsun+ ds * np.sin( b - np.arctan( Zsun / Rsun) )
    rh0=density_function(rd, zd,h=h )
    val=np.trapz(rh0*(ds**2), x=ds)
    return val

def easy_non_integral_volume(l, b, dmin, dmax, h):
    #this avoids integrals, to be used in likelhood methods like pymc
    nsamp=100
    ds = np.linspace(dmin,dmax,nsamp)
    rd=np.sqrt( (ds * np.cos( b ) )**2 + Rsun * (Rsun - 2 * ds * np.cos( b ) * np.cos( l ) ) )
    zd=Zsun+ ds * np.sin( b - np.arctan( Zsun / Rsun) )
    rh0=density_function(rd, zd,h=h )
    return np.nansum(rh0*ds**3)
    

def interpolated_lf(spt):
    sortedindex=np.argsort(lf['spt'])
    xs=lf['spt']
    ys=lf['phi']
    val=np.interp(spt,  xs[sortedindex], ys[sortedindex])
    return val

def bin_in_five_subtypes(sp_types, number):
    ranges=[[17, 20], [20, 25], [25, 30], [30, 35], [35, 42]]
    numbers=[]
    for r in ranges:
        idx= np.logical_and((r[0]<=sp_types), (r[1]>sp_types))
        numbers.append(np.nansum(number[idx]))
    return numbers

def compute_volume(h):
    vs=[]
    for s in wispsim.SPGRID:
        dmax, dmin=DLIMITS[s]
        vls=[]
        for l, b in lbs:
            vls.append(easy_non_integral_volume(l,b,dmin, dmax, h))
        vs.append(np.nansum(vls))
    return vs


    

In [5]:

    

interpolated_lf_values=interpolated_lf(wispsim.SPGRID)


    

In [6]:

    

pnts=wisps.OBSERVED_POINTINGS
lbs=[[pnt.coord.galactic.l.radian,pnt.coord.galactic.b.radian] for pnt in pnts]


    

In [7]:

    

DLIMITS=pnts[0].dist_limits


    

In [8]:

    

NOBS=np.array([148, 26, 3, 12, 2])


    

In [9]:

    

#@as_op(itypes=[tt.dscalar], otypes=[tt.dscalar])
def lnlike(h):
    vol=np.array(compute_volume(h))*domega
    nsim=np.array(bin_in_five_subtypes(wispsim.SPGRID, (vol*interpolated_lf_values)))
    #compute chi-square
    nobs_err=np.sqrt(NOBS)
    chi_square=np.nansum((nsim-NOBS)**2/(nobs_err**2))
    return -chi_square*0.5

def lnprior(h):
    if 10.0 < h < 1000.:
        return 0.0
    return -np.inf

def lnprob(theta):
    lp = lnprior(theta)
    if not np.isfinite(lp):
        return -np.inf
    return lp + lnlike(theta)


    

In [10]:

    

#forget pymc, doesn't work well with integrals
with pm.Model() as model:
    h=pm.Uniform('h', lower=100, upper=1000)
    like = pm.DensityDist('like', lnlike, observed={'h': h})
    trace = pm.sample(draws=100, cores=2, step=pm.Metropolis())


    

    




---------------------------------------------------------------------------
KeyboardInterrupt                         Traceback (most recent call last)
<ipython-input-10-2f61fa542c74> in <module>
      2 with pm.Model() as model:
      3     h=pm.Uniform('h', lower=100, upper=1000)
----> 4     like = pm.DensityDist('like', lnlike, observed={'h': h})
      5     trace = pm.sample(draws=100, cores=2, step=pm.Metropolis())

~/anaconda3/lib/python3.7/site-packages/pymc3/distributions/distribution.py in __new__(cls, name, *args, **kwargs)
     44             total_size = kwargs.pop('total_size', None)
     45             dist = cls.dist(*args, **kwargs)
---> 46             return model.Var(name, dist, data, total_size)
     47         else:
     48             raise TypeError("Name needs to be a string but got: {}".format(name))

~/anaconda3/lib/python3.7/site-packages/pymc3/model.py in Var(self, name, dist, data, total_size)
    843             with self:
    844                 var = MultiObservedRV(name=name, data=data, distribution=dist,
--> 845                                       total_size=total_size, model=self)
    846             self.observed_RVs.append(var)
    847             if var.missing_values:

~/anaconda3/lib/python3.7/site-packages/pymc3/model.py in __init__(self, name, data, distribution, total_size, model)
   1442         self.missing_values = [datum.missing_values for datum in self.data.values()
   1443                                if datum.missing_values is not None]
-> 1444         self.logp_elemwiset = distribution.logp(**self.data)
   1445         # The logp might need scaling in minibatches.
   1446         # This is done in `Factor`.

<ipython-input-9-fa5ba6237de9> in lnlike(h)
      1 #@as_op(itypes=[tt.dscalar], otypes=[tt.dscalar])
      2 def lnlike(h):
----> 3     vol=np.array(compute_volume(h))*domega
      4     nsim=np.array(bin_in_five_subtypes(wispsim.SPGRID, (vol*interpolated_lf_values)))
      5     #compute chi-square

<ipython-input-4-a338fa65093b> in compute_volume(h)
     53         vls=[]
     54         for l, b in lbs:
---> 55             vls.append(easy_non_integral_volume(l,b,dmin, dmax, h))
     56         vs.append(np.nansum(vls))
     57     return vs

<ipython-input-4-a338fa65093b> in easy_non_integral_volume(l, b, dmin, dmax, h)
     28     rd=np.sqrt( (ds * np.cos( b ) )**2 + Rsun * (Rsun - 2 * ds * np.cos( b ) * np.cos( l ) ) )
     29     zd=Zsun+ ds * np.sin( b - np.arctan( Zsun / Rsun) )
---> 30     rh0=density_function(rd, zd,h=h )
     31     return np.nansum(rh0*ds**3)
     32 

<ipython-input-4-a338fa65093b> in density_function(r, z, h)
     11     zpart=np.exp(-abs(z-Zsun)/h)
     12     rpart=np.exp(-(r-Rsun)/l)
---> 13     return zpart*rpart
     14 
     15 def custom_volume(l,b,dmin, dmax, h):

~/anaconda3/lib/python3.7/site-packages/theano/tensor/var.py in __mul__(self, other)
    153         # and the return value in that case
    154         try:
--> 155             return theano.tensor.mul(self, other)
    156         except (NotImplementedError, AsTensorError):
    157             return NotImplemented

~/anaconda3/lib/python3.7/site-packages/theano/gof/op.py in __call__(self, *inputs, **kwargs)
    668                 # compute output value once with test inputs to validate graph
    669                 thunk = node.op.make_thunk(node, storage_map, compute_map,
--> 670                                            no_recycling=[])
    671                 thunk.inputs = [storage_map[v] for v in node.inputs]
    672                 thunk.outputs = [storage_map[v] for v in node.outputs]

~/anaconda3/lib/python3.7/site-packages/theano/gof/op.py in make_thunk(self, node, storage_map, compute_map, no_recycling, impl)
    953             try:
    954                 return self.make_c_thunk(node, storage_map, compute_map,
--> 955                                          no_recycling)
    956             except (NotImplementedError, utils.MethodNotDefined):
    957                 # We requested the c code, so don't catch the error.

~/anaconda3/lib/python3.7/site-packages/theano/gof/op.py in make_c_thunk(self, node, storage_map, compute_map, no_recycling)
    856         _logger.debug('Trying CLinker.make_thunk')
    857         outputs = cl.make_thunk(input_storage=node_input_storage,
--> 858                                 output_storage=node_output_storage)
    859         thunk, node_input_filters, node_output_filters = outputs
    860 

~/anaconda3/lib/python3.7/site-packages/theano/gof/cc.py in make_thunk(self, input_storage, output_storage, storage_map, keep_lock)
   1215         cthunk, module, in_storage, out_storage, error_storage = self.__compile__(
   1216             input_storage, output_storage, storage_map,
-> 1217             keep_lock=keep_lock)
   1218 
   1219         res = _CThunk(cthunk, init_tasks, tasks, error_storage, module)

~/anaconda3/lib/python3.7/site-packages/theano/gof/cc.py in __compile__(self, input_storage, output_storage, storage_map, keep_lock)
   1155                                             output_storage,
   1156                                             storage_map,
-> 1157                                             keep_lock=keep_lock)
   1158         return (thunk,
   1159                 module,

~/anaconda3/lib/python3.7/site-packages/theano/gof/cc.py in cthunk_factory(self, error_storage, in_storage, out_storage, storage_map, keep_lock)
   1622                 node.op.prepare_node(node, storage_map, None, 'c')
   1623             module = get_module_cache().module_from_key(
-> 1624                 key=key, lnk=self, keep_lock=keep_lock)
   1625 
   1626         vars = self.inputs + self.outputs + self.orphans

~/anaconda3/lib/python3.7/site-packages/theano/gof/cmodule.py in module_from_key(self, key, lnk, keep_lock)
   1153         # Is the source code already in the cache?
   1154         module_hash = get_module_hash(src_code, key)
-> 1155         module = self._get_from_hash(module_hash, key, keep_lock=keep_lock)
   1156         if module is not None:
   1157             return module

~/anaconda3/lib/python3.7/site-packages/theano/gof/cmodule.py in _get_from_hash(self, module_hash, key, keep_lock)
   1063                 if (key[0] and not key_broken and
   1064                         self.check_for_broken_eq):
-> 1065                     self.check_key(key, key_data.key_pkl)
   1066             self._update_mappings(key, key_data, module.__file__, check_in_keys=not key_broken)
   1067             return module

~/anaconda3/lib/python3.7/site-packages/theano/gof/cmodule.py in check_key(self, key, key_pkl)
   1267         # part of the key is not broken.
   1268         for other in self.similar_keys.get(get_safe_part(key), []):
-> 1269             if other is not key and other == key and hash(other) != hash(key):
   1270                 raise AssertionError(
   1271                     "Found two keys that are equal but have a different hash. "

~/anaconda3/lib/python3.7/site-packages/theano/gof/utils.py in __eq__(self, other)
    196                     return (type(self) == type(other) and
    197                             tuple(getattr(self, a) for a in props) ==
--> 198                             tuple(getattr(other, a) for a in props))
    199                 dct['__eq__'] = __eq__
    200 

KeyboardInterrupt: 

In [12]:

    

ndim, nwalkers = 1, 100
pos = [ 1000*np.random.randn(ndim) for i in range(nwalkers)]
import emcee
sampler = emcee.EnsembleSampler(nwalkers, ndim, lnprob)


    

In [ ]:

    

#ugh very slow probably won't converge


    

In [ ]:

    

samples=sampler.run_mcmc(pos, 1000)


    

In [23]:

    

sampl = sampler.chain.flatten()


    

In [25]:

    

plt.hist(sampl)


    

    
Out[25]:





(array([  600.,   800.,  1054.,  1928.,   219., 14472.,   338.,    75.,
          387.,   127.]),
 array([-2659.69274494, -2148.18375678, -1636.67476862, -1125.16578046,
         -613.65679231,  -102.14780415,   409.36118401,   920.87017217,
         1432.37916032,  1943.88814848,  2455.39713664]),
 <a list of 10 Patch objects>)






    

In [ ]: