robclewley/Scenario-local-linear.yml

## readme.md

      
    Raw
  

              readme.md
            
          
    Literate modeling first demo

Start by reading the setup at https://gist.github.com/robclewley/0e54e7becf3bb017ea61#file-scenario-local-linear-yml

  
## index.html
<!doctype html>
<html>
  <head>
  </head>
  <body>
    <iframe src="http://nbviewer.ipython.org/gist/robclewley/0e54e7becf3bb017ea61/studycontext1.ipynb" width="600" height="800"></iframe>
  </body>
</html>

## model.py
"""
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-local-linear.yml
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


## setup.yml
pars:
    eps: 0.1
    a: 0.5
target_dom:
    x: [1, 1.3]
    y: [-0.9, -0.75]
error_tol: 0.002

## studycontext1.ipynb

      
Display the source blob

    
Display the rendered blob

    
    Raw
  

              studycontext1.ipynb
            
          
      Loading

      Sorry, something went wrong. Reload?
      Sorry, we cannot display this file.
      Sorry, this file is invalid so it cannot be displayed.
      
          Viewer requires iframe.
      
    
## studycontext1.py
"""
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')

## viz_errors1.png

      
    Raw
  

              viz_errors1.png
            
          
## viz_VF_overlay1.png

      
    Raw
  

              viz_VF_overlay1.png
	<!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 - (xxx/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
	"""
	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_half1.1, 'y': ydom_half1.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')