mktranstrum/MBAM

Newer

Older

Raw Normal view History

100644 51 lines (39 sloc) 1.61 KB

Feb 9, 2016 geodesic.py script and an example of usage.
	1	# An example of calculating a geodesic
	2	# model takes the form y(t,x) = e^(-x[0]t) + e^(-x[1]t) for time points t = 0.5, 2.0
	3	# We enforce x[i] > 0 by going to log parameters:
	4	# y(t,x) = e^(-exp(x[0])t) + e^(-exp(x[1])t)
	5
	6	from geodesic import geodesic, InitialVelocity
	7	import numpy as np
	8	exp = np.exp
	9
	10	# Model predictions
	11	def r(x):
	12	return np.array([ exp(-exp(x[0])0.5) + exp(-exp(x[1])0.5), exp(-exp(x[0])2) + exp(-exp(x[1])2)])
	13
	14
	15	# Jacobian
	16	def j(x):
	17	return np.array([ [ -0.5exp(x[0])exp(-x[0]0.5), -2.0exp(x[0])exp(-x[0]2)],
	18	[ -0.5exp(x[1])exp(-x[1]0.5), -2.0exp(x[1])exp(-x[1]2)] ] )
	19
	20	# Directional second derivative
	21	def Avv(x,v):
	22	h = 1e-4
	23	return (r(x + hv) + r(x - hv) - 2*r(x))/h/h
	24
	25
	26
	27	# Choose starting parameters
	28	x = np.log([1.0, 2.0])
	29	v = InitialVelocity(x, j, Avv)
	30
	31	# Callback function used to monitor the geodesic after each step
	32	def callback(geo):
	33	# Integrate until the norm of the velocity has grown by a factor of 10
	34	# and print out some diagnotistic along the way
	35	print "Iteration: %i, tau: %f, \|v\| = %f" %(len(geo.vs), geo.ts[-1], np.linalg.norm(geo.vs[-1]))
	36	return np.linalg.norm(geo.vs[-1]) < 10.0
	37
	38	# Construct the geodesic
	39	# It is usually not necessary to be very accurate here, so we set small tolerances
	40	geo = geodesic(r, j, Avv, 2, 2, x, v, atol = 1e-2, rtol = 1e-2, callback = callback)
	41
	42	# Integrate
	43	geo.integrate(25.0)
	44	import pylab
	45	# plot the geodesic path to find the limit
	46	# This should show the singularity at the "fold line" x[0] = x[1]
	47	pylab.plot(geo.ts, geo.xs)
	48	pylab.xlabel("tau")
	49	pylab.ylabel("Parameter Values")
	50	pylab.show()