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 = ' *** Random towards kinetics ***'
tex2 = '(Sep-2023, MC)'; tex3 = '29_18:55'
print(tex1+(81-len(tex1)-len(tex2))*' '+tex2+'\n'+73*' '+tex3)
#=======================================================================
tfinal, rates, concmu, concsigma, incumb, tol, ntrials, irepeat, \
  iseed, npoints, ishow = sys.argv[1:]

v = [tfinal, tol]; [tfinal, tol] = [float(i) for i in v]
v = [incumb, irepeat, iseed, npoints, ishow]
[incumb, irepeat, iseed, npoints, ishow] = [int(i) for i in v]
ntrials = float(ntrials); noyes = ['No', 'Yes']
print (' {:<24}{:>12}    |'.format('t: final,', tfinal))
rates     = np.array(np.matrix(rates)).ravel()
nrates = len(rates)
print (' {:<30}{:>6}    | {:>2} rates'.
  format('Kinetic rate constants:', '', nrates))
print (3*' ', *rates)
concmu    = np.array(np.matrix(concmu)).ravel()
concsigma = np.array(np.matrix(concsigma)).ravel()
nreac  = len(concmu)
print (' {:<30}{:>6}    | {:>2} reactants'.
  format('Concentrations: mu, sigma', '', nreac))
print (3*' ', *concmu)
print (3*' ', *concsigma)
if len(concsigma)-nreac: print(' %Err Size mismatch.'); exit ()
incumb -= 1
print (' {:<24}{:>6}{:>6.1%}    | (tolerance)'.\
  format('Incumbent component,', incumb, tol))
npoints += 1
form = ' {:<14}{:>10.2e}{:>12}    | Repeatable: {:<3}{:6}Plot points,{:>5}'
print(form.format('Trials, seed,', ntrials, iseed, noyes[irepeat], \
  '', npoints))
print(' {:<24}{:>12}    |'.format('Show,', ishow))
print(' '+40*'-'+'+'+39*'-')
#=================================Input End=============================
def rhs(t, y):
    [aa, bb, cc] = y; zk1, zk2 = rates
    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===============================
times = np.linspace(0, tfinal, num=npoints)

y_initial = concmu
soln = solve_ivp(rhs, [0, tfinal], y_initial, t_eval=times)
t_v = soln.t
y_v = [soln.y[0], soln.y[1], soln.y[2]]
if mytest: print('%Tst len(y_v),', len(y_v))
a_mat = np.vstack([soln.y[0], soln.y[1], soln.y[2]]).T

x_v = t_v
for i in range (npoints):
    print(x_v[i], *a_mat[i,:], file=tmpfilexy_wra)

print (' {:<30}{:>6}    |'.format('Final concentrations:', ''))
y_final = a_mat[-1,:]
print(("{:>9.4g}"*nreac).format(*y_final))

tmpfilexy_wra.seek(0)
init = ', '.join(str(h) for h in concmu)

gnuplot_load = \
f'''wid = 620; set term pngcairo font 'Times, 14' size wid, wid/1.618
set output '{imgfile_nam}'
set title 'Kinetics: from {init}' 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 object 1 circle at 0, 0 size 1.0 fc rgb "navy" lw 2
set style line 12 lc rgb 'red'    lw 4 pt 19 ps 1
set style line 13 lc rgb 'blue'   lw 4 pt 3 ps 1
set style line 14 lc rgb 'orange' lw 4 pt 4 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 'A' w lines ls 12, \
 '' usi 1:3 title 'B' w lines ls 13, \
 '' usi 1:4 title 'C' 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)
###s = os.system('/usr/bin/wc '+imgfile)

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

#=============================Simulation================================
tex = 'Simulation'
print (' '*(30-len(tex)//2) + '=== === === {} === === ==='.format(tex))
target = y_final[incumb]
print(' {:<24}{:>12.4g}    | (from previous)'.\
  format('Target value,', target))

#if iseed and irepeat==1: np.random.seed(iseed)
if irepeat==1: np.random.seed(iseed)

y_summary = np.empty((int(ntrials),nreac))
if mytest: print ('%Tst np.shape(y_summary):', np.shape(y_summary))

spec_low, spec_high = (1 - tol) * target, (1 + tol) * target
good_val = 0
for irun in range(int(ntrials)):
    y_initial = np.random.default_rng().normal(concmu, concsigma)
    soln = solve_ivp(rhs, [0, tfinal], y_initial, t_eval=times)
    a_mat = np.vstack([soln.y[0], soln.y[1], soln.y[2]]).T
    y_summary[irun] = a_mat[-1,:]
    x = y_summary[irun,incumb]
    if (x - spec_low) * (spec_high - x) > 0: good_val += 1
good_fraction = good_val / ntrials
print(' {:<24}{:>12.2%}    | within {:<6.1%}'.\
  format('Acceptable fraction,', good_fraction, tol))

if mytest: print ('%Tst last y_summary:', y_summary[-1,:])
nbins = npoints
hist, bin_edges = np.histogram(y_summary[:,incumb], bins=nbins, \
  density=True)

if mytest:
    print ('%Tst bin_edges:', end='')
    for x in bin_edges: print(f'{x:8.3g}', end='')
    print ("\n%Tst hist:     ", end='')
    for x in hist: print(f'{x:8.3g}', end='')

tmpfilexy_wra.close()
tmpfilexy_wra = open(tmpfilexy_nam, 'w')
for i in range (nbins):
    print(bin_edges[i], hist[i], file=tmpfilexy_wra)

max_y = np.max(hist)
if mytest: print ('%Tst max hist,', max_y)
max_y *= 1.05

xleft, xright = 0.6*y_final[incumb], (2-0.6)*y_final[incumb]
tmpfilexy_wra.seek(0)
gnuplot_load = \
f'''wid = 620
set term pngcairo dashed font 'Times, 14' size wid, wid/1.618
set output '{imgfile_nam}'
set title 'Incumbent component: {incumb+1}' offset 0, -0.3
set key left top Right
set timestamp font 'Helvetica-Bold, 9' offset 0, 0
set xlabel '{{/Times-Italic c}}' offset +10., 0.75
set ylabel 'freq.' offset 1.5, 0
set xrange [{xleft}:{xright}]
set yrange [0:{max_y}]
#stats '{tmpfilexy_nam}' using 1:2 nooutput; max_y = STATS_max_y
#set object 1 circle at 0, 0 size 1.0 fc rgb "navy" lw 2
#dashtype 2 - -  3 . . 4 -- . 5 -- .. 7=1 8=3 9=4 etc.
set arrow from {spec_low},  0 to {spec_low}, {max_y} nohead \
  linecolor "black" lw 2 dashtype 4
set arrow from {spec_high}, 0 to {spec_high}, {max_y} nohead \
  linecolor "black" lw 2 dashtype 4
set style line 12 lc rgb '#006400'  lw 2 pt 6 ps 0.5
set style line 13 lc rgb 'blue'   lw 4 pt 3 ps 1
set style line 14 lc rgb 'orange' lw 4 pt 4 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 'f' w linespoints ls 12
'''

if mytest: print("%Tst gnuplot_load:\n"+gnuplot_load)

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/gnu2.png")
os.system("/bin/cp -f " + imgfile_nam + \
  " /afs/ist.utl.pt/users/3/8/ist11038/web/tmp/gnu2.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 a_mat[i,:]: print(f"{x:12.5g}", end="")
        print() #print(f"{n:4d}", end="")
print(' End of program')