Commit b229a013 authored by Jan Fousek's avatar Jan Fousek

Cleanup and notes

parent 91447409
Pipeline #2466 failed with stages
# CodeJam 9 code generators
Here we explore the code synthesis without the loopy framework. The basic idea is to stick with the simple structured model representation and use a combination of AST transformations and active code templates to target various scenarios and architectures. Ideally, the templates will be derived from existing hand-written implementations.
Relevant files:
* `models.py` : example of model structured definition
* `dsl2python.py` : simple templated synthesis of a TVB-like Python class with a phase-plane interactive plot
* `template.py` active template for Python model class
* `ppi.py` : phase-plane interactive plot
* `dsl2C.py` : example of AST transformations for C-like languages
* `template.cu` active template for Cuda dfun function
## Model description
This was adapted from the loopy implementation. Note, that this is more a data-structure (string drifts), than a class with executable functions. Could be as well a JSON object...
While it would be better to have a functional model class, it is difficult to create one without introducing redundancy in the notation -- especially in the drift function definitions. So let's just stick with this:
```Python
class G2DO:
"Generic nonlinear 2-D (phase plane) oscillator."
state = 'W', 'V'
limit = (-5, 5), (-5, 5)
input = 'c_0'
param = 'a'
const = {'tau': 1.0, 'I': 0.0, 'a': -2.0, 'b': -10.0, 'c': 0.0, 'd': 0.02,
'e': 3.0, 'f': 1.0, 'g': 0.0, 'alpha': 1.0, 'beta': 1.0,
'gamma': 1.0}
drift = (
'd * tau * (alpha*W - f*V**3 + e*V**2 + g*V + gamma*I + gamma*c_0)',
'd * (a + b*V + c*V**2 - beta*W) / tau'
)
diffs = 1e-3, 1e-3
obsrv = 'W', 'V'
```
## Code synthesis
For both C and Python we use the `mako` template engine to fill templates preparing the boilerplate for the `drift` expressions.
The advantage of the `drift` having Python syntax is, that we can use the built-in `ast` library for parsing, and the `ctree` library for basic transformations and code synthesis. For the demo, the simple operator conversion had to be extended to cope with `**`.
In future, we can use the `ctree` framework for more complex transformations, such as array index arithmetic etc.
A demo Python code synthesized for the phase-plane interactive plot:
```Python
class Generic2D:
def __init__(self,tau=1.0, I=0.0, a=-2.0, b=-10.0, c=0.0, d=0.02, e=3.0, f=1.0, g=0.0, alpha=1.0, beta=1.0, gamma=1.0):
self.limit = ((-5, 5), (-5, 5))
self.tau = tau
self.I = I
self.a = a
self.b = b
self.c = c
self.d = d
self.e = e
self.f = f
self.g = g
self.alpha = alpha
self.beta = beta
self.gamma = gamma
def dfun(self,state_variables, *args, **kwargs):
W, V = state_variables
tau = self.tau
I = self.I
a = self.a
b = self.b
c = self.c
d = self.d
e = self.e
f = self.f
g = self.g
alpha = self.alpha
beta = self.beta
gamma = self.gamma
c_0 = 0.0
dW = d * tau * (alpha*W - f*V**3 + e*V**2 + g*V + gamma*I + gamma*c_0)
dV = d * (a + b*V + c*V**2 - beta*W) / tau
return [W, V]
```
A demo code synthesized to fit into an old code generator for TVB.
```C
void model_dfun(
float * _dx, float *_x, float *mmpr, float*input
)
{
float tau = mmpr[n_thr*0];
float I = mmpr[n_thr*1];
float a = mmpr[n_thr*2];
float b = mmpr[n_thr*3];
float c = mmpr[n_thr*4];
float d = mmpr[n_thr*5];
float e = mmpr[n_thr*6];
float f = mmpr[n_thr*7];
float g = mmpr[n_thr*8];
float alpha = mmpr[n_thr*9];
float beta = mmpr[n_thr*10];
float gamma = mmpr[n_thr*11];
float W = _x[n_thr*0];
float V = _x[n_thr*1];
float c_0 = input[0];
_dx[n_thr*0] = d * tau * (alpha * W - f * pow(V, 3) + e * pow(V, 2) + g * V + gamma * I + gamma * c_0);
_dx[n_thr*1] = d * tau * (alpha * W - f * pow(V, 3) + e * pow(V, 2) + g * V + gamma * I + gamma * c_0);
}
```
......@@ -46,10 +46,10 @@ def dfuns_to_c(dfuns):
if __name__=='__main__':
from ppi import G2DO
from models import G2DO
from mako.template import Template
template = Template(filename='cuda_model.cu')
template = Template(filename='template.cu')
print(template.render( svar=G2DO.state,const=G2DO.const, cvar=G2DO.input,
dfuns=dfuns_to_c(G2DO.drift)))
......
from ppi import G2DO
from models import G2DO
from mako.template import Template
import argparse
parser = argparse.ArgumentParser(description='Demo Python code generator')
parser.add_argument('--plot', default=False, action='store_true',
help='launch a phase plane interactive plot')
args = parser.parse_args()
template = Template(filename='template.py')
print(template.render( name='Generic2D',
model_str = template.render( name='Generic2D',
const=G2DO.const,
limit=G2DO.limit,
sv=G2DO.state,
drift=G2DO.drift,
input=G2DO.input.split() ) )
input=G2DO.input.split() )
print(model_str)
if args.plot:
from ppi import PhasePlaneInteractive
exec(model_str)
model = Generic2D()
ppi = PhasePlaneInteractive(model)
ppi()
class G2DO:
"Generic nonlinear 2-D (phase plane) oscillator."
state = 'W', 'V'
limit = (-5, 5), (-5, 5)
input = 'c_0'
param = 'a'
const = {'tau': 1.0, 'I': 0.0, 'a': -2.0, 'b': -10.0, 'c': 0.0, 'd': 0.02,
'e': 3.0, 'f': 1.0, 'g': 0.0, 'alpha': 1.0, 'beta': 1.0,
'gamma': 1.0}
drift = (
'd * tau * (alpha*W - f*V**3 + e*V**2 + g*V + gamma*I + gamma*c_0)',
'd * (a + b*V + c*V**2 - beta*W) / tau'
)
diffs = 1e-3, 1e-3
obsrv = 'W', 'V'
......@@ -6,22 +6,6 @@ from matplotlib.widgets import Button
TRAJ_STEPS = 4096
class G2DO:
"Generic nonlinear 2-D (phase plane) oscillator."
state = 'W', 'V'
limit = (-5, 5), (-5, 5)
input = 'c_0'
param = 'a'
const = {'tau': 1.0, 'I': 0.0, 'a': -2.0, 'b': -10.0, 'c': 0.0, 'd': 0.02,
'e': 3.0, 'f': 1.0, 'g': 0.0, 'alpha': 1.0, 'beta': 1.0,
'gamma': 1.0}
drift = (
'd * tau * (alpha*W - f*V**3 + e*V**2 + g*V + gamma*I + gamma*c_0)',
'd * (a + b*V + c*V**2 - beta*W) / tau'
)
diffs = 1e-3, 1e-3
obsrv = 'W', 'V'
class Oscillator:
def __init__(self, eta=0.07674, gamma=1.21, epsilon=12.3):
......@@ -30,7 +14,7 @@ class Oscillator:
self.gamma = gamma
self.epsilon = epsilon
def dfun(self, state_variables):
def dfun(self, state_variables, *args, **kwargs):
psi1, psi2 = state_variables
......@@ -39,72 +23,72 @@ class Oscillator:
return [dpsi1, dpsi2]
class ODEintAdapter:
class PhasePlaneInteractive:
def __init__(self, model):
self.model = model
def dfun(self, state, t):
return model.dfun(state)
def plot_trajectory(self,x0):
tspan = np.linspace(0, 200, TRAJ_STEPS)
ys = odeint(self.model.dfun, x0, tspan)
self.ax.plot(ys[:,0], ys[:,1], 'b-') # path
self.ax.plot([ys[0,0]], [ys[0,1]], 'o') # start
self.ax.plot([ys[-1,0]], [ys[-1,1]], 's') # end
def onclick(self, event):
self.plot_trajectory([event.xdata,event.ydata])
plt.draw()
if __name__ == '__main__':
model = Oscillator()
model_ode = ODEintAdapter(model)
def __call__(self):
y1 = np.linspace(-2.0, 2.0, 20)
y2 = np.linspace(-2.0, 2.0, 20)
y1 = np.linspace(-2.0, 2.0, 20)
y2 = np.linspace(-2.0, 2.0, 20)
Y1, Y2 = np.meshgrid(y1, y2)
Y1, Y2 = np.meshgrid(y1, y2)
u, v = np.zeros(Y1.shape), np.zeros(Y2.shape)
u, v = np.zeros(Y1.shape), np.zeros(Y2.shape)
NI, NJ = Y1.shape
NI, NJ = Y1.shape
for i in range(NI):
for j in range(NJ):
x = Y1[i, j]
y = Y2[i, j]
yprime = self.model.dfun([x, y])
u[i,j] = yprime[0]
v[i,j] = yprime[1]
fig, self.ax = plt.subplots()
plt.subplots_adjust(bottom=0.2)
for i in range(NI):
for j in range(NJ):
x = Y1[i, j]
y = Y2[i, j]
yprime = model.dfun([x, y])
u[i,j] = yprime[0]
v[i,j] = yprime[1]
fig, ax = plt.subplots()
plt.subplots_adjust(bottom=0.2)
def clear(event):
self.ax.clear()
Q = self.ax.quiver(Y1, Y2, u, v, color='r')
self.ax.set_xlabel('$y_1$')
self.ax.set_ylabel('$y_2$')
self.ax.set_xlim(-2.0,2.0)
self.ax.set_ylim(-2.0,2.0)
plt.draw()
def clear(event):
ax.clear()
Q = ax.quiver(Y1, Y2, u, v, color='r')
ax.set_xlabel('$y_1$')
ax.set_ylabel('$y_2$')
ax.set_xlim(-2.0,2.0)
ax.set_ylim(-2.0,2.0)
plt.draw()
def plot_trajectory(x0, model):
tspan = np.linspace(0, 200, TRAJ_STEPS)
ys = odeint(model.dfun, x0, tspan)
ax.plot(ys[:,0], ys[:,1], 'b-') # path
ax.plot([ys[0,0]], [ys[0,1]], 'o') # start
ax.plot([ys[-1,0]], [ys[-1,1]], 's') # end
clear([0.0, 0.0])
cid = fig.canvas.mpl_connect('button_press_event', self.onclick)
plt.xlabel('$y_1$')
plt.ylabel('$y_2$')
plt.xlim([-2, 2])
plt.ylim([-2, 2])
def onclick(event):
plot_trajectory([event.xdata,event.ydata], model_ode)
plt.draw()
clear([0.0, 0.0])
cid = fig.canvas.mpl_connect('button_press_event', onclick)
plt.xlabel('$y_1$')
plt.ylabel('$y_2$')
plt.xlim([-2, 2])
plt.ylim([-2, 2])
axclear = plt.axes([0.8, 0.02, 0.1, 0.075])
bclear = Button(axclear, 'Clear')
bclear.on_clicked(clear)
plt.show()
axclear = plt.axes([0.8, 0.02, 0.1, 0.075])
bclear = Button(axclear, 'Clear')
bclear.on_clicked(clear)
plt.show()
if __name__ == '__main__':
model = Oscillator()
ppi = PhasePlaneInteractive(model)
ppi()
......@@ -10,7 +10,7 @@ ${c}=${val}${'' if loop.last else ', '}\
self.${c} = ${c}
% endfor
def dfun(self,state_variables):
def dfun(self,state_variables, *args, **kwargs):
<%
sv_csl = ", ".join(sv)
%>
......
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment