import numpy

class ufo():

    location = 2
    height = 1

    def EmitPulse(self):
        import random
        import math
        dir = random.uniform(0, 2 * 3.14159265)
        #
        #        _U
        #       /d|h
        # +---+l--+
        # 0   1   2
        #
        # tan(d) = (2 - l) / h
        # h * tan(d) = 2 - l
        # l = 2 - h * tan(d)

        h = 2 - self.height * math.tan(dir)
        return dir, h


def logPosteriorPosition(pos, priorLocation, hs):
    # Sivia-Skilling Data Analysis Bayesian Tutorial 2.38
    priors2, hs2 = numpy.meshgrid(priorLocation, hs)
    ph2 = numpy.log(pos ** 2 + (hs2 - priors2) ** 2)
    po3 = - ph2.sum(axis=0)

    #postss = []
    #for i, p in enumerate(priorLocation):
        #po = -sum(numpy.log(height ** 2 + (hs - p) ** 2))
        #print(i,p,po)
        #postss.append(po)

    #print(postss)
    print(po3)
    #print(len(postss), len(po3))

    return po3


l = ufo()

priorsh = numpy.linspace(-2, 10, 100)

metans = []
metahs = []
metaposts = []
for n in range(6):
    hs = []
    for i in range(10**n):
        d, h = l.EmitPulse()
        hs.append(h)

    h_mean = sum(hs) / len(hs)
    logPost = logPosteriorPosition(l.height, priorsh, numpy.array(hs))
    logPostmax = numpy.argmax(logPost)
    h_post = priorsh[logPostmax]
    metans.append(n)
    metahs.append(h_mean)
    metaposts.append(h_post)
    print(n, h_mean, h_post)

import pickle
f = open("detections.pickle", 'wb')
pickle.dump(hs, f)

import matplotlib
matplotlib.interactive(True)
from matplotlib import pylab

pylab.figure()
pylab.scatter(metans, metahs)
pylab.scatter(metans, metaposts)
pylab.axhline(l.location)
pylab.show()

pylab.figure()
pylab.hist(hs, range=(0, 5), bins=100)
pylab.show()