ode integrator method plots

plot_n_link_pendulum.py

#!/usr/bin/env python
# -*- coding: utf-8 -*-
import os
import math
from itertools import groupby, count
import numpy as np
import matplotlib.pylab as plt
import seaborn as sns


def plot_results():
    methods = ['lambdify', 'theano', 'cython']

    fig = None
    axes = None
    error = None
    rows = None
    cols = None
    for method in methods:
        datafile = os.path.join(os.path.dirname(os.path.abspath(__file__)),
                                '.'.join(['n_link_pendulum_on_cart_regression',
                                          method, 'npy']))
        x = np.load(datafile)
        t = np.linspace(0, x.shape[0] - 1, x.shape[0])

        if fig is None:
            cols = 2
            rows = math.ceil(x.shape[1] / 2)
            fig, axes = plt.subplots(rows, cols, sharex=True)
            axes = axes.ravel()
            error = np.zeros(x.shape)

        if method == 'lambdify':
            error = x
        elif method == 'cython':
            error -= x

        for i in range(x.shape[1]):
            ax = axes[i]
            ax.plot(t, x[:, i], label=method)
            if method == methods[-1]:
                ax.legend()
    fig.suptitle('n-link pendulum states vs time')
    fig.savefig('benchmark-states.png')

    fig, axes = plt.subplots(rows, cols, sharex=True)
    axes = axes.ravel()
    for i in range(error.shape[1]):
        ax = axes[i]
        ax.plot(error[:, i])
        ax.title
    fig.suptitle('n-link pendulum state difference (lambdify - cython) vs time')
    fig.savefig('benchmark-error.png')

    eps = np.finfo(x.dtype).eps
    tol = 1e10*eps
    tol_str = '1e10*eps'
    for i in range(error.shape[1]):
        e = error[:, i] > tol
        if e.any():
            ind = np.nonzero(e)[0]
            ind = [ind[0], ind[-1]]
            print('for state {} simulation exceeds {} '.format(i, tol_str) +
                  'difference between indices: {}'.format(ind))

    print('{} = {}'.format(tol_str, tol))
    plt.show()


if __name__ == "__main__":
    plot_results()