In [82]:

    

%pylab inline
%config InlineBackend.figure_format = 'retina'
import exocartographer.gp_emission_map as gem
import exocartographer.util as u
import healpy
import IPython.display as disp
import matplotlib.animation as anim
from matplotlib.colors import LogNorm
import scipy.optimize as so
import scipy.signal as ss


    

    




Populating the interactive namespace from numpy and matplotlib






    




WARNING: pylab import has clobbered these variables: ['disp']
`%matplotlib` prevents importing * from pylab and numpy

In [7]:

    

lc1 = loadtxt('epic211098454.txt')
errorbar(lc1[:,0], lc1[:,1], lc1[:,2])


    

    
Out[7]:





<Container object of 3 artists>






    

In [88]:

    

fs, psd = ss.welch(lc1[:,1], fs=1.0/diff(lc1[:,0])[0], nperseg=2048)
loglog(fs, psd)


    

    
Out[88]:





[<matplotlib.lines.Line2D at 0x1155d3f90>]






    

In [175]:

    

P0 = 1.129 # From README


    

In [189]:

    

reload(gem)


    

    
Out[189]:





<module 'exocartographer.gp_emission_map' from '/Users/farr/Documents/code/exocartographer/exocartographer/gp_emission_map.py'>

In [190]:

    

mappost = gem.EmissionMapPosterior(lc1[:,0], log(lc1[:,1]), lc1[:,1]/lc1[:,2])


    

In [191]:

    

gppbest = so.fmin_powell(lambda x: -mappost.gp_marginalised_posterior(x), zeros(6))


    

    




---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
<ipython-input-191-86c9a289e5c8> in <module>()
----> 1 gppbest = so.fmin_powell(lambda x: -mappost.gp_marginalised_posterior(x), zeros(6))

/Users/farr/Library/Enthought/Canopy_64bit/User/lib/python2.7/site-packages/scipy/optimize/optimize.pyc in fmin_powell(func, x0, args, xtol, ftol, maxiter, maxfun, full_output, disp, retall, callback, direc)
   2280             'return_all': retall}
   2281 
-> 2282     res = _minimize_powell(func, x0, args, callback=callback, **opts)
   2283 
   2284     if full_output:

/Users/farr/Library/Enthought/Canopy_64bit/User/lib/python2.7/site-packages/scipy/optimize/optimize.pyc in _minimize_powell(func, x0, args, callback, xtol, ftol, maxiter, maxfev, disp, direc, return_all, **unknown_options)
   2339         direc = asarray(direc, dtype=float)
   2340 
-> 2341     fval = squeeze(func(x))
   2342     x1 = x.copy()
   2343     iter = 0

/Users/farr/Library/Enthought/Canopy_64bit/User/lib/python2.7/site-packages/scipy/optimize/optimize.pyc in function_wrapper(*wrapper_args)
    287     def function_wrapper(*wrapper_args):
    288         ncalls[0] += 1
--> 289         return function(*(wrapper_args + args))
    290 
    291     return ncalls, function_wrapper

<ipython-input-191-86c9a289e5c8> in <lambda>(x)
----> 1 gppbest = so.fmin_powell(lambda x: -mappost.gp_marginalised_posterior(x), zeros(6))

/Users/farr/Documents/code/exocartographer/exocartographer/gp_emission_map.py in gp_marginalised_posterior(self, p)
    187         SVsV = np.dot(SVs, V)
    188 
--> 189         SVsVI = SVsV + np.eye(self.intensity.shape[0])
    190         ldSVsVI = np.linalg.slogdet(SVsVI)[1]
    191 

ValueError: operands could not be broadcast together with shapes (192,192) (3204,3204) 

In [177]:

    

p0 = mappost.to_params(zeros(mappost.nparams))
p0['logit_wn_rel_amp'] = u.logit(0.1, mappost.wn_low, mappost.wn_high)
p0['logit_spatial_scale'] = u.logit(pi/4.0, mappost.spatial_scale_low, mappost.spatial_scale_high)
p0['log_period'] = log(P0)
p0['logit_cos_theta'] = u.logit(0.5)
p0['log_intensity_map'] = randn(mappost.npix) # initial guess
p0['log_sigma'] = log(std(log(lc1[:,1]))/std(mappost.intensity_series(p0)))
p0['log_intensity_map'] = exp(p0['log_sigma'])*randn(mappost.npix)
p0['mu'] = mean(p0['log_intensity_map']) + mean(log(lc1[:,1])) - mean(mappost.intensity_series(p0))
p0['log_intensity_map'] += p0['mu']


    

In [178]:

    

def all_but_period_posterior(x):
    xx = insert(x, 4, P0)
    return mappost(xx)


    

In [179]:

    

plot(lc1[:,0], log(lc1[:,1]), '-k')
plot(lc1[:,0], mappost.intensity_series(p0))


    

    
Out[179]:





[<matplotlib.lines.Line2D at 0x121b02190>]






    

In [180]:

    

pp = p0.reshape((1,)).view(float)
cb_best = concatenate((pp[:4], pp[5:]))
def cb(x):
    global cb_best
    cb_best = x
    p = mappost.to_params(insert(x, 4, P0))
    m = exp(p['log_intensity_map'])
    disp.clear_output(wait=True)
    figure(1)
    clf()
    healpy.mollview(m, cmap='viridis', fig=1, norm=LogNorm())
    disp.display(gcf())
    figure(2)
    clf()
    plot(lc1[:,0], log(lc1[:,1]), '-k')
    plot(lc1[:,0], mappost.intensity_series(p), '-b')
    disp.display(gcf())


    

In [182]:

    

pbest = so.fmin_powell(lambda x: -all_but_period_posterior(x), cb_best, callback=cb)


    

    













    













    




---------------------------------------------------------------------------
KeyboardInterrupt                         Traceback (most recent call last)
<ipython-input-182-33a2cdd3f8bf> in <module>()
----> 1 pbest = so.fmin_powell(lambda x: -all_but_period_posterior(x), cb_best, callback=cb)

/Users/farr/Library/Enthought/Canopy_64bit/User/lib/python2.7/site-packages/scipy/optimize/optimize.pyc in fmin_powell(func, x0, args, xtol, ftol, maxiter, maxfun, full_output, disp, retall, callback, direc)
   2280             'return_all': retall}
   2281 
-> 2282     res = _minimize_powell(func, x0, args, callback=callback, **opts)
   2283 
   2284     if full_output:

/Users/farr/Library/Enthought/Canopy_64bit/User/lib/python2.7/site-packages/scipy/optimize/optimize.pyc in _minimize_powell(func, x0, args, callback, xtol, ftol, maxiter, maxfev, disp, direc, return_all, **unknown_options)
   2351             fx2 = fval
   2352             fval, x, direc1 = _linesearch_powell(func, x, direc1,
-> 2353                                                  tol=xtol * 100)
   2354             if (fx2 - fval) > delta:
   2355                 delta = fx2 - fval

/Users/farr/Library/Enthought/Canopy_64bit/User/lib/python2.7/site-packages/scipy/optimize/optimize.pyc in _linesearch_powell(func, p, xi, tol)
   2173     def myfunc(alpha):
   2174         return func(p + alpha*xi)
-> 2175     alpha_min, fret, iter, num = brent(myfunc, full_output=1, tol=tol)
   2176     xi = alpha_min*xi
   2177     return squeeze(fret), p + xi, xi

/Users/farr/Library/Enthought/Canopy_64bit/User/lib/python2.7/site-packages/scipy/optimize/optimize.pyc in brent(func, args, brack, tol, full_output, maxiter)
   1912     options = {'xtol': tol,
   1913                'maxiter': maxiter}
-> 1914     res = _minimize_scalar_brent(func, brack, args, **options)
   1915     if full_output:
   1916         return res['x'], res['fun'], res['nit'], res['nfev']

/Users/farr/Library/Enthought/Canopy_64bit/User/lib/python2.7/site-packages/scipy/optimize/optimize.pyc in _minimize_scalar_brent(func, brack, args, xtol, maxiter, **unknown_options)
   1944                   full_output=True, maxiter=maxiter)
   1945     brent.set_bracket(brack)
-> 1946     brent.optimize()
   1947     x, fval, nit, nfev = brent.get_result(full_output=True)
   1948     return OptimizeResult(fun=fval, x=x, nit=nit, nfev=nfev,

/Users/farr/Library/Enthought/Canopy_64bit/User/lib/python2.7/site-packages/scipy/optimize/optimize.pyc in optimize(self)
   1750         # set up for optimization
   1751         func = self.func
-> 1752         xa, xb, xc, fa, fb, fc, funcalls = self.get_bracket_info()
   1753         _mintol = self._mintol
   1754         _cg = self._cg

/Users/farr/Library/Enthought/Canopy_64bit/User/lib/python2.7/site-packages/scipy/optimize/optimize.pyc in get_bracket_info(self)
   1724         ### carefully DOCUMENT any CHANGES in core ##
   1725         if brack is None:
-> 1726             xa, xb, xc, fa, fb, fc, funcalls = bracket(func, args=args)
   1727         elif len(brack) == 2:
   1728             xa, xb, xc, fa, fb, fc, funcalls = bracket(func, xa=brack[0],

/Users/farr/Library/Enthought/Canopy_64bit/User/lib/python2.7/site-packages/scipy/optimize/optimize.pyc in bracket(func, xa, xb, args, grow_limit, maxiter)
   2099     _gold = 1.618034
   2100     _verysmall_num = 1e-21
-> 2101     fa = func(*(xa,) + args)
   2102     fb = func(*(xb,) + args)
   2103     if (fa < fb):                      # Switch so fa > fb

/Users/farr/Library/Enthought/Canopy_64bit/User/lib/python2.7/site-packages/scipy/optimize/optimize.pyc in myfunc(alpha)
   2172     """
   2173     def myfunc(alpha):
-> 2174         return func(p + alpha*xi)
   2175     alpha_min, fret, iter, num = brent(myfunc, full_output=1, tol=tol)
   2176     xi = alpha_min*xi

/Users/farr/Library/Enthought/Canopy_64bit/User/lib/python2.7/site-packages/scipy/optimize/optimize.pyc in function_wrapper(*wrapper_args)
    287     def function_wrapper(*wrapper_args):
    288         ncalls[0] += 1
--> 289         return function(*(wrapper_args + args))
    290 
    291     return ncalls, function_wrapper

<ipython-input-182-33a2cdd3f8bf> in <lambda>(x)
----> 1 pbest = so.fmin_powell(lambda x: -all_but_period_posterior(x), cb_best, callback=cb)

<ipython-input-178-b03173b49712> in all_but_period_posterior(x)
      1 def all_but_period_posterior(x):
      2     xx = insert(x, 4, P0)
----> 3     return mappost(xx)

/Users/farr/Documents/code/exocartographer/exocartographer/gp_emission_map.py in __call__(self, p)

/Users/farr/Documents/code/exocartographer/exocartographer/gp_emission_map.py in logpdata(self, p)
    151         p = self.to_params(p)
    152 
--> 153         lp = 0.0
    154 
    155         lp += flat_logit_log_prior(p['logit_wn_rel_amp'],

/Users/farr/Documents/code/exocartographer/exocartographer/gp_emission_map.py in intensity_series(self, p)
    137         sigma = self.sigma(p)
    138         wn_rel_amp = inv_logit(p['logit_wn_rel_amp'], low=self.wn_low, high=self.wn_high)
--> 139         sp_scale = self.spatial_scale(p)
    140 
    141         return gm.map_logprior(p['log_intensity_map'], p['mu'], sigma, wn_rel_amp, sp_scale)

/Users/farr/Documents/code/exocartographer/exocartographer/gp_emission_map.py in spatial_intensity_series(self, p)
    134     def logmapprior(self, p):
    135         p = self.to_params(p)
--> 136 
    137         sigma = self.sigma(p)
    138         wn_rel_amp = inv_logit(p['logit_wn_rel_amp'], low=self.wn_low, high=self.wn_high)

KeyboardInterrupt: 





    













    

In [183]:

    

reload(gem)


    

    
Out[183]:





<module 'exocartographer.gp_emission_map' from '/Users/farr/Documents/code/exocartographer/exocartographer/gp_emission_map.py'>

In [ ]: