#!/usr/bin/python -Wall

# ================================================================
# John Kerl
# kerl at math dot arizona dot edu
# 2008-05-12
# ================================================================

from   __future__ import division # 1/2 = 0.5, not 0.
import random
from   math import *
import copy
import sys
import re

# ----------------------------------------------------------------
def usage():
	print >> sys.stderr, "Usage: %s [options]" % (sys.argv[0])
	print >> sys.stderr, "Options:"
	print >> sys.stderr, "  T=[number]"
	print >> sys.stderr, "  nt=[number]"
	print >> sys.stderr, "  X0=[number]"
	print >> sys.stderr, "  nX=[number]"
	print >> sys.stderr, "  alpha=[number]"
	print >> sys.stderr, "  beta=[number]"
	sys.exit(1)

# ----------------------------------------------------------------
# SDE:      dX_t  = alpha X_t dt + beta X_t db_t
# Solution:  X_t  = X_0 exp((alpha - beta^2/2)t) exp(beta b_t).
# Mean:    E[X_t] = X_0 exp(alpha t).

T  = 2
nt = 2000
X0 = 1
nX = 20
alpha = 2
beta  = 3

argc = len(sys.argv)
for argi in range(1, argc):
	if   re.match(r'^T=',     sys.argv[argi]): T     = float(sys.argv[argi][2:])
	elif re.match(r'^nt=',    sys.argv[argi]): nt    = int  (sys.argv[argi][3:])
	elif re.match(r'^X0=',    sys.argv[argi]): X0    = float(sys.argv[argi][3:])
	elif re.match(r'^nX=',    sys.argv[argi]): nX    = int  (sys.argv[argi][3:])
	elif re.match(r'^alpha=', sys.argv[argi]): alpha = float(sys.argv[argi][6:])
	elif re.match(r'^beta=',  sys.argv[argi]): beta  = float(sys.argv[argi][5:])
	else: usage()

t  =  0
dt = T/nt
exact  = [X0]  * nX
approx = [X0]  * nX
Bt     = [0.0] * nX
sqrtdt = sqrt(dt)

print >> sys.stderr, "# alpha     = %11.7f" % (alpha)
print >> sys.stderr, "# beta**2/2 = %11.7f" % (beta**2/2)
print >> sys.stderr, "# diff      = %11.7f" % (alpha - beta**2/2)

while (t <= T):
	for k in range(0, nX):
		exact[k] = X0 * exp((alpha - beta**2/2) * t) * exp(beta * Bt[k])

	print "%11.7f" % (t),
	for k in range(0, nX):
		#print "%11.7f %11.7f %9.3e" % \
		#	(exact[k], approx[k], exact[k]-approx[k]),
		print "%11.7f" % (exact[k]),
	print "%11.7f" % (X0 * exp(alpha * t)) # This is E[X_t].

	for k in range(0, nX):
		dB  = random.gauss(0, sqrtdt)
		approx[k] += alpha * approx[k] * dt + beta * approx[k] * dB
		Bt[k] += dB
	t  += dt