import sys
import numpy as np
import pandas as pd
from scipy.optimize import minimize
from scipy.integrate import solve_ivp
import subprocess
import base64
import os

tst = 0
#=================================Input Begin===========================
if tst: print('%Test len(sys.argv),', str(len(sys.argv))+'; '
  "sys.argv:\n", ' ', *sys.argv)

tex1 = " *** D'Ambrosio ODEs ***"
tex2 = '(Feb-2024, MC)'; tex3 = '21_21:49'
print(tex1+(81-len(tex1)-len(tex2))*' '+tex2+'\n'+73*' '+tex3)
user_name = os.getenv("USER")
internetprefix = "file://" if user_name is None else ''
if tst: print ('%Test user_name:', user_name)
#=======================================================================
tini, tfin, rates, params, inivals, meth, npoints, ishow = sys.argv[1:]

v = [npoints, ishow]; [npoints, ishow] = [int(i) for i in v]
v = [tini, tfin]; [tini, tfin] = [float(i) for i in v]

def strToMatrix (h_str):
   r_v = h_str.split('\n'); amat = [[float(x) for x in
     row.strip().split()] for row in r_v]; return np.array(amat)

rates  = rates.replace(',', ' '); params = params.replace(',', ' ')
inivals = inivals.replace(',', ' ')
rates  = strToMatrix(rates.strip())
params = strToMatrix(params.strip())
inivals = strToMatrix(inivals.strip())

print (' {:<12}{:>12}{:>12}    |'.format('Time in', tini, tfin))

rates_v  = np.array(np.matrix(rates)).ravel()
params_v = np.array(np.matrix(params)).ravel()
npar = params_v.size
print (' {:<36}    | {} values'.\
  format('Parameters:', params_v.size))
print (3*' ', *params_v)

ynivals_v = np.array(np.matrix(inivals)).ravel()
nval = ynivals_v.size
print (' {:<36}    | {} values'.\
  format('Init vals:', ynivals_v.size))
print (3*' ', *ynivals_v)

#if guess_vec.size != npar: print ('%Err Size mismatch.'); exit ()
print(' {:<24}{:>12}    |'.format('Method,', meth))
print(' {:<24}{:>12}    |'.format('Points,', npoints))
print(' {:<24}{:>12}    |'.format('Show,', ('No', 'Yes')[ishow]))
print(' '+40*'-'+'+'+39*'-')
#=================================Input End=============================

#==d e f================================================================
def ode_system (t, y_v, aa, ff):
    xx, yy, zz = y_v; zk1, zk2, zk3, zk4, zk5 = rates_v
    dx = zk1*aa*yy -zk2*xx*yy + zk3*aa*xx -2*zk4*xx*xx
    dy =-zk1*aa*yy -zk2*xx*yy + ff*zk5*aa*zz
    dz = zk3*aa*xx -zk5*aa*zz
    return dx, dy, dz

aa, ff = params_v
t_v = np.linspace(tini, tfin, npoints)
t_span = [t_v[0], t_v[-1]]
print (' {:<20}{:>7.3g}{:>9.3g}    | (integration interval)'.\
  format("'SolveIVP' 't_span',", *t_span))
print (" 'SolveIVP' 'y0':"+' '*24+'| Initial value')
print (3*' ', *ynivals_v)

sol_ivp = solve_ivp(ode_system, t_span, ynivals_v,
  method=meth, t_eval=t_v,
  args=(params_v[0], params_v[1])) #PARTICULAR

#============================Prepare plot===============================
ncols = 3
y_calc_mat = np.zeros((npoints,ncols))
for j in range(ncols):
    y_calc_mat[:,j] = sol_ivp.y[j]

#================================PLOT===================================
temporary = '/afs/ist.utl.pt/users/3/8/ist11038/web/tmp' # = '/tmp'
basedir = temporary + '/tx' + str(os.getpid())
if tst: print('%Test FILES:')
tmpfilexy_nam = basedir+'tmpfxy.dat'
tmpfilexy_wra = open(tmpfilexy_nam, 'w')
if tst: print("%Test FILE tmpfilexy_nam:\n", '  '+tmpfilexy_nam)
in_file_nam = basedir+'tmp.dat'
in_file_wra = open(in_file_nam, 'w')

imgfile_nam = basedir+'tmp.png'
if tst: print("%Test FILE imgfile_nam:\n", '  '+imgfile_nam)
imgfile_wra = open(imgfile_nam, 'w')

#===============================Display=================================
nrows = npoints
fort = ' {:>10.5g}'*(1+ncols)
for i in range(nrows):
    print (fort.format(t_v[i], *y_calc_mat[i,:]), \
      file=tmpfilexy_wra)

tmpfilexy_wra.seek(0)
gnuplot_load = \
f'''wid = 620; set term png font 'Times, 14' size wid, wid/1.618
set output '{imgfile_nam}'
set title 'Solved ODEs' offset 0, -0.4
set key center right Left font ', 11'
set timestamp font 'Helvetica-Bold, 9' offset 0, -0.25
set xlabel '{{/Times-Italic t}}' offset +9., 0.75
set ylabel '{{/Times-Italic F}}' rotate by 0 offset 1.5, 0
set xrange [*:*]; set yrange [0:*]
set tics nomirror
# + x +x, (clos/ope) 4 5 sq, 6 7 o, 8 9 tri, 10 11 tri inv
set style line 12 lc rgb 'red'   lw 1 pt 7 ps 0.5
set style line 13 lc rgb '#A36C' lw 1 pt 7 ps 0.5 # jade
set style line 14 lc rgb 'blue'  lw 1 pt 7 ps 0.5 # pt 11 ps 2
set style line 15 lc rgb 'red'  lw 3 pt 4 ps 1
set style line 16 lc rgb '#A36C' lw 3 pt 10 ps 1
set style line 17 lc rgb 'blue' lw 3 pt 3 ps 1
plot '{tmpfilexy_nam}' usi 1:2 title 'X' w linespoints ls 12, \
 '' usi 1:3 title 'Y' w linespoints ls 13, \
 '' usi 1:4 title 'Z' w linespoints ls 14
'''

if tst: print("\n%Test gnuplot_load:\n"+gnuplot_load)

in_file_wra.seek(0)
with open(in_file_nam, 'w') as pli:
    lines = gnuplot_load.splitlines(True)
    pli.writelines(lines); pli.flush()
    subprocess.call(['/usr/bin/gnuplot', in_file_nam], shell=False)

imgfile_wra.seek(0)
if tst: print ("%Test To copy ... /bin/cp -f\n  " + imgfile_nam + \
  "\n  /afs/ist.utl.pt/users/3/8/ist11038/web/tmp/gnu.png")
os.system("/bin/cp -f " + imgfile_nam + \
  " /afs/ist.utl.pt/users/3/8/ist11038/web/tmp/gnu.png")
if tst: print ("%Test ... From copy\n")

with open(imgfile_nam, 'rb') as image_file:
    encoded = base64.b64encode(image_file.read())
encoded = str(encoded)[2:-1]
os.remove(in_file_nam)
os.remove(imgfile_nam); os.remove(tmpfilexy_nam)

print('<center><!--PLOT:-->')
print('<img alt="Plot" src="data:image/png; base64, '+encoded+'" />')
print('</center>')

#===============================Display=================================
if ishow==0:
    print(' End of program'); exit ()

from openpyxl import Workbook
workbook = Workbook(); workbook.remove(workbook['Sheet'])#it's keyword !

import random
mysheet = ''.join(random.sample('BCDFGHJKLMNPQRSTVWXZ', 3)) + '_' + \
          ''.join(random.sample('01234567890', 2))
sht = workbook.create_sheet(title=mysheet); sht = workbook.active

print(' Show_i'+'_'*12+'t'+('_'*6+'y_exp.')*3+('_'*5+'y_calc.')*3)
fort = ' {:>5})' + '{:>12.4g}'*(1+2*ncols)
for i in range(nrows):
    print (fort.format(i, t_v[i], *y_calc_mat[i,:]))
    row =             (i, t_v[i], *y_calc_mat[i,:])
    sht.append(row)
excel_file = '/afs/ist.utl.pt/users/3/8/ist11038/web/tmp/output.xlsx'
workbook.save(excel_file)

print ('\n <a href=\n' + \
  ' "http://web.tecnico.ulisboa.pt/ist11038/tmp/output.xlsx"\n>' + \
  '<b>output.xlsx</b></a> -- Your sheetname: "' + mysheet + \
  '" (Verify)')

print('\n End of program')