{%- macro render_model(model) %}
class {{ model.name }}(BaseModel):

    def __init__(self):

        # declare constants
        {%- for const in model.c %}
        {{ const.name }} = {{ render_expression(expr=const.value) }};
        {%- endfor %}

        # declare parameters
        {%- for param in model.p %}
        {{ param.name }} = sympy.symbols('{{ param.name }}');
        {%- endfor %}
        self.p = sympy.Matrix([
        {%- for param in model.p %}
            {{ param.name }}{% if not loop.last %},{% endif %}
        {%- endfor %}]);
        self.p0 = {
        {%- for param in model.p %}
            '{{ param.name }}': {{ render_expression(expr=param.value) -}}
                {% if not loop.last %},{% endif %}
        {%- endfor %}
        };

        # declare  discrete states
        {%- for var in model.z %}
        {{ var.name }} = sympy.symbols('{{ var.name }}');
        {%- endfor %}
        self.z = sympy.Matrix([
        {%- for var in model.z %}
            {{ var.name }}{% if not loop.last %},{% endif %}
        {%- endfor %}]);
        self.z0 = {
        {%- for var in model.z %}
            '{{ var.name }}': 0{% if not loop.last %},{% endif %}
        {%- endfor %}
        };

        # declare inputs
        {%- for var in model.u %}
        {{ var.name }} = sympy.symbols('{{ var.name }}');
        {%- endfor %}
        self.u = sympy.Matrix([
        {%- for var in model.u %}
            {{ var.name }}{% if not loop.last %},{% endif %}
        {%- endfor %}]);

        # declare outputs
        {%- for var in model.y %}
        {{ var.name }} = sympy.symbols('{{ var.name }}');
        {%- endfor %}
        self.y = sympy.Matrix([
        {%- for var in model.y %}
            {{ var.name }}{% if not loop.last %},{% endif %}
        {%- endfor %}]);

        # declare continuous states
        {%- for var in model.x %}
        {{ var.name }} = sympy.symbols('{{ var.name }}');
        {%- endfor %}
        self.x = sympy.Matrix([
        {%- for var in model.x %}
            {{ var.name }}{% if not loop.last %},{% endif %}
        {%- endfor %}]);
        self.x0 = {
        {%- for var in model.x %}
            '{{ var.name }}': 0{% if not loop.last %},{% endif %}
        {%- endfor %}
        };

        # ode
        self.ode = sympy.Matrix([
        {%- for expr in model.ode %}
            {{ render_expression(expr=expr) }}{% if not loop.last %},{% endif %}
        {%- endfor %}]);

        # alg
        self.alg = sympy.Matrix([])

{%- endmacro %}

{%- macro binary(op, expr) %}
    {{- render_expression(expr=expr.lhs) }}{{ op }}{{ render_expression(expr=expr.rhs) -}}
{%- endmacro -%}

{%- macro render_expression(expr) %}
    {%- if "UnsignedReal" in expr %}
        {{- expr.UnsignedReal }}
    {%- elif "UnsignedInteger" in expr %}
        {{- expr.UnsignedInteger }}
    {%- elif "Ref" in expr %}
        {{- expr.Ref.comp.name }}
    {%- elif "Add" in expr %}
        {{- binary(op=" + ", expr=expr.Add) -}}
    {%- elif "Sub" in expr %}
        {{- binary(op=" - ", expr=expr.Sub) -}}
    {%- elif "Mul" in expr %}
        {{- binary(op="*", expr=expr.Mul) -}}
    {%- elif "Div" in expr %}
        {{- binary(op="/", expr=expr.Div) -}}
    {%- elif "Exp" in expr %}
        {{- binary(op="^", expr=expr.Exp) -}}
    {%- elif "Or" in expr -%}
        ca.logic_or(
            {{- render_expression(expr=expr.Or.lhs) -}},
            {{- render_expression(expr=expr.Or.rhs) -}})
    {%- elif "And" in expr -%}
        ca.logic_and(
            {{- render_expression(expr=expr.And.lhs) -}},
            {{- render_expression(expr=expr.And.rhs) -}})
    {%- elif "Equal" in expr -%}
            {{- render_expression(expr=expr.Equal.lhs) -}} ==
            {{- render_expression(expr=expr.Equal.rhs) -}}
    {%- endif %}
{%- endmacro -%}

{#{{ model | pprint }}#}
import sympy
import numpy as np
import scipy.integrate

class BaseModel:

    def __repr__(self):
        return repr(self.__dict__)

    def simulate(self, t=None, u=None):
        if t is None:
            t = np.arange(0, 1, 0.01)
        if u is None:
            u = np.zeros(len(t))
        if np.shape(u) == ():
            u = u*np.ones(len(t))

        p0 = np.array([self.p0[k] for k in self.p0.keys()]);
        z0 = np.array([self.z0[k] for k in self.z0.keys()]);
        x0 = np.array([self.x0[k] for k in self.x0.keys()]);
        f = sympy.lambdify([self.x, self.u, self.z, self.p], list(self.ode))
        res = scipy.integrate.solve_ivp(
            fun=lambda ti, x: f(x, [np.interp(ti, t, u)], z0, p0),
            t_span=[t[0], t[-1]],
            y0=x0,
            t_eval=t)
        return {
            't': t,
            'x': res['y'].T
        }

{{ render_model(model=model) }}