# -*- coding: utf-8 -*-
"""
Created on Tue Apr 28 12:12:36 2020

@author: hsauro
"""

import tellurium as te
import numpy as np
import math
from scipy.interpolate import CubicSpline
import matplotlib.pyplot as plt

def convertNumber (x):
    ax = math.fabs (x)
    if x >= 0:
        return ' + ' + format (ax, "f")
    else:
        return ' - ' + format (ax, "f") 
    
xm = np.array ([0.0, 1, 2, 3, 4, 5, 6, 8])
ym = np.array ([1, 0.25, 1.2, 0.65, 1, 0.85, 0.75, 0.8])

nData = len (xm);

cs = CubicSpline (xm, ym)

# Check it works, plot the spline on its own
x = []; y = []
xvalue = 0
for i in range (160):
    x.append (xvalue) 
    y.append (cs (xvalue))
    xvalue = xvalue + 0.05
plt.plot (x, y)

strequ = 'Xo := piecewise (';
t1 = ''
for k in range (nData-1):
    hstep = (xm[k+1] - xm[k])/20;
    # Optimized cubic equation
    delay = '(time-' + str (xm[k]) + ')';

    coef0s = convertNumber (cs.c[3][k]);
    coef1s = convertNumber (cs.c[2][k]);
    coef2s = convertNumber (cs.c[1][k]);

    t2 = '(('  + format (cs.c[0][k], 'f') + '*' + delay + coef2s + ')*' + delay + coef1s + ')*' + delay + coef0s;
    t1 = t1 + t2 + ', (time >=' + str (xm[k])+ ') && ' + '(time <= ' + str (xm[k+1]) + ')';
    # Continuation character
    if k != nData - 2:
       t1 = t1 + ',\ ' + '\n'
strequ = strequ + t1 + ');'
print (strequ)

r = te.loada ('''

    $Xo -> S1; Xo;
    S1 -> $X1; S1;      

Xo := piecewise (((-0.835435*(time-0.0) + 3.356306)*(time-0.0) - 3.270871)*(time-0.0) + 1.000000, (time >=0.0) && (time <= 1.0),\ 
((-0.835435*(time-1.0) + 0.850000)*(time-1.0) + 0.935435)*(time-1.0) + 0.250000, (time >=1.0) && (time <= 2.0),\ 
((0.977177*(time-2.0) - 1.656306)*(time-2.0) + 0.129129)*(time-2.0) + 1.200000, (time >=2.0) && (time <= 3.0),\ 
((-0.673272*(time-3.0) + 1.275224)*(time-3.0) - 0.251953)*(time-3.0) + 0.650000, (time >=3.0) && (time <= 4.0),\ 
((0.315910*(time-4.0) - 0.744591)*(time-4.0) + 0.278681)*(time-4.0) + 1.000000, (time >=4.0) && (time <= 5.0),\ 
((-0.040368*(time-5.0) + 0.203139)*(time-5.0) - 0.262771)*(time-5.0) + 0.850000, (time >=5.0) && (time <= 6.0),\ 
((-0.040368*(time-6.0) + 0.082035)*(time-6.0) + 0.022403)*(time-6.0) + 0.750000, (time >=6.0) && (time <= 8.0));

''')

m = r.simulate (0., 8, 100, ['time', 'Xo', 'S1'])
r.plot()