import sys
import numpy as np
import subprocess
import base64 #import datetime
import os
from scipy.integrate import solve_ivp

mytest = 0
#=================================Input Begin===========================
if mytest: print('%Tst len(sys.argv),', str(len(sys.argv))+'; ' \
      "sys.argv:\n", ' ', sys.argv)
tex1 = ' *** Kinetic setting ***'
tex2 = '(Sep-2023, MC)'; tex3 = '12_00:25'
print(tex1+(81-len(tex1)-len(tex2))*' '+tex2+'\n'+73*' '+tex3)
#=======================================================================
file_size, file_tmp, max_file_size, tfinal, table, rates, \
  concinit, incumb, npoints, ishow = sys.argv[1:]
mode = 'Ver' # not used

v = [file_size, max_file_size, incumb, npoints, ishow]
[file_size, max_file_size, incumb, npoints, ishow] = [int(i) for i in v]
v = [tfinal]; [tfinal] = [float(i) for i in v]

print (' {:<17}{:>12}{:>7}    | File_tmp, {:<32}'.\
  format('File_size, max,', file_size, max_file_size, file_tmp))
print (' {:<24}{:>12}    |'.format('t (final),', tfinal))

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

print (' {:<36}    | Kinetic data:'.format('time____concentrations'))
table = strToMatrix(table.strip())
nrows, ncols = table[:,0].size, table[0,:].size

t_v = table[:,0]; a_m = table[:,1:]
for i in range(nrows):
    print(f" {t_v[i]:4.3g}", end="")
    for x in a_m[i,0:]: print(f"{x:8.3g}", end="")
    print()
rates = np.array(np.matrix(rates)).ravel()
nrates = len(rates)
print (' {:<36}    | {:>2} rates'.\
  format('Kinetic rate constants:', nrates))
print (3*' ', *rates)

concinit = np.array(np.matrix(concinit)).ravel()
ncomps = len(concinit)
print (' {:<36}    | {:>2} compounds'.\
  format('Concentrations -- initial:', ncomps))
print (3*' ', *concinit)
###if len(concmu) != len(concsigma): print(' %Err Size mismatch.'); exit ()

# print (' {:<24}{:>12}    |'.format('Incumbent compound,', incumb))
# dict = {"Ver": "Verify", "Opt": "Optimize"}
# print (' {:<24}{:>12}    |'.format('Computing mode,', dict[mode]))
noyes = ['No', 'Yes']
npoints += 1
print(' {:<24}{:>12}    |'.format('Show data,', noyes[ishow]))
print(' '+40*'-'+'+'+39*'-')
#=================================Input End=============================
def ode_system(t, y, par):
    aa, bb, cc = y; zk1, zk2 = par
    dadt = -zk1 * aa
    dbdt = zk1 * aa - zk2 * bb
    dcdt = zk2 * bb
    return [dadt, dbdt, dcdt]

#=============================Preparation===============================
user_name = os.getenv("USER")
if mytest: print ('%Tst User:', user_name)
internetprefix = "file://" if user_name is None else '' # not used

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

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

#=============================Integration===============================
# t_vec = np.linspace(0, tfinal, num=npoints)
t_vec = np.zeros(nrows); ncomps = ncols - 1
exp_m = np.zeros((nrows, ncomps)); calc_m = np.zeros((nrows, ncomps))
t_vec = table[:,0]
for j in range(ncomps): exp_m[:,j] = table[:,j+1]

if mode=='Opt':
    print ('%Err Not implemented.'); exit ()
    import mod_optimize
    rates = mod_optimize.adjust (t_vec, exp_m, rates)

t_last = t_vec[-1]
y_init = concinit
soln = solve_ivp(ode_system, [0, t_last], y_init, t_eval=t_vec, \
  args=(rates, ))
y_vec = [soln.y[0], soln.y[1], soln.y[2]]
if mytest: print('%Tst len(y_vec),', len(y_vec))
a_m = np.vstack([soln.y[0], soln.y[1], soln.y[2]]).T
print (' {:<36}    |'.format('Final concentrations:'))
y_fin = a_m[-1,:]
print ('  '+("{:>9.4g}"*ncomps).format(*y_fin))

for j in range(ncomps): calc_m[:,j] = soln.y[j]

if mytest: print ('%Tst table a_m .shape,', table.shape, a_m.shape)
discrep = 0
for i in range(nrows):
    for j in range(ncomps): discrep += (calc_m[i,j] - exp_m[i,j])**2
print ('discrep,', discrep)

x_v = t_vec
for i in range (nrows):
#   print (x_v[i], *a_m[i,:], file=tmpfilexy_wra)
    print (x_v[i], *exp_m[i,:], *calc_m[i,:], file=tmpfilexy_wra)

tmpfilexy_wra.seek(0)
hinit = ', '.join(str(h) for h in concinit)

gnuplot_load = \
f'''wid = 620; set term pngcairo font 'Times, 14' size wid, wid/1.618
set output '{imgfile_nam}'
set title 'Kinetics: from {hinit}' offset 0, -0.3
set key left top Right
set timestamp font 'Helvetica-Bold, 9' offset 0, 0
set xlabel '{{/Times-Italic t}}' offset +10., 0.75
set ylabel 'concentrations' offset 1.5, 0
set xrange [0:*]; set yrange [0:*]; set xtics nomirror
#set object 1 circle at 0, 0 size 1.0 fc rgb "navy" lw 2
set style line 12 lc rgb 'orange' lw 2 ps 1
set style line 13 lc rgb 'blue'   lw 2 ps 1
set style line 14 lc rgb 'black'  lw 3 ps 1
set style line 19 lc rgb 'red' lw 4 pt 19 ps 2 # pt 4 ps 3

plot '{tmpfilexy_nam}' \
    usi 1:2 title 'Ae' w linespoints ls 12 pt 19, \
 '' usi 1:5 title 'Ac' w lines ls 12, \
 '' usi 1:3 title 'Be' w linespoints ls 13 pt 3, \
 '' usi 1:6 title 'Bc' w lines ls 13, \
 '' usi 1:4 title 'Ce' w linespoints ls 14 pt 4, \
 '' usi 1:7 title 'Cc' w lines ls 14
'''
if mytest: print("%Tst 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 mytest: print ("%Tst To copy...  /bin/cp -f\n  " + imgfile_nam + \
  "\n  /afs/ist.utl.pt/users/3/8/ist11038/web/tmp/gnu1.png")
os.system("/bin/cp -f " + imgfile_nam + \
  " /afs/ist.utl.pt/users/3/8/ist11038/web/tmp/gnu1.png")
if mytest: print ('%Tst ...From copy')

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>')

if ishow:
    print(' Show_i___________x__________yy...')
    for i in range(npoints):
        print(f" {i:6d}{t_v[i]:12.5g}", end="")
        for x in calc_m[i,:]: print(f"{x:12.5g}", end="")
        print() #print(f"{n:4d}", end="")
print(' End of program')