Gaussian2d implementation and test

Gaussian2d.py

#!/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)