Skip to content

HTTPS clone URL

Subversion checkout URL

You can clone with
or
.
Download ZIP
Browse files

added some examples: QR decomposition and explicit euler

  • Loading branch information...
commit cc8812378c6886ca77bbc6e7d327b7b3674fc25d 1 parent 52ad8ce
@b45ch1 authored
View
4 adolc/linalg.py
@@ -16,7 +16,9 @@ def qr(in_A):
"""
QR decomposition of A
- Q,R = qr(A)
+ Q.T,R = qr(A)
+
+ I.e. returns the Q.T instead of Q!
"""
# input checks
View
69 examples/explicit_euler.py
@@ -0,0 +1,69 @@
+#!/usr/bin/env python
+
+from numpy import *
+from pylab import *
+import adolc
+
+
+
+def explicit_euler(x0,f,ts,p,q):
+ N = size(ts)
+ if isinstance(p[0],adolc._adolc.adouble):
+ x = array([adolc.adouble(0) for m in range(Nm)])
+ else:
+ x = zeros(N)
+ x[0] = x0
+ for n in range(1,N):
+ h = ts[n] - ts[n-1]
+ x[n]= x[n-1] + h*f(ts[n-1],x[n-1],p,q)
+ return x
+
+
+def f(t,x,p,q):
+ """ rhs of the ODE"""
+ return p[1] + q[0]*x
+
+# analytic solutions
+def phi(t,p,q):
+ return ((p[1] + q[0]*p[0]) * exp(t/q[0]) - p[1])/q[0]
+
+def phip0(t,p,q):
+ return exp(q[0]*t)
+
+def phip1(t,p,q):
+ return (exp(q[0]*t)-1.)/q[0]
+
+
+p = array([10.,2.])
+q = array([-1.])
+v = concatenate((p,q))
+
+Np = size(p)
+Nq = size(q)
+Nv = size(v)
+Nm = 100
+ts = linspace(0,10,Nm)
+
+# TAPING THE INTEGRATOR
+adolc.trace_on(1)
+av = adolc.adouble(v)
+adolc.independent(av)
+ax = explicit_euler(av[0],f,ts,av[:Np],av[Np:])
+adolc.dependent(ax)
+adolc.trace_off()
+
+# COMPUTING FUNCTION AND JACOBIAN FROM THE TAPE
+y = adolc.zos_forward(1,v,0)
+J = adolc.jacobian(1,v)
+
+x_plot = plot(ts, y,'b')
+x_analytical_plot = plot(ts,phi(ts,p,q),'b.')
+
+xp0_plot = plot(ts, J[:,0], 'g')
+xp0_analytical_plot = plot(ts, phip0(ts,p,q), 'g.')
+
+xp1_plot = plot(ts, J[:,1], 'r')
+xp1_analytical_plot = plot(ts, phip1(ts,p,q), 'r.')
+
+show()
+print J
View
30 examples/prettyplotting.py
@@ -0,0 +1,30 @@
+
+import pylab as pyl
+
+#SETUP PLOTTING PARAMETERS
+fig_width_pt = 350.0 # Get this from LaTeX using \showthe\columnwidth
+inches_per_pt = 1.0/72.25 # Convert pt to inches
+golden_mean = (pyl.sqrt(5)-1.0)/2.0 # Aesthetic ratio
+fig_width = fig_width_pt*inches_per_pt # width in inches
+fig_height =fig_width*golden_mean # height in inches
+fig_size = [fig_width,fig_height]
+legend_padding = 0.05
+params = {
+ 'backend': 'ps',
+ 'ps.usedistiller': 'xpdf',
+ 'font.family' : 'sans-serif',
+ 'font.style' : 'normal',
+ 'font.variant' : 'normal',
+ 'font.weight' : 'normal', #bold
+ 'font.stretch' : 'normal',
+ 'font.size' : 'normal', #large
+ 'axes.labelsize': 11,
+ 'text.fontsize': 10,
+ 'title.fontsize':10,
+ 'legend.fontsize':10,
+ 'xtick.labelsize': 11,
+ 'ytick.labelsize': 11,
+ 'lines.markersize':5,
+ 'text.usetex': True,
+ 'figure.figsize': fig_size}
+pyl.rcParams.update(params)
View
29 examples/qr_decomposition.py
@@ -0,0 +1,29 @@
+from adolc import *
+import adolc
+import adolc.linalg
+from numpy import *
+
+N = 4
+A = array([[ r*N + c for c in range(N) ] for r in range(N)],dtype=float)
+
+aA = adouble(A)
+trace_on(1)
+independent(aA)
+(aQT, aR) = adolc.linalg.qr(aA)
+dependent(aQT)
+dependent(aR)
+trace_off()
+
+V = zeros((N**2,1,1),dtype=float)
+V[:,0,0] = 1.
+(result,Z) = hov_forward(1, ravel(A), V)
+QT = result[:N**2].reshape((N,N))
+R = result[N**2:].reshape((N,N))
+
+#(QT,R) = myqr(A)
+#print dot(QT.T,R) - A
+
+QTdot = Z[:N**2,0,0].reshape((N,N))
+Rdot = Z[N**2:,0,0].reshape((N,N))
+
+print QTdot.T + dot(dot(QT.T,QTdot), QT.T)
Please sign in to comment.
Something went wrong with that request. Please try again.