Skip to content

Instantly share code, notes, and snippets.

@grisu48
Last active September 4, 2019 08:41
Show Gist options
  • Save grisu48/0172f003c4e751a5a6873b03b1facaff to your computer and use it in GitHub Desktop.
Save grisu48/0172f003c4e751a5a6873b03b1facaff to your computer and use it in GitHub Desktop.
Gaussian2d implementation and test
#!/usr/bin/env python3
# -*- coding: utf-8 -*-
import random
import numpy as np
from math import *
import time
import sys
class Gauss2d :
'''
2D Gaussian:
f(x,y) = A * exp( - (x-x0)^2/(2s_x^2) - (y-y0)^2/(2s_y^2) )
'''
def __init__(self, A, x0, s_x, y0, s_y) :
self._a = A
self._x0 = x0
self._s_x = s_x
self._y0 = y0
self._s_y = s_y
def get(self, x,y) :
a1 = (x-self._x0)**2/(2*self._s_x**2)
a2 = (y-self._y0)**2/(2*self._s_y**2)
return self._a * exp(-(a1+a2))
def volume(self) :
'''
Analytical volume under the gaussian
'''
return 2*pi*self._a*self._s_x*self._s_y
def __str__(self) :
return "Gauss2d(%e, [%f, %f], [%f, %f])" % (self._a, self._x0, self._s_x, self._y0, self._s_y)
def gen_range(xb, d) :
x_min, x_max = xb
n = int(ceil((x_max-x_min / d)))
ret = [x_min + d*i for i in range(n)]
return np.array(ret)
def randfloat(f_min = 0.0, f_max = 1.0) :
f = random.randint(0, 100000) / 100000.0
return f_min + f*(f_max-f_min)
if __name__ == "__main__" :
nTests = 1000
tol = 1e-3 # Allowed tolerance
# Bounds
xb = [-50, 50]
yb = [-50, 50]
dx, dy = 0.1, 0.1
X = gen_range(xb, dx)
Y = gen_range(yb, dy)
#print(len(X), X)
#print(len(Y), Y)
def area(gauss) :
ret = 0.0
for x in X :
for y in Y :
ret+=gauss.get(x,y)
return ret*dx*dy
print("Testing %d gaussians ... " % nTests)
print("---------------------------------------")
print(" Hit CTRL+C if you had enough ... ")
print("---------------------------------------")
random.seed(time.time())
maxDev = 0.0
try :
for i in range(nTests) :
A = randfloat(1,10)
x0 = randfloat(-1,1)
sx = max(dx,dy)+randfloat(0,1)
y0 = randfloat(-1,1)
sy = max(dx,dy)+randfloat(0,1)
gauss = Gauss2d(A,x0,sx,y0,sy)
sys.stdout.write(" [%4d/%4d] " % (i+1,nTests))
sys.stdout.write("%-80s ... " % (str(gauss)))
sys.stdout.flush()
got, expect = area(gauss), gauss.volume()
maxDev = max(maxDev, (abs(got/expect - 1.0)))
sys.stdout.write("(%8.5f/%8.5f) = %.5f\n" % (got, expect, (got/expect)))
except KeyboardInterrupt :
sys.stdout.write("\n")
sys.stdout.flush()
sys.stderr.write("User termination\n")
sys.stderr.flush()
if maxDev > tol :
print("Uh oh: Max ratio: %e" % (maxDev+1.0))
sys.exit(1)
else :
print("All good. Max ratio: %e" % (maxDev+1.0))
sys.exit(0)
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment