Start by reading the setup at https://gist.github.com/robclewley/0e54e7becf3bb017ea61#file-scenario-local-linear-yml
Last active
August 29, 2015 14:18
-
-
Save robclewley/0e54e7becf3bb017ea61 to your computer and use it in GitHub Desktop.
Literate modeling first development ideas - continued on github
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| <!doctype html> | |
| <html> | |
| <head> | |
| </head> | |
| <body> | |
| <iframe src="http://nbviewer.ipython.org/gist/robclewley/0e54e7becf3bb017ea61/studycontext1.ipynb" width="600" height="800"></iframe> | |
| </body> | |
| </html> |
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| """ | |
| 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() |
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| 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 | |
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| pars: | |
| eps: 0.1 | |
| a: 0.5 | |
| target_dom: | |
| x: [1, 1.3] | |
| y: [-0.9, -0.75] | |
| error_tol: 0.002 |
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| """ | |
| 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') |
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
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?
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
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