#!/usr/bin/env python

import numpy as np
import matplotlib.pyplot as plt
import matplotlib.ticker as mtick
from tkinter import Tk
from tkinter.filedialog import askopenfilename, asksaveasfile, askdirectory
import pandas as pd
import tkinter as tk
from scipy.signal import find_peaks


def error_plot(folder,t_step,r_criteria,save):
    # Load iterations and residual error
    histor = folder + '/histor.dat'
    input = open(histor, 'r')
    output = open(folder + "/histor_better.dat", 'w')
    output.writelines(line.strip() +'\n' for line in input)
    input.close()
    output.close()
    error_info = pd.read_csv(folder + "/histor_better.dat", sep=' ', header=None, usecols=(0,1,2))

    # Select only last iteration of residual error
    error=[] 
    for n in range(1,error_info.shape[0]):
        if (error_info[0][n]>error_info[0][n-1]):
            error.append(error_info[2][n-1])

    time = np.linspace(start = t_step, stop = len(error)*t_step, num = len(error))
    
    # Plots of interest
    # Liniar Scale
    plt.figure()  
    plt.plot(time,error)
    plt.plot(time,r_criteria*np.ones(len(error)),'r')
    plt.ylabel('Residual error')
    plt.xlabel('Time steps')
    plt.title('Last nonlinear residual error for each time step')
    plt.grid(True)
    if save: plt.savefig(plt_folder + case + '_Last_nonlin_res_error.pdf')

    # Semilog scale
    plt.figure()  
    plt.semilogy(time,error)
    plt.semilogy(time,r_criteria*np.ones(len(error)),'r')
    plt.ylabel('Residual error')
    plt.xlabel('Time steps')
    plt.title('Log - Last nonlinear residual error for each time step')
    plt.grid(True)
    if save: plt.savefig(plt_folder + case + '_Log_Last_nonlin_res_error.pdf')
    

def periodicity(project,folder,dt,T_cyc,n_cyc):
    pressure = np.loadtxt(folder+'/PHistRCR.dat',skiprows=2,)
    time = np.linspace(0,T_cyc*n_cyc,round(T_cyc/dt*n_cyc))
    peak_P = []
    peak_P_pos = []
    Nc = round(T_cyc/dt)
    for i in range(0,n_cyc):
        peak_P.append(np.amax(pressure[i*Nc:Nc*(i+1),-1])/1333.22)
        peak_P_pos.append(np.where(pressure[i*Nc:Nc*(i+1),-1] == np.amax(pressure[i*Nc:Nc*(i+1),-1]))[0][0]+Nc*i)
    
    peak_Pdiff = [peak_P[n]-peak_P[n-1] for n in range(1,len(peak_P))]
    peak_Pdiff = list(map(abs, peak_Pdiff)) 
    
    fig, ax = plt.subplots()
    ax.plot(time,pressure[:,-1]/1333.22)
    ax.plot(time[peak_P_pos], peak_P,'ro',label='Cylce pike')
    ax.set(xlabel='time [s]', ylabel='Pressure [mmHg]',
        title='Pressure @ last outlet')
    ax.spines['right'].set_visible(False)
    ax.spines['top'].set_visible(False)
    ax.legend(loc=0)
    plt.show()
    
    if (peak_Pdiff[-1]<=1): print('The numerical simulation \'{0}\' has achieve periodicity!\nSystolic Blood Pressure (SBP):\nsecond-last cycle = {1:.2f} mmHg,\nlast cycle = {2:.2f} mmHg,\n\u0394mmHg = {3:.2f} mmHg'.format(project,peak_P[-2],peak_P[-1],peak_Pdiff[-1]))  


def pressure(folder,N_ts,T_cyc,dt,n_cyc):
    pressure = np.loadtxt(folder+'/PHistRCR.dat',skiprows=2,)
    Nc = round(T_cyc/dt)
    time = np.linspace(0,T_cyc,Nc)
    fig, ax = plt.subplots()
    for i in range(0,pressure.shape[1]):
        ax.plot(time,pressure[N_ts-Nc:N_ts,-i]/1333.22,label='ROI-'+str(i+2))
        
    ax.set(xlabel='time [s]', ylabel='Pressure [mmHg]',
        title='Pressure @ each outlet')
    ax.spines['right'].set_visible(False)
    ax.spines['top'].set_visible(False)
    ax.legend(loc=0)
    plt.show()