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')
@robclewley
Copy link
Author

Meta Project overview

First pass at a structured document style that defines a computational modeling scenario, to be transformed into an executable StudyContext object.

Note: UAC = user acceptance criteria

For the purposes of a minimal demonstration, this contrived mock-up pretends that it's difficult to match a linear model to some vector field data.

This whole git history is not the complete story of the "literate model". The version of the interactive computational environment used by these scripts also needs to be captured, e.g. with Sumatra.

Ideally, the concatenation of all triple quoted string segments would form a YAML or ArchieML description of the tasks. The ipython notebook here is a start of that, but the YAML is not displayed properly, and some other preprocessor has to extract all the YAML into one document anyway. To weave the document otherwise (in plain .py form), a preprocessor would have to switch literal code with comment text.

Would like to track the git history of the file to connect dependencies of studycontext content objects

The continuation of this project only makes sense with sub-directories, thus it is now a regular github project

@robclewley
Copy link
Author

Meta notes about project activity timeline:

  • Started a new Wing IDE project
  • Initialized git based on a clone of this gist
  • Created the black box model library (with a public interface) [details of model invisible to investigator]
  • Made parameters yaml file, setup.yml [details of model invisible to investigator]
  • Public interface = access to the scenario "givens" only
  • Created study context python script based on iterative edits of the the original Scenario specification yaml file
  • Perform (meta) sanity checks on setup and commit "empty" script that contributes no steps to the solution

@robclewley
Copy link
Author

Problems with format so far:

  • comment structure is annoying: want triple-quote docstrings inside the functions, and not repeated outside for the studycontext metadata; don't want to keep indenting exactly the same amount for next note
  • May need a directive at start of each comment to signify to preprocessor to parse this differently
  • what to do if functions defined during development of studycontext1.py need to be pushed to a common library? how does the provenance markup showing the order of sub-tasks in that file get affected?
  • each executable 'step' should be an individually executable script, so that different sub-tasks can easily be repeated and their setup iterated without having to comment out unwanted other steps
  • Need scripts to automate git add/commit/push with automated message part including "step #" of the process in the commit message
  • How to add math markup in a way that will be transformed to readable math in the browser?

@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