First, we are implementng a powerlaw parameterization for planet radius and perido/insolation using the PyStan No Uturn MCMC sampler, in order to calculate occurrence rates based off of work by Chris Burke and the open science blog post by Dan Foreman-Mackey: The original exopop.ipynb from Dan Foreman-Mackey, is found at http://dan.iel.fm/downloads/notebooks/exopop.ipynb, and the accompanying blog post is at http://dan.iel.fm/posts/exopop/ Some parts of this code were developed at the exoplanet population hackathon by Ravi and Joe Catanzarite.

8/15/2016 This PyStan model version uses pre-computed completeness as input, for testing and developement.


In [227]:
#%%=========================================================
# Intialize
# Had to put __future__ before other imports
from __future__ import division, print_function
from matplotlib import rcParams
rcParams["savefig.dpi"] = 100
rcParams["font.size"] = 20
import os
import requests
import pandas as pd
from cStringIO import StringIO

# The following commands are for use in the notebook environment
# Comment them out if this code is to be run in a python or ipython shell
%matplotlib inline
%config InlineBackend.figure_format = "retina"

In [228]:
#%%=========================================================
# This function downloads and caches a dataset from the exoplanet archive
# into a pandas frame
def get_catalog(name, basepath="data"):
    fn = os.path.join(basepath, "{0}.h5".format(name))
    if os.path.exists(fn):
        return pd.read_hdf(fn, name)
    if not os.path.exists(basepath):
        os.makedirs(basepath)
    print("Downloading {0}...".format(name))
    url = ("http://exoplanetarchive.ipac.caltech.edu/cgi-bin/nstedAPI/"
           "nph-nstedAPI?table={0}&select=*").format(name)
    r = requests.get(url)
    if r.status_code != requests.codes.ok:
        r.raise_for_status()
    fh = StringIO(r.content)
    df = pd.read_csv(fh)
    df.to_hdf(fn, name, format="t")
    return df

In [229]:
#%%=========================================================
# Get the stellar catalog and make selection cuts   
import numpy as np

#!!!!! Q1-Q17 DR 24 (9.2 pipeline)
stlr = get_catalog("q1_q17_dr24_stellar")

#!!!!! Q1-Q16 (9.1 pipeline)
# stlr = get_catalog("q1_q16_stellar")

# !!!!! Select G dwarfs.
# m = (5300 <= stlr.teff) & (stlr.teff <= 6000)

# !!!!! Select K dwarfs.
m = (3900 <= stlr.teff) & (stlr.teff <= 5300)

# !!!!! Select M dwarfs.
# m = (2400 <= stlr.teff) & (stlr.teff <= 3900)

# stellar radius cut
#m &= stlr.radius <= 1.15

# Only include stars with sufficient data coverage:

# Minimum dataspan of 2 years
m &= stlr.dataspan > 365.25*2.

# Minimum dutycycle of 1/2
m &= stlr.dutycycle > 0.50

# Minimum data coverage of 2 years
m &= stlr.dutycycle*stlr.dataspan > 365.25*2.

# minimum rms cdpp at 7.5 hour pulse
#m &= stlr.rrmscdpp07p5 <= 1000.

# Only select stars with mass estimates.
m &= np.isfinite(stlr.mass)

# put selected data into pandas data frame
base_stlr = pd.DataFrame(stlr)
stlr = pd.DataFrame(stlr[m])

print("Selected {0} targets after cuts".format(len(stlr)))


Selected 40685 targets after cuts

In [230]:
#%%=========================================================
# Plot an HR diagram of the selected targets
import matplotlib.pyplot as pl

pl.plot(base_stlr.teff, base_stlr.logg, ".k", ms=3, alpha=0.5)
pl.plot(stlr.teff, stlr.logg, ".r", ms=3, alpha=0.5)
pl.xlim(9500, 2500)
pl.ylim(5.5, 3)
pl.ylabel("$\log g$");
pl.xlabel("$T_\mathrm{eff}$");



In [231]:
#%%=========================================================

# !!!!! Get the planet catalog and make selection cuts 

# !!!!! Q1-Q17 DR24 (9.2 pipeline)
kois = get_catalog("q1_q17_dr24_koi")

# !!!!! Q1-Q16 planet catalog (9.1 pipeline)
# kois = get_catalog("q1_q16_koi")

# Set insolation and planet radius ranges
# period_rng = (20, 320)
rp_rng = (0.75, 2.5)
insolation_rng = (0.2, 20)

# Join on the stellar list.
kois = pd.merge(kois, stlr[["kepid"]], on="kepid", how="inner")

# Only select the KOIs in the relevant part of parameter space.
m = kois.koi_pdisposition == "CANDIDATE"
base_kois = pd.DataFrame(kois[m])

# Select insolation range instead of period range
m &= (insolation_rng[0] <= kois.koi_insol) & (kois.koi_insol <= insolation_rng[1])
# m &= (period_rng[0] <= kois.koi_period) & (kois.koi_period <= period_rng[1])

m &= np.isfinite(kois.koi_prad) & (rp_rng[0] <= kois.koi_prad) & (kois.koi_prad <= rp_rng[1])

# !!!!! Only include PCs with MES > 15 for 9.2 (Q1-Q17 Dr24)
m &= kois.koi_max_mult_ev > 15
# !!!!! Comment out above statement for 9.1 (Q1-Q16) -- Note all max_mult_ev seem to be NaNs anyyway

# Panda data frame for selected kois
# Note that kois now contains only the selected planet candidates
kois = pd.DataFrame(kois[m])

#print("min insolation = {0} ".format(np.min((kois.koi_insol))))
#print("max insolation = {0} ".format(np.max((kois.koi_insol))))
#print(kois.koi_insol)


#print("min period = {0} ".format(np.min(kois.koi_period)))
#print("max period = {0} ".format(np.max(kois.koi_period)))
#print("Selected {0} KOIs after cuts".format(len(kois)))

In [232]:
#%%=========================================================
# Plot the measurements with error bars, in insolation-radius parameter space
yerr = np.abs(np.array(base_kois[["koi_prad_err2", "koi_prad_err1"]])).T
pl.errorbar(base_kois.koi_insol, base_kois.koi_prad, yerr=yerr, fmt=".k", ms=4,
            capsize=0, alpha=0.3)
pl.plot(kois.koi_insol, kois.koi_prad, ".k", ms=6)
pl.fill_between(insolation_rng, [rp_rng[1], rp_rng[1]], [rp_rng[0], rp_rng[0]], color="r", alpha=0.2)
pl.xlim(insolation_rng + 0.1 * np.array([-1, 1]))
pl.ylim(rp_rng + 0.5 * np.array([-1, 1]))
pl.xlabel("insolation [Earth units]")
pl.ylabel("$R_p \, [R_\oplus]$");



In [233]:
#%%=========================================================
# Completeness model helper functions (radius and orbital period)

from scipy.stats import gamma

def get_duration(period, aor, e):
    """
    Equation (1) from Burke et al. This estimates the transit
    duration in the same units as the input period. There is a
    typo in the paper (24/4 = 6 != 4).
    
    :param period: the period in any units of your choosing
    :param aor:    the dimensionless semi-major axis (scaled
                   by the stellar radius)
    :param e:      the eccentricity of the orbit
    
    """
    return 0.25 * period * np.sqrt(1 - e**2) / aor

def get_a(period, mstar, Go4pi=2945.4625385377644/(4*np.pi*np.pi)):
    """
    Compute the semi-major axis of an orbit in Solar radii.
    
    :param period: the period in days
    :param mstar:  the stellar mass in Solar masses
    
    """
    return (Go4pi*period*period*mstar) ** (1./3)

def get_delta(k, c=1.0874, s=1.0187):
    """
    Estimate the approximate expected transit depth as a function
    of radius ratio. There might be a typo here. In the paper it
    uses c + s*k but in the public code, it is c - s*k:
    https://github.com/christopherburke/KeplerPORTs
    
    :param k: the dimensionless radius ratio between the planet and
              the star
    
    """
    delta_max = k*k * (c + s*k)
    
    # !!!!! For Q1-Q16 (9.1 pipeline) DFM used a multiplier 0f 0.84 instead of 1.0 in the equation below
    # return 0.84*delta_max
    
    # !!!!! for Q1-Q17 DR24 (9.2 pipeline), use a multiplier of 1 ????? Ask Chris
    return 1.0* delta_max

# 14 pulse durations
cdpp_cols = [k for k in stlr.keys() if k.startswith("rrmscdpp")]
cdpp_vals = np.array([k[-4:].replace("p", ".") for k in cdpp_cols], dtype=float)

def get_mes(star, period, rp, tau, re=0.009171):
    """
    Estimate the multiple event statistic value for a transit.
    
    :param star:   a pandas row giving the stellar properties
    :param period: the period in days
    :param rp:     the planet radius in Earth radii
    :param tau:    the transit duration in hours
    
    """
    # Interpolate the RMS CDPP corresponding to the transit duration.
    cdpp = np.array(star[cdpp_cols], dtype=float)
    sigma = np.interp(tau, cdpp_vals, cdpp)

    # Compute the radius ratio and estimate the S/N.
    k = rp * re / star.radius
    snr = get_delta(k) * 1e6 / sigma
    
    # Scale by the estimated number of transits.
    ntrn = star.dataspan * star.dutycycle / period 
    return snr * np.sqrt(ntrn)

# Pre-compute and freeze the gamma function from Equation (5) in
# Burke et al.

# !!!!! Q1-Q16 (9.1 pipeline): DFM used the parameters below
# pgam = gamma(4.65, loc=0., scale=0.98)
# !!!!! But the parameters for FGK are (4.35, 0, 1.05) according to astro-ph 1507.05097 (Christiansen)
# pgam = gamma(4.35, loc=0., scale=1.05)

# !!!!! Parameters for 9.2 pipeline Q1-Q17 DR24 are from Jessie Christiansen
# !!!!! Note that these parameters do not apply to M stars!
pgam = gamma(103.0113, loc=0., scale=0.10583)

# mesthres_cols are column names for the 14 pulse durations
mesthres_cols = [k for k in stlr.keys() if k.startswith("mesthres")]

# pulse_durations_obs are the 14 pulse durations
pulse_durations_obs = np.array([k[-4:].replace("p", ".") for k in mesthres_cols],
                         dtype=float)
def get_pdet(star, aor, period, rp, e):
    """
    Equation (5) from Burke et al. Estimate the detection efficiency
    for a transit.
    
    :param star:   a pandas row giving the stellar properties
    :param aor:    the dimensionless semi-major axis (scaled
                   by the stellar radius)
    :param period: the period in days
    :param rp:     the planet radius in Earth radii
    :param e:      the orbital eccentricity
    
    """
    # mest is the interpolated MES threshold corresponding to the transit duration
    # tau is the pulse duration
    # pulse_durations_obs are the 14 pulse durations
    # np.array(star[mesthres_cols],dtype=float) are the coresponding MES thresholds
    tau = get_duration(period, aor, e) * 24.
    mes = get_mes(star, period, rp, tau)
    mest = np.interp(tau, pulse_durations_obs,
                     np.array(star[mesthres_cols],
                              dtype=float))
    x = mes - 4.1 - (mest - 7.1)
    
    # !!!!! DFM originally used no multiplier in the equation below;
    # The multiplier of 0.78442 must be a 'plateau factor' from Jessie, for 9.2
    # for the 9.2 Q1-Q17 DR24 data, provided by Jessie
    return 0.78442*pgam.cdf(x)

    # !!!!! For the 9.1 data
    # return pgam.cdf(x)

def get_pwin(star, period):
    """
    Equation (6) from Burke et al. Estimates the window function
    using a binomial distribution.
    
    :param star:   a pandas row giving the stellar properties
    :param period: the period in days
    
    """
    M = star.dataspan / period
    f = star.dutycycle
    omf = 1.0 - f
    pw = 1 - omf**M - M*f*omf**(M-1) - 0.5*M*(M-1)*f*f*omf**(M-2)
    msk = (pw >= 0.0) * (M >= 2.0)
    return pw * msk

def get_pgeom(aor, e):
    """
    The geometric transit probability.
    
    See e.g. Kipping (2014) for the eccentricity factor
    http://arxiv.org/abs/1408.1393
    
    :param aor: the dimensionless semi-major axis (scaled
                by the stellar radius)
    :param e:   the orbital eccentricity

    """
    return 1. / (aor * (1 - e*e)) * (aor > 1.0)

In [234]:
#================================================================================================
#%% This cell contains functions and code for calculating completeness in the
#       parameter space of [ insolation , planet radius ] 
    
def get_completeness(star, period, rp, e, with_geom=True):
    """
    A helper function to combine all the completeness effects.
    
    :param star:      a pandas row giving the stellar properties (each row is for one star)
    :param period:    the period in days 
    :param rp:        the planet radius in Earth radii
    :param e:         the orbital eccentricity
    :param with_geom: include the geometric transit probability?
    
    """
    aor = get_a(period, star.mass) / star.radius
    pdet = get_pdet(star, aor, period, rp, e)
    pwin = get_pwin(star, period)
    if not with_geom:
        return pdet * pwin
    pgeom = get_pgeom(aor, e)
    return pdet * pwin * pgeom
    print(len(rp))
# Construct grid for planet radius
rp2 = np.linspace(rp_rng[0], rp_rng[1], 61)

# Construct grid for insolation
insolation = np.linspace(insolation_rng[0], insolation_rng[1], 57)
#insolation_grid, rp_grid2 = np.meshgrid(insolation, rp2, indexing="xy")
insolation_grid, rp_grid2 = np.meshgrid(insolation, rp2, indexing="ij")


    
#def get_completeness_from_insolation(star,  rp_grid2, insolation_grid, e, with_geom=True):
def get_completeness_from_insolation(star, insolation_grid, rp_grid2, e, with_geom=True):
    
    # compute the periods corresponding to an insolation grid
    insolation = insolation_grid
    period_grid = get_period_from_insolation( star , insolation )
    
    # completeness 
    completeness = get_completeness(star, period_grid, rp_grid2, e, with_geom=True)
    
    return completeness
    

# Add a function to compute insolation on the period grid, for a given star 
def get_insolation_from_period( star , period ):
    
    # Get needed stellar parameters
    teffStar = star.teff
    teffSun = 5777
    rStar = star.radius
    mStar = star.mass
    
    # Semimajor axis of planet in AU
    aPlanet = mStar**(1.0/3.0) * (period/365.25)**(2.0/3.0) 
    
    # Compute insolation
    insolation = (teffStar/teffSun)**4*(rStar/1)**2*(1/aPlanet)**2
    
    return insolation
    
def get_period_from_insolation( star , insolation ):
    
    # Get needed stellar parameters
    teffStar = star.teff
    teffSun = 5777
    rStar = star.radius
    mStar = star.mass
    
    # Get semimajor axis from star properties and insolation, using
    # insolation = ( teffStar / teffSun )**4 * ( rStar / 1)**2 * ( 1 / aPlanet )**2
    aPlanet = ( ( teffStar / teffSun )**4 * ( rStar / 1)**2 / insolation )**(0.5)
    
    # Get orbit period in days from semimajor axis of planet in AU and start properties, using
    # aPlanet = mStar**(1.0/3.0) * (period/365.25)**(2.0/3.0)
    period = 365.25 * ( aPlanet/( mStar**(1.0/3.0) ) )**(3.0/2.0)
    
    return period

In [235]:
#%%==============================================================================================
# Test: compute completeness grid in [ insolation , planet radius ] 
# parameter space for the first star in the catalog
# new_completeness_grid_single_star = get_completeness_from_insolation(stlr.iloc[0], insolation_grid, rp_grid2, 0.0, with_geom=True)

# Marginalize detection contours over all selected targets
# including the geometric factor. This takes a few minutes.
new_completeness = np.zeros_like(insolation_grid)
for _, star in stlr.iterrows():
    new_completeness += get_completeness_from_insolation(star, insolation_grid, rp_grid2, 0.0, with_geom=True)
    #new_completeness += get_completeness_from_insolation(star, rp_grid2, insolation_grid, 0.0, with_geom=True)

In [236]:
#%%=========================================================
# Reproducing Figure 1 from Burke paper, which is the
# completeness contour (not including geometric effect) for an example target
 
# Repeating for new insolation completeness.  
    
# Choose the first star
# star = stlr[stlr.kepid == 10593626].iloc[0]
# First M star
#star = stlr[stlr.kepid == 10031707].iloc[0]
star = stlr[stlr.kepid == stlr.kepid.iloc[0]].iloc[0]

# Compute the completeness map on a grid.
X, Y = insolation_grid, rp_grid2, 
Z = get_completeness_from_insolation(star, X, Y, 0.0, with_geom=False)

# Plot with the same contour levels as the figure. Add some contours at low end of completeness.
c = pl.contour(X, Y, Z, [0.0005, 0.001, 0.005, 0.01, 0.05, 0.1, 0.25, 0.5, 0.75, 0.9, 0.95, 0.99], colors="k")
pl.clabel(c, fontsize=12, inline=1, fmt="%.4f") 
pl.xlabel("insolation [Earth units]")
pl.ylabel("$R_p \, [R_\oplus]$")
pl.title("det. eff. for KIC {0}".format(np.min(stlr.kepid.iloc[0])));



In [237]:
#%%=========================================================
# Plot the average new_completeness contour (radius-insolation)
# Include the geometric effect
pl.pcolor(insolation_grid, rp_grid2, new_completeness, cmap="Blues")
c = pl.contour(insolation_grid, rp_grid2, new_completeness / len(stlr),
               colors="k", alpha=0.8)
pl.clabel(c, fontsize=12, inline=1, fmt="%.3f")
pl.title("mean pipeline detection efficiency")
pl.xlabel("insolation [Earth units]")
pl.ylabel("$R_p \, [R_\oplus]$");



In [238]:
#%%=========================================================
# Population inference with an independent power law model
# Using modified code above that computes 
# completeness in the parameter space of [ insolation , planet radius ], 

# A double power law model for the population.
def population_model_insolation(theta, insolation, rp):
    # Parameters 
    # lnf0 is normalization, 
    # beta is exponent of insolation power law,
    # alpha is exponent of radius power law
    lnf0, beta, alpha = theta
    v = np.exp(lnf0) * np.ones_like(insolation)
    for x, rng, n in zip((insolation, rp),
                         (insolation_rng, rp_rng),
                         (beta, alpha)):
        n1 = n + 1
        v *= x**n*n1 / (rng[1]**n1-rng[0]**n1)
    return v

# The ln-likelihood function given at the top of this post.
# change to insolation from planet catalog

# Insolation and radius for planets in catalog
# koi_periods = np.array(kois.koi_period)
koi_insolation = np.array(kois.koi_insol)
koi_rps = np.array(kois.koi_prad)

# Parameter space volume in each bin of [insolation, radius] grid
# Note the bins are not uniformly spaced in insolation
vol = np.diff(insolation_grid, axis=0)[:, :-1] * np.diff(rp_grid2, axis=1)[:-1, :]
def lnlike(theta):
    pop = population_model_insolation(theta, insolation_grid, rp_grid2) * new_completeness
    pop = 0.5 * (pop[:-1, :-1] + pop[1:, 1:])
    norm = np.sum(pop * vol)
    #print(norm)
    ll = np.sum(np.log(population_model_insolation(theta, koi_insolation, koi_rps))) - norm
    return ll if np.isfinite(ll) else -np.inf
# The ln-probability function is just proportional to the ln-likelihood
# since we're assuming uniform priors.
bounds = [(-5, 5), (-5, 5), (-5, 5)]
def lnprob(theta):
    # Broad uniform priors.
    for t, rng in zip(theta, bounds):
        if not rng[0] < t < rng[1]:
            return -np.inf
    return lnlike(theta)

# The negative ln-likelihood is useful for optimization.
# Optimizers want to *minimize* your function.
def nll(theta):
    ll = lnlike(theta)
    return -ll if np.isfinite(ll) else 1e15

In [239]:
#%%=========================================================
# Maximum likelihood solution by minimizing negative log-likelihood
from scipy.optimize import minimize
# Initial guess for logF, beta, and alpha
theta_0 = np.array([1, 0.66, -1.5])

r = minimize(nll, theta_0, method="L-BFGS-B", bounds=bounds)
print(r)
# r.x is the vector of parameters (logF, beta, and alpha) from the maximum likelihood solution


      fun: 532.02713091922919
 hess_inv: <3x3 LbfgsInvHessProduct with dtype=float64>
      jac: array([ -5.68434189e-05,  -1.13686838e-04,  -4.54747351e-05])
  message: 'CONVERGENCE: REL_REDUCTION_OF_F_<=_FACTR*EPSMCH'
     nfev: 52
      nit: 10
   status: 0
  success: True
        x: array([ 0.77440222, -1.18780177, -0.97610962])

In [240]:
#%%=========================================================
# Plot the maximum likelihood solution

# We'll reuse these functions to plot all of our results.

# This function plots the density of samples of a double power law,
# as a function of the x0 input, marginalized over the y input
def make_plot(pop_comp, x0, x, y, ax):
    # pop_comp is a two-dimensional array of completeness values
    # pop_comp.shape is 57(insolation) x 61(radius)  
    # x0 is: bin edges of second variable (radius), for model fitting
    # x is rebinning of x0 into a coarser grid for the plots
    # y is: bin edges of first variable (insolation), for model fitting
    # Mid-bin values of the 2D array pop_comp, along the first dimension -- insolation
    pop = 0.5 * (pop_comp[:, 1:] + pop_comp[:, :-1])
    # on first call, pop is 56 x 61

    # Integrate completeness over the first variable, y, insolation
    # np.diff(y)[None, :, None] is the parameter space interval in y
    pop = np.sum(pop * np.diff(y)[None, :, None], axis=1)
    # After above command, pop is 1D, collapsed onto the radius dimension
    
    # Credible regions in x: radius
    # x is used only to get the parameter space interval dx in radius for the plot
    # Note: Assumes bin spacing np.diff(x) is uniform
    a, b, c, d, e = np.percentile(pop * np.diff(x)[0], [2.5, 16, 50, 84, 97.5], axis=0)
    # print(c)

    ax.fill_between(x0, a, e, color="k", alpha=0.1, edgecolor="none")
    ax.fill_between(x0, b, d, color="k", alpha=0.3, edgecolor="none")
    # c is the median value of the distribution over 
    ax.plot(x0, c, "k", lw=1)

def plot_results(samples):
    # Loop through the samples and compute the list of population models.
    samples = np.atleast_2d(samples)
    
    # print(samples.shape)
    # print(len(samples))
    # len(samples) is the length of the MCMC chain
    pop = np.empty((len(samples), insolation_grid.shape[0], insolation_grid.shape[1]))
    gamma_earth = np.empty((len(samples)))
    for i, p in enumerate(samples):
        # insolation_grid and rp_grid2 are meshgrids, 57x61
        # print(p)
        # print(rp_grid2.shape)
        # print(insolation_grid.shape)
        # print(i)
        # power law planet density on insolation, radius grid
        pop[i] = population_model_insolation(p, insolation_grid, rp_grid2)
        # planet density at the point corresponding to earth (insolation = 1 and radius = 1)
        gamma_earth[i] = population_model_insolation(p, 1.0, 1.0) * 1.0
        # print(gamma_earth[i])
        
    # Set up 4x4 grid of plots
    fig, axes = pl.subplots(2, 2, figsize=(10, 8))
    fig.subplots_adjust(wspace=0.4, hspace=0.4)
    
    # Histogram of planet radius over a new grid
    # Using a coarser grid for the plot
    dx = 5*(rp2[1] - rp2[0])
    x = np.arange(rp_rng[0], rp_rng[1] + dx, dx)
    n, _ = np.histogram(koi_rps, x)
    
    # Plot the predicted radius distribution against the observed radius distribution
    ax = axes[0, 0]
    
    # Predicted radius distribution
    # rp2 is of length 61
    # nsolation is of length 57
    # pop is 57 x 61 -- so it's insolation x radius
    # new_completeness is 1 x 57 x 61
    make_plot(pop * new_completeness[None, :, :], rp2, x, insolation, ax)
    
    # Observed radius distribution, with Poisson errors
    ax.errorbar(0.5*(x[:-1]+x[1:]), n, yerr=np.sqrt(n), fmt=".k",
                capsize=0)
    ax.set_xlim(rp_rng[0], rp_rng[1])
    ax.set_xlabel("$R_p\,[R_\oplus]$")
    ax.set_ylabel("# of detected planets")
    
    # Plot the true radius distribution.
    ax = axes[0, 1]
    make_plot(pop, rp2, x, insolation, ax)
    ax.set_xlim(rp_rng[0], rp_rng[1])
    # ax.set_ylim(0, 0.37)
    ax.set_xlabel("$R_p\,[R_\oplus]$")
    ax.set_ylabel("$\mathrm{d}N / \mathrm{d}R$; $\Delta R = 0.25\,R_\oplus$")
    
    # Histogram of insolation over a new grid.
    # Using a coarser grid for the plot
    dx = 5*(insolation[1] - insolation[0])
    
    x = np.arange(insolation_rng[0], insolation_rng[1] + dx, dx)
    n, _ = np.histogram(koi_insolation, x)
    
    # Plot the predicted insolation distribution against the observed insolation distribution
    ax = axes[1, 0]
    
    # Predicted insolation distribution
    make_plot(np.swapaxes(pop * new_completeness[None, :, :], 1, 2), insolation, x, rp2, ax)
    
    # Observed insolation distribution, with Poisson errors
    ax.errorbar(0.5*(x[:-1]+x[1:]), n, yerr=np.sqrt(n), fmt=".k",
                capsize=0)
    # ax.set_xlim(insolation_rng[0], insolation_rng[1])
    ax.set_xlim(0.5, insolation_rng[1])
    ax.set_ylim(0, 25)
    ax.set_xlabel("insolation, [Earth units]")
    ax.set_ylabel("# of detected planets")
    
    # Plot the true insolation distribution.
    ax = axes[1, 1]
    make_plot(np.swapaxes(pop, 1, 2), insolation, x, rp2, ax)
    ax.set_xlim(insolation_rng[0], insolation_rng[1])
    ax.set_xlabel("insolation, [Earth units]")
    ax.set_ylabel("$\mathrm{d}N / \mathrm{d}I$; $\Delta I = $")
    
    return gamma_earth

In [241]:
#%%=========================================================

# This line calls all the plotting machinery above, and returns the maximum likelihood value of gamma_earth
# Note that r.x is the set of parameters [ lnf0, beta, alpha ] returned by the maximum likelihood fit
print(plot_results(r.x));

# Or, try your own values for the parameters
# thetaTry = np.array([.01, -1.2, -1.1])
# print(plot_results(thetaTry));


[ 0.42880752]

In [242]:
#%%=========================================================
# Sample from the posterior probability distribution for the population parameters using emcee
import emcee

ndim, nwalkers = len(r.x), 16
pos = [r.x + 1e-5 * np.random.randn(ndim) for i in range(nwalkers)]
sampler = emcee.EnsembleSampler(nwalkers, ndim, lnprob)

# Burn in.
pos, _, _ = sampler.run_mcmc(pos, 1000)
sampler.reset()

# Production.
pos, _, _ = sampler.run_mcmc(pos, 4000)

In [243]:
#%%=========================================================
# Triangle plot of PDFs using DFMs corner package
import corner
corner.corner(sampler.flatchain, labels=[r"$\ln F$", r"$\beta$", r"$\alpha$"]);



In [244]:
#%%=========================================================
# Plot marginalized posteriors of N and dN/dR
# For Rp, marginalize over P
# For P, marginalize over Rp
# For N, plot also the data and the error bars

# Problem -- plotting all the chains maxes out my PCs memory.    
# Solution -- plot only the last 4000 chains 
gamma_earth = plot_results(sampler.flatchain[60000:63999,:])



In [245]:
#%%=========================================================
# Plot the PDF of gamma_earth
pl.hist(np.log10(gamma_earth), 50, histtype="step", color="k")
pl.gca().set_yticklabels([])
pl.title("the rate of Earth analogs")
pl.xlabel(r"$\log_{10}\Gamma_\oplus $");
#pl.xlabel(r"$\log_{10}\Gamma_\oplus = \left. \log_{10}\mathrm{d}^2 N / \mathrm{d} Insolation \, \mathrm{d} R_p \right |_\oplus$");


Out[245]:
<matplotlib.text.Text at 0x118161cd0>

In [246]:
#%%=========================================================
# Integrate the planet density over a given range in insolation and radius
#   to get the exoplanet occurrence rate predicted by the power law in that region
def integrated_gamma(theta,insolation1,insolation2,radius1,radius2):
    lnf0, beta, alpha = theta
    
    # Parameter Space Boundaries for our model
    insol_rng = (0.2, 20)
    radius_rng = (0.75, 2.5)
    
    # Compute exoplanet occurrence rate integrated over chosen region of [ insolation , radius] parameter space
    integral_over_insolation = (insolation2**(beta + 1) - insolation1**(beta + 1))/(insol_rng[1]**(beta + 1) - insol_rng[0]**(beta + 1))
    integral_over_radius = (radius2**(alpha + 1) - radius1**(alpha + 1))/(radius_rng[1]**(alpha + 1) - radius_rng[0]**(alpha + 1))
    eta = integral_over_insolation*integral_over_radius*np.exp(lnf0)
    
    return eta
    

#%% Compute and plot the posterior PDF for the exoplanet occurence rate in a 
#   desired region of period, radius parameter space

# Initialize   
int_gamma_samples = np.empty(len(sampler.flatchain))

# !!!!! Choose insolation limits according to stellar type

# Select G dwarfs
# insolation1 = 0.295
# insolation2 = 1.824

# Select K dwarfs
insolation1 = 0.235
insolation2 = 1.681

# Select M dwarfs
# insolation1 = 0.205
# insolation2 = 1.514

# !!!!! Choose radius limits
radius1 = 1.0
radius2 = 1.6

for i, theta in enumerate(sampler.flatchain):
    int_gamma_samples[i] = integrated_gamma(theta,insolation1,insolation2,radius1,radius2)

# print result for this parameter space region
print("16, 50, 84 percentile range for integrated gamma = {0} ".format(np.percentile(int_gamma_samples,[16,50,84])))

             
# Plot the posterior of int gamma_samples over selected planet radius and period range             
pl.hist(int_gamma_samples, 50, histtype="step", color="k")
pl.gca().set_yticklabels([])
pl.title("Integrated ocurrence rate over radius and insolation")
pl.xlabel(r"integrated  $\Gamma$ in parameter space region [insolation1, insolation2, radius1, radius2] = {0}".format([insolation1,insolation2,radius1,radius2]));
# pl.xlabel(r"$\log_{10}\Gamma_\oplus = \left. \log_{10}\mathrm{d}N / \mathrm{d}\ln P \, \mathrm{d}\ln R_p \right |_\oplus$");
pl.xticks


16, 50, 84 percentile range for integrated gamma = [ 0.31885863  0.43170725  0.57174896] 
Out[246]:
<function matplotlib.pyplot.xticks>

In [247]:
import pystan
import numpy as np
import gc
np.set_printoptions(threshold='nan')
import matplotlib as plt

#%%===============================================================================================
# Pre-compute some values to be used in the Stan model

# change to column major for PySta version (overwrite old version here)
insolation_grid, rp_grid2 = np.meshgrid(insolation, rp2, indexing="xy")

print(insolation_grid, rp_grid2)
# pre compute completeness for testing PyStan Model.

# Column names for RMS CDPP
cdpp_cols = [k for k in stlr.keys() if k.startswith("rrmscdpp")]

# Values of RMS CDPP corresponding the the 14 pulse durations for all selected stars (column)
#cdpp_obs = np.array(star[cdpp_cols], dtype=float) # <- for one row (one star) change star for stlr.
cdpp_obs_all = np.array(stlr[cdpp_cols], dtype=float)

# Values for the 14 pulse durations derived from the column names
pulse_durations_obs = np.array([k[-4:].replace("p", ".") for k in mesthres_cols],dtype=float)
nPulses = len(pulse_durations_obs)

# Column names for the 14 MES thresholds
mesthres_cols = [k for k in stlr.keys() if k.startswith("mesthres")]

# MES threshold values corresponding to the 14 pulse durations
mesthres_obs_all = np.array(stlr[mesthres_cols],dtype=float)
#print(mesthres_obs_all)                        
# mest <- np.interp(tau, pulse_durations_obs, mesthres_obs, dtype=float)
# x <- mes - 4.1 - (mest - 7.1);

# Period grid
period_rng = (20, 320)
period_grid = np.linspace(period_rng[0], period_rng[1], 57)
nPeriodGrid = len(period_grid)
#print(period_grid)

# Planet radius grid
planet_radius_rng = (0.75, 2.5)
planet_radius_grid = np.linspace(planet_radius_rng[0], planet_radius_rng[1], 61)
nRadiusGrid = len(planet_radius_grid)

# Insolation grid
insolation_rng = (0.2, 20)
insolation_grid = np.linspace(insolation_rng[0], insolation_rng[1], 57)
nInsolationGrid = len(insolation_grid)

# Make 2D meshgrids for bin edges.
planet_radius_grid_vol, insolation_grid_vol = np.meshgrid(planet_radius_grid, insolation_grid, indexing="xy")

# 2D grid of insolation-radius parameter space bin volumes at bin centers
vol = np.diff(insolation_grid_vol, axis=0)[:, :-1] * np.diff(planet_radius_grid_vol, axis=1)[:-1, :]
volShape = list(vol.shape)

# Stellar parameters (for all stars in the stellar catalog)
teffStar = stlr.teff
teffSun = 5777
rStar = stlr.radius
mStar = stlr.mass

# number of stars in the stellar catalog:
nStar = len(stlr)

# Insolation and radius for planets in the planet catalog
# koi_periods = np.array(kois.koi_period)
koi_insolation = np.array(kois.koi_insol)
koi_rps = np.array(kois.koi_prad)

# Number of planets selected from the KOI catalog
nKois = len(kois)

# Eccentricity is set to zero for now.
e = 0.0

# Dataspan is length of time between first and last observations of each star
star_dataspan = stlr.dataspan
# Dutycycle is the fraction of cadences with valid observations
star_dutycycle = stlr.dutycycle

# print(dataspan)
# print(dutycycle)
# print(koi_rps)

# Garbage collection
gc.isenabled()
#gc.get_objects()
gc.collect()


[[  0.2          0.55357143   0.90714286   1.26071429   1.61428571
    1.96785714   2.32142857   2.675        3.02857143   3.38214286
    3.73571429   4.08928571   4.44285714   4.79642857   5.15         5.50357143
    5.85714286   6.21071429   6.56428571   6.91785714   7.27142857   7.625
    7.97857143   8.33214286   8.68571429   9.03928571   9.39285714
    9.74642857  10.1         10.45357143  10.80714286  11.16071429
   11.51428571  11.86785714  12.22142857  12.575       12.92857143
   13.28214286  13.63571429  13.98928571  14.34285714  14.69642857  15.05
   15.40357143  15.75714286  16.11071429  16.46428571  16.81785714
   17.17142857  17.525       17.87857143  18.23214286  18.58571429
   18.93928571  19.29285714  19.64642857  20.        ]
 [  0.2          0.55357143   0.90714286   1.26071429   1.61428571
    1.96785714   2.32142857   2.675        3.02857143   3.38214286
    3.73571429   4.08928571   4.44285714   4.79642857   5.15         5.50357143
    5.85714286   6.21071429   6.56428571   6.91785714   7.27142857   7.625
    7.97857143   8.33214286   8.68571429   9.03928571   9.39285714
    9.74642857  10.1         10.45357143  10.80714286  11.16071429
   11.51428571  11.86785714  12.22142857  12.575       12.92857143
   13.28214286  13.63571429  13.98928571  14.34285714  14.69642857  15.05
   15.40357143  15.75714286  16.11071429  16.46428571  16.81785714
   17.17142857  17.525       17.87857143  18.23214286  18.58571429
   18.93928571  19.29285714  19.64642857  20.        ]
 [  0.2          0.55357143   0.90714286   1.26071429   1.61428571
    1.96785714   2.32142857   2.675        3.02857143   3.38214286
    3.73571429   4.08928571   4.44285714   4.79642857   5.15         5.50357143
    5.85714286   6.21071429   6.56428571   6.91785714   7.27142857   7.625
    7.97857143   8.33214286   8.68571429   9.03928571   9.39285714
    9.74642857  10.1         10.45357143  10.80714286  11.16071429
   11.51428571  11.86785714  12.22142857  12.575       12.92857143
   13.28214286  13.63571429  13.98928571  14.34285714  14.69642857  15.05
   15.40357143  15.75714286  16.11071429  16.46428571  16.81785714
   17.17142857  17.525       17.87857143  18.23214286  18.58571429
   18.93928571  19.29285714  19.64642857  20.        ]
 [  0.2          0.55357143   0.90714286   1.26071429   1.61428571
    1.96785714   2.32142857   2.675        3.02857143   3.38214286
    3.73571429   4.08928571   4.44285714   4.79642857   5.15         5.50357143
    5.85714286   6.21071429   6.56428571   6.91785714   7.27142857   7.625
    7.97857143   8.33214286   8.68571429   9.03928571   9.39285714
    9.74642857  10.1         10.45357143  10.80714286  11.16071429
   11.51428571  11.86785714  12.22142857  12.575       12.92857143
   13.28214286  13.63571429  13.98928571  14.34285714  14.69642857  15.05
   15.40357143  15.75714286  16.11071429  16.46428571  16.81785714
   17.17142857  17.525       17.87857143  18.23214286  18.58571429
   18.93928571  19.29285714  19.64642857  20.        ]
 [  0.2          0.55357143   0.90714286   1.26071429   1.61428571
    1.96785714   2.32142857   2.675        3.02857143   3.38214286
    3.73571429   4.08928571   4.44285714   4.79642857   5.15         5.50357143
    5.85714286   6.21071429   6.56428571   6.91785714   7.27142857   7.625
    7.97857143   8.33214286   8.68571429   9.03928571   9.39285714
    9.74642857  10.1         10.45357143  10.80714286  11.16071429
   11.51428571  11.86785714  12.22142857  12.575       12.92857143
   13.28214286  13.63571429  13.98928571  14.34285714  14.69642857  15.05
   15.40357143  15.75714286  16.11071429  16.46428571  16.81785714
   17.17142857  17.525       17.87857143  18.23214286  18.58571429
   18.93928571  19.29285714  19.64642857  20.        ]
 [  0.2          0.55357143   0.90714286   1.26071429   1.61428571
    1.96785714   2.32142857   2.675        3.02857143   3.38214286
    3.73571429   4.08928571   4.44285714   4.79642857   5.15         5.50357143
    5.85714286   6.21071429   6.56428571   6.91785714   7.27142857   7.625
    7.97857143   8.33214286   8.68571429   9.03928571   9.39285714
    9.74642857  10.1         10.45357143  10.80714286  11.16071429
   11.51428571  11.86785714  12.22142857  12.575       12.92857143
   13.28214286  13.63571429  13.98928571  14.34285714  14.69642857  15.05
   15.40357143  15.75714286  16.11071429  16.46428571  16.81785714
   17.17142857  17.525       17.87857143  18.23214286  18.58571429
   18.93928571  19.29285714  19.64642857  20.        ]
 [  0.2          0.55357143   0.90714286   1.26071429   1.61428571
    1.96785714   2.32142857   2.675        3.02857143   3.38214286
    3.73571429   4.08928571   4.44285714   4.79642857   5.15         5.50357143
    5.85714286   6.21071429   6.56428571   6.91785714   7.27142857   7.625
    7.97857143   8.33214286   8.68571429   9.03928571   9.39285714
    9.74642857  10.1         10.45357143  10.80714286  11.16071429
   11.51428571  11.86785714  12.22142857  12.575       12.92857143
   13.28214286  13.63571429  13.98928571  14.34285714  14.69642857  15.05
   15.40357143  15.75714286  16.11071429  16.46428571  16.81785714
   17.17142857  17.525       17.87857143  18.23214286  18.58571429
   18.93928571  19.29285714  19.64642857  20.        ]
 [  0.2          0.55357143   0.90714286   1.26071429   1.61428571
    1.96785714   2.32142857   2.675        3.02857143   3.38214286
    3.73571429   4.08928571   4.44285714   4.79642857   5.15         5.50357143
    5.85714286   6.21071429   6.56428571   6.91785714   7.27142857   7.625
    7.97857143   8.33214286   8.68571429   9.03928571   9.39285714
    9.74642857  10.1         10.45357143  10.80714286  11.16071429
   11.51428571  11.86785714  12.22142857  12.575       12.92857143
   13.28214286  13.63571429  13.98928571  14.34285714  14.69642857  15.05
   15.40357143  15.75714286  16.11071429  16.46428571  16.81785714
   17.17142857  17.525       17.87857143  18.23214286  18.58571429
   18.93928571  19.29285714  19.64642857  20.        ]
 [  0.2          0.55357143   0.90714286   1.26071429   1.61428571
    1.96785714   2.32142857   2.675        3.02857143   3.38214286
    3.73571429   4.08928571   4.44285714   4.79642857   5.15         5.50357143
    5.85714286   6.21071429   6.56428571   6.91785714   7.27142857   7.625
    7.97857143   8.33214286   8.68571429   9.03928571   9.39285714
    9.74642857  10.1         10.45357143  10.80714286  11.16071429
   11.51428571  11.86785714  12.22142857  12.575       12.92857143
   13.28214286  13.63571429  13.98928571  14.34285714  14.69642857  15.05
   15.40357143  15.75714286  16.11071429  16.46428571  16.81785714
   17.17142857  17.525       17.87857143  18.23214286  18.58571429
   18.93928571  19.29285714  19.64642857  20.        ]
 [  0.2          0.55357143   0.90714286   1.26071429   1.61428571
    1.96785714   2.32142857   2.675        3.02857143   3.38214286
    3.73571429   4.08928571   4.44285714   4.79642857   5.15         5.50357143
    5.85714286   6.21071429   6.56428571   6.91785714   7.27142857   7.625
    7.97857143   8.33214286   8.68571429   9.03928571   9.39285714
    9.74642857  10.1         10.45357143  10.80714286  11.16071429
   11.51428571  11.86785714  12.22142857  12.575       12.92857143
   13.28214286  13.63571429  13.98928571  14.34285714  14.69642857  15.05
   15.40357143  15.75714286  16.11071429  16.46428571  16.81785714
   17.17142857  17.525       17.87857143  18.23214286  18.58571429
   18.93928571  19.29285714  19.64642857  20.        ]
 [  0.2          0.55357143   0.90714286   1.26071429   1.61428571
    1.96785714   2.32142857   2.675        3.02857143   3.38214286
    3.73571429   4.08928571   4.44285714   4.79642857   5.15         5.50357143
    5.85714286   6.21071429   6.56428571   6.91785714   7.27142857   7.625
    7.97857143   8.33214286   8.68571429   9.03928571   9.39285714
    9.74642857  10.1         10.45357143  10.80714286  11.16071429
   11.51428571  11.86785714  12.22142857  12.575       12.92857143
   13.28214286  13.63571429  13.98928571  14.34285714  14.69642857  15.05
   15.40357143  15.75714286  16.11071429  16.46428571  16.81785714
   17.17142857  17.525       17.87857143  18.23214286  18.58571429
   18.93928571  19.29285714  19.64642857  20.        ]
 [  0.2          0.55357143   0.90714286   1.26071429   1.61428571
    1.96785714   2.32142857   2.675        3.02857143   3.38214286
    3.73571429   4.08928571   4.44285714   4.79642857   5.15         5.50357143
    5.85714286   6.21071429   6.56428571   6.91785714   7.27142857   7.625
    7.97857143   8.33214286   8.68571429   9.03928571   9.39285714
    9.74642857  10.1         10.45357143  10.80714286  11.16071429
   11.51428571  11.86785714  12.22142857  12.575       12.92857143
   13.28214286  13.63571429  13.98928571  14.34285714  14.69642857  15.05
   15.40357143  15.75714286  16.11071429  16.46428571  16.81785714
   17.17142857  17.525       17.87857143  18.23214286  18.58571429
   18.93928571  19.29285714  19.64642857  20.        ]
 [  0.2          0.55357143   0.90714286   1.26071429   1.61428571
    1.96785714   2.32142857   2.675        3.02857143   3.38214286
    3.73571429   4.08928571   4.44285714   4.79642857   5.15         5.50357143
    5.85714286   6.21071429   6.56428571   6.91785714   7.27142857   7.625
    7.97857143   8.33214286   8.68571429   9.03928571   9.39285714
    9.74642857  10.1         10.45357143  10.80714286  11.16071429
   11.51428571  11.86785714  12.22142857  12.575       12.92857143
   13.28214286  13.63571429  13.98928571  14.34285714  14.69642857  15.05
   15.40357143  15.75714286  16.11071429  16.46428571  16.81785714
   17.17142857  17.525       17.87857143  18.23214286  18.58571429
   18.93928571  19.29285714  19.64642857  20.        ]
 [  0.2          0.55357143   0.90714286   1.26071429   1.61428571
    1.96785714   2.32142857   2.675        3.02857143   3.38214286
    3.73571429   4.08928571   4.44285714   4.79642857   5.15         5.50357143
    5.85714286   6.21071429   6.56428571   6.91785714   7.27142857   7.625
    7.97857143   8.33214286   8.68571429   9.03928571   9.39285714
    9.74642857  10.1         10.45357143  10.80714286  11.16071429
   11.51428571  11.86785714  12.22142857  12.575       12.92857143
   13.28214286  13.63571429  13.98928571  14.34285714  14.69642857  15.05
   15.40357143  15.75714286  16.11071429  16.46428571  16.81785714
   17.17142857  17.525       17.87857143  18.23214286  18.58571429
   18.93928571  19.29285714  19.64642857  20.        ]
 [  0.2          0.55357143   0.90714286   1.26071429   1.61428571
    1.96785714   2.32142857   2.675        3.02857143   3.38214286
    3.73571429   4.08928571   4.44285714   4.79642857   5.15         5.50357143
    5.85714286   6.21071429   6.56428571   6.91785714   7.27142857   7.625
    7.97857143   8.33214286   8.68571429   9.03928571   9.39285714
    9.74642857  10.1         10.45357143  10.80714286  11.16071429
   11.51428571  11.86785714  12.22142857  12.575       12.92857143
   13.28214286  13.63571429  13.98928571  14.34285714  14.69642857  15.05
   15.40357143  15.75714286  16.11071429  16.46428571  16.81785714
   17.17142857  17.525       17.87857143  18.23214286  18.58571429
   18.93928571  19.29285714  19.64642857  20.        ]
 [  0.2          0.55357143   0.90714286   1.26071429   1.61428571
    1.96785714   2.32142857   2.675        3.02857143   3.38214286
    3.73571429   4.08928571   4.44285714   4.79642857   5.15         5.50357143
    5.85714286   6.21071429   6.56428571   6.91785714   7.27142857   7.625
    7.97857143   8.33214286   8.68571429   9.03928571   9.39285714
    9.74642857  10.1         10.45357143  10.80714286  11.16071429
   11.51428571  11.86785714  12.22142857  12.575       12.92857143
   13.28214286  13.63571429  13.98928571  14.34285714  14.69642857  15.05
   15.40357143  15.75714286  16.11071429  16.46428571  16.81785714
   17.17142857  17.525       17.87857143  18.23214286  18.58571429
   18.93928571  19.29285714  19.64642857  20.        ]
 [  0.2          0.55357143   0.90714286   1.26071429   1.61428571
    1.96785714   2.32142857   2.675        3.02857143   3.38214286
    3.73571429   4.08928571   4.44285714   4.79642857   5.15         5.50357143
    5.85714286   6.21071429   6.56428571   6.91785714   7.27142857   7.625
    7.97857143   8.33214286   8.68571429   9.03928571   9.39285714
    9.74642857  10.1         10.45357143  10.80714286  11.16071429
   11.51428571  11.86785714  12.22142857  12.575       12.92857143
   13.28214286  13.63571429  13.98928571  14.34285714  14.69642857  15.05
   15.40357143  15.75714286  16.11071429  16.46428571  16.81785714
   17.17142857  17.525       17.87857143  18.23214286  18.58571429
   18.93928571  19.29285714  19.64642857  20.        ]
 [  0.2          0.55357143   0.90714286   1.26071429   1.61428571
    1.96785714   2.32142857   2.675        3.02857143   3.38214286
    3.73571429   4.08928571   4.44285714   4.79642857   5.15         5.50357143
    5.85714286   6.21071429   6.56428571   6.91785714   7.27142857   7.625
    7.97857143   8.33214286   8.68571429   9.03928571   9.39285714
    9.74642857  10.1         10.45357143  10.80714286  11.16071429
   11.51428571  11.86785714  12.22142857  12.575       12.92857143
   13.28214286  13.63571429  13.98928571  14.34285714  14.69642857  15.05
   15.40357143  15.75714286  16.11071429  16.46428571  16.81785714
   17.17142857  17.525       17.87857143  18.23214286  18.58571429
   18.93928571  19.29285714  19.64642857  20.        ]
 [  0.2          0.55357143   0.90714286   1.26071429   1.61428571
    1.96785714   2.32142857   2.675        3.02857143   3.38214286
    3.73571429   4.08928571   4.44285714   4.79642857   5.15         5.50357143
    5.85714286   6.21071429   6.56428571   6.91785714   7.27142857   7.625
    7.97857143   8.33214286   8.68571429   9.03928571   9.39285714
    9.74642857  10.1         10.45357143  10.80714286  11.16071429
   11.51428571  11.86785714  12.22142857  12.575       12.92857143
   13.28214286  13.63571429  13.98928571  14.34285714  14.69642857  15.05
   15.40357143  15.75714286  16.11071429  16.46428571  16.81785714
   17.17142857  17.525       17.87857143  18.23214286  18.58571429
   18.93928571  19.29285714  19.64642857  20.        ]
 [  0.2          0.55357143   0.90714286   1.26071429   1.61428571
    1.96785714   2.32142857   2.675        3.02857143   3.38214286
    3.73571429   4.08928571   4.44285714   4.79642857   5.15         5.50357143
    5.85714286   6.21071429   6.56428571   6.91785714   7.27142857   7.625
    7.97857143   8.33214286   8.68571429   9.03928571   9.39285714
    9.74642857  10.1         10.45357143  10.80714286  11.16071429
   11.51428571  11.86785714  12.22142857  12.575       12.92857143
   13.28214286  13.63571429  13.98928571  14.34285714  14.69642857  15.05
   15.40357143  15.75714286  16.11071429  16.46428571  16.81785714
   17.17142857  17.525       17.87857143  18.23214286  18.58571429
   18.93928571  19.29285714  19.64642857  20.        ]
 [  0.2          0.55357143   0.90714286   1.26071429   1.61428571
    1.96785714   2.32142857   2.675        3.02857143   3.38214286
    3.73571429   4.08928571   4.44285714   4.79642857   5.15         5.50357143
    5.85714286   6.21071429   6.56428571   6.91785714   7.27142857   7.625
    7.97857143   8.33214286   8.68571429   9.03928571   9.39285714
    9.74642857  10.1         10.45357143  10.80714286  11.16071429
   11.51428571  11.86785714  12.22142857  12.575       12.92857143
   13.28214286  13.63571429  13.98928571  14.34285714  14.69642857  15.05
   15.40357143  15.75714286  16.11071429  16.46428571  16.81785714
   17.17142857  17.525       17.87857143  18.23214286  18.58571429
   18.93928571  19.29285714  19.64642857  20.        ]
 [  0.2          0.55357143   0.90714286   1.26071429   1.61428571
    1.96785714   2.32142857   2.675        3.02857143   3.38214286
    3.73571429   4.08928571   4.44285714   4.79642857   5.15         5.50357143
    5.85714286   6.21071429   6.56428571   6.91785714   7.27142857   7.625
    7.97857143   8.33214286   8.68571429   9.03928571   9.39285714
    9.74642857  10.1         10.45357143  10.80714286  11.16071429
   11.51428571  11.86785714  12.22142857  12.575       12.92857143
   13.28214286  13.63571429  13.98928571  14.34285714  14.69642857  15.05
   15.40357143  15.75714286  16.11071429  16.46428571  16.81785714
   17.17142857  17.525       17.87857143  18.23214286  18.58571429
   18.93928571  19.29285714  19.64642857  20.        ]
 [  0.2          0.55357143   0.90714286   1.26071429   1.61428571
    1.96785714   2.32142857   2.675        3.02857143   3.38214286
    3.73571429   4.08928571   4.44285714   4.79642857   5.15         5.50357143
    5.85714286   6.21071429   6.56428571   6.91785714   7.27142857   7.625
    7.97857143   8.33214286   8.68571429   9.03928571   9.39285714
    9.74642857  10.1         10.45357143  10.80714286  11.16071429
   11.51428571  11.86785714  12.22142857  12.575       12.92857143
   13.28214286  13.63571429  13.98928571  14.34285714  14.69642857  15.05
   15.40357143  15.75714286  16.11071429  16.46428571  16.81785714
   17.17142857  17.525       17.87857143  18.23214286  18.58571429
   18.93928571  19.29285714  19.64642857  20.        ]
 [  0.2          0.55357143   0.90714286   1.26071429   1.61428571
    1.96785714   2.32142857   2.675        3.02857143   3.38214286
    3.73571429   4.08928571   4.44285714   4.79642857   5.15         5.50357143
    5.85714286   6.21071429   6.56428571   6.91785714   7.27142857   7.625
    7.97857143   8.33214286   8.68571429   9.03928571   9.39285714
    9.74642857  10.1         10.45357143  10.80714286  11.16071429
   11.51428571  11.86785714  12.22142857  12.575       12.92857143
   13.28214286  13.63571429  13.98928571  14.34285714  14.69642857  15.05
   15.40357143  15.75714286  16.11071429  16.46428571  16.81785714
   17.17142857  17.525       17.87857143  18.23214286  18.58571429
   18.93928571  19.29285714  19.64642857  20.        ]
 [  0.2          0.55357143   0.90714286   1.26071429   1.61428571
    1.96785714   2.32142857   2.675        3.02857143   3.38214286
    3.73571429   4.08928571   4.44285714   4.79642857   5.15         5.50357143
    5.85714286   6.21071429   6.56428571   6.91785714   7.27142857   7.625
    7.97857143   8.33214286   8.68571429   9.03928571   9.39285714
    9.74642857  10.1         10.45357143  10.80714286  11.16071429
   11.51428571  11.86785714  12.22142857  12.575       12.92857143
   13.28214286  13.63571429  13.98928571  14.34285714  14.69642857  15.05
   15.40357143  15.75714286  16.11071429  16.46428571  16.81785714
   17.17142857  17.525       17.87857143  18.23214286  18.58571429
   18.93928571  19.29285714  19.64642857  20.        ]
 [  0.2          0.55357143   0.90714286   1.26071429   1.61428571
    1.96785714   2.32142857   2.675        3.02857143   3.38214286
    3.73571429   4.08928571   4.44285714   4.79642857   5.15         5.50357143
    5.85714286   6.21071429   6.56428571   6.91785714   7.27142857   7.625
    7.97857143   8.33214286   8.68571429   9.03928571   9.39285714
    9.74642857  10.1         10.45357143  10.80714286  11.16071429
   11.51428571  11.86785714  12.22142857  12.575       12.92857143
   13.28214286  13.63571429  13.98928571  14.34285714  14.69642857  15.05
   15.40357143  15.75714286  16.11071429  16.46428571  16.81785714
   17.17142857  17.525       17.87857143  18.23214286  18.58571429
   18.93928571  19.29285714  19.64642857  20.        ]
 [  0.2          0.55357143   0.90714286   1.26071429   1.61428571
    1.96785714   2.32142857   2.675        3.02857143   3.38214286
    3.73571429   4.08928571   4.44285714   4.79642857   5.15         5.50357143
    5.85714286   6.21071429   6.56428571   6.91785714   7.27142857   7.625
    7.97857143   8.33214286   8.68571429   9.03928571   9.39285714
    9.74642857  10.1         10.45357143  10.80714286  11.16071429
   11.51428571  11.86785714  12.22142857  12.575       12.92857143
   13.28214286  13.63571429  13.98928571  14.34285714  14.69642857  15.05
   15.40357143  15.75714286  16.11071429  16.46428571  16.81785714
   17.17142857  17.525       17.87857143  18.23214286  18.58571429
   18.93928571  19.29285714  19.64642857  20.        ]
 [  0.2          0.55357143   0.90714286   1.26071429   1.61428571
    1.96785714   2.32142857   2.675        3.02857143   3.38214286
    3.73571429   4.08928571   4.44285714   4.79642857   5.15         5.50357143
    5.85714286   6.21071429   6.56428571   6.91785714   7.27142857   7.625
    7.97857143   8.33214286   8.68571429   9.03928571   9.39285714
    9.74642857  10.1         10.45357143  10.80714286  11.16071429
   11.51428571  11.86785714  12.22142857  12.575       12.92857143
   13.28214286  13.63571429  13.98928571  14.34285714  14.69642857  15.05
   15.40357143  15.75714286  16.11071429  16.46428571  16.81785714
   17.17142857  17.525       17.87857143  18.23214286  18.58571429
   18.93928571  19.29285714  19.64642857  20.        ]
 [  0.2          0.55357143   0.90714286   1.26071429   1.61428571
    1.96785714   2.32142857   2.675        3.02857143   3.38214286
    3.73571429   4.08928571   4.44285714   4.79642857   5.15         5.50357143
    5.85714286   6.21071429   6.56428571   6.91785714   7.27142857   7.625
    7.97857143   8.33214286   8.68571429   9.03928571   9.39285714
    9.74642857  10.1         10.45357143  10.80714286  11.16071429
   11.51428571  11.86785714  12.22142857  12.575       12.92857143
   13.28214286  13.63571429  13.98928571  14.34285714  14.69642857  15.05
   15.40357143  15.75714286  16.11071429  16.46428571  16.81785714
   17.17142857  17.525       17.87857143  18.23214286  18.58571429
   18.93928571  19.29285714  19.64642857  20.        ]
 [  0.2          0.55357143   0.90714286   1.26071429   1.61428571
    1.96785714   2.32142857   2.675        3.02857143   3.38214286
    3.73571429   4.08928571   4.44285714   4.79642857   5.15         5.50357143
    5.85714286   6.21071429   6.56428571   6.91785714   7.27142857   7.625
    7.97857143   8.33214286   8.68571429   9.03928571   9.39285714
    9.74642857  10.1         10.45357143  10.80714286  11.16071429
   11.51428571  11.86785714  12.22142857  12.575       12.92857143
   13.28214286  13.63571429  13.98928571  14.34285714  14.69642857  15.05
   15.40357143  15.75714286  16.11071429  16.46428571  16.81785714
   17.17142857  17.525       17.87857143  18.23214286  18.58571429
   18.93928571  19.29285714  19.64642857  20.        ]
 [  0.2          0.55357143   0.90714286   1.26071429   1.61428571
    1.96785714   2.32142857   2.675        3.02857143   3.38214286
    3.73571429   4.08928571   4.44285714   4.79642857   5.15         5.50357143
    5.85714286   6.21071429   6.56428571   6.91785714   7.27142857   7.625
    7.97857143   8.33214286   8.68571429   9.03928571   9.39285714
    9.74642857  10.1         10.45357143  10.80714286  11.16071429
   11.51428571  11.86785714  12.22142857  12.575       12.92857143
   13.28214286  13.63571429  13.98928571  14.34285714  14.69642857  15.05
   15.40357143  15.75714286  16.11071429  16.46428571  16.81785714
   17.17142857  17.525       17.87857143  18.23214286  18.58571429
   18.93928571  19.29285714  19.64642857  20.        ]
 [  0.2          0.55357143   0.90714286   1.26071429   1.61428571
    1.96785714   2.32142857   2.675        3.02857143   3.38214286
    3.73571429   4.08928571   4.44285714   4.79642857   5.15         5.50357143
    5.85714286   6.21071429   6.56428571   6.91785714   7.27142857   7.625
    7.97857143   8.33214286   8.68571429   9.03928571   9.39285714
    9.74642857  10.1         10.45357143  10.80714286  11.16071429
   11.51428571  11.86785714  12.22142857  12.575       12.92857143
   13.28214286  13.63571429  13.98928571  14.34285714  14.69642857  15.05
   15.40357143  15.75714286  16.11071429  16.46428571  16.81785714
   17.17142857  17.525       17.87857143  18.23214286  18.58571429
   18.93928571  19.29285714  19.64642857  20.        ]
 [  0.2          0.55357143   0.90714286   1.26071429   1.61428571
    1.96785714   2.32142857   2.675        3.02857143   3.38214286
    3.73571429   4.08928571   4.44285714   4.79642857   5.15         5.50357143
    5.85714286   6.21071429   6.56428571   6.91785714   7.27142857   7.625
    7.97857143   8.33214286   8.68571429   9.03928571   9.39285714
    9.74642857  10.1         10.45357143  10.80714286  11.16071429
   11.51428571  11.86785714  12.22142857  12.575       12.92857143
   13.28214286  13.63571429  13.98928571  14.34285714  14.69642857  15.05
   15.40357143  15.75714286  16.11071429  16.46428571  16.81785714
   17.17142857  17.525       17.87857143  18.23214286  18.58571429
   18.93928571  19.29285714  19.64642857  20.        ]
 [  0.2          0.55357143   0.90714286   1.26071429   1.61428571
    1.96785714   2.32142857   2.675        3.02857143   3.38214286
    3.73571429   4.08928571   4.44285714   4.79642857   5.15         5.50357143
    5.85714286   6.21071429   6.56428571   6.91785714   7.27142857   7.625
    7.97857143   8.33214286   8.68571429   9.03928571   9.39285714
    9.74642857  10.1         10.45357143  10.80714286  11.16071429
   11.51428571  11.86785714  12.22142857  12.575       12.92857143
   13.28214286  13.63571429  13.98928571  14.34285714  14.69642857  15.05
   15.40357143  15.75714286  16.11071429  16.46428571  16.81785714
   17.17142857  17.525       17.87857143  18.23214286  18.58571429
   18.93928571  19.29285714  19.64642857  20.        ]
 [  0.2          0.55357143   0.90714286   1.26071429   1.61428571
    1.96785714   2.32142857   2.675        3.02857143   3.38214286
    3.73571429   4.08928571   4.44285714   4.79642857   5.15         5.50357143
    5.85714286   6.21071429   6.56428571   6.91785714   7.27142857   7.625
    7.97857143   8.33214286   8.68571429   9.03928571   9.39285714
    9.74642857  10.1         10.45357143  10.80714286  11.16071429
   11.51428571  11.86785714  12.22142857  12.575       12.92857143
   13.28214286  13.63571429  13.98928571  14.34285714  14.69642857  15.05
   15.40357143  15.75714286  16.11071429  16.46428571  16.81785714
   17.17142857  17.525       17.87857143  18.23214286  18.58571429
   18.93928571  19.29285714  19.64642857  20.        ]
 [  0.2          0.55357143   0.90714286   1.26071429   1.61428571
    1.96785714   2.32142857   2.675        3.02857143   3.38214286
    3.73571429   4.08928571   4.44285714   4.79642857   5.15         5.50357143
    5.85714286   6.21071429   6.56428571   6.91785714   7.27142857   7.625
    7.97857143   8.33214286   8.68571429   9.03928571   9.39285714
    9.74642857  10.1         10.45357143  10.80714286  11.16071429
   11.51428571  11.86785714  12.22142857  12.575       12.92857143
   13.28214286  13.63571429  13.98928571  14.34285714  14.69642857  15.05
   15.40357143  15.75714286  16.11071429  16.46428571  16.81785714
   17.17142857  17.525       17.87857143  18.23214286  18.58571429
   18.93928571  19.29285714  19.64642857  20.        ]
 [  0.2          0.55357143   0.90714286   1.26071429   1.61428571
    1.96785714   2.32142857   2.675        3.02857143   3.38214286
    3.73571429   4.08928571   4.44285714   4.79642857   5.15         5.50357143
    5.85714286   6.21071429   6.56428571   6.91785714   7.27142857   7.625
    7.97857143   8.33214286   8.68571429   9.03928571   9.39285714
    9.74642857  10.1         10.45357143  10.80714286  11.16071429
   11.51428571  11.86785714  12.22142857  12.575       12.92857143
   13.28214286  13.63571429  13.98928571  14.34285714  14.69642857  15.05
   15.40357143  15.75714286  16.11071429  16.46428571  16.81785714
   17.17142857  17.525       17.87857143  18.23214286  18.58571429
   18.93928571  19.29285714  19.64642857  20.        ]
 [  0.2          0.55357143   0.90714286   1.26071429   1.61428571
    1.96785714   2.32142857   2.675        3.02857143   3.38214286
    3.73571429   4.08928571   4.44285714   4.79642857   5.15         5.50357143
    5.85714286   6.21071429   6.56428571   6.91785714   7.27142857   7.625
    7.97857143   8.33214286   8.68571429   9.03928571   9.39285714
    9.74642857  10.1         10.45357143  10.80714286  11.16071429
   11.51428571  11.86785714  12.22142857  12.575       12.92857143
   13.28214286  13.63571429  13.98928571  14.34285714  14.69642857  15.05
   15.40357143  15.75714286  16.11071429  16.46428571  16.81785714
   17.17142857  17.525       17.87857143  18.23214286  18.58571429
   18.93928571  19.29285714  19.64642857  20.        ]
 [  0.2          0.55357143   0.90714286   1.26071429   1.61428571
    1.96785714   2.32142857   2.675        3.02857143   3.38214286
    3.73571429   4.08928571   4.44285714   4.79642857   5.15         5.50357143
    5.85714286   6.21071429   6.56428571   6.91785714   7.27142857   7.625
    7.97857143   8.33214286   8.68571429   9.03928571   9.39285714
    9.74642857  10.1         10.45357143  10.80714286  11.16071429
   11.51428571  11.86785714  12.22142857  12.575       12.92857143
   13.28214286  13.63571429  13.98928571  14.34285714  14.69642857  15.05
   15.40357143  15.75714286  16.11071429  16.46428571  16.81785714
   17.17142857  17.525       17.87857143  18.23214286  18.58571429
   18.93928571  19.29285714  19.64642857  20.        ]
 [  0.2          0.55357143   0.90714286   1.26071429   1.61428571
    1.96785714   2.32142857   2.675        3.02857143   3.38214286
    3.73571429   4.08928571   4.44285714   4.79642857   5.15         5.50357143
    5.85714286   6.21071429   6.56428571   6.91785714   7.27142857   7.625
    7.97857143   8.33214286   8.68571429   9.03928571   9.39285714
    9.74642857  10.1         10.45357143  10.80714286  11.16071429
   11.51428571  11.86785714  12.22142857  12.575       12.92857143
   13.28214286  13.63571429  13.98928571  14.34285714  14.69642857  15.05
   15.40357143  15.75714286  16.11071429  16.46428571  16.81785714
   17.17142857  17.525       17.87857143  18.23214286  18.58571429
   18.93928571  19.29285714  19.64642857  20.        ]
 [  0.2          0.55357143   0.90714286   1.26071429   1.61428571
    1.96785714   2.32142857   2.675        3.02857143   3.38214286
    3.73571429   4.08928571   4.44285714   4.79642857   5.15         5.50357143
    5.85714286   6.21071429   6.56428571   6.91785714   7.27142857   7.625
    7.97857143   8.33214286   8.68571429   9.03928571   9.39285714
    9.74642857  10.1         10.45357143  10.80714286  11.16071429
   11.51428571  11.86785714  12.22142857  12.575       12.92857143
   13.28214286  13.63571429  13.98928571  14.34285714  14.69642857  15.05
   15.40357143  15.75714286  16.11071429  16.46428571  16.81785714
   17.17142857  17.525       17.87857143  18.23214286  18.58571429
   18.93928571  19.29285714  19.64642857  20.        ]
 [  0.2          0.55357143   0.90714286   1.26071429   1.61428571
    1.96785714   2.32142857   2.675        3.02857143   3.38214286
    3.73571429   4.08928571   4.44285714   4.79642857   5.15         5.50357143
    5.85714286   6.21071429   6.56428571   6.91785714   7.27142857   7.625
    7.97857143   8.33214286   8.68571429   9.03928571   9.39285714
    9.74642857  10.1         10.45357143  10.80714286  11.16071429
   11.51428571  11.86785714  12.22142857  12.575       12.92857143
   13.28214286  13.63571429  13.98928571  14.34285714  14.69642857  15.05
   15.40357143  15.75714286  16.11071429  16.46428571  16.81785714
   17.17142857  17.525       17.87857143  18.23214286  18.58571429
   18.93928571  19.29285714  19.64642857  20.        ]
 [  0.2          0.55357143   0.90714286   1.26071429   1.61428571
    1.96785714   2.32142857   2.675        3.02857143   3.38214286
    3.73571429   4.08928571   4.44285714   4.79642857   5.15         5.50357143
    5.85714286   6.21071429   6.56428571   6.91785714   7.27142857   7.625
    7.97857143   8.33214286   8.68571429   9.03928571   9.39285714
    9.74642857  10.1         10.45357143  10.80714286  11.16071429
   11.51428571  11.86785714  12.22142857  12.575       12.92857143
   13.28214286  13.63571429  13.98928571  14.34285714  14.69642857  15.05
   15.40357143  15.75714286  16.11071429  16.46428571  16.81785714
   17.17142857  17.525       17.87857143  18.23214286  18.58571429
   18.93928571  19.29285714  19.64642857  20.        ]
 [  0.2          0.55357143   0.90714286   1.26071429   1.61428571
    1.96785714   2.32142857   2.675        3.02857143   3.38214286
    3.73571429   4.08928571   4.44285714   4.79642857   5.15         5.50357143
    5.85714286   6.21071429   6.56428571   6.91785714   7.27142857   7.625
    7.97857143   8.33214286   8.68571429   9.03928571   9.39285714
    9.74642857  10.1         10.45357143  10.80714286  11.16071429
   11.51428571  11.86785714  12.22142857  12.575       12.92857143
   13.28214286  13.63571429  13.98928571  14.34285714  14.69642857  15.05
   15.40357143  15.75714286  16.11071429  16.46428571  16.81785714
   17.17142857  17.525       17.87857143  18.23214286  18.58571429
   18.93928571  19.29285714  19.64642857  20.        ]
 [  0.2          0.55357143   0.90714286   1.26071429   1.61428571
    1.96785714   2.32142857   2.675        3.02857143   3.38214286
    3.73571429   4.08928571   4.44285714   4.79642857   5.15         5.50357143
    5.85714286   6.21071429   6.56428571   6.91785714   7.27142857   7.625
    7.97857143   8.33214286   8.68571429   9.03928571   9.39285714
    9.74642857  10.1         10.45357143  10.80714286  11.16071429
   11.51428571  11.86785714  12.22142857  12.575       12.92857143
   13.28214286  13.63571429  13.98928571  14.34285714  14.69642857  15.05
   15.40357143  15.75714286  16.11071429  16.46428571  16.81785714
   17.17142857  17.525       17.87857143  18.23214286  18.58571429
   18.93928571  19.29285714  19.64642857  20.        ]
 [  0.2          0.55357143   0.90714286   1.26071429   1.61428571
    1.96785714   2.32142857   2.675        3.02857143   3.38214286
    3.73571429   4.08928571   4.44285714   4.79642857   5.15         5.50357143
    5.85714286   6.21071429   6.56428571   6.91785714   7.27142857   7.625
    7.97857143   8.33214286   8.68571429   9.03928571   9.39285714
    9.74642857  10.1         10.45357143  10.80714286  11.16071429
   11.51428571  11.86785714  12.22142857  12.575       12.92857143
   13.28214286  13.63571429  13.98928571  14.34285714  14.69642857  15.05
   15.40357143  15.75714286  16.11071429  16.46428571  16.81785714
   17.17142857  17.525       17.87857143  18.23214286  18.58571429
   18.93928571  19.29285714  19.64642857  20.        ]
 [  0.2          0.55357143   0.90714286   1.26071429   1.61428571
    1.96785714   2.32142857   2.675        3.02857143   3.38214286
    3.73571429   4.08928571   4.44285714   4.79642857   5.15         5.50357143
    5.85714286   6.21071429   6.56428571   6.91785714   7.27142857   7.625
    7.97857143   8.33214286   8.68571429   9.03928571   9.39285714
    9.74642857  10.1         10.45357143  10.80714286  11.16071429
   11.51428571  11.86785714  12.22142857  12.575       12.92857143
   13.28214286  13.63571429  13.98928571  14.34285714  14.69642857  15.05
   15.40357143  15.75714286  16.11071429  16.46428571  16.81785714
   17.17142857  17.525       17.87857143  18.23214286  18.58571429
   18.93928571  19.29285714  19.64642857  20.        ]
 [  0.2          0.55357143   0.90714286   1.26071429   1.61428571
    1.96785714   2.32142857   2.675        3.02857143   3.38214286
    3.73571429   4.08928571   4.44285714   4.79642857   5.15         5.50357143
    5.85714286   6.21071429   6.56428571   6.91785714   7.27142857   7.625
    7.97857143   8.33214286   8.68571429   9.03928571   9.39285714
    9.74642857  10.1         10.45357143  10.80714286  11.16071429
   11.51428571  11.86785714  12.22142857  12.575       12.92857143
   13.28214286  13.63571429  13.98928571  14.34285714  14.69642857  15.05
   15.40357143  15.75714286  16.11071429  16.46428571  16.81785714
   17.17142857  17.525       17.87857143  18.23214286  18.58571429
   18.93928571  19.29285714  19.64642857  20.        ]
 [  0.2          0.55357143   0.90714286   1.26071429   1.61428571
    1.96785714   2.32142857   2.675        3.02857143   3.38214286
    3.73571429   4.08928571   4.44285714   4.79642857   5.15         5.50357143
    5.85714286   6.21071429   6.56428571   6.91785714   7.27142857   7.625
    7.97857143   8.33214286   8.68571429   9.03928571   9.39285714
    9.74642857  10.1         10.45357143  10.80714286  11.16071429
   11.51428571  11.86785714  12.22142857  12.575       12.92857143
   13.28214286  13.63571429  13.98928571  14.34285714  14.69642857  15.05
   15.40357143  15.75714286  16.11071429  16.46428571  16.81785714
   17.17142857  17.525       17.87857143  18.23214286  18.58571429
   18.93928571  19.29285714  19.64642857  20.        ]
 [  0.2          0.55357143   0.90714286   1.26071429   1.61428571
    1.96785714   2.32142857   2.675        3.02857143   3.38214286
    3.73571429   4.08928571   4.44285714   4.79642857   5.15         5.50357143
    5.85714286   6.21071429   6.56428571   6.91785714   7.27142857   7.625
    7.97857143   8.33214286   8.68571429   9.03928571   9.39285714
    9.74642857  10.1         10.45357143  10.80714286  11.16071429
   11.51428571  11.86785714  12.22142857  12.575       12.92857143
   13.28214286  13.63571429  13.98928571  14.34285714  14.69642857  15.05
   15.40357143  15.75714286  16.11071429  16.46428571  16.81785714
   17.17142857  17.525       17.87857143  18.23214286  18.58571429
   18.93928571  19.29285714  19.64642857  20.        ]
 [  0.2          0.55357143   0.90714286   1.26071429   1.61428571
    1.96785714   2.32142857   2.675        3.02857143   3.38214286
    3.73571429   4.08928571   4.44285714   4.79642857   5.15         5.50357143
    5.85714286   6.21071429   6.56428571   6.91785714   7.27142857   7.625
    7.97857143   8.33214286   8.68571429   9.03928571   9.39285714
    9.74642857  10.1         10.45357143  10.80714286  11.16071429
   11.51428571  11.86785714  12.22142857  12.575       12.92857143
   13.28214286  13.63571429  13.98928571  14.34285714  14.69642857  15.05
   15.40357143  15.75714286  16.11071429  16.46428571  16.81785714
   17.17142857  17.525       17.87857143  18.23214286  18.58571429
   18.93928571  19.29285714  19.64642857  20.        ]
 [  0.2          0.55357143   0.90714286   1.26071429   1.61428571
    1.96785714   2.32142857   2.675        3.02857143   3.38214286
    3.73571429   4.08928571   4.44285714   4.79642857   5.15         5.50357143
    5.85714286   6.21071429   6.56428571   6.91785714   7.27142857   7.625
    7.97857143   8.33214286   8.68571429   9.03928571   9.39285714
    9.74642857  10.1         10.45357143  10.80714286  11.16071429
   11.51428571  11.86785714  12.22142857  12.575       12.92857143
   13.28214286  13.63571429  13.98928571  14.34285714  14.69642857  15.05
   15.40357143  15.75714286  16.11071429  16.46428571  16.81785714
   17.17142857  17.525       17.87857143  18.23214286  18.58571429
   18.93928571  19.29285714  19.64642857  20.        ]
 [  0.2          0.55357143   0.90714286   1.26071429   1.61428571
    1.96785714   2.32142857   2.675        3.02857143   3.38214286
    3.73571429   4.08928571   4.44285714   4.79642857   5.15         5.50357143
    5.85714286   6.21071429   6.56428571   6.91785714   7.27142857   7.625
    7.97857143   8.33214286   8.68571429   9.03928571   9.39285714
    9.74642857  10.1         10.45357143  10.80714286  11.16071429
   11.51428571  11.86785714  12.22142857  12.575       12.92857143
   13.28214286  13.63571429  13.98928571  14.34285714  14.69642857  15.05
   15.40357143  15.75714286  16.11071429  16.46428571  16.81785714
   17.17142857  17.525       17.87857143  18.23214286  18.58571429
   18.93928571  19.29285714  19.64642857  20.        ]
 [  0.2          0.55357143   0.90714286   1.26071429   1.61428571
    1.96785714   2.32142857   2.675        3.02857143   3.38214286
    3.73571429   4.08928571   4.44285714   4.79642857   5.15         5.50357143
    5.85714286   6.21071429   6.56428571   6.91785714   7.27142857   7.625
    7.97857143   8.33214286   8.68571429   9.03928571   9.39285714
    9.74642857  10.1         10.45357143  10.80714286  11.16071429
   11.51428571  11.86785714  12.22142857  12.575       12.92857143
   13.28214286  13.63571429  13.98928571  14.34285714  14.69642857  15.05
   15.40357143  15.75714286  16.11071429  16.46428571  16.81785714
   17.17142857  17.525       17.87857143  18.23214286  18.58571429
   18.93928571  19.29285714  19.64642857  20.        ]
 [  0.2          0.55357143   0.90714286   1.26071429   1.61428571
    1.96785714   2.32142857   2.675        3.02857143   3.38214286
    3.73571429   4.08928571   4.44285714   4.79642857   5.15         5.50357143
    5.85714286   6.21071429   6.56428571   6.91785714   7.27142857   7.625
    7.97857143   8.33214286   8.68571429   9.03928571   9.39285714
    9.74642857  10.1         10.45357143  10.80714286  11.16071429
   11.51428571  11.86785714  12.22142857  12.575       12.92857143
   13.28214286  13.63571429  13.98928571  14.34285714  14.69642857  15.05
   15.40357143  15.75714286  16.11071429  16.46428571  16.81785714
   17.17142857  17.525       17.87857143  18.23214286  18.58571429
   18.93928571  19.29285714  19.64642857  20.        ]
 [  0.2          0.55357143   0.90714286   1.26071429   1.61428571
    1.96785714   2.32142857   2.675        3.02857143   3.38214286
    3.73571429   4.08928571   4.44285714   4.79642857   5.15         5.50357143
    5.85714286   6.21071429   6.56428571   6.91785714   7.27142857   7.625
    7.97857143   8.33214286   8.68571429   9.03928571   9.39285714
    9.74642857  10.1         10.45357143  10.80714286  11.16071429
   11.51428571  11.86785714  12.22142857  12.575       12.92857143
   13.28214286  13.63571429  13.98928571  14.34285714  14.69642857  15.05
   15.40357143  15.75714286  16.11071429  16.46428571  16.81785714
   17.17142857  17.525       17.87857143  18.23214286  18.58571429
   18.93928571  19.29285714  19.64642857  20.        ]
 [  0.2          0.55357143   0.90714286   1.26071429   1.61428571
    1.96785714   2.32142857   2.675        3.02857143   3.38214286
    3.73571429   4.08928571   4.44285714   4.79642857   5.15         5.50357143
    5.85714286   6.21071429   6.56428571   6.91785714   7.27142857   7.625
    7.97857143   8.33214286   8.68571429   9.03928571   9.39285714
    9.74642857  10.1         10.45357143  10.80714286  11.16071429
   11.51428571  11.86785714  12.22142857  12.575       12.92857143
   13.28214286  13.63571429  13.98928571  14.34285714  14.69642857  15.05
   15.40357143  15.75714286  16.11071429  16.46428571  16.81785714
   17.17142857  17.525       17.87857143  18.23214286  18.58571429
   18.93928571  19.29285714  19.64642857  20.        ]
 [  0.2          0.55357143   0.90714286   1.26071429   1.61428571
    1.96785714   2.32142857   2.675        3.02857143   3.38214286
    3.73571429   4.08928571   4.44285714   4.79642857   5.15         5.50357143
    5.85714286   6.21071429   6.56428571   6.91785714   7.27142857   7.625
    7.97857143   8.33214286   8.68571429   9.03928571   9.39285714
    9.74642857  10.1         10.45357143  10.80714286  11.16071429
   11.51428571  11.86785714  12.22142857  12.575       12.92857143
   13.28214286  13.63571429  13.98928571  14.34285714  14.69642857  15.05
   15.40357143  15.75714286  16.11071429  16.46428571  16.81785714
   17.17142857  17.525       17.87857143  18.23214286  18.58571429
   18.93928571  19.29285714  19.64642857  20.        ]
 [  0.2          0.55357143   0.90714286   1.26071429   1.61428571
    1.96785714   2.32142857   2.675        3.02857143   3.38214286
    3.73571429   4.08928571   4.44285714   4.79642857   5.15         5.50357143
    5.85714286   6.21071429   6.56428571   6.91785714   7.27142857   7.625
    7.97857143   8.33214286   8.68571429   9.03928571   9.39285714
    9.74642857  10.1         10.45357143  10.80714286  11.16071429
   11.51428571  11.86785714  12.22142857  12.575       12.92857143
   13.28214286  13.63571429  13.98928571  14.34285714  14.69642857  15.05
   15.40357143  15.75714286  16.11071429  16.46428571  16.81785714
   17.17142857  17.525       17.87857143  18.23214286  18.58571429
   18.93928571  19.29285714  19.64642857  20.        ]
 [  0.2          0.55357143   0.90714286   1.26071429   1.61428571
    1.96785714   2.32142857   2.675        3.02857143   3.38214286
    3.73571429   4.08928571   4.44285714   4.79642857   5.15         5.50357143
    5.85714286   6.21071429   6.56428571   6.91785714   7.27142857   7.625
    7.97857143   8.33214286   8.68571429   9.03928571   9.39285714
    9.74642857  10.1         10.45357143  10.80714286  11.16071429
   11.51428571  11.86785714  12.22142857  12.575       12.92857143
   13.28214286  13.63571429  13.98928571  14.34285714  14.69642857  15.05
   15.40357143  15.75714286  16.11071429  16.46428571  16.81785714
   17.17142857  17.525       17.87857143  18.23214286  18.58571429
   18.93928571  19.29285714  19.64642857  20.        ]
 [  0.2          0.55357143   0.90714286   1.26071429   1.61428571
    1.96785714   2.32142857   2.675        3.02857143   3.38214286
    3.73571429   4.08928571   4.44285714   4.79642857   5.15         5.50357143
    5.85714286   6.21071429   6.56428571   6.91785714   7.27142857   7.625
    7.97857143   8.33214286   8.68571429   9.03928571   9.39285714
    9.74642857  10.1         10.45357143  10.80714286  11.16071429
   11.51428571  11.86785714  12.22142857  12.575       12.92857143
   13.28214286  13.63571429  13.98928571  14.34285714  14.69642857  15.05
   15.40357143  15.75714286  16.11071429  16.46428571  16.81785714
   17.17142857  17.525       17.87857143  18.23214286  18.58571429
   18.93928571  19.29285714  19.64642857  20.        ]] [[ 0.75        0.75        0.75        0.75        0.75        0.75        0.75
   0.75        0.75        0.75        0.75        0.75        0.75        0.75
   0.75        0.75        0.75        0.75        0.75        0.75        0.75
   0.75        0.75        0.75        0.75        0.75        0.75        0.75
   0.75        0.75        0.75        0.75        0.75        0.75        0.75
   0.75        0.75        0.75        0.75        0.75        0.75        0.75
   0.75        0.75        0.75        0.75        0.75        0.75        0.75
   0.75        0.75        0.75        0.75        0.75        0.75        0.75
   0.75      ]
 [ 0.77916667  0.77916667  0.77916667  0.77916667  0.77916667  0.77916667
   0.77916667  0.77916667  0.77916667  0.77916667  0.77916667  0.77916667
   0.77916667  0.77916667  0.77916667  0.77916667  0.77916667  0.77916667
   0.77916667  0.77916667  0.77916667  0.77916667  0.77916667  0.77916667
   0.77916667  0.77916667  0.77916667  0.77916667  0.77916667  0.77916667
   0.77916667  0.77916667  0.77916667  0.77916667  0.77916667  0.77916667
   0.77916667  0.77916667  0.77916667  0.77916667  0.77916667  0.77916667
   0.77916667  0.77916667  0.77916667  0.77916667  0.77916667  0.77916667
   0.77916667  0.77916667  0.77916667  0.77916667  0.77916667  0.77916667
   0.77916667  0.77916667  0.77916667]
 [ 0.80833333  0.80833333  0.80833333  0.80833333  0.80833333  0.80833333
   0.80833333  0.80833333  0.80833333  0.80833333  0.80833333  0.80833333
   0.80833333  0.80833333  0.80833333  0.80833333  0.80833333  0.80833333
   0.80833333  0.80833333  0.80833333  0.80833333  0.80833333  0.80833333
   0.80833333  0.80833333  0.80833333  0.80833333  0.80833333  0.80833333
   0.80833333  0.80833333  0.80833333  0.80833333  0.80833333  0.80833333
   0.80833333  0.80833333  0.80833333  0.80833333  0.80833333  0.80833333
   0.80833333  0.80833333  0.80833333  0.80833333  0.80833333  0.80833333
   0.80833333  0.80833333  0.80833333  0.80833333  0.80833333  0.80833333
   0.80833333  0.80833333  0.80833333]
 [ 0.8375      0.8375      0.8375      0.8375      0.8375      0.8375
   0.8375      0.8375      0.8375      0.8375      0.8375      0.8375
   0.8375      0.8375      0.8375      0.8375      0.8375      0.8375
   0.8375      0.8375      0.8375      0.8375      0.8375      0.8375
   0.8375      0.8375      0.8375      0.8375      0.8375      0.8375
   0.8375      0.8375      0.8375      0.8375      0.8375      0.8375
   0.8375      0.8375      0.8375      0.8375      0.8375      0.8375
   0.8375      0.8375      0.8375      0.8375      0.8375      0.8375
   0.8375      0.8375      0.8375      0.8375      0.8375      0.8375
   0.8375      0.8375      0.8375    ]
 [ 0.86666667  0.86666667  0.86666667  0.86666667  0.86666667  0.86666667
   0.86666667  0.86666667  0.86666667  0.86666667  0.86666667  0.86666667
   0.86666667  0.86666667  0.86666667  0.86666667  0.86666667  0.86666667
   0.86666667  0.86666667  0.86666667  0.86666667  0.86666667  0.86666667
   0.86666667  0.86666667  0.86666667  0.86666667  0.86666667  0.86666667
   0.86666667  0.86666667  0.86666667  0.86666667  0.86666667  0.86666667
   0.86666667  0.86666667  0.86666667  0.86666667  0.86666667  0.86666667
   0.86666667  0.86666667  0.86666667  0.86666667  0.86666667  0.86666667
   0.86666667  0.86666667  0.86666667  0.86666667  0.86666667  0.86666667
   0.86666667  0.86666667  0.86666667]
 [ 0.89583333  0.89583333  0.89583333  0.89583333  0.89583333  0.89583333
   0.89583333  0.89583333  0.89583333  0.89583333  0.89583333  0.89583333
   0.89583333  0.89583333  0.89583333  0.89583333  0.89583333  0.89583333
   0.89583333  0.89583333  0.89583333  0.89583333  0.89583333  0.89583333
   0.89583333  0.89583333  0.89583333  0.89583333  0.89583333  0.89583333
   0.89583333  0.89583333  0.89583333  0.89583333  0.89583333  0.89583333
   0.89583333  0.89583333  0.89583333  0.89583333  0.89583333  0.89583333
   0.89583333  0.89583333  0.89583333  0.89583333  0.89583333  0.89583333
   0.89583333  0.89583333  0.89583333  0.89583333  0.89583333  0.89583333
   0.89583333  0.89583333  0.89583333]
 [ 0.925       0.925       0.925       0.925       0.925       0.925       0.925
   0.925       0.925       0.925       0.925       0.925       0.925       0.925
   0.925       0.925       0.925       0.925       0.925       0.925       0.925
   0.925       0.925       0.925       0.925       0.925       0.925       0.925
   0.925       0.925       0.925       0.925       0.925       0.925       0.925
   0.925       0.925       0.925       0.925       0.925       0.925       0.925
   0.925       0.925       0.925       0.925       0.925       0.925       0.925
   0.925       0.925       0.925       0.925       0.925       0.925       0.925
   0.925     ]
 [ 0.95416667  0.95416667  0.95416667  0.95416667  0.95416667  0.95416667
   0.95416667  0.95416667  0.95416667  0.95416667  0.95416667  0.95416667
   0.95416667  0.95416667  0.95416667  0.95416667  0.95416667  0.95416667
   0.95416667  0.95416667  0.95416667  0.95416667  0.95416667  0.95416667
   0.95416667  0.95416667  0.95416667  0.95416667  0.95416667  0.95416667
   0.95416667  0.95416667  0.95416667  0.95416667  0.95416667  0.95416667
   0.95416667  0.95416667  0.95416667  0.95416667  0.95416667  0.95416667
   0.95416667  0.95416667  0.95416667  0.95416667  0.95416667  0.95416667
   0.95416667  0.95416667  0.95416667  0.95416667  0.95416667  0.95416667
   0.95416667  0.95416667  0.95416667]
 [ 0.98333333  0.98333333  0.98333333  0.98333333  0.98333333  0.98333333
   0.98333333  0.98333333  0.98333333  0.98333333  0.98333333  0.98333333
   0.98333333  0.98333333  0.98333333  0.98333333  0.98333333  0.98333333
   0.98333333  0.98333333  0.98333333  0.98333333  0.98333333  0.98333333
   0.98333333  0.98333333  0.98333333  0.98333333  0.98333333  0.98333333
   0.98333333  0.98333333  0.98333333  0.98333333  0.98333333  0.98333333
   0.98333333  0.98333333  0.98333333  0.98333333  0.98333333  0.98333333
   0.98333333  0.98333333  0.98333333  0.98333333  0.98333333  0.98333333
   0.98333333  0.98333333  0.98333333  0.98333333  0.98333333  0.98333333
   0.98333333  0.98333333  0.98333333]
 [ 1.0125      1.0125      1.0125      1.0125      1.0125      1.0125
   1.0125      1.0125      1.0125      1.0125      1.0125      1.0125
   1.0125      1.0125      1.0125      1.0125      1.0125      1.0125
   1.0125      1.0125      1.0125      1.0125      1.0125      1.0125
   1.0125      1.0125      1.0125      1.0125      1.0125      1.0125
   1.0125      1.0125      1.0125      1.0125      1.0125      1.0125
   1.0125      1.0125      1.0125      1.0125      1.0125      1.0125
   1.0125      1.0125      1.0125      1.0125      1.0125      1.0125
   1.0125      1.0125      1.0125      1.0125      1.0125      1.0125
   1.0125      1.0125      1.0125    ]
 [ 1.04166667  1.04166667  1.04166667  1.04166667  1.04166667  1.04166667
   1.04166667  1.04166667  1.04166667  1.04166667  1.04166667  1.04166667
   1.04166667  1.04166667  1.04166667  1.04166667  1.04166667  1.04166667
   1.04166667  1.04166667  1.04166667  1.04166667  1.04166667  1.04166667
   1.04166667  1.04166667  1.04166667  1.04166667  1.04166667  1.04166667
   1.04166667  1.04166667  1.04166667  1.04166667  1.04166667  1.04166667
   1.04166667  1.04166667  1.04166667  1.04166667  1.04166667  1.04166667
   1.04166667  1.04166667  1.04166667  1.04166667  1.04166667  1.04166667
   1.04166667  1.04166667  1.04166667  1.04166667  1.04166667  1.04166667
   1.04166667  1.04166667  1.04166667]
 [ 1.07083333  1.07083333  1.07083333  1.07083333  1.07083333  1.07083333
   1.07083333  1.07083333  1.07083333  1.07083333  1.07083333  1.07083333
   1.07083333  1.07083333  1.07083333  1.07083333  1.07083333  1.07083333
   1.07083333  1.07083333  1.07083333  1.07083333  1.07083333  1.07083333
   1.07083333  1.07083333  1.07083333  1.07083333  1.07083333  1.07083333
   1.07083333  1.07083333  1.07083333  1.07083333  1.07083333  1.07083333
   1.07083333  1.07083333  1.07083333  1.07083333  1.07083333  1.07083333
   1.07083333  1.07083333  1.07083333  1.07083333  1.07083333  1.07083333
   1.07083333  1.07083333  1.07083333  1.07083333  1.07083333  1.07083333
   1.07083333  1.07083333  1.07083333]
 [ 1.1         1.1         1.1         1.1         1.1         1.1         1.1
   1.1         1.1         1.1         1.1         1.1         1.1         1.1
   1.1         1.1         1.1         1.1         1.1         1.1         1.1
   1.1         1.1         1.1         1.1         1.1         1.1         1.1
   1.1         1.1         1.1         1.1         1.1         1.1         1.1
   1.1         1.1         1.1         1.1         1.1         1.1         1.1
   1.1         1.1         1.1         1.1         1.1         1.1         1.1
   1.1         1.1         1.1         1.1         1.1         1.1         1.1
   1.1       ]
 [ 1.12916667  1.12916667  1.12916667  1.12916667  1.12916667  1.12916667
   1.12916667  1.12916667  1.12916667  1.12916667  1.12916667  1.12916667
   1.12916667  1.12916667  1.12916667  1.12916667  1.12916667  1.12916667
   1.12916667  1.12916667  1.12916667  1.12916667  1.12916667  1.12916667
   1.12916667  1.12916667  1.12916667  1.12916667  1.12916667  1.12916667
   1.12916667  1.12916667  1.12916667  1.12916667  1.12916667  1.12916667
   1.12916667  1.12916667  1.12916667  1.12916667  1.12916667  1.12916667
   1.12916667  1.12916667  1.12916667  1.12916667  1.12916667  1.12916667
   1.12916667  1.12916667  1.12916667  1.12916667  1.12916667  1.12916667
   1.12916667  1.12916667  1.12916667]
 [ 1.15833333  1.15833333  1.15833333  1.15833333  1.15833333  1.15833333
   1.15833333  1.15833333  1.15833333  1.15833333  1.15833333  1.15833333
   1.15833333  1.15833333  1.15833333  1.15833333  1.15833333  1.15833333
   1.15833333  1.15833333  1.15833333  1.15833333  1.15833333  1.15833333
   1.15833333  1.15833333  1.15833333  1.15833333  1.15833333  1.15833333
   1.15833333  1.15833333  1.15833333  1.15833333  1.15833333  1.15833333
   1.15833333  1.15833333  1.15833333  1.15833333  1.15833333  1.15833333
   1.15833333  1.15833333  1.15833333  1.15833333  1.15833333  1.15833333
   1.15833333  1.15833333  1.15833333  1.15833333  1.15833333  1.15833333
   1.15833333  1.15833333  1.15833333]
 [ 1.1875      1.1875      1.1875      1.1875      1.1875      1.1875
   1.1875      1.1875      1.1875      1.1875      1.1875      1.1875
   1.1875      1.1875      1.1875      1.1875      1.1875      1.1875
   1.1875      1.1875      1.1875      1.1875      1.1875      1.1875
   1.1875      1.1875      1.1875      1.1875      1.1875      1.1875
   1.1875      1.1875      1.1875      1.1875      1.1875      1.1875
   1.1875      1.1875      1.1875      1.1875      1.1875      1.1875
   1.1875      1.1875      1.1875      1.1875      1.1875      1.1875
   1.1875      1.1875      1.1875      1.1875      1.1875      1.1875
   1.1875      1.1875      1.1875    ]
 [ 1.21666667  1.21666667  1.21666667  1.21666667  1.21666667  1.21666667
   1.21666667  1.21666667  1.21666667  1.21666667  1.21666667  1.21666667
   1.21666667  1.21666667  1.21666667  1.21666667  1.21666667  1.21666667
   1.21666667  1.21666667  1.21666667  1.21666667  1.21666667  1.21666667
   1.21666667  1.21666667  1.21666667  1.21666667  1.21666667  1.21666667
   1.21666667  1.21666667  1.21666667  1.21666667  1.21666667  1.21666667
   1.21666667  1.21666667  1.21666667  1.21666667  1.21666667  1.21666667
   1.21666667  1.21666667  1.21666667  1.21666667  1.21666667  1.21666667
   1.21666667  1.21666667  1.21666667  1.21666667  1.21666667  1.21666667
   1.21666667  1.21666667  1.21666667]
 [ 1.24583333  1.24583333  1.24583333  1.24583333  1.24583333  1.24583333
   1.24583333  1.24583333  1.24583333  1.24583333  1.24583333  1.24583333
   1.24583333  1.24583333  1.24583333  1.24583333  1.24583333  1.24583333
   1.24583333  1.24583333  1.24583333  1.24583333  1.24583333  1.24583333
   1.24583333  1.24583333  1.24583333  1.24583333  1.24583333  1.24583333
   1.24583333  1.24583333  1.24583333  1.24583333  1.24583333  1.24583333
   1.24583333  1.24583333  1.24583333  1.24583333  1.24583333  1.24583333
   1.24583333  1.24583333  1.24583333  1.24583333  1.24583333  1.24583333
   1.24583333  1.24583333  1.24583333  1.24583333  1.24583333  1.24583333
   1.24583333  1.24583333  1.24583333]
 [ 1.275       1.275       1.275       1.275       1.275       1.275       1.275
   1.275       1.275       1.275       1.275       1.275       1.275       1.275
   1.275       1.275       1.275       1.275       1.275       1.275       1.275
   1.275       1.275       1.275       1.275       1.275       1.275       1.275
   1.275       1.275       1.275       1.275       1.275       1.275       1.275
   1.275       1.275       1.275       1.275       1.275       1.275       1.275
   1.275       1.275       1.275       1.275       1.275       1.275       1.275
   1.275       1.275       1.275       1.275       1.275       1.275       1.275
   1.275     ]
 [ 1.30416667  1.30416667  1.30416667  1.30416667  1.30416667  1.30416667
   1.30416667  1.30416667  1.30416667  1.30416667  1.30416667  1.30416667
   1.30416667  1.30416667  1.30416667  1.30416667  1.30416667  1.30416667
   1.30416667  1.30416667  1.30416667  1.30416667  1.30416667  1.30416667
   1.30416667  1.30416667  1.30416667  1.30416667  1.30416667  1.30416667
   1.30416667  1.30416667  1.30416667  1.30416667  1.30416667  1.30416667
   1.30416667  1.30416667  1.30416667  1.30416667  1.30416667  1.30416667
   1.30416667  1.30416667  1.30416667  1.30416667  1.30416667  1.30416667
   1.30416667  1.30416667  1.30416667  1.30416667  1.30416667  1.30416667
   1.30416667  1.30416667  1.30416667]
 [ 1.33333333  1.33333333  1.33333333  1.33333333  1.33333333  1.33333333
   1.33333333  1.33333333  1.33333333  1.33333333  1.33333333  1.33333333
   1.33333333  1.33333333  1.33333333  1.33333333  1.33333333  1.33333333
   1.33333333  1.33333333  1.33333333  1.33333333  1.33333333  1.33333333
   1.33333333  1.33333333  1.33333333  1.33333333  1.33333333  1.33333333
   1.33333333  1.33333333  1.33333333  1.33333333  1.33333333  1.33333333
   1.33333333  1.33333333  1.33333333  1.33333333  1.33333333  1.33333333
   1.33333333  1.33333333  1.33333333  1.33333333  1.33333333  1.33333333
   1.33333333  1.33333333  1.33333333  1.33333333  1.33333333  1.33333333
   1.33333333  1.33333333  1.33333333]
 [ 1.3625      1.3625      1.3625      1.3625      1.3625      1.3625
   1.3625      1.3625      1.3625      1.3625      1.3625      1.3625
   1.3625      1.3625      1.3625      1.3625      1.3625      1.3625
   1.3625      1.3625      1.3625      1.3625      1.3625      1.3625
   1.3625      1.3625      1.3625      1.3625      1.3625      1.3625
   1.3625      1.3625      1.3625      1.3625      1.3625      1.3625
   1.3625      1.3625      1.3625      1.3625      1.3625      1.3625
   1.3625      1.3625      1.3625      1.3625      1.3625      1.3625
   1.3625      1.3625      1.3625      1.3625      1.3625      1.3625
   1.3625      1.3625      1.3625    ]
 [ 1.39166667  1.39166667  1.39166667  1.39166667  1.39166667  1.39166667
   1.39166667  1.39166667  1.39166667  1.39166667  1.39166667  1.39166667
   1.39166667  1.39166667  1.39166667  1.39166667  1.39166667  1.39166667
   1.39166667  1.39166667  1.39166667  1.39166667  1.39166667  1.39166667
   1.39166667  1.39166667  1.39166667  1.39166667  1.39166667  1.39166667
   1.39166667  1.39166667  1.39166667  1.39166667  1.39166667  1.39166667
   1.39166667  1.39166667  1.39166667  1.39166667  1.39166667  1.39166667
   1.39166667  1.39166667  1.39166667  1.39166667  1.39166667  1.39166667
   1.39166667  1.39166667  1.39166667  1.39166667  1.39166667  1.39166667
   1.39166667  1.39166667  1.39166667]
 [ 1.42083333  1.42083333  1.42083333  1.42083333  1.42083333  1.42083333
   1.42083333  1.42083333  1.42083333  1.42083333  1.42083333  1.42083333
   1.42083333  1.42083333  1.42083333  1.42083333  1.42083333  1.42083333
   1.42083333  1.42083333  1.42083333  1.42083333  1.42083333  1.42083333
   1.42083333  1.42083333  1.42083333  1.42083333  1.42083333  1.42083333
   1.42083333  1.42083333  1.42083333  1.42083333  1.42083333  1.42083333
   1.42083333  1.42083333  1.42083333  1.42083333  1.42083333  1.42083333
   1.42083333  1.42083333  1.42083333  1.42083333  1.42083333  1.42083333
   1.42083333  1.42083333  1.42083333  1.42083333  1.42083333  1.42083333
   1.42083333  1.42083333  1.42083333]
 [ 1.45        1.45        1.45        1.45        1.45        1.45        1.45
   1.45        1.45        1.45        1.45        1.45        1.45        1.45
   1.45        1.45        1.45        1.45        1.45        1.45        1.45
   1.45        1.45        1.45        1.45        1.45        1.45        1.45
   1.45        1.45        1.45        1.45        1.45        1.45        1.45
   1.45        1.45        1.45        1.45        1.45        1.45        1.45
   1.45        1.45        1.45        1.45        1.45        1.45        1.45
   1.45        1.45        1.45        1.45        1.45        1.45        1.45
   1.45      ]
 [ 1.47916667  1.47916667  1.47916667  1.47916667  1.47916667  1.47916667
   1.47916667  1.47916667  1.47916667  1.47916667  1.47916667  1.47916667
   1.47916667  1.47916667  1.47916667  1.47916667  1.47916667  1.47916667
   1.47916667  1.47916667  1.47916667  1.47916667  1.47916667  1.47916667
   1.47916667  1.47916667  1.47916667  1.47916667  1.47916667  1.47916667
   1.47916667  1.47916667  1.47916667  1.47916667  1.47916667  1.47916667
   1.47916667  1.47916667  1.47916667  1.47916667  1.47916667  1.47916667
   1.47916667  1.47916667  1.47916667  1.47916667  1.47916667  1.47916667
   1.47916667  1.47916667  1.47916667  1.47916667  1.47916667  1.47916667
   1.47916667  1.47916667  1.47916667]
 [ 1.50833333  1.50833333  1.50833333  1.50833333  1.50833333  1.50833333
   1.50833333  1.50833333  1.50833333  1.50833333  1.50833333  1.50833333
   1.50833333  1.50833333  1.50833333  1.50833333  1.50833333  1.50833333
   1.50833333  1.50833333  1.50833333  1.50833333  1.50833333  1.50833333
   1.50833333  1.50833333  1.50833333  1.50833333  1.50833333  1.50833333
   1.50833333  1.50833333  1.50833333  1.50833333  1.50833333  1.50833333
   1.50833333  1.50833333  1.50833333  1.50833333  1.50833333  1.50833333
   1.50833333  1.50833333  1.50833333  1.50833333  1.50833333  1.50833333
   1.50833333  1.50833333  1.50833333  1.50833333  1.50833333  1.50833333
   1.50833333  1.50833333  1.50833333]
 [ 1.5375      1.5375      1.5375      1.5375      1.5375      1.5375
   1.5375      1.5375      1.5375      1.5375      1.5375      1.5375
   1.5375      1.5375      1.5375      1.5375      1.5375      1.5375
   1.5375      1.5375      1.5375      1.5375      1.5375      1.5375
   1.5375      1.5375      1.5375      1.5375      1.5375      1.5375
   1.5375      1.5375      1.5375      1.5375      1.5375      1.5375
   1.5375      1.5375      1.5375      1.5375      1.5375      1.5375
   1.5375      1.5375      1.5375      1.5375      1.5375      1.5375
   1.5375      1.5375      1.5375      1.5375      1.5375      1.5375
   1.5375      1.5375      1.5375    ]
 [ 1.56666667  1.56666667  1.56666667  1.56666667  1.56666667  1.56666667
   1.56666667  1.56666667  1.56666667  1.56666667  1.56666667  1.56666667
   1.56666667  1.56666667  1.56666667  1.56666667  1.56666667  1.56666667
   1.56666667  1.56666667  1.56666667  1.56666667  1.56666667  1.56666667
   1.56666667  1.56666667  1.56666667  1.56666667  1.56666667  1.56666667
   1.56666667  1.56666667  1.56666667  1.56666667  1.56666667  1.56666667
   1.56666667  1.56666667  1.56666667  1.56666667  1.56666667  1.56666667
   1.56666667  1.56666667  1.56666667  1.56666667  1.56666667  1.56666667
   1.56666667  1.56666667  1.56666667  1.56666667  1.56666667  1.56666667
   1.56666667  1.56666667  1.56666667]
 [ 1.59583333  1.59583333  1.59583333  1.59583333  1.59583333  1.59583333
   1.59583333  1.59583333  1.59583333  1.59583333  1.59583333  1.59583333
   1.59583333  1.59583333  1.59583333  1.59583333  1.59583333  1.59583333
   1.59583333  1.59583333  1.59583333  1.59583333  1.59583333  1.59583333
   1.59583333  1.59583333  1.59583333  1.59583333  1.59583333  1.59583333
   1.59583333  1.59583333  1.59583333  1.59583333  1.59583333  1.59583333
   1.59583333  1.59583333  1.59583333  1.59583333  1.59583333  1.59583333
   1.59583333  1.59583333  1.59583333  1.59583333  1.59583333  1.59583333
   1.59583333  1.59583333  1.59583333  1.59583333  1.59583333  1.59583333
   1.59583333  1.59583333  1.59583333]
 [ 1.625       1.625       1.625       1.625       1.625       1.625       1.625
   1.625       1.625       1.625       1.625       1.625       1.625       1.625
   1.625       1.625       1.625       1.625       1.625       1.625       1.625
   1.625       1.625       1.625       1.625       1.625       1.625       1.625
   1.625       1.625       1.625       1.625       1.625       1.625       1.625
   1.625       1.625       1.625       1.625       1.625       1.625       1.625
   1.625       1.625       1.625       1.625       1.625       1.625       1.625
   1.625       1.625       1.625       1.625       1.625       1.625       1.625
   1.625     ]
 [ 1.65416667  1.65416667  1.65416667  1.65416667  1.65416667  1.65416667
   1.65416667  1.65416667  1.65416667  1.65416667  1.65416667  1.65416667
   1.65416667  1.65416667  1.65416667  1.65416667  1.65416667  1.65416667
   1.65416667  1.65416667  1.65416667  1.65416667  1.65416667  1.65416667
   1.65416667  1.65416667  1.65416667  1.65416667  1.65416667  1.65416667
   1.65416667  1.65416667  1.65416667  1.65416667  1.65416667  1.65416667
   1.65416667  1.65416667  1.65416667  1.65416667  1.65416667  1.65416667
   1.65416667  1.65416667  1.65416667  1.65416667  1.65416667  1.65416667
   1.65416667  1.65416667  1.65416667  1.65416667  1.65416667  1.65416667
   1.65416667  1.65416667  1.65416667]
 [ 1.68333333  1.68333333  1.68333333  1.68333333  1.68333333  1.68333333
   1.68333333  1.68333333  1.68333333  1.68333333  1.68333333  1.68333333
   1.68333333  1.68333333  1.68333333  1.68333333  1.68333333  1.68333333
   1.68333333  1.68333333  1.68333333  1.68333333  1.68333333  1.68333333
   1.68333333  1.68333333  1.68333333  1.68333333  1.68333333  1.68333333
   1.68333333  1.68333333  1.68333333  1.68333333  1.68333333  1.68333333
   1.68333333  1.68333333  1.68333333  1.68333333  1.68333333  1.68333333
   1.68333333  1.68333333  1.68333333  1.68333333  1.68333333  1.68333333
   1.68333333  1.68333333  1.68333333  1.68333333  1.68333333  1.68333333
   1.68333333  1.68333333  1.68333333]
 [ 1.7125      1.7125      1.7125      1.7125      1.7125      1.7125
   1.7125      1.7125      1.7125      1.7125      1.7125      1.7125
   1.7125      1.7125      1.7125      1.7125      1.7125      1.7125
   1.7125      1.7125      1.7125      1.7125      1.7125      1.7125
   1.7125      1.7125      1.7125      1.7125      1.7125      1.7125
   1.7125      1.7125      1.7125      1.7125      1.7125      1.7125
   1.7125      1.7125      1.7125      1.7125      1.7125      1.7125
   1.7125      1.7125      1.7125      1.7125      1.7125      1.7125
   1.7125      1.7125      1.7125      1.7125      1.7125      1.7125
   1.7125      1.7125      1.7125    ]
 [ 1.74166667  1.74166667  1.74166667  1.74166667  1.74166667  1.74166667
   1.74166667  1.74166667  1.74166667  1.74166667  1.74166667  1.74166667
   1.74166667  1.74166667  1.74166667  1.74166667  1.74166667  1.74166667
   1.74166667  1.74166667  1.74166667  1.74166667  1.74166667  1.74166667
   1.74166667  1.74166667  1.74166667  1.74166667  1.74166667  1.74166667
   1.74166667  1.74166667  1.74166667  1.74166667  1.74166667  1.74166667
   1.74166667  1.74166667  1.74166667  1.74166667  1.74166667  1.74166667
   1.74166667  1.74166667  1.74166667  1.74166667  1.74166667  1.74166667
   1.74166667  1.74166667  1.74166667  1.74166667  1.74166667  1.74166667
   1.74166667  1.74166667  1.74166667]
 [ 1.77083333  1.77083333  1.77083333  1.77083333  1.77083333  1.77083333
   1.77083333  1.77083333  1.77083333  1.77083333  1.77083333  1.77083333
   1.77083333  1.77083333  1.77083333  1.77083333  1.77083333  1.77083333
   1.77083333  1.77083333  1.77083333  1.77083333  1.77083333  1.77083333
   1.77083333  1.77083333  1.77083333  1.77083333  1.77083333  1.77083333
   1.77083333  1.77083333  1.77083333  1.77083333  1.77083333  1.77083333
   1.77083333  1.77083333  1.77083333  1.77083333  1.77083333  1.77083333
   1.77083333  1.77083333  1.77083333  1.77083333  1.77083333  1.77083333
   1.77083333  1.77083333  1.77083333  1.77083333  1.77083333  1.77083333
   1.77083333  1.77083333  1.77083333]
 [ 1.8         1.8         1.8         1.8         1.8         1.8         1.8
   1.8         1.8         1.8         1.8         1.8         1.8         1.8
   1.8         1.8         1.8         1.8         1.8         1.8         1.8
   1.8         1.8         1.8         1.8         1.8         1.8         1.8
   1.8         1.8         1.8         1.8         1.8         1.8         1.8
   1.8         1.8         1.8         1.8         1.8         1.8         1.8
   1.8         1.8         1.8         1.8         1.8         1.8         1.8
   1.8         1.8         1.8         1.8         1.8         1.8         1.8
   1.8       ]
 [ 1.82916667  1.82916667  1.82916667  1.82916667  1.82916667  1.82916667
   1.82916667  1.82916667  1.82916667  1.82916667  1.82916667  1.82916667
   1.82916667  1.82916667  1.82916667  1.82916667  1.82916667  1.82916667
   1.82916667  1.82916667  1.82916667  1.82916667  1.82916667  1.82916667
   1.82916667  1.82916667  1.82916667  1.82916667  1.82916667  1.82916667
   1.82916667  1.82916667  1.82916667  1.82916667  1.82916667  1.82916667
   1.82916667  1.82916667  1.82916667  1.82916667  1.82916667  1.82916667
   1.82916667  1.82916667  1.82916667  1.82916667  1.82916667  1.82916667
   1.82916667  1.82916667  1.82916667  1.82916667  1.82916667  1.82916667
   1.82916667  1.82916667  1.82916667]
 [ 1.85833333  1.85833333  1.85833333  1.85833333  1.85833333  1.85833333
   1.85833333  1.85833333  1.85833333  1.85833333  1.85833333  1.85833333
   1.85833333  1.85833333  1.85833333  1.85833333  1.85833333  1.85833333
   1.85833333  1.85833333  1.85833333  1.85833333  1.85833333  1.85833333
   1.85833333  1.85833333  1.85833333  1.85833333  1.85833333  1.85833333
   1.85833333  1.85833333  1.85833333  1.85833333  1.85833333  1.85833333
   1.85833333  1.85833333  1.85833333  1.85833333  1.85833333  1.85833333
   1.85833333  1.85833333  1.85833333  1.85833333  1.85833333  1.85833333
   1.85833333  1.85833333  1.85833333  1.85833333  1.85833333  1.85833333
   1.85833333  1.85833333  1.85833333]
 [ 1.8875      1.8875      1.8875      1.8875      1.8875      1.8875
   1.8875      1.8875      1.8875      1.8875      1.8875      1.8875
   1.8875      1.8875      1.8875      1.8875      1.8875      1.8875
   1.8875      1.8875      1.8875      1.8875      1.8875      1.8875
   1.8875      1.8875      1.8875      1.8875      1.8875      1.8875
   1.8875      1.8875      1.8875      1.8875      1.8875      1.8875
   1.8875      1.8875      1.8875      1.8875      1.8875      1.8875
   1.8875      1.8875      1.8875      1.8875      1.8875      1.8875
   1.8875      1.8875      1.8875      1.8875      1.8875      1.8875
   1.8875      1.8875      1.8875    ]
 [ 1.91666667  1.91666667  1.91666667  1.91666667  1.91666667  1.91666667
   1.91666667  1.91666667  1.91666667  1.91666667  1.91666667  1.91666667
   1.91666667  1.91666667  1.91666667  1.91666667  1.91666667  1.91666667
   1.91666667  1.91666667  1.91666667  1.91666667  1.91666667  1.91666667
   1.91666667  1.91666667  1.91666667  1.91666667  1.91666667  1.91666667
   1.91666667  1.91666667  1.91666667  1.91666667  1.91666667  1.91666667
   1.91666667  1.91666667  1.91666667  1.91666667  1.91666667  1.91666667
   1.91666667  1.91666667  1.91666667  1.91666667  1.91666667  1.91666667
   1.91666667  1.91666667  1.91666667  1.91666667  1.91666667  1.91666667
   1.91666667  1.91666667  1.91666667]
 [ 1.94583333  1.94583333  1.94583333  1.94583333  1.94583333  1.94583333
   1.94583333  1.94583333  1.94583333  1.94583333  1.94583333  1.94583333
   1.94583333  1.94583333  1.94583333  1.94583333  1.94583333  1.94583333
   1.94583333  1.94583333  1.94583333  1.94583333  1.94583333  1.94583333
   1.94583333  1.94583333  1.94583333  1.94583333  1.94583333  1.94583333
   1.94583333  1.94583333  1.94583333  1.94583333  1.94583333  1.94583333
   1.94583333  1.94583333  1.94583333  1.94583333  1.94583333  1.94583333
   1.94583333  1.94583333  1.94583333  1.94583333  1.94583333  1.94583333
   1.94583333  1.94583333  1.94583333  1.94583333  1.94583333  1.94583333
   1.94583333  1.94583333  1.94583333]
 [ 1.975       1.975       1.975       1.975       1.975       1.975       1.975
   1.975       1.975       1.975       1.975       1.975       1.975       1.975
   1.975       1.975       1.975       1.975       1.975       1.975       1.975
   1.975       1.975       1.975       1.975       1.975       1.975       1.975
   1.975       1.975       1.975       1.975       1.975       1.975       1.975
   1.975       1.975       1.975       1.975       1.975       1.975       1.975
   1.975       1.975       1.975       1.975       1.975       1.975       1.975
   1.975       1.975       1.975       1.975       1.975       1.975       1.975
   1.975     ]
 [ 2.00416667  2.00416667  2.00416667  2.00416667  2.00416667  2.00416667
   2.00416667  2.00416667  2.00416667  2.00416667  2.00416667  2.00416667
   2.00416667  2.00416667  2.00416667  2.00416667  2.00416667  2.00416667
   2.00416667  2.00416667  2.00416667  2.00416667  2.00416667  2.00416667
   2.00416667  2.00416667  2.00416667  2.00416667  2.00416667  2.00416667
   2.00416667  2.00416667  2.00416667  2.00416667  2.00416667  2.00416667
   2.00416667  2.00416667  2.00416667  2.00416667  2.00416667  2.00416667
   2.00416667  2.00416667  2.00416667  2.00416667  2.00416667  2.00416667
   2.00416667  2.00416667  2.00416667  2.00416667  2.00416667  2.00416667
   2.00416667  2.00416667  2.00416667]
 [ 2.03333333  2.03333333  2.03333333  2.03333333  2.03333333  2.03333333
   2.03333333  2.03333333  2.03333333  2.03333333  2.03333333  2.03333333
   2.03333333  2.03333333  2.03333333  2.03333333  2.03333333  2.03333333
   2.03333333  2.03333333  2.03333333  2.03333333  2.03333333  2.03333333
   2.03333333  2.03333333  2.03333333  2.03333333  2.03333333  2.03333333
   2.03333333  2.03333333  2.03333333  2.03333333  2.03333333  2.03333333
   2.03333333  2.03333333  2.03333333  2.03333333  2.03333333  2.03333333
   2.03333333  2.03333333  2.03333333  2.03333333  2.03333333  2.03333333
   2.03333333  2.03333333  2.03333333  2.03333333  2.03333333  2.03333333
   2.03333333  2.03333333  2.03333333]
 [ 2.0625      2.0625      2.0625      2.0625      2.0625      2.0625
   2.0625      2.0625      2.0625      2.0625      2.0625      2.0625
   2.0625      2.0625      2.0625      2.0625      2.0625      2.0625
   2.0625      2.0625      2.0625      2.0625      2.0625      2.0625
   2.0625      2.0625      2.0625      2.0625      2.0625      2.0625
   2.0625      2.0625      2.0625      2.0625      2.0625      2.0625
   2.0625      2.0625      2.0625      2.0625      2.0625      2.0625
   2.0625      2.0625      2.0625      2.0625      2.0625      2.0625
   2.0625      2.0625      2.0625      2.0625      2.0625      2.0625
   2.0625      2.0625      2.0625    ]
 [ 2.09166667  2.09166667  2.09166667  2.09166667  2.09166667  2.09166667
   2.09166667  2.09166667  2.09166667  2.09166667  2.09166667  2.09166667
   2.09166667  2.09166667  2.09166667  2.09166667  2.09166667  2.09166667
   2.09166667  2.09166667  2.09166667  2.09166667  2.09166667  2.09166667
   2.09166667  2.09166667  2.09166667  2.09166667  2.09166667  2.09166667
   2.09166667  2.09166667  2.09166667  2.09166667  2.09166667  2.09166667
   2.09166667  2.09166667  2.09166667  2.09166667  2.09166667  2.09166667
   2.09166667  2.09166667  2.09166667  2.09166667  2.09166667  2.09166667
   2.09166667  2.09166667  2.09166667  2.09166667  2.09166667  2.09166667
   2.09166667  2.09166667  2.09166667]
 [ 2.12083333  2.12083333  2.12083333  2.12083333  2.12083333  2.12083333
   2.12083333  2.12083333  2.12083333  2.12083333  2.12083333  2.12083333
   2.12083333  2.12083333  2.12083333  2.12083333  2.12083333  2.12083333
   2.12083333  2.12083333  2.12083333  2.12083333  2.12083333  2.12083333
   2.12083333  2.12083333  2.12083333  2.12083333  2.12083333  2.12083333
   2.12083333  2.12083333  2.12083333  2.12083333  2.12083333  2.12083333
   2.12083333  2.12083333  2.12083333  2.12083333  2.12083333  2.12083333
   2.12083333  2.12083333  2.12083333  2.12083333  2.12083333  2.12083333
   2.12083333  2.12083333  2.12083333  2.12083333  2.12083333  2.12083333
   2.12083333  2.12083333  2.12083333]
 [ 2.15        2.15        2.15        2.15        2.15        2.15        2.15
   2.15        2.15        2.15        2.15        2.15        2.15        2.15
   2.15        2.15        2.15        2.15        2.15        2.15        2.15
   2.15        2.15        2.15        2.15        2.15        2.15        2.15
   2.15        2.15        2.15        2.15        2.15        2.15        2.15
   2.15        2.15        2.15        2.15        2.15        2.15        2.15
   2.15        2.15        2.15        2.15        2.15        2.15        2.15
   2.15        2.15        2.15        2.15        2.15        2.15        2.15
   2.15      ]
 [ 2.17916667  2.17916667  2.17916667  2.17916667  2.17916667  2.17916667
   2.17916667  2.17916667  2.17916667  2.17916667  2.17916667  2.17916667
   2.17916667  2.17916667  2.17916667  2.17916667  2.17916667  2.17916667
   2.17916667  2.17916667  2.17916667  2.17916667  2.17916667  2.17916667
   2.17916667  2.17916667  2.17916667  2.17916667  2.17916667  2.17916667
   2.17916667  2.17916667  2.17916667  2.17916667  2.17916667  2.17916667
   2.17916667  2.17916667  2.17916667  2.17916667  2.17916667  2.17916667
   2.17916667  2.17916667  2.17916667  2.17916667  2.17916667  2.17916667
   2.17916667  2.17916667  2.17916667  2.17916667  2.17916667  2.17916667
   2.17916667  2.17916667  2.17916667]
 [ 2.20833333  2.20833333  2.20833333  2.20833333  2.20833333  2.20833333
   2.20833333  2.20833333  2.20833333  2.20833333  2.20833333  2.20833333
   2.20833333  2.20833333  2.20833333  2.20833333  2.20833333  2.20833333
   2.20833333  2.20833333  2.20833333  2.20833333  2.20833333  2.20833333
   2.20833333  2.20833333  2.20833333  2.20833333  2.20833333  2.20833333
   2.20833333  2.20833333  2.20833333  2.20833333  2.20833333  2.20833333
   2.20833333  2.20833333  2.20833333  2.20833333  2.20833333  2.20833333
   2.20833333  2.20833333  2.20833333  2.20833333  2.20833333  2.20833333
   2.20833333  2.20833333  2.20833333  2.20833333  2.20833333  2.20833333
   2.20833333  2.20833333  2.20833333]
 [ 2.2375      2.2375      2.2375      2.2375      2.2375      2.2375
   2.2375      2.2375      2.2375      2.2375      2.2375      2.2375
   2.2375      2.2375      2.2375      2.2375      2.2375      2.2375
   2.2375      2.2375      2.2375      2.2375      2.2375      2.2375
   2.2375      2.2375      2.2375      2.2375      2.2375      2.2375
   2.2375      2.2375      2.2375      2.2375      2.2375      2.2375
   2.2375      2.2375      2.2375      2.2375      2.2375      2.2375
   2.2375      2.2375      2.2375      2.2375      2.2375      2.2375
   2.2375      2.2375      2.2375      2.2375      2.2375      2.2375
   2.2375      2.2375      2.2375    ]
 [ 2.26666667  2.26666667  2.26666667  2.26666667  2.26666667  2.26666667
   2.26666667  2.26666667  2.26666667  2.26666667  2.26666667  2.26666667
   2.26666667  2.26666667  2.26666667  2.26666667  2.26666667  2.26666667
   2.26666667  2.26666667  2.26666667  2.26666667  2.26666667  2.26666667
   2.26666667  2.26666667  2.26666667  2.26666667  2.26666667  2.26666667
   2.26666667  2.26666667  2.26666667  2.26666667  2.26666667  2.26666667
   2.26666667  2.26666667  2.26666667  2.26666667  2.26666667  2.26666667
   2.26666667  2.26666667  2.26666667  2.26666667  2.26666667  2.26666667
   2.26666667  2.26666667  2.26666667  2.26666667  2.26666667  2.26666667
   2.26666667  2.26666667  2.26666667]
 [ 2.29583333  2.29583333  2.29583333  2.29583333  2.29583333  2.29583333
   2.29583333  2.29583333  2.29583333  2.29583333  2.29583333  2.29583333
   2.29583333  2.29583333  2.29583333  2.29583333  2.29583333  2.29583333
   2.29583333  2.29583333  2.29583333  2.29583333  2.29583333  2.29583333
   2.29583333  2.29583333  2.29583333  2.29583333  2.29583333  2.29583333
   2.29583333  2.29583333  2.29583333  2.29583333  2.29583333  2.29583333
   2.29583333  2.29583333  2.29583333  2.29583333  2.29583333  2.29583333
   2.29583333  2.29583333  2.29583333  2.29583333  2.29583333  2.29583333
   2.29583333  2.29583333  2.29583333  2.29583333  2.29583333  2.29583333
   2.29583333  2.29583333  2.29583333]
 [ 2.325       2.325       2.325       2.325       2.325       2.325       2.325
   2.325       2.325       2.325       2.325       2.325       2.325       2.325
   2.325       2.325       2.325       2.325       2.325       2.325       2.325
   2.325       2.325       2.325       2.325       2.325       2.325       2.325
   2.325       2.325       2.325       2.325       2.325       2.325       2.325
   2.325       2.325       2.325       2.325       2.325       2.325       2.325
   2.325       2.325       2.325       2.325       2.325       2.325       2.325
   2.325       2.325       2.325       2.325       2.325       2.325       2.325
   2.325     ]
 [ 2.35416667  2.35416667  2.35416667  2.35416667  2.35416667  2.35416667
   2.35416667  2.35416667  2.35416667  2.35416667  2.35416667  2.35416667
   2.35416667  2.35416667  2.35416667  2.35416667  2.35416667  2.35416667
   2.35416667  2.35416667  2.35416667  2.35416667  2.35416667  2.35416667
   2.35416667  2.35416667  2.35416667  2.35416667  2.35416667  2.35416667
   2.35416667  2.35416667  2.35416667  2.35416667  2.35416667  2.35416667
   2.35416667  2.35416667  2.35416667  2.35416667  2.35416667  2.35416667
   2.35416667  2.35416667  2.35416667  2.35416667  2.35416667  2.35416667
   2.35416667  2.35416667  2.35416667  2.35416667  2.35416667  2.35416667
   2.35416667  2.35416667  2.35416667]
 [ 2.38333333  2.38333333  2.38333333  2.38333333  2.38333333  2.38333333
   2.38333333  2.38333333  2.38333333  2.38333333  2.38333333  2.38333333
   2.38333333  2.38333333  2.38333333  2.38333333  2.38333333  2.38333333
   2.38333333  2.38333333  2.38333333  2.38333333  2.38333333  2.38333333
   2.38333333  2.38333333  2.38333333  2.38333333  2.38333333  2.38333333
   2.38333333  2.38333333  2.38333333  2.38333333  2.38333333  2.38333333
   2.38333333  2.38333333  2.38333333  2.38333333  2.38333333  2.38333333
   2.38333333  2.38333333  2.38333333  2.38333333  2.38333333  2.38333333
   2.38333333  2.38333333  2.38333333  2.38333333  2.38333333  2.38333333
   2.38333333  2.38333333  2.38333333]
 [ 2.4125      2.4125      2.4125      2.4125      2.4125      2.4125
   2.4125      2.4125      2.4125      2.4125      2.4125      2.4125
   2.4125      2.4125      2.4125      2.4125      2.4125      2.4125
   2.4125      2.4125      2.4125      2.4125      2.4125      2.4125
   2.4125      2.4125      2.4125      2.4125      2.4125      2.4125
   2.4125      2.4125      2.4125      2.4125      2.4125      2.4125
   2.4125      2.4125      2.4125      2.4125      2.4125      2.4125
   2.4125      2.4125      2.4125      2.4125      2.4125      2.4125
   2.4125      2.4125      2.4125      2.4125      2.4125      2.4125
   2.4125      2.4125      2.4125    ]
 [ 2.44166667  2.44166667  2.44166667  2.44166667  2.44166667  2.44166667
   2.44166667  2.44166667  2.44166667  2.44166667  2.44166667  2.44166667
   2.44166667  2.44166667  2.44166667  2.44166667  2.44166667  2.44166667
   2.44166667  2.44166667  2.44166667  2.44166667  2.44166667  2.44166667
   2.44166667  2.44166667  2.44166667  2.44166667  2.44166667  2.44166667
   2.44166667  2.44166667  2.44166667  2.44166667  2.44166667  2.44166667
   2.44166667  2.44166667  2.44166667  2.44166667  2.44166667  2.44166667
   2.44166667  2.44166667  2.44166667  2.44166667  2.44166667  2.44166667
   2.44166667  2.44166667  2.44166667  2.44166667  2.44166667  2.44166667
   2.44166667  2.44166667  2.44166667]
 [ 2.47083333  2.47083333  2.47083333  2.47083333  2.47083333  2.47083333
   2.47083333  2.47083333  2.47083333  2.47083333  2.47083333  2.47083333
   2.47083333  2.47083333  2.47083333  2.47083333  2.47083333  2.47083333
   2.47083333  2.47083333  2.47083333  2.47083333  2.47083333  2.47083333
   2.47083333  2.47083333  2.47083333  2.47083333  2.47083333  2.47083333
   2.47083333  2.47083333  2.47083333  2.47083333  2.47083333  2.47083333
   2.47083333  2.47083333  2.47083333  2.47083333  2.47083333  2.47083333
   2.47083333  2.47083333  2.47083333  2.47083333  2.47083333  2.47083333
   2.47083333  2.47083333  2.47083333  2.47083333  2.47083333  2.47083333
   2.47083333  2.47083333  2.47083333]
 [ 2.5         2.5         2.5         2.5         2.5         2.5         2.5
   2.5         2.5         2.5         2.5         2.5         2.5         2.5
   2.5         2.5         2.5         2.5         2.5         2.5         2.5
   2.5         2.5         2.5         2.5         2.5         2.5         2.5
   2.5         2.5         2.5         2.5         2.5         2.5         2.5
   2.5         2.5         2.5         2.5         2.5         2.5         2.5
   2.5         2.5         2.5         2.5         2.5         2.5         2.5
   2.5         2.5         2.5         2.5         2.5         2.5         2.5
   2.5       ]]
Out[247]:
4052

In [248]:
# Check that everything is there
print(new_completeness.shape)
print(cdpp_cols)
print(len(cdpp_obs_all))
print(len(pulse_durations_obs))
print(mesthres_cols)
print(len(mesthres_obs_all))
print(nPeriodGrid)
print(nRadiusGrid)
print(nInsolationGrid)
print(planet_radius_grid_vol.shape)
print(insolation_grid_vol.shape)
print(vol.shape)
print(nStar)
print(len(koi_insolation))
print(len(koi_rps))
print(nKois)
print(planet_radius_rng)
print(insolation_rng)
print(insolation_grid.shape)


(57, 61)
['rrmscdpp01p5', 'rrmscdpp02p0', 'rrmscdpp02p5', 'rrmscdpp03p0', 'rrmscdpp03p5', 'rrmscdpp04p5', 'rrmscdpp05p0', 'rrmscdpp06p0', 'rrmscdpp07p5', 'rrmscdpp09p0', 'rrmscdpp10p5', 'rrmscdpp12p0', 'rrmscdpp12p5', 'rrmscdpp15p0']
40685
14
['mesthres01p5', 'mesthres02p0', 'mesthres02p5', 'mesthres03p0', 'mesthres03p5', 'mesthres04p5', 'mesthres05p0', 'mesthres06p0', 'mesthres07p5', 'mesthres09p0', 'mesthres10p5', 'mesthres12p0', 'mesthres12p5', 'mesthres15p0']
40685
57
61
57
(57, 61)
(57, 61)
(56, 60)
40685
117
117
117
(0.75, 2.5)
(0.2, 20)
(57,)

In [255]:
from pystan import StanModel

stan_model_occ = """

functions {


// Interpolation Functions: these were copied from the Stan users group

int intFloor(int leftStart, int rightStart, real iReal) {
  // This is absurd. Use bisection algorithm to find int floor.
  int left;
  int right;
  // int mid; !!!!! But from below, it looks as though mid won't always be int!!!!!
  int mid;

  left <- leftStart;
  right <- rightStart;

  while((left + 1) < right) {
    // int mid; -- moved this to the declarations block at the top
    // print("left, right, mid, i, ", left, ", ", right, ", ", mid, ", ", iReal);
    mid <- left + (right - left) / 2; // is this meant to be 2 (forcing integer division), or 2.?
    if(iReal < mid) {
      right <- mid;
    }
    else {
      left <- mid;
    }
  }
  return left;
}

// Interpolate arr using a non-integral index i
// Note: 1 <= i <= length(arr)

real interpolateLinear(real[] arr, real i) {
  int iLeft;
  real valLeft;
  int iRight;
  real valRight;

  print("interpolating ", i);

  // Get i, value at left. If exact time match, then return value.
  iLeft <- intFloor(1, size(arr), i);
  valLeft <- arr[iLeft];
  if(iLeft == i) {
    return valLeft;
  }

  // Get i, value at right.
  iRight <- iLeft + 1;
  valRight <- arr[iRight];

  // Linearly interpolate between values at left and right.
  return valLeft + (valRight - valLeft) * (i - iLeft);
}

real get_duration(real period, real aor, real e){
/*
    Equation (1) from Burke et al. This estimates the transit
    duration in the same units as the input period. There is a
    typo in the paper (24/4 = 6 != 4).
    
    :param period: the period in any units of your choosing
    :param aor:    the dimensionless semi-major axis (scaled
                  by the stellar radius)
    :param e:      the eccentricity of the orbit
*/

    //return 0.25 * period * sqrt(1 - exp(2)) / aor; !!!!! 7/5/2016 fixed typo !!!!!
    return 0.25 * period * sqrt(1 - e^2) / aor;
}
        
real get_a(real period, real mstar) {
/*
   Compute the semi-major axis of an orbit in Solar radii.
    
    :param period: the period in days
    :param mstar:  the stellar mass in Solar masses
   
    
*/
    real pi;
    real Go4pi;
    
    pi <- 3.14159265359;
    Go4pi <- 2945.4625385377644/(4*pi*pi);
    return (Go4pi*period*period*mstar) ^ (1.0/3.0);
}

real get_delta(real k) {
/*
    Estimate the approximate expected transit depth as a function
    of radius ratio. There might be a typo here. In the paper it
    uses c + s*k but in the public code, it is c - s*k:
    https://github.com/christopherburke/KeplerPORTs
    
    :param k: the dimensionless radius ratio between the planet and
              the star
    
*/ 
    real c;
    real s;
    real delta_max;
    
    c <- 1.0874; 
    s <- 1.0187;
    delta_max <- k*k * (c + s*k);
    
/*  # !!!!! For Q1-Q16 (9.1 pipeline) DFM used a multiplier 0f 0.84 instead of 1.0 in the equation below
    # return 0.84*delta_max
    
    # !!!!! for Q1-Q17 DR24 (9.2 pipeline), use a multiplier of 1 
*/
    return 1.0* delta_max;
}



real get_mes(real rStar, real star_dataspan, real star_dutycycle, real period, real rp, real tau, real[] pulse_durations_obs, real[] cdpp_table) {

/*
    Estimate the multiple event statistic value for a transit. 
    This function operates on one star at a time.  
    
    :param rStar:   stellar properties 
    :param star_dataspan:   stellar properties
    :param star_dutycycle:   stellar properties

    :param period: the period in days
    :param rp:     the planet radius in Earth radii
    :param tau:    the transit duration in hours
    :param sigma:  the rms cdpp in ppm
    
    This needs to eventually account for MES smearing effect / impact parameter. 
    
*/    

    real frac_index;
    real sigma; // interpolated CDPP at pulse duration of this transit.
    real re;  // radius of the Earth
    real k; // radius ratios 
    real snr;
    real ntrn; // number of transits 
    int nPulses;
    
    nPulses <- 14;
    
//  Interpolate to the RMS CDPP corresponding to the transit duration.

//  Fractional index for interpolation functions.
    
    frac_index <- ((tau - pulse_durations_obs[1])/(pulse_durations_obs[14] - pulse_durations_obs[1])) * nPulses;
    
// sigma is estimated RMS CDPP in parts per million (ppm)

    sigma <- interpolateLinear( cdpp_table, frac_index ); 

// Compute the radius ratio and estimate the S/N for a single transit

    re <- 0.009171; //Stellar units
    k <- rp * re / rStar;  // planet to star radius ratio

// Signal to noise of a single transit at this target. 
// Note transit depth is converted to ppm to be consistent with sigma
// TBD: snr needs to be corrected for impact parameter.  Multiply by some scale factor to account 
// for the expected impact parameter. 

    snr <- get_delta(k) * 1e6 / sigma;  
    
// Scale by the estimated number of observed transits to estimate the MES
    
    ntrn <- star_dataspan * star_dutycycle / period ;
    return snr * sqrt(ntrn) ;

// TBD: import and use the numerical window function and the one-sigma depth function instead.
 
}



real get_pdet(real[] cdpp_table, real star_dataspan, real star_dutycycle, real rStar, real aor, real period, real rp, real e, real[] pulse_durations_obs, real[] mesthres_table) {
/*
    Equation (5) from Burke et al. Estimate the detection efficiency
    for a transit at a given star.
    
    :param star:   a pandas row giving the stellar properties !!!!! not needed here
   
    :param aor:    the dimensionless semi-major axis (scaled
                   by the stellar radius)
    :param period: the period in days
    :param rp:     the planet radius in Earth radii
    :param e:      the eccentricity
    :param pulse_durations_obs: a vector of 14 trial pulse durations 
    :param mesthres_table: a vector of 14 MES thresholds corresponding to the 14 trial pulse durations   
    :param cdpp_table: a vector of 14 MES thresholds corresponding to the 14 trial pulse durations   

    
*/

    // tau is the transit duration in hours
    real tau;
    real mes;
    real frac_index;
    real mest;
    real x;
   
    
    tau <- get_duration(period, aor, e) * 24.;
    
    // mes is the Multiple Event Statistic corresponding to this signal; it is like SNR 
    mes <- get_mes(rStar, star_dataspan, star_dutycycle, period, rp, tau, pulse_durations_obs, cdpp_table );
   
    // Fractional index for interpolation functions.
    frac_index <- ((tau - pulse_durations_obs[0])/(pulse_durations_obs[size(pulse_durations_obs)-1] - pulse_durations_obs[0])) * size(pulse_durations_obs);
    
    // Interpolate MES threshold corresponding to transit duration
    mest <- interpolateLinear( mesthres_table, frac_index );
    
    // Argument for the gamma CDF model
    x <- mes - 4.1 - (mest - 7.1);
    
/*    # !!!!! DFM originally used no multiplier in the equation below;
      # The multiplier of 0.78442 must be a 'plateau factor' from Jessie, for 9.2
      # for the 9.2 Q1-Q17 DR24 data, provided by Jessie
*/
    return 0.78442*gamma_cdf(x,103.0113,0.10583);
}


real get_pwin(real star_dataspan, real star_dutycycle, real period) {
/*
    Equation (6) from Burke et al. Estimates the window function
    using a binomial distribution. 
    
    This is over a specific scalar period (period_grid using the jth grid element),
    and for a given star
    
    :param star:   a pandas row giving the stellar properties
    :param period: the period in days
    
*/  
    
    real M; // dataspan divided by period
    real f; // star dutycycle
    real omf; //1 - star dutycycle
    real pw; // probability of the window function
    real msk; // mas to select the probability bigger than zero and data span greater than two years
    
    
    M <- star_dataspan / period;
    f <- star_dutycycle;
    omf <- 1.0 - f;
    pw <- 1 - omf^M - M*f*omf^(M-1) - 0.5*M*(M-1)*f*f*omf^(M-2);
    msk <- (pw >= 0.0) * (M >= 2.0);
    return pw * msk;


}

real get_pgeom(real aor, real e) {
/*
    The geometric transit probability.
    
    See e.g. Kipping (2014) for the eccentricity factor
    http://arxiv.org/abs/1408.1393
    
    :param aor: the dimensionless semi-major axis (scaled
                by the stellar radius)
    :param e:   the orbital eccentricity

*/    
    return 1. / (aor * (1 - e*e)) * (aor>1.0);
}

real[,] get_completeness(int nRadiusGrid, int nPeriodGrid, real[] mesthres_table, real[] pulse_durations_obs, real[] cdpp_table, real mStar, real rStar, real[] period_grid_0, real[] planet_radius_grid, real e, real star_dataspan, real star_dutycycle) {

/*
    A helper function to combine all the completeness effects for one star at a time.
    The inputs are 1D period and radius grids here. 
    
    :param mStar:     mass of star in solar units (one star)
    :param rStar:     radius of star in solar units (one star)
    :param star_dataspan:     length of observation data set in days (one star)
    :param star_dutycycle:    fraction of dataspan with valid data (one star)
 
    :param period_grid:    the period in days  1D grid (over set of periods for defined grid)
    :param planet_radius_grid:   the planet radius in Earth radii   1D grid (over set of radii for defined grid)
    :param e:  the orbital eccentricity (zero for now, over every planet)
    
*/  
    real aor[nPeriodGrid];
    real pdet[nRadiusGrid, nPeriodGrid]; 
    real pwin[nPeriodGrid];
    real pgeom[nPeriodGrid];
    real completeness[nRadiusGrid, nPeriodGrid];
    
    for (i in 1:nRadiusGrid) 
        for (j in 1:nPeriodGrid) {
            aor[j] <- get_a(period_grid_0[j], mStar) / rStar;
            pdet[i,j] <- get_pdet(cdpp_table, star_dataspan, star_dutycycle, rStar, aor[j], period_grid_0[j], planet_radius_grid[i], e, pulse_durations_obs, mesthres_table);
            pwin[j] <- get_pwin(star_dataspan, star_dutycycle, period_grid_0[j]);
            pgeom[j] <- get_pgeom(aor[j], e);
            completeness[i,j] <- pdet[i,j] * pwin[j] * pgeom[j];
    }
    return completeness;

}   


real[] get_period_from_insolation(int nInsolationGrid, real teffStar , real teffSun, real rStar, real mStar, real[] insolation ) {
/*   Convert 1D insolation grid to 1D period grid, given star properties
     :param teffStar:
     :param teffSun:
     :param rStar:
     :param mStar:
     :param insolation:

      Get semimajor axis from star properties and insolation, using
      insolation <- ( teffStar / teffSun )^4 * ( rStar / 1.0)^2 * ( 1 / aPlanet )^2
      
      Get orbit period in days from semimajor axis of planet in AU and start properties, using
      aPlanet <- mStar^(1.0/3.0) * (period/365.25)^(2.0/3.0)
*/
      
      real aPlanet[nInsolationGrid];
      real period[nInsolationGrid];
      
      for (j in 1:nInsolationGrid) {
        aPlanet[j] <- ( ( teffStar / teffSun )^4 * ( rStar / 1.0)^2 / insolation[j] )^(0.5);
      // for (j in 1:nInsolationGrid) 
        period[j] <- 365.25 * ( aPlanet[j]/( mStar^(1.0/3.0) ) )^(3.0/2.0);
        }
      
      return period;
}

//
//real get_insolation_from_period( real teffStar , real teffSun, real rStar, real mStar, real period ) {
//
// /*
//     :param teffStar:
//     :param teffSun:
//     :param rStar:
//     :param mStar:
//     :param period:
//        Leaving this in here for now in order to adapt from insolation to period and vice versa.
// */
//        
//      real aPlanet;
//      real insolation;
//      
//      // Semimajor axis of planet in AU
//      aPlanet <- (mStar^(1.0/3.0)) * ((period/365.25)^(2.0/3.0)); 
//      
//    
//      // Insolation
//      insolation <- ((teffStar/teffSun)^4)*((rStar/1)^2)*((1/aPlanet)^2);
//      
//      return insolation;
// }



real[,] get_completeness_from_insolation(real star_dataspan, real star_dutycycle, real[] cdpp_table, real[] pulse_durations_obs, int nStar, int nRadiusGrid, int nInsolationGrid, real[] mesthres_table, real teffStar, real teffSun, real rStar, real mStar, real[] insolation_grid, real[] planet_radius_grid, real e) {

/*
    This function really means "Get completeness using the period computed from the 
    input 1D insolation grid.
    
    nInsolationGrid equals nPeriodGrid by construction to minimize confusion etc. 
    
    :param teffStar:
    :param teffSun:
    :param rStar:
    :param mStar:
    :param insolation_grid:
    :param planet_radius_grid:
    :param e:
*/  
    
    real completeness[nRadiusGrid,nInsolationGrid]; 
    real period_grid_0[nInsolationGrid];
    
    
    period_grid_0 <- get_period_from_insolation(nInsolationGrid, teffStar, teffSun, rStar, mStar, insolation_grid);

    completeness <- get_completeness(nRadiusGrid, nInsolationGrid, mesthres_table, pulse_durations_obs, cdpp_table, mStar, rStar, period_grid_0, planet_radius_grid, e, star_dataspan, star_dutycycle);

    return completeness;
}    


real[,] sum_completeness(real[,] mesthres_obs_all, real[,] cdpp_obs_all, real[] pulse_durations_obs, real[] star_dataspan, real[] star_dutycycle, int nStar, int nRadiusGrid, int nInsolationGrid, real[] teffStar, real teffSun, real[] rStar, real[] mStar, real[] insolation_grid, real[] planet_radius_grid, real e) {

/* Marginalize detection contours over all selected targets
  including the geometric factor. This takes a few minutes.
  
  This will be a time sink in Stan likely.  Future work: figure 
  clever way to optimize this.  
*/

    real sum_completeness_temp[nRadiusGrid,nInsolationGrid];
    real temp_completeness[nRadiusGrid,nInsolationGrid];
    real cdpp_table[14];
    real mesthres_table[14];



    for (i in 1:nRadiusGrid)
        for (j in 1:nInsolationGrid) {
            sum_completeness_temp[i,j] <- 0.0;
        }
            
    for (k in 1:nStar) {
        for (m in 1:14) {
            cdpp_table[m] <- cdpp_obs_all[k,m];
            mesthres_table[m] <- mesthres_obs_all[k,m];
        }
        temp_completeness <- get_completeness_from_insolation(star_dataspan[k], star_dutycycle[k], cdpp_table, pulse_durations_obs, nStar, nRadiusGrid, nInsolationGrid, mesthres_table, teffStar[k], teffSun, rStar[k], mStar[k], insolation_grid, planet_radius_grid, e);
        for (i in 1:nRadiusGrid)
            for (j in 1:nInsolationGrid) {
                sum_completeness_temp[i,j] <- sum_completeness_temp[i,j] + temp_completeness[i,j];
            }
    }
    return sum_completeness_temp;
}


real[,] population_model_insolation(int nInsolationGrid, int nRadiusGrid, 
real lnf0, real alpha, real beta, real[] insolation_grid, real[] planet_radius_grid, 
real[] planet_radius_rng, real[] insolation_rng) {

/* Population inference with an independent power law model
   Using modified code above that computes 
   completeness in the parameter space of [ insolation , planet radius ], 

   A double power law model for the population.

   :param nInsolationGrid: size of insolation grid
   :param nRadiusGrid: size of radius grid
   :param lnf0: normalization, 
   :param beta: exponent of insolation power law,
   :param alpha: exponent of radius power law
   :param insolation_grid:       //1D grid
   :param planet_radius_grid:    //1D grid
   :param insolation_rng: // lower and upper boundary values of insolation grid
   :param radius_rng:     // lower and upper boundary values of radius grid
   
*/    

// Declare internal variables
   real powerLawNumberDensity[nRadiusGrid,nInsolationGrid];

// Loop over phase space bins, compute number density in each bin, 
//      from the power law using the input lnf0, alpha and beta parameters
//   for (i in 1:nRadiusGrid) 
//     for (j in 1:nInsolationGrid) {
     for (j in 1:nInsolationGrid) 
       for (i in 1:nRadiusGrid) {

       powerLawNumberDensity[i,j] <- ( exp( lnf0 ) * ( alpha+1 ) * ( beta+1 ) * ( planet_radius_grid[i]^alpha ) * ( insolation_grid[j]^beta ) ) / ( ( planet_radius_rng[2]^( alpha+1 ) - ( planet_radius_rng[1]^( alpha+1 ) ) ) * ( ( insolation_rng[2]^( beta+1 ) ) - ( insolation_rng[1]^( beta+1 ) ) ) );
     }
   return powerLawNumberDensity;
}

real lnlike(real alpha, real beta, real lnf0, int nKois, int nStar, real[] planet_radius_rng, 
real[] insolation_rng, real[,] mesthres_obs_all, real[,] cdpp_obs_all, real[] pulse_durations_obs, 
real[] star_dataspan, real[] star_dutycycle, real[,] vol, real[] insolation_grid, real[] planet_radius_grid, 
real[] koi_insolation, real[] koi_rps, int nRadiusGrid, int nInsolationGrid, real[] teffStar, 
real teffSun, real[] rStar, real[] mStar, real e, real[,] new_completeness ) {

/* log likelihood function

    :param alpha:
    :param beta:
    :param lnf0:
    :param nKois
    :param nStar
    :param planet_radius_rng
    :param insolation_rng
    :param mesthres_obs_all
    :param cdpp_obs_all
    :param pulse_durations_obs
    :param star_dataspan
    :param star_dutycycle
    :param vol:
    :param insolation_grid:
    :param planet_radius_grid:
    :param koi_insolation:
    :param koi_rps:
    :param nRadiusGrid
    :param nInsolationGrid
    :param teffStar
    :param teffSun
    :param rStar
    :param mStar
    :param e
    :param new_completeness
        
*/

// Declare internal variables

    real pop[nRadiusGrid,nInsolationGrid];
    real pop1[nRadiusGrid,nInsolationGrid];
    real pop2[nRadiusGrid,nInsolationGrid];
    real populationDensity[nRadiusGrid-1,nInsolationGrid-1];
    real sumExpectedCounts;
    real logPop[nKois];
    real lnLikelihood;
    real sumLogPop;
    real lnLikelihood0;
    real popModelInsol0[nRadiusGrid,nInsolationGrid];
    real popModelInsol[nKois,nKois];
    
     
    // Compute the population model for a grid over [ insolation, radius ] parameter space
    popModelInsol0 <- population_model_insolation(nInsolationGrid, nRadiusGrid, lnf0, alpha, beta, insolation_grid, planet_radius_grid, planet_radius_rng, insolation_rng);

    // Compute the expected density of detections per bin in the [ insolation, radius ] parameter space,
    // based on the population model and the completeness
        for (i in 1:nRadiusGrid)
            for (j in 1: nInsolationGrid) {
//        for (j in 1: nInsolationGrid) 
//            for (i in 1:nRadiusGrid) {

            //sum_completeness_0[i,j] <- new_completeness[j,i];
             pop[i,j] <- popModelInsol0[i,j] * new_completeness[j,i];
            // pop[i,j] <- popModelInsol0[i,j] * new_completeness[i,j];

        }
        
    // Initialize
    sumExpectedCounts <- 0.0; 
    sumLogPop <-0.0;
    
    for (i in 1:nRadiusGrid-1)
      for (j in 1: nInsolationGrid-1) {
//      for (j in 1: nInsolationGrid-1) 
//          for (i in 1:nRadiusGrid-1) {

      
         // Compute the number density of detections on a grid of bin centers
         
         // Figured out the proper syntax for translation of the Python line below into Stan
         // pop <- 0.5*(pop[:-1, :-1] + pop[1:, 1:]) is python syntax from DFM notebook.
         // Note that indices start with 1 in Stan
         pop1[i,j] <- pop[i,j];
         pop2[i,j] <- pop[i+1,j+1];
         populationDensity[i,j] <- 0.5*(pop1[i,j]+pop2[i,j]);
         
         // Integrate the planet counts over the phase space volume
         sumExpectedCounts <- sumExpectedCounts + (populationDensity[i,j]*vol[j,i]);
      }
   
   // Population model evaluated over the list of detected planets
   
   popModelInsol <- population_model_insolation(nKois, nKois, lnf0, alpha, beta, koi_insolation, koi_rps, planet_radius_rng, insolation_rng);

    for (k in 1:nKois) {
        logPop[k] <- log( popModelInsol[k,k] );
        sumLogPop <- sumLogPop + logPop[k];
    }
    
    // Combine the two terms to form the log-likelihood function for the Inhomogeneous Poisson Process 
    lnLikelihood <- sumLogPop - sumExpectedCounts;
    
    // Catch error: if lnLikelihood has bad value, return a very small number
    lnLikelihood0 <- if_else( is_inf(lnLikelihood) || is_nan(lnLikelihood) , -1e15, lnLikelihood );
    
   return lnLikelihood0;
   
} 

}
// end function block


data {

// Declare precomputed variables as data

int<lower=1> nKois;

int<lower=1> nPulses;

int<lower=1> nRadiusGrid;

int<lower=1> nInsolationGrid;

int<lower=1> nStar;

real star_dataspan[nStar];

real star_dutycycle[nStar];

real cdpp_obs_all[nStar, nPulses];

real mesthres_obs_all[nStar, nPulses];

real pulse_durations_obs[nPulses];

real insolation_rng[2]; 

real planet_radius_rng[2];

real insolation_grid[nInsolationGrid];

real planet_radius_grid[nRadiusGrid];

real teffStar[nStar];

real teffSun;

real rStar[nStar];

real mStar[nStar];

// real vol[nRadiusGrid-1,nInsolationGrid-1];
 real vol[nInsolationGrid-1,nRadiusGrid-1];

real koi_insolation[nKois];

real koi_rps[nKois];

real e;

real new_completeness[nInsolationGrid,nRadiusGrid];
//real new_completeness[nRadiusGrid,nInsolationGrid];


}


parameters {

// For v1: only compute population level parameter posteriors.
real<lower=-5,upper=5> alpha; 
real<lower=-5,upper=5> beta;
real<lower=-5,upper=5> lnf0;

}

model {

    // Priors
    alpha ~ uniform(-5,5);
    beta ~ uniform(-5,5);
    lnf0 ~ uniform(-5,5);

    // Sample using lnLikelihood function
    increment_log_prob(lnlike(alpha, beta, lnf0, nKois, nStar, planet_radius_rng, insolation_rng, mesthres_obs_all, cdpp_obs_all, pulse_durations_obs, star_dataspan, star_dutycycle, vol, insolation_grid, planet_radius_grid, koi_insolation, koi_rps, nRadiusGrid, nInsolationGrid, teffStar, teffSun, rStar, mStar, e, new_completeness));
}


"""

data = {'new_completeness':new_completeness, 'star_dutycycle':star_dutycycle, 'star_dataspan':star_dataspan, 
        'nStar':nStar, 'nKois':nKois,'nPulses':nPulses,'nRadiusGrid':nRadiusGrid, 'nInsolationGrid':nInsolationGrid, 
        'cdpp_obs_all':cdpp_obs_all, 'mesthres_obs_all':mesthres_obs_all, 'pulse_durations_obs':pulse_durations_obs, 
        'insolation_rng':insolation_rng, 'planet_radius_rng':planet_radius_rng,
        'insolation_grid':insolation_grid, 'planet_radius_grid':planet_radius_grid, 
        'vol':vol, 'teffStar':teffStar,'teffSun':teffSun, 'rStar':rStar, 'mStar':mStar,
        'koi_insolation':koi_insolation, 'koi_rps':koi_rps, 'e':e}

sm = StanModel(model_code=stan_model_occ)

#init = [{'lnf0':0.82,'alpha':-0.97,'beta':-1.26}]

#fit = sm.sampling(data=data, iter=10, chains=2, init=init, n_jobs=-1)
fit = pystan.stan(model_code=stan_model_occ, data=data, iter=100, chains=5, n_jobs=-1, verbose=True);
#fit = pystan.stan(model_code=stan_model_occ, data=data, iter=100, chains=1, init=init, n_jobs=-1, verbose=True);


#get_inits(fit)

la = fit.extract(permuted=True)  # return a dictionary of arrays
alpha = la['alpha']
beta = la['beta']
lnf0 = la['lnf0']
print(alpha)
print(beta)
print(lnf0)

a = fit.extract(permuted=False)
print(a)
print(fit)
print(alpha)

#plt.histogram(alpha)
#fit.plot()


Compiling /var/folders/v0/fmjb82wj6tzbyvl5rpgkyt600000gp/T/tmpjRtdjS/pystan/stanfit4anon_model_9e5045edf8799b6e16462f3517df64f2_a72157596188ec1daae6f6858b197bc8.pyx because it changed.
[1/1] Cythonizing /var/folders/v0/fmjb82wj6tzbyvl5rpgkyt600000gp/T/tmpjRtdjS/pystan/stanfit4anon_model_9e5045edf8799b6e16462f3517df64f2_a72157596188ec1daae6f6858b197bc8.pyx
building 'stanfit4anon_model_9e5045edf8799b6e16462f3517df64f2_a72157596188ec1daae6f6858b197bc8' extension
creating /var/folders/v0/fmjb82wj6tzbyvl5rpgkyt600000gp/T/tmpjRtdjS/pystan/var
creating /var/folders/v0/fmjb82wj6tzbyvl5rpgkyt600000gp/T/tmpjRtdjS/pystan/var/folders
creating /var/folders/v0/fmjb82wj6tzbyvl5rpgkyt600000gp/T/tmpjRtdjS/pystan/var/folders/v0
creating /var/folders/v0/fmjb82wj6tzbyvl5rpgkyt600000gp/T/tmpjRtdjS/pystan/var/folders/v0/fmjb82wj6tzbyvl5rpgkyt600000gp
creating /var/folders/v0/fmjb82wj6tzbyvl5rpgkyt600000gp/T/tmpjRtdjS/pystan/var/folders/v0/fmjb82wj6tzbyvl5rpgkyt600000gp/T
creating /var/folders/v0/fmjb82wj6tzbyvl5rpgkyt600000gp/T/tmpjRtdjS/pystan/var/folders/v0/fmjb82wj6tzbyvl5rpgkyt600000gp/T/tmpjRtdjS
creating /var/folders/v0/fmjb82wj6tzbyvl5rpgkyt600000gp/T/tmpjRtdjS/pystan/var/folders/v0/fmjb82wj6tzbyvl5rpgkyt600000gp/T/tmpjRtdjS/pystan
gcc -fno-strict-aliasing -I//anaconda/include -arch x86_64 -DNDEBUG -g -fwrapv -O3 -Wall -Wstrict-prototypes -DBOOST_RESULT_OF_USE_TR1 -DBOOST_NO_DECLTYPE -DBOOST_DISABLE_ASSERTS -DEIGEN_NO_DEBUG -I/var/folders/v0/fmjb82wj6tzbyvl5rpgkyt600000gp/T/tmpjRtdjS/pystan -I//anaconda/lib/python2.7/site-packages/pystan -I//anaconda/lib/python2.7/site-packages/pystan/stan/src -I//anaconda/lib/python2.7/site-packages/pystan/stan/lib/stan_math_2.8.0 -I//anaconda/lib/python2.7/site-packages/pystan/stan/lib/stan_math_2.8.0/lib/eigen_3.2.4 -I//anaconda/lib/python2.7/site-packages/pystan/stan/lib/stan_math_2.8.0/lib/boost_1.58.0 -I//anaconda/lib/python2.7/site-packages/numpy/core/include -I//anaconda/include/python2.7 -c /var/folders/v0/fmjb82wj6tzbyvl5rpgkyt600000gp/T/tmpjRtdjS/pystan/stanfit4anon_model_9e5045edf8799b6e16462f3517df64f2_a72157596188ec1daae6f6858b197bc8.cpp -o /var/folders/v0/fmjb82wj6tzbyvl5rpgkyt600000gp/T/tmpjRtdjS/pystan/var/folders/v0/fmjb82wj6tzbyvl5rpgkyt600000gp/T/tmpjRtdjS/pystan/stanfit4anon_model_9e5045edf8799b6e16462f3517df64f2_a72157596188ec1daae6f6858b197bc8.o -O0 -ftemplate-depth-256 -Wno-unused-function -Wno-uninitialized
g++ -bundle -undefined dynamic_lookup -L//anaconda/lib -arch x86_64 -arch x86_64 /var/folders/v0/fmjb82wj6tzbyvl5rpgkyt600000gp/T/tmpjRtdjS/pystan/var/folders/v0/fmjb82wj6tzbyvl5rpgkyt600000gp/T/tmpjRtdjS/pystan/stanfit4anon_model_9e5045edf8799b6e16462f3517df64f2_a72157596188ec1daae6f6858b197bc8.o -L//anaconda/lib -o /var/folders/v0/fmjb82wj6tzbyvl5rpgkyt600000gp/T/tmpjRtdjS/pystan/stanfit4anon_model_9e5045edf8799b6e16462f3517df64f2_a72157596188ec1daae6f6858b197bc8.so
//anaconda/lib/python2.7/multiprocessing/queues.py:392: UserWarning: Pickling fit objects is an experimental feature!
The relevant StanModel instance must be pickled along with this fit object.
When unpickling the StanModel must be unpickled first.
  return send(obj)
//anaconda/lib/python2.7/multiprocessing/queues.py:392: UserWarning: Pickling fit objects is an experimental feature!
The relevant StanModel instance must be pickled along with this fit object.
When unpickling the StanModel must be unpickled first.
  return send(obj)
//anaconda/lib/python2.7/multiprocessing/queues.py:392: UserWarning: Pickling fit objects is an experimental feature!
The relevant StanModel instance must be pickled along with this fit object.
When unpickling the StanModel must be unpickled first.
  return send(obj)
//anaconda/lib/python2.7/multiprocessing/queues.py:392: UserWarning: Pickling fit objects is an experimental feature!
The relevant StanModel instance must be pickled along with this fit object.
When unpickling the StanModel must be unpickled first.
  return send(obj)
[-0.68000748 -0.75385242 -1.26590258 -0.64541516 -1.0469596  -1.46868077
 -0.67274615 -1.00623015 -0.39783321 -0.67916089 -0.92132741 -0.8969744
 -0.98029966 -0.73893268 -0.82678056 -1.28306831 -0.67664594 -1.09535881
 -0.49533788 -0.40730999 -1.35502688 -0.35376848 -0.81643857 -0.79558197
 -0.55123448 -0.40253118 -0.83203483 -1.542966   -0.9108598  -0.84630263
 -0.94499063 -0.80476288 -0.92239249 -0.96969388 -1.3178869  -0.80360112
 -0.2424555  -1.06614426 -0.84717531 -1.0651934  -1.01683801 -1.3890711
 -0.86952874 -0.92831297 -0.91925952 -1.30212707 -1.49542206 -1.0892086
 -0.8680106  -0.71214767 -1.30533379 -1.50471282 -1.37532898  0.1283055
 -0.62767677 -0.04436449 -1.51053084 -1.15063749 -1.25537695 -0.88882539
 -0.18814142 -1.42417073 -0.46526133 -0.7798144  -0.6081195   0.11655109
 -1.35117018 -1.5441118   0.45234498 -1.11556017 -0.22483225 -1.27898474
 -1.02729275 -1.30129528 -0.72664392 -1.58268297 -1.14204438 -0.85241702
  0.3946147  -0.57858245 -0.91372525 -1.30143117 -0.59026606 -1.44385761
 -0.13592646 -1.13460961 -0.19573153 -0.8698606  -0.13669311 -1.24440612
  0.44674547 -0.44944977 -1.32023584 -0.25128905 -1.58916074 -0.85849686
 -0.20416894 -0.75398236 -0.1122601  -0.85502016 -0.86853149 -1.51097897
 -1.01626219 -1.31629134 -1.28280611 -1.37718423 -1.38393365 -1.05604269
 -0.62527023 -1.22152791 -1.01227634 -0.98565246 -1.12318991 -1.01919252
 -1.10677458 -0.85664993 -1.13228028 -1.30157838 -1.29338378 -0.5403608
 -0.81122288 -1.22979685 -1.29585685 -0.71063734 -1.3503698  -1.06909684
 -1.10390433 -0.22192091 -1.08125012 -1.12427512 -0.71688546 -1.1085952
 -1.18346893 -1.08445223 -1.25832959 -1.19968247 -0.79021949 -0.99707397
 -1.28779107 -0.52770089 -1.42404185  0.01687237 -0.77292183 -1.1052459
 -1.11571549 -1.27154062 -1.11380834 -0.16187578 -1.23084779 -1.22574523
 -0.70518482 -0.49906001 -0.7615905  -0.41812537 -0.73317902 -1.68142585
 -0.81071626 -0.23488113 -0.34822746 -1.10727835 -0.42260074 -1.6546409
 -1.49261527 -0.90127569 -1.49261527 -1.66000934 -1.43324552 -1.65233464
 -1.85062283 -2.09576626 -0.44010007 -0.45888953 -0.77245503 -1.65601458
 -1.12686292 -1.04074301 -0.92326733 -0.97476455 -0.84675893 -1.80742848
 -1.40186706 -1.04094079 -0.70979814 -1.60096215 -0.41856511 -0.72610862
 -1.4562446  -1.0665244  -0.46706505 -1.38198375 -1.25760671 -1.0808657
 -1.59556154 -0.32189281 -0.79597663 -1.93255621 -1.59207192 -2.07910113
 -1.07180095 -1.38608077 -0.82031434 -0.30511276 -1.16681897 -0.8772653
 -0.93650205 -0.76794997 -1.15759309 -1.1121371  -0.71193485 -1.1868807
 -0.5005378  -0.47215389 -1.01214534 -0.46620318 -1.53544579 -0.55705495
 -1.05277741 -0.35849828 -1.14055328 -1.30270878 -0.74878048 -0.99471398
 -1.21451568 -1.13348879 -0.74782102 -0.57042687 -0.28628111 -0.69752296
 -0.81037884 -1.09036425 -0.94776831 -1.12827957 -0.77564526 -1.08667526
 -0.69415821 -1.01947868 -1.59297464 -1.00680347 -0.90627221 -0.4519196
 -0.90744999 -1.73242357 -1.09805708 -0.57042687 -0.5683554  -1.17054105
 -0.5330147  -1.13867337 -1.05079947 -1.09936987]
[-1.20092628 -1.23484342 -1.0553343  -1.21413307 -1.1528624  -1.00247398
 -1.13331387 -1.24345616 -1.21570161 -1.13095937 -1.23368786 -1.11386924
 -1.26298557 -1.21795743 -1.23656883 -1.03312909 -1.09512173 -1.10083257
 -1.09894226 -1.18897223 -1.24783055 -1.16689321 -1.21972399 -1.09167573
 -1.21296507 -1.21821422 -1.20199536 -1.16780663 -1.15706038 -1.2195977
 -1.08218701 -1.29842305 -1.13823666 -1.11599108 -1.34050913 -1.11299431
 -1.2323085  -1.21540174 -1.13068884 -1.18872004 -1.18055062 -1.30965146
 -1.09032991 -1.16132574 -1.32109129 -0.98795414 -1.00907567 -1.19272873
 -1.1243346  -1.12361037 -1.10821663 -1.17738443 -1.22455557 -1.03310741
 -1.03597123 -1.24436703 -1.17917602 -1.18928852 -1.15233162 -1.18152638
 -1.22713112 -1.06361461 -1.1475578  -1.04564411 -0.98622308 -1.086325
 -1.22917036 -1.12602747 -1.15349318 -1.22585288 -1.3127462  -1.29455011
 -1.08991444 -1.13774731 -1.27981386 -1.2008288  -1.08321584 -1.23440782
 -1.19507477 -1.19750999 -1.15246855 -1.2044199  -1.16520623 -1.10892883
 -1.22120122 -1.16012457 -1.10525201 -1.18396092 -1.25360111 -1.10211965
 -1.15104914 -1.12522267 -1.2322023  -1.1014059  -1.19116059 -1.18822333
 -1.145244   -1.0364668  -1.28027864 -1.10860701 -1.21225012 -1.1928545
 -1.26509939 -1.25829871 -0.95756949 -1.27551775 -1.27539814 -1.14053721
 -1.23456153 -1.01055271 -1.0770077  -1.11709501 -1.2726455  -1.16908333
 -1.19390213 -1.08788532 -1.17159992 -1.30830324 -1.19195912 -1.12764486
 -1.30260041 -1.30193192 -1.27267262 -1.27678998 -1.13640839 -1.13440044
 -1.16164917 -1.07805082 -1.12276964 -1.30776689 -1.165147   -1.25262209
 -1.27523623 -1.38605015 -1.25347327 -1.20418529 -1.22878436 -1.16261971
 -1.25643285 -1.24557851 -1.3892268  -1.07847293 -1.22336656 -1.29170231
 -1.27180404 -1.24292159 -1.15145668 -1.15623749 -1.07398888 -1.28773654
 -1.08794301 -1.32779081 -1.08486844 -1.18468844 -1.22618301 -1.12779714
 -1.11032419 -1.11917    -1.21915602 -1.10139078 -1.11648006 -1.22715693
 -1.3198319  -1.20239781 -1.3198319  -1.11661024 -1.1384323  -1.26116302
 -1.15684649 -1.20413546 -1.20688997 -1.22399145 -1.26585195 -1.13601548
 -1.21827423 -1.1982366  -1.28075756 -1.16285525 -1.18920911 -1.19930631
 -1.15716772 -1.23542501 -1.23852793 -1.26337771 -1.16599833 -1.12988091
 -1.21621193 -1.09729244 -1.17935524 -1.32860811 -1.32526828 -1.11137244
 -1.14289986 -1.2116045  -1.14760806 -1.23793989 -1.24630911 -1.2607567
 -1.11308524 -1.17064721 -1.19952055 -1.0096089  -0.99786084 -1.20965535
 -1.34817507 -1.02277358 -0.99683187 -1.06760198 -1.23102129 -1.14638215
 -1.16258498 -1.27908825 -1.02180321 -1.27318833 -1.19845001 -1.33538316
 -1.15537536 -1.17682787 -0.94998105 -1.13873793 -1.20379094 -1.16436511
 -1.1921816  -0.98805246 -1.08919005 -1.14850926 -1.01728872 -1.14529813
 -1.14725873 -1.19511491 -1.15144866 -1.1942263  -1.25040596 -1.18167979
 -1.28739293 -1.07590459 -1.11958914 -1.30166296 -1.07110617 -1.13555696
 -1.08229934 -1.25057443 -1.18613849 -1.14850926 -1.1848607  -1.14147193
 -1.18121657 -1.07562931 -1.05984524 -1.05035022]
[ 0.73086913  0.64810963  0.56878067  0.72041081  0.5937074   0.65193444
  0.55219517  0.85356084  0.73727791  0.59120398  0.76451307  0.7537889
  0.88653464  0.64159234  0.76994339  0.56192836  0.5987182   0.50117929
  0.40987228  0.63661995  1.03836848  0.54520465  0.87351098  0.59630923
  0.72589104  0.75621683  0.77082884  0.8805167   0.65069078  0.74122479
  0.65679031  0.83154767  0.54314094  0.78305424  0.97438479  0.6153257
  0.49965195  0.89045414  0.55709116  1.01770818  0.6209848   0.99172581
  0.67508393  0.5023004   0.87552588  0.56328046  0.69670449  0.84520896
  0.65940491  0.5846734   0.79611668  0.86700505  0.84606381  0.26044636
  0.48182564  0.58429466  0.81706642  0.60401714  0.74775775  0.8266128
  0.61009965  0.83138507  0.60681409  0.54522293  0.436969    0.22657292
  0.86938165  0.84089776  0.36830528  0.5773203   0.69488115  1.05925651
  0.8320463   0.78362558  0.77273369  1.02741521  0.69205713  0.85065279
  0.36320177  0.54460597  0.78738525  0.88181797  0.41735064  0.88010211
  0.69670233  0.86753463  0.47016681  0.8580018   0.6444744   0.70860396
  0.33656242  0.61561393  0.91664626  0.45728012  0.84999163  0.59417025
  0.33340564  0.58278662  0.60877127  0.66585106  0.71734343  0.83738293
  0.85411359  0.82125968  0.4371124   0.9900744   1.15274213  0.79423327
  0.76465043  0.49059322  0.7634129   0.77844011  0.86939105  0.54338167
  0.90698032  0.74802361  0.76617492  0.84271752  1.01971943  0.62284163
  0.93625723  0.85260307  0.82841856  0.85326376  0.80779797  0.78200198
  0.84410814  0.41123592  0.79151351  1.00514863  0.64691362  1.00585668
  0.89123638  1.17852515  1.10898668  0.78025156  0.75417698  0.82719382
  0.98702615  0.87535245  1.2105958   0.4165459   0.78819116  1.16130393
  0.81168032  1.05782349  0.87799436  0.39733827  0.88862384  1.03750381
  0.52522008  0.89050472  0.60446892  0.46849854  0.80359401  0.94087793
  0.65562887  0.52676982  0.50178843  0.62441533  0.39872294  0.88555622
  1.0322837   0.84859994  1.0322837   1.00553182  0.99372765  0.90751845
  1.05761421  1.25549017  0.57663827  0.57905584  0.7152327   0.93120308
  0.84834005  0.85763917  0.91152795  0.67203839  0.77912509  0.86080019
  0.98484949  0.88356495  0.90452753  1.02246438  0.53221169  0.52767401
  0.83932087  0.72644958  0.6710936   1.05176971  1.00405497  0.72289568
  1.07151644  0.50541191  0.67577905  1.18518612  0.99680704  1.16392687
  0.69165574  0.87459053  0.71220881  0.41719332  0.62246208  0.76229497
  0.98023645  0.50823119  0.61535978  0.46568292  0.71991862  0.77408198
  0.68246394  0.88060749  0.64440348  0.86577391  1.01276028  0.99007515
  0.69913173  0.68588993  0.46885892  0.806739    0.62003297  0.65033595
  0.82628709  0.41255323  0.55611729  0.56142137  0.44124459  0.71009267
  0.79050294  0.94142907  0.56614859  0.85905918  0.90660957  0.71862665
  0.82473775  0.64360147  0.90631321  0.88499765  0.65526825  0.59065397
  0.67881442  1.04132942  0.92011807  0.56142137  0.47937445  0.78777449
  0.6891703   0.49264677  0.64878596  0.383757  ]
[[[ -9.21327412e-01  -1.23368786e+00   7.64513073e-01  -5.29748681e+02]
  [ -7.26643915e-01  -1.27981386e+00   7.72733690e-01  -5.30441981e+02]
  [ -1.18346893e+00  -1.27523623e+00   8.91236384e-01  -5.30283233e+02]
  [ -1.10727835e+00  -1.10139078e+00   6.24415335e-01  -5.30229250e+02]
  [ -8.20314342e-01  -1.19952055e+00   7.12208807e-01  -5.29545742e+02]]

 [[ -8.16438571e-01  -1.21972399e+00   8.73510982e-01  -5.30185030e+02]
  [ -1.30129528e+00  -1.13774731e+00   7.83625584e-01  -5.30042352e+02]
  [ -1.22574523e+00  -1.28773654e+00   1.03750381e+00  -5.30300941e+02]
  [ -9.74764550e-01  -1.16285525e+00   6.72038389e-01  -5.29686810e+02]
  [ -8.10378843e-01  -1.14725873e+00   7.90502937e-01  -5.30436156e+02]]

 [[ -8.04762875e-01  -1.29842305e+00   8.31547668e-01  -5.30514419e+02]
  [ -1.88141417e-01  -1.22713112e+00   6.10099650e-01  -5.31553698e+02]
  [ -1.08445223e+00  -1.38605015e+00   1.17852515e+00  -5.32406376e+02]
  [ -1.06652440e+00  -1.09729244e+00   7.26449584e-01  -5.30087700e+02]
  [ -5.68355403e-01  -1.18486070e+00   4.79374451e-01  -5.31587820e+02]]

 [[ -1.01683801e+00  -1.18055062e+00   6.20984803e-01  -5.30746572e+02]
  [ -7.53982363e-01  -1.03646680e+00   5.82786622e-01  -5.31069666e+02]
  [ -1.10524590e+00  -1.29170231e+00   1.16130393e+00  -5.32603964e+02]
  [ -1.08086570e+00  -1.11137244e+00   7.22895679e-01  -5.29846987e+02]
  [ -2.86281106e-01  -1.01728872e+00   4.41244585e-01  -5.32511626e+02]]

 [[ -9.69693875e-01  -1.11599108e+00   7.83054235e-01  -5.30426809e+02]
  [ -8.52417021e-01  -1.23440782e+00   8.50652790e-01  -5.29741905e+02]
  [ -1.11571549e+00  -1.27180404e+00   8.11680321e-01  -5.30766412e+02]
  [ -9.01275688e-01  -1.20239781e+00   8.48599940e-01  -5.29781748e+02]
  [ -3.05112761e-01  -1.00960890e+00   4.17193321e-01  -5.32406221e+02]]

 [[ -9.22392486e-01  -1.13823666e+00   5.43140942e-01  -5.30725470e+02]
  [ -8.88825389e-01  -1.18152638e+00   8.26612802e-01  -5.29874678e+02]
  [ -1.10390433e+00  -1.16164917e+00   8.44108143e-01  -5.29798742e+02]
  [ -1.12686292e+00  -1.21827423e+00   8.48340055e-01  -5.29549761e+02]
  [ -6.94158212e-01  -1.28739293e+00   8.24737749e-01  -5.30415929e+02]]

 [[ -8.68010595e-01  -1.12433460e+00   6.59404908e-01  -5.29632713e+02]
  [ -6.27676767e-01  -1.03597123e+00   4.81825640e-01  -5.30893499e+02]
  [ -1.13228028e+00  -1.17159992e+00   7.66174924e-01  -5.29584267e+02]
  [ -1.04094079e+00  -1.23542501e+00   8.83564947e-01  -5.29588805e+02]
  [ -3.58498281e-01  -1.17682787e+00   6.85889925e-01  -5.31156268e+02]]

 [[ -9.10859795e-01  -1.15706038e+00   6.50690785e-01  -5.29665628e+02]
  [ -7.79814400e-01  -1.04564411e+00   5.45222927e-01  -5.30580576e+02]
  [ -1.22152791e+00  -1.01055271e+00   4.90593222e-01  -5.32397943e+02]
  [ -1.07180095e+00  -1.11308524e+00   6.91655744e-01  -5.29794112e+02]
  [ -4.51919600e-01  -1.13555696e+00   5.90653969e-01  -5.30378399e+02]]

 [[ -9.44990627e-01  -1.08218701e+00   6.56790314e-01  -5.30104679e+02]
  [ -1.14204438e+00  -1.08321584e+00   6.92057131e-01  -5.30218635e+02]
  [ -8.56649930e-01  -1.08788532e+00   7.48023612e-01  -5.31262022e+02]
  [ -4.18125365e-01  -1.18468844e+00   4.68498536e-01  -5.31535100e+02]
  [ -4.72153893e-01  -1.27908825e+00   8.80607489e-01  -5.31802558e+02]]

 [[ -1.09535881e+00  -1.10083257e+00   5.01179292e-01  -5.31903034e+02]
  [ -1.15063749e+00  -1.18928852e+00   6.04017140e-01  -5.32271462e+02]
  [ -1.42404185e+00  -1.38922680e+00   1.21059580e+00  -5.32457528e+02]
  [ -1.49261527e+00  -1.31983190e+00   1.03228370e+00  -5.31768691e+02]
  [ -4.66203177e-01  -1.27318833e+00   8.65773913e-01  -5.31672436e+02]]

 [[ -1.31788690e+00  -1.34050913e+00   9.74384795e-01  -5.32063786e+02]
  [ -1.11556017e+00  -1.22585288e+00   5.77320301e-01  -5.34112347e+02]
  [ -1.38393365e+00  -1.27539814e+00   1.15274213e+00  -5.31345295e+02]
  [ -1.49261527e+00  -1.31983190e+00   1.03228370e+00  -5.31768691e+02]
  [ -5.57054955e-01  -1.33538316e+00   9.90075146e-01  -5.32600234e+02]]

 [[ -1.35502688e+00  -1.24783055e+00   1.03836848e+00  -5.30309533e+02]
  [ -1.02729275e+00  -1.08991444e+00   8.32046299e-01  -5.31927202e+02]
  [ -1.51097897e+00  -1.19285450e+00   8.37382929e-01  -5.31177946e+02]
  [ -1.59556154e+00  -1.14289986e+00   1.07151644e+00  -5.32780439e+02]
  [ -7.47821025e-01  -1.08919005e+00   5.56117291e-01  -5.30041896e+02]]

 [[ -6.79160890e-01  -1.13095937e+00   5.91203979e-01  -5.29805305e+02]
  [ -1.58916074e+00  -1.19116059e+00   8.49991634e-01  -5.31660577e+02]
  [ -1.25832959e+00  -1.25347327e+00   1.10898668e+00  -5.31384425e+02]
  [ -1.80742848e+00  -1.19930631e+00   8.60800195e-01  -5.34138313e+02]
  [ -9.94713982e-01  -1.16436511e+00   6.50335948e-01  -5.29960066e+02]]

 [[ -7.12147670e-01  -1.12361037e+00   5.84673400e-01  -5.29812932e+02]
  [ -1.51053084e+00  -1.17917602e+00   8.17066424e-01  -5.31217919e+02]
  [ -1.28280611e+00  -9.57569488e-01   4.37112404e-01  -5.33889973e+02]
  [ -1.68142585e+00  -1.12779714e+00   9.40877928e-01  -5.31726721e+02]
  [ -5.00537799e-01  -1.16258498e+00   6.82463939e-01  -5.30451395e+02]]

 [[ -1.06614426e+00  -1.21540174e+00   8.90454138e-01  -5.29630517e+02]
  [ -1.44385761e+00  -1.10892883e+00   8.80102106e-01  -5.31040013e+02]
  [ -1.12427512e+00  -1.30776689e+00   1.00514863e+00  -5.30370868e+02]
  [ -1.65601458e+00  -1.13601548e+00   9.31203078e-01  -5.31476692e+02]
  [ -5.33014702e-01  -1.18121657e+00   6.89170305e-01  -5.30147207e+02]]

 [[ -1.26590258e+00  -1.05533430e+00   5.68780671e-01  -5.31609066e+02]
  [ -1.37532898e+00  -1.22455557e+00   8.46063814e-01  -5.30646321e+02]
  [ -1.27154062e+00  -1.24292159e+00   1.05782349e+00  -5.30664628e+02]
  [ -1.65464090e+00  -1.22715693e+00   8.85556216e-01  -5.32563789e+02]
  [ -1.05277741e+00  -1.15537536e+00   6.99131732e-01  -5.29654576e+02]]

 [[ -1.28306831e+00  -1.03312909e+00   5.61928358e-01  -5.31846304e+02]
  [ -1.24440612e+00  -1.10211965e+00   7.08603962e-01  -5.30246258e+02]
  [ -1.28779107e+00  -1.25643285e+00   9.87026148e-01  -5.30011303e+02]
  [ -1.65233464e+00  -1.26116302e+00   9.07518449e-01  -5.33198301e+02]
  [ -1.01947868e+00  -1.07590459e+00   6.43601466e-01  -5.30153580e+02]]

 [[ -9.28312967e-01  -1.16132574e+00   5.02300395e-01  -5.32083616e+02]
  [ -1.30533379e+00  -1.10821663e+00   7.96116680e-01  -5.30347807e+02]
  [ -1.01626219e+00  -1.26509939e+00   8.54113588e-01  -5.29890627e+02]
  [ -1.59207192e+00  -1.24630911e+00   9.96807039e-01  -5.31083837e+02]
  [ -9.36502054e-01  -1.34817507e+00   9.80236453e-01  -5.31149769e+02]]

 [[ -1.04695960e+00  -1.15286240e+00   5.93707400e-01  -5.30800355e+02]
  [ -1.25537695e+00  -1.15233162e+00   7.47757750e-01  -5.30035229e+02]
  [ -7.10637339e-01  -1.27678998e+00   8.53263761e-01  -5.30284825e+02]
  [ -1.85062283e+00  -1.15684649e+00   1.05761421e+00  -5.32628105e+02]
  [ -1.08667526e+00  -1.18167979e+00   7.18626652e-01  -5.29821350e+02]]

 [[ -6.76645942e-01  -1.09512173e+00   5.98718203e-01  -5.30139134e+02]
  [ -1.30143117e+00  -1.20441990e+00   8.81817971e-01  -5.29818809e+02]
  [ -8.11222882e-01  -1.30260041e+00   9.36257226e-01  -5.30554968e+02]
  [ -1.93255621e+00  -1.23793989e+00   1.18518612e+00  -5.32905318e+02]
  [ -1.53544579e+00  -1.19845001e+00   1.01276028e+00  -5.30706817e+02]]

 [[ -9.19259519e-01  -1.32109129e+00   8.75525876e-01  -5.30917296e+02]
  [ -1.54411180e+00  -1.12602747e+00   8.40897758e-01  -5.31031032e+02]
  [ -7.72921834e-01  -1.22336656e+00   7.88191159e-01  -5.29667701e+02]
  [ -2.07910113e+00  -1.26075670e+00   1.16392687e+00  -5.34438856e+02]
  [ -1.30270878e+00  -1.13873793e+00   8.06738997e-01  -5.30014959e+02]]

 [[ -8.96974400e-01  -1.11386924e+00   7.53788897e-01  -5.30369168e+02]
  [ -1.13460961e+00  -1.16012457e+00   8.67534628e-01  -5.29991819e+02]
  [ -1.12318991e+00  -1.27264550e+00   8.69391051e-01  -5.30192089e+02]
  [ -2.09576626e+00  -1.20413546e+00   1.25549017e+00  -5.34688404e+02]
  [ -5.70426867e-01  -1.14850926e+00   5.61421368e-01  -5.30037847e+02]]

 [[ -1.38907110e+00  -1.30965146e+00   9.91725806e-01  -5.31230580e+02]
  [ -1.42417073e+00  -1.06361461e+00   8.31385073e-01  -5.31750537e+02]
  [ -1.10859520e+00  -1.25262209e+00   1.00585668e+00  -5.30313170e+02]
  [ -1.38198375e+00  -1.32860811e+00   1.05176971e+00  -5.31237265e+02]
  [ -5.70426867e-01  -1.14850926e+00   5.61421368e-01  -5.30037847e+02]]

 [[ -7.95581966e-01  -1.09167573e+00   5.96309231e-01  -5.29957421e+02]
  [ -1.35117018e+00  -1.22917036e+00   8.69381646e-01  -5.30345052e+02]
  [ -1.23084779e+00  -1.07398888e+00   8.88623841e-01  -5.32837468e+02]
  [ -1.43324552e+00  -1.13843230e+00   9.93727651e-01  -5.31752350e+02]
  [ -7.67949968e-01  -1.02277358e+00   5.08231191e-01  -5.30956314e+02]]

 [[ -8.03601116e-01  -1.11299431e+00   6.15325700e-01  -5.29755286e+02]
  [ -1.32023584e+00  -1.23220230e+00   9.16646260e-01  -5.29960563e+02]
  [ -1.30157838e+00  -1.30830324e+00   8.42717523e-01  -5.32904585e+02]
  [ -1.40186706e+00  -1.15716772e+00   9.84849494e-01  -5.31103588e+02]
  [ -1.11213710e+00  -1.06760198e+00   4.65682916e-01  -5.32264994e+02]]

 [[ -7.53852421e-01  -1.23484342e+00   6.48109629e-01  -5.30646230e+02]
  [ -1.27898474e+00  -1.29455011e+00   1.05925651e+00  -5.30444613e+02]
  [ -1.29585685e+00  -1.27267262e+00   8.28418559e-01  -5.31671283e+02]
  [ -1.45624460e+00  -1.21621193e+00   8.39320869e-01  -5.31149059e+02]
  [ -1.16681897e+00  -9.97860837e-01   6.22462078e-01  -5.31841307e+02]]

 [[ -7.38932683e-01  -1.21795743e+00   6.41592343e-01  -5.30289686e+02]
  [ -8.58496864e-01  -1.18822333e+00   5.94170251e-01  -5.30619122e+02]
  [ -1.31629134e+00  -1.25829871e+00   8.21259679e-01  -5.31462867e+02]
  [ -1.60096215e+00  -1.26337771e+00   1.02246438e+00  -5.31237606e+02]
  [ -1.15759309e+00  -9.96831868e-01   6.15359775e-01  -5.31827952e+02]]

 [[ -8.69528741e-01  -1.09032991e+00   6.75083931e-01  -5.30189355e+02]
  [ -5.90266056e-01  -1.16520623e+00   4.17350637e-01  -5.32376017e+02]
  [ -1.29338378e+00  -1.19195912e+00   1.01971943e+00  -5.31029734e+02]
  [ -1.66000934e+00  -1.11661024e+00   1.00553182e+00  -5.32473320e+02]
  [ -1.13348879e+00  -9.88052463e-01   4.12553235e-01  -5.33011969e+02]]

 [[ -8.26780564e-01  -1.23656883e+00   7.69943390e-01  -5.29698259e+02]
  [ -4.65261330e-01  -1.14755780e+00   6.06814092e-01  -5.30282218e+02]
  [ -1.37718423e+00  -1.27551775e+00   9.90074396e-01  -5.30430077e+02]
  [ -1.38608077e+00  -1.17064721e+00   8.74590527e-01  -5.30086925e+02]
  [ -1.14055328e+00  -9.49981054e-01   4.68858916e-01  -5.32909046e+02]]

 [[ -8.46302628e-01  -1.21959770e+00   7.41224790e-01  -5.29619300e+02]
  [ -5.78582452e-01  -1.19750999e+00   5.44605972e-01  -5.30850208e+02]
  [ -1.35036980e+00  -1.13640839e+00   8.07797973e-01  -5.30173477e+02]
  [ -8.46758930e-01  -1.18920911e+00   7.79125091e-01  -5.29545956e+02]
  [ -9.47768313e-01  -1.15144866e+00   5.66148593e-01  -5.30702672e+02]]

 [[ -2.42455499e-01  -1.23230850e+00   4.99651949e-01  -5.32357352e+02]
  [ -6.08119502e-01  -9.86223076e-01   4.36968996e-01  -5.31895343e+02]
  [ -1.05604269e+00  -1.14053721e+00   7.94233274e-01  -5.29801097e+02]
  [ -1.25760671e+00  -1.32526828e+00   1.00405497e+00  -5.30957579e+02]
  [ -1.01214534e+00  -1.02180321e+00   6.44403482e-01  -5.31404404e+02]]

 [[ -4.07309994e-01  -1.18897223e+00   6.36619953e-01  -5.30434182e+02]
  [ -1.35926455e-01  -1.22120122e+00   6.96702331e-01  -5.32423813e+02]
  [ -7.16885462e-01  -1.16514700e+00   6.46913619e-01  -5.29628413e+02]
  [ -4.22600741e-01  -1.11648006e+00   3.98722941e-01  -5.31349311e+02]
  [ -1.05079947e+00  -1.05984524e+00   6.48785956e-01  -5.30442719e+02]]

 [[ -4.02531178e-01  -1.21821422e+00   7.56216833e-01  -5.31165119e+02]
  [ -2.24832252e-01  -1.31274620e+00   6.94881148e-01  -5.32710973e+02]
  [ -6.25270235e-01  -1.23456153e+00   7.64650434e-01  -5.30029374e+02]
  [ -2.34881133e-01  -1.11917000e+00   5.26769822e-01  -5.31328220e+02]
  [ -1.09936987e+00  -1.05035022e+00   3.83757004e-01  -5.33752187e+02]]

 [[ -3.97833214e-01  -1.21570161e+00   7.37277908e-01  -5.31004941e+02]
  [ -1.12260102e-01  -1.28027864e+00   6.08771274e-01  -5.32807815e+02]
  [ -5.40360802e-01  -1.12764486e+00   6.22841630e-01  -5.30271949e+02]
  [ -7.05184821e-01  -1.08794301e+00   5.25220083e-01  -5.30158196e+02]
  [ -1.13867337e+00  -1.07562931e+00   4.92646768e-01  -5.32035645e+02]]

 [[ -4.95337877e-01  -1.09894226e+00   4.09872277e-01  -5.31066264e+02]
  [ -1.95731534e-01  -1.10525201e+00   4.70166807e-01  -5.31467530e+02]
  [ -7.90219492e-01  -1.22878436e+00   7.54176981e-01  -5.29673980e+02]
  [ -7.95976628e-01  -1.14760806e+00   6.75779053e-01  -5.29559014e+02]
  [ -9.07449989e-01  -1.08229934e+00   6.78814415e-01  -5.30303998e+02]]

 [[ -6.80007479e-01  -1.20092628e+00   7.30869131e-01  -5.29728324e+02]
  [ -4.43644934e-02  -1.24436703e+00   5.84294658e-01  -5.32516526e+02]
  [ -2.21920910e-01  -1.07805082e+00   4.11235923e-01  -5.31561139e+02]
  [ -7.33179015e-01  -1.22618301e+00   8.03594005e-01  -5.29840043e+02]
  [ -1.09036425e+00  -1.19511491e+00   9.41429068e-01  -5.30377743e+02]]

 [[ -6.45415163e-01  -1.21413307e+00   7.20410814e-01  -5.29818642e+02]
  [ -1.36693108e-01  -1.25360111e+00   6.44474400e-01  -5.32088326e+02]
  [ -1.61875785e-01  -1.15623749e+00   3.97338270e-01  -5.32129003e+02]
  [ -7.72455028e-01  -1.26585195e+00   7.15232702e-01  -5.30670267e+02]
  [ -1.09805708e+00  -1.18613849e+00   9.20118065e-01  -5.30223891e+02]]

 [[ -3.53768479e-01  -1.16689321e+00   5.45204651e-01  -5.30621535e+02]
  [  1.28305500e-01  -1.03310741e+00   2.60446358e-01  -5.33928420e+02]
  [  1.68723718e-02  -1.07847293e+00   4.16545898e-01  -5.32889980e+02]
  [ -8.10716263e-01  -1.11032419e+00   6.55628874e-01  -5.29847607e+02]
  [ -1.12827957e+00  -1.19422630e+00   8.59059183e-01  -5.29548900e+02]]

 [[ -5.51234478e-01  -1.21296507e+00   7.25891044e-01  -5.30129800e+02]
  [  1.16551093e-01  -1.08632500e+00   2.26572920e-01  -5.34141885e+02]
  [ -5.27700887e-01  -1.24557851e+00   8.75352451e-01  -5.31687420e+02]
  [ -1.04074301e+00  -1.19823660e+00   8.57639172e-01  -5.29573344e+02]
  [ -8.77265300e-01  -1.20965535e+00   7.62294968e-01  -5.29475923e+02]]

 [[ -1.00623015e+00  -1.24345616e+00   8.53560841e-01  -5.29607700e+02]
  [  3.94614698e-01  -1.19507477e+00   3.63201768e-01  -5.35406524e+02]
  [ -8.68531495e-01  -1.21225012e+00   7.17343432e-01  -5.29682218e+02]
  [ -9.23267329e-01  -1.28075756e+00   9.11527947e-01  -5.30018435e+02]
  [ -1.73242357e+00  -1.25057443e+00   1.04132942e+00  -5.31865130e+02]]

 [[ -8.32034827e-01  -1.20199536e+00   7.70828840e-01  -5.29501515e+02]
  [  4.52344982e-01  -1.15349318e+00   3.68305276e-01  -5.35403773e+02]
  [ -9.85652455e-01  -1.11709501e+00   7.78440109e-01  -5.30284527e+02]
  [ -7.09798139e-01  -1.23852793e+00   9.04527529e-01  -5.30977699e+02]
  [ -1.59297464e+00  -1.11958914e+00   9.06313207e-01  -5.31358493e+02]]

 [[ -1.30212707e+00  -9.87954144e-01   5.63280462e-01  -5.32423960e+02]
  [  4.46745472e-01  -1.15104914e+00   3.36562416e-01  -5.35412325e+02]
  [ -1.01227634e+00  -1.07700770e+00   7.63412903e-01  -5.31158010e+02]
  [ -4.99060007e-01  -1.32779081e+00   8.90504722e-01  -5.31847790e+02]
  [ -7.11934852e-01  -1.23102129e+00   7.19918624e-01  -5.29865681e+02]]

 [[ -1.54296600e+00  -1.16780663e+00   8.80516704e-01  -5.30799235e+02]
  [ -2.51289055e-01  -1.10140590e+00   4.57280125e-01  -5.31230045e+02]
  [ -1.01919252e+00  -1.16908333e+00   5.43381668e-01  -5.32002869e+02]
  [ -7.26108623e-01  -1.12988091e+00   5.27674006e-01  -5.30222861e+02]
  [ -7.75645263e-01  -1.25040596e+00   9.06609574e-01  -5.30466118e+02]]

 [[ -1.49542206e+00  -1.00907567e+00   6.96704491e-01  -5.32609249e+02]
  [ -2.04168944e-01  -1.14524400e+00   3.33405638e-01  -5.32865529e+02]
  [ -1.22979685e+00  -1.30193192e+00   8.52603066e-01  -5.31863514e+02]
  [ -4.18565111e-01  -1.16599833e+00   5.32211687e-01  -5.30523482e+02]
  [ -1.17054105e+00  -1.14147193e+00   7.87774489e-01  -5.29732568e+02]]

 [[ -1.46868077e+00  -1.00247398e+00   6.51934444e-01  -5.32669815e+02]
  [ -8.69860598e-01  -1.18396092e+00   8.58001804e-01  -5.30315093e+02]
  [ -1.11380834e+00  -1.15145668e+00   8.77994365e-01  -5.30352879e+02]
  [ -4.67065053e-01  -1.17935524e+00   6.71093602e-01  -5.30353437e+02]
  [ -1.18688070e+00  -1.14638215e+00   7.74081976e-01  -5.29711751e+02]]

 [[ -1.06519340e+00  -1.18872004e+00   1.01770818e+00  -5.32367103e+02]
  [ -9.13725253e-01  -1.15246855e+00   7.87385249e-01  -5.29861320e+02]
  [ -1.10677458e+00  -1.19390213e+00   9.06980324e-01  -5.29908107e+02]
  [ -4.40100066e-01  -1.20688997e+00   5.76638270e-01  -5.30668931e+02]
  [ -1.21451568e+00  -1.19218160e+00   8.26287087e-01  -5.29661867e+02]]

 [[ -9.80299663e-01  -1.26298557e+00   8.86534637e-01  -5.29771793e+02]
  [ -8.55020162e-01  -1.10860701e+00   6.65851056e-01  -5.29836536e+02]
  [ -1.19968247e+00  -1.20418529e+00   7.80251556e-01  -5.30011838e+02]
  [ -4.58889534e-01  -1.22399145e+00   5.79055841e-01  -5.30961160e+02]
  [ -6.97522962e-01  -1.14529813e+00   7.10092674e-01  -5.30027769e+02]]

 [[ -6.72746145e-01  -1.13331387e+00   5.52195173e-01  -5.29985418e+02]
  [ -4.49449772e-01  -1.12522267e+00   6.15613931e-01  -5.30696954e+02]
  [ -1.08125012e+00  -1.12276964e+00   7.91513514e-01  -5.30064593e+02]
  [ -3.48227459e-01  -1.21915602e+00   5.01788432e-01  -5.31891136e+02]
  [ -7.48780483e-01  -1.20379094e+00   6.20032973e-01  -5.30270027e+02]]

 [[ -8.47175307e-01  -1.13068884e+00   5.57091163e-01  -5.30168251e+02]
  [ -1.50471282e+00  -1.17738443e+00   8.67005051e-01  -5.30682233e+02]
  [ -1.06909684e+00  -1.13440044e+00   7.82001981e-01  -5.29796652e+02]
  [ -3.21892814e-01  -1.21160450e+00   5.05411914e-01  -5.31651139e+02]
  [ -1.00680347e+00  -1.30166296e+00   8.84997651e-01  -5.30461230e+02]]

 [[ -1.08920860e+00  -1.19272873e+00   8.45208960e-01  -5.29509858e+02]
  [ -1.58268297e+00  -1.20082880e+00   1.02741521e+00  -5.30887178e+02]
  [ -9.97073967e-01  -1.16261971e+00   8.27193819e-01  -5.29845396e+02]
  [ -7.61590503e-01  -1.08486844e+00   6.04468923e-01  -5.30120463e+02]
  [ -9.06272206e-01  -1.07110617e+00   6.55268254e-01  -5.30376732e+02]]]
Inference for Stan model: anon_model_9e5045edf8799b6e16462f3517df64f2.
5 chains, each with iter=100; warmup=50; thin=1; 
post-warmup draws per chain=50, total post-warmup draws=250.

        mean se_mean     sd   2.5%    25%    50%    75%  97.5%  n_eff   Rhat
alpha  -0.94    0.03   0.44   -1.7  -1.24  -0.99   -0.7  -0.02  250.0   1.25
beta   -1.18  5.3e-3   0.08  -1.33  -1.23  -1.18  -1.12   -1.0  250.0   1.04
lnf0    0.74    0.01    0.2   0.38    0.6   0.76   0.88   1.16  250.0   1.14
lp__  -530.9    0.08   1.23 -534.1 -531.7 -530.6 -530.0 -529.5  250.0   1.06

Samples were drawn using NUTS(diag_e) at Thu Aug 18 15:59:19 2016.
For each parameter, n_eff is a crude measure of effective sample size,
and Rhat is the potential scale reduction factor on split chains (at 
convergence, Rhat=1).
[-0.68000748 -0.75385242 -1.26590258 -0.64541516 -1.0469596  -1.46868077
 -0.67274615 -1.00623015 -0.39783321 -0.67916089 -0.92132741 -0.8969744
 -0.98029966 -0.73893268 -0.82678056 -1.28306831 -0.67664594 -1.09535881
 -0.49533788 -0.40730999 -1.35502688 -0.35376848 -0.81643857 -0.79558197
 -0.55123448 -0.40253118 -0.83203483 -1.542966   -0.9108598  -0.84630263
 -0.94499063 -0.80476288 -0.92239249 -0.96969388 -1.3178869  -0.80360112
 -0.2424555  -1.06614426 -0.84717531 -1.0651934  -1.01683801 -1.3890711
 -0.86952874 -0.92831297 -0.91925952 -1.30212707 -1.49542206 -1.0892086
 -0.8680106  -0.71214767 -1.30533379 -1.50471282 -1.37532898  0.1283055
 -0.62767677 -0.04436449 -1.51053084 -1.15063749 -1.25537695 -0.88882539
 -0.18814142 -1.42417073 -0.46526133 -0.7798144  -0.6081195   0.11655109
 -1.35117018 -1.5441118   0.45234498 -1.11556017 -0.22483225 -1.27898474
 -1.02729275 -1.30129528 -0.72664392 -1.58268297 -1.14204438 -0.85241702
  0.3946147  -0.57858245 -0.91372525 -1.30143117 -0.59026606 -1.44385761
 -0.13592646 -1.13460961 -0.19573153 -0.8698606  -0.13669311 -1.24440612
  0.44674547 -0.44944977 -1.32023584 -0.25128905 -1.58916074 -0.85849686
 -0.20416894 -0.75398236 -0.1122601  -0.85502016 -0.86853149 -1.51097897
 -1.01626219 -1.31629134 -1.28280611 -1.37718423 -1.38393365 -1.05604269
 -0.62527023 -1.22152791 -1.01227634 -0.98565246 -1.12318991 -1.01919252
 -1.10677458 -0.85664993 -1.13228028 -1.30157838 -1.29338378 -0.5403608
 -0.81122288 -1.22979685 -1.29585685 -0.71063734 -1.3503698  -1.06909684
 -1.10390433 -0.22192091 -1.08125012 -1.12427512 -0.71688546 -1.1085952
 -1.18346893 -1.08445223 -1.25832959 -1.19968247 -0.79021949 -0.99707397
 -1.28779107 -0.52770089 -1.42404185  0.01687237 -0.77292183 -1.1052459
 -1.11571549 -1.27154062 -1.11380834 -0.16187578 -1.23084779 -1.22574523
 -0.70518482 -0.49906001 -0.7615905  -0.41812537 -0.73317902 -1.68142585
 -0.81071626 -0.23488113 -0.34822746 -1.10727835 -0.42260074 -1.6546409
 -1.49261527 -0.90127569 -1.49261527 -1.66000934 -1.43324552 -1.65233464
 -1.85062283 -2.09576626 -0.44010007 -0.45888953 -0.77245503 -1.65601458
 -1.12686292 -1.04074301 -0.92326733 -0.97476455 -0.84675893 -1.80742848
 -1.40186706 -1.04094079 -0.70979814 -1.60096215 -0.41856511 -0.72610862
 -1.4562446  -1.0665244  -0.46706505 -1.38198375 -1.25760671 -1.0808657
 -1.59556154 -0.32189281 -0.79597663 -1.93255621 -1.59207192 -2.07910113
 -1.07180095 -1.38608077 -0.82031434 -0.30511276 -1.16681897 -0.8772653
 -0.93650205 -0.76794997 -1.15759309 -1.1121371  -0.71193485 -1.1868807
 -0.5005378  -0.47215389 -1.01214534 -0.46620318 -1.53544579 -0.55705495
 -1.05277741 -0.35849828 -1.14055328 -1.30270878 -0.74878048 -0.99471398
 -1.21451568 -1.13348879 -0.74782102 -0.57042687 -0.28628111 -0.69752296
 -0.81037884 -1.09036425 -0.94776831 -1.12827957 -0.77564526 -1.08667526
 -0.69415821 -1.01947868 -1.59297464 -1.00680347 -0.90627221 -0.4519196
 -0.90744999 -1.73242357 -1.09805708 -0.57042687 -0.5683554  -1.17054105
 -0.5330147  -1.13867337 -1.05079947 -1.09936987]

In [256]:
corner.corner(np.hstack((lnf0.reshape(-1,1), beta.reshape(-1,1), alpha.reshape(-1,1))), labels=[r"$\ln F$", r"$\beta$", r"$\alpha$"]);
#corner.corner(np.hstack((beta.reshape(-1,1), alpha.reshape(-1,1))), labels=[r"$\beta$", r"$\alpha$"]);


print(alpha.shape)


(250,)

In [ ]:


In [ ]: