ORTHOGONAL POLYNOMIAL DENSITY ESTIMATION

In [1]:

    

import numpy as np
import matplotlib.pyplot as plt
import seaborn as sbn
import pandas as pd
from uuid import uuid4

from lpde.geometry import WidthOf, Window, PointAt, BoundingBox, Mapper, Grid
from lpde.estimators import ParallelEstimator
from lpde.estimators.datatypes import Event, Degree, Action, Scalings, Signal
from lpde.producers import MockParams, MockProducer
from lpde.producers.distributions import gaussian


    

In [2]:

    

%matplotlib notebook


    

In [5]:

    

legendre_width = WidthOf(1.8)

center = PointAt(0, 0)
window = Window(1.8, 1.8)
bounds = BoundingBox(center, window)

mapper = Mapper(bounds, legendre_width)

degree = Degree(20, 20)
params = MockParams(20, 100, 1, gaussian)
density = ParallelEstimator(degree, mapper, params)

action = Action.ADD
point = PointAt(0.5, 0.5)
event = Event(uuid4(), action, point)

grid = Grid(100, 100)


    

In [ ]:

    

def gaussian():
    x, y = np.random.multivariate_normal((0,0), ((0.1,0), (0,0.1)))
    if (-0.9 <= x <= 0.9) and (-0.9 <= y <= 0.9):
        return x, y
    else:
        return gaussian()

def uniform():
    return np.random.uniform(low=-0.9, high=0.9, size=2)

locations = []

def new_event(dist):
    location = dist()
    locations.append(location)
    point = PointAt(*location)
    return Event(uuid4(), Action(1), point)

def random_event(dist):
    event_type = np.random.randint(low=-1, high=2)
    if event_type == 1:
        location = dist()
        point = PointAt(*location)
        return Event(uuid4(), Action(1), point)
    elif event_type == 0:
        location = dist()
        point = PointAt(*location)
        column = density._phi.sample(1, axis=1).columns.values[0]
        return Event(column, Action(0), point)
    column = density._phi.sample(1, axis=1).columns.values[0]
    return Event(column, Action(-1))


    

In [6]:

    

density.controller.start(1, 1.0)


    

    




/Users/georg/Documents/Python/LPDE/lpde/estimators/parallel/minimizer.py:78: RuntimeWarning: overflow encountered in double_scalars
  return self.__neg_log_l(c[1:]) + c[0]*self.__norm(c[1:])
/Users/georg/Documents/Python/LPDE/lpde/estimators/parallel/minimizer.py:86: RuntimeWarning: overflow encountered in square
  return -log(square(c.dot(self.__phi_ijn))).sum()
/Users/georg/Documents/Python/LPDE/lpde/estimators/parallel/minimizer.py:82: RuntimeWarning: overflow encountered in multiply
  self.__grad_c[1:] = self.__grad_neg_log_l(c[1:]) + 2.0*c[0]*c[1:]

In [8]:

    

density.producer.start()


    

In [13]:

    

density.at(point)


    

    
Out[13]:





3.8338478920782548

In [34]:

    

density.controller.stop()


    

In [ ]:

    

#%%time
for i in range(100):
    density.update_with(new_event(gaussian))
#while not density.controller.phi_queue.empty():
#    pass


    

In [ ]:

    

%%time
for i in range(1000):
    density.update_with(random_event(gaussian))
while not density._phi_queue_empty:
    pass


    

In [ ]:

    

%%time
for i in range(1000):
    density.at(new_event(uniform).location)
while not density._phi_queue_empty:
    pass


    

In [33]:

    

n_hat = density.on(grid)

fig, ax = plt.subplots()
ax.set(xlabel=r'$x$', ylabel=r'$y$')
#ax.scatter(*zip(*locations), s=5, c='k')
contour = ax.imshow(n_hat,
                    cmap='viridis',
                    alpha=0.9,
                    extent=bounds.x_range+bounds.y_range,
                    origin='lower',
                    interpolation='bicubic')
cbar = plt.colorbar(contour, ax=ax, label=r'$n(x)$')
fig.tight_layout()


    

In [ ]:

    

from bokeh.io import output_file, show
from bokeh.plotting import figure
from bokeh.tile_providers import STAMEN_TONER, STAMEN_TERRAIN

bound = 20000000 # meters
fig = figure(tools='pan, wheel_zoom', x_range=(-bound, bound), y_range=(-bound, bound))
fig.axis.visible = False
fig.add_tile(STAMEN_TERRAIN)
output_file("stamen_toner_plot.html")
show(fig)


    

In [ ]:

    

from bokeh.plotting import figure, show, output_file

p = figure(x_range=bounds.x_range, y_range=bounds.y_range)

# must give a vector of image data for image parameter
p.image(image=[n_hat],
        x=bounds.x_range[0],
        y=bounds.y_range[0],
        dw=bounds.window[0],
        dh=bounds.window[1],
        palette="Spectral11",
        alpha=0.8)

output_file("image.html", title="image.py example")

show(p)  # open a browser


    

In [ ]:

    

def f():
    pass


    

In [ ]:

    

type(f)


    

In [ ]:

    

callable(function)


    

In [ ]: