Skip to content

Instantly share code, notes, and snippets.

@robclewley
Last active August 29, 2015 14:18
Show Gist options
  • Save robclewley/0e54e7becf3bb017ea61 to your computer and use it in GitHub Desktop.
Save robclewley/0e54e7becf3bb017ea61 to your computer and use it in GitHub Desktop.
Literate modeling first development ideas - continued on github
<!doctype html>
<html>
<head>
</head>
<body>
<iframe src="http://nbviewer.ipython.org/gist/robclewley/0e54e7becf3bb017ea61/studycontext1.ipynb" width="600" height="800"></iframe>
</body>
</html>
"""
Van der Pol equations to generate a nonlinear 2D vector field
and save a black box object
"""
from PyDSTool import *
from PyDSTool.Toolbox import phaseplane as pp
import sys
# public exports
__all__ = ['F', 'target_dom', 'compute', 'L2_tol']
# ----
import yaml
with open('setup.yml') as f:
setup = yaml.load(f)
L2_tol = setup['error_tol']
#pars = {'eps': 1e-1, 'a': 0.5}
pars = setup['pars']
# target_dom = {'x': [-2.5, 2.5], 'y': [-2, 2]}
target_dom = setup['target_dom']
# Convenience variables
xdom = target_dom['x']
ydom = target_dom['y']
xInterval = Interval('xdom', float, xdom, abseps=1e-3)
yInterval = Interval('ydom', float, ydom, abseps=1e-3)
def build():
icdict = {'x': pars['a'],
'y': pars['a'] - pars['a']**3/3}
xstr = '(y - (x*x*x/3 - x))/eps'
ystr = 'a - x'
#DSargs.fnspecs = {'Jacobian': (['t','x','y'],
# """[[(1-x*x)/eps, 1/eps ],
# [ -1, 0 ]]""")}
algparams = {'max_pts': 3000, 'init_step': 0.01, 'stiff': True}
compatGens = findGenSubClasses('ODEsystem')
class ODEComponent(LeafComponent):
compatibleGens=compatGens
targetLangs=targetLangs
# compatibleSubcomponents=(sys,)
vdpC = ODEComponent('vdp_gen')
vdpC.add([Var(xstr, 'x', domain=target_dom['x'], specType='RHSfuncSpec'),
Var(ystr, 'y', domain=target_dom['y'], specType='RHSfuncSpec'),
Par(pars['a'], 'a'), Par(pars['eps'], 'eps')])
vdpC.flattenSpec()
MC = ModelConstructor('vdp',
generatorspecs={'vdp_gen':
{'modelspec': vdpC,
'target': 'Vode_ODEsystem',
'algparams': algparams}
},
indepvar=('t',[0,10]),
eventtol=1e-5, activateAllBounds=True,
withStdEvts={'vdp_gen': True},
stdEvtArgs={'term': True,
'precise': True,
'eventtol': 1e-5},
abseps=1e-7,
withJac={'vdp_gen': True})
# withJac option for automatic calc of Jacobian
return MC.getModel()
try:
vdp = loadObjects('model.sav')[0]
except:
vdp = build()
saveObjects(vdp, 'model.sav')
def F(xarray):
x, y = xarray
assert x in xInterval
assert y in yInterval
return vdp.Rhs(0, {'x':x, 'y':y}, asarray=True)
def compute(xarray, trajname='test'):
x, y = xarray
assert x in xInterval
assert y in yInterval
icdict = {'x': x,
'y': y}
vdp.set(ics=icdict)
vdp.compute(trajname, tdata=[0,5])
pts = vdp.sample(trajname)
return pts
def sanity_check():
print target_dom
x1 = sum(xdom)/2
y1 = sum(ydom)/2
print F((x1,y1))
pts = compute((x1,y1), 'test')
print len(pts)
plt.plot(pts['x'], pts['y'], 'k-o')
if __name__ == '__main__':
vdp_gen = vdp.registry.values()[0]
pp.plot_PP_vf(vdp_gen, 'x', 'y', subdomain=target_dom,
N=20, scale_exp=-1)
sanity_check()
plt.show()
scenario-metadata:
name: local linear fit
difficulty: easy
time: short
tags: graphical, ODE
Given:
- 2D, weakly nonlinear, determininstic vector field given as a black box function of (x,y) -> (f_x(x,y), f_y(x,y))
- finite, rectangular domain [x0,x1] X [y0,y1]
- Euclidean error tolerance of forward trajectories (chosen to be feasible for the given function f, unless task scenario is open to being proven inconsistent)
Goal:
- Derive a UAC-satisfying local linear ODE approximation in the domain (where nonlinearity is fairly weak)
UAC:
- Predict all forward trajectories up to given max tolerance of Euclidean error
pars:
eps: 0.1
a: 0.5
target_dom:
x: [1, 1.3]
y: [-0.9, -0.75]
error_tol: 0.002
Display the source blob
Display the rendered blob
Raw
Loading
Sorry, something went wrong. Reload?
Sorry, we cannot display this file.
Sorry, this file is invalid so it cannot be displayed.
"""
studycontext-metadata:
ID: 1
tag: solve scenario
"""
import PyDSTool as dst
from PyDSTool.Toolbox import phaseplane as pp
import numpy as np
from matplotlib import pyplot as plt
"""
studycontext-givens:
tag: import public interface to the target model
"""
from model import F, target_dom, compute, L2_tol
# Convenience variables
xdom = target_dom['x']
ydom = target_dom['y']
xdom_half = sum(xdom)/2
ydom_half = sum(ydom)/2
xdom_width = xdom[1]-xdom[0]
ydom_width = ydom[1]-ydom[0]
xInterval = dst.Interval('xdom', float, xdom, abseps=1e-3)
yInterval = dst.Interval('ydom', float, ydom, abseps=1e-3)
# Convenience function for orientation of new users
# (not explicit part of study context)
def demo_sim():
print target_dom
print F((xdom_half,ydom_half))
pts = compute((xdom_half,ydom_half), 'test')
print len(pts)
plt.plot(pts['x'], pts['y'])
plt.show()
"""
studycontext-step:
tag: build linear model
notes:
- define build function
"""
def build_lin():
# make local linear system spec
DSargs = dst.args(name='lin')
xfn_str = '(x0+yfx*y - x)/taux'
yfn_str = '(y0+xfy*x - y)/tauy'
DSargs.varspecs = {'x': xfn_str, 'y': yfn_str}
DSargs.xdomain = {'x': xdom, 'y': ydom}
DSargs.pars = {'x0': xdom_half, 'y0': ydom_half,
'xfy': 1, 'yfx': 1,
'taux': 1, 'tauy': 1}
DSargs.algparams = {'init_step':0.001,
'max_step': 0.001,
'max_pts': 10000}
DSargs.checklevel = 0
DSargs.tdata = [0, 10]
DSargs.ics = {'x': xdom_half*1.1, 'y': ydom_half*1.1}
DSargs.fnspecs = {'Jacobian': (['t', 'x', 'y'],
"""[[-1/taux, -yfx/taux],
[-xfy/tauy, -1/tauy]]""")}
return dst.Generator.Vode_ODEsystem(DSargs)
"""
- instantiate local model
"""
lin = build_lin()
"""
studycontext-uac:
tag: define uac
notes:
- test sample 30 x 30 grid of ICs in domain
"""
def make_mesh_pts(N=30):
xsamples = xInterval.uniformSample(xdom_width/N, avoidendpoints=True)
ysamples = yInterval.uniformSample(ydom_width/N, avoidendpoints=True)
xmesh, ymesh = np.meshgrid(xsamples,ysamples)
return xmesh, ymesh, np.dstack((xmesh,ymesh)).reshape(N*N,2)
xmesh5, ymesh5, mesh_pts_test5 = make_mesh_pts(5)
xmesh30, ymesh30, mesh_pts_test30 = make_mesh_pts(30)
"""
- mirror signature of F for linear model
"""
def LF(pt):
x, y = pt
return lin.Rhs(0, {'x':x, 'y':y})
"""
- L2 metric between vector fields
"""
def metric(pt):
# could usefully vectorize this
return pp.dist(LF(pt), F(pt))
def condition(m, tol):
return m < tol
"""
studycontext-goal:
tag: define goal
notes: functions with richer diagnostic feedback about spatial errors
"""
def error_pts(mesh_pts):
return np.array([metric(pt) for pt in mesh_pts])
def test_goal(mesh_pts, goal_tol=L2_tol):
errors_array = error_pts(mesh_pts)
result = condition(np.max(errors_array), goal_tol)
return dst.args(result=result,
errors=errors_array)
"""
studycontext-test:
tag: initial test that goal fails with a small sample
"""
test = test_goal(mesh_pts_test5)
assert not test.result
"""
studycontext-step:
tag: visualize mesh of errors
notes:
- figure 1
"""
def viz_errors(mesh_pts, errors, fignum, scaling=2):
fig = plt.figure(fignum)
fig.clf()
for pt, err in zip(mesh_pts, errors):
plt.plot(pt[0], pt[1], 'ko', markersize=scaling*err)
plt.show()
test30 = test_goal(mesh_pts_test30)
viz_errors(mesh_pts_test30, test30.errors, 1)
fig1 = plt.figure(1)
fig1.savefig('viz_errors1.png')
"""
studycontext-step:
tag: visualize mesh of linear vectors overlaid with actual vectors
notes:
- figure 2
"""
def get_all_Fs(F, mesh_pts):
return np.array([F(pt) for pt in mesh_pts])
# store globally for later convenience
all_Fs = get_all_Fs(F, mesh_pts_test30)
all_LFs = get_all_Fs(LF, mesh_pts_test30)
def Fmesh(xmesh, ymesh, F):
dxs, dys = np.zeros(xmesh.shape, float), np.zeros(ymesh.shape, float)
for xi, x in enumerate(xmesh[0]):
for yi, y in enumerate(ymesh.T[0]):
dx, dy = F((x,y))
# note order of indices
dxs[yi,xi] = dx
dys[yi,xi] = dy
return dxs, dys
Fmeshes5 = Fmesh(xmesh5, ymesh5, F)
LFmeshes5 = Fmesh(xmesh5, ymesh5, LF)
Fmeshes30 = Fmesh(xmesh30, ymesh30, F)
LFmeshes30 = Fmesh(xmesh30, ymesh30, LF)
def viz_VF(Fmeshes, meshes, fignum, col):
plt.figure(fignum)
Fxmesh, Fymesh = Fmeshes
xmesh, ymesh = meshes
plt.quiver(xmesh, ymesh, Fxmesh, Fymesh, angles='xy', pivot='middle',
scale=10, lw=0.5, width=0.005*max(xdom_width,ydom_width),
minshaft=2, minlength=0.1,
units='inches', color=col)
# test case
viz_VF(Fmeshes5, (xmesh5, ymesh5), 2, 'r')
viz_VF(LFmeshes5, (xmesh5, ymesh5), 2, 'k')
#viz_VF(Fmeshes30, (xmesh30, ymesh30), 2, 'r')
fig2 = plt.figure(2)
fig2.savefig('viz_VF_overlay1.png')
@tonyfast
Copy link

tonyfast commented Apr 2, 2015

So right now, the names in your files only have a meaning in the context of PyDSTool. The information is hard to liberate and parse.

Python files are meant for code and documentation. The semantic information is orthogonal to Python. You're structured document is in your language.

It should be an Ipython Notebook. The cell structure allows you to use mixed languages, python and yaml orthogonally. So magic can parse the yaml on execute.

I really think the first and easiest thing to change maybe the presentation layer of the python files are a notebook. Then people can see it.

Similarly, we can create a bl.ocks view of it and add readme.md to describe the larger context.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment