minor change

functions.py

@@ -5,18 +5,49 @@ Created on Thu Jul 16 15:08:27 2020
 @author: Aloma Blanch
 """
 
-
+import matplotlib
 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 statistics import mean 
-from scipy.signal import find_peaks
+import matplotlib.pyplot as plt
 import matplotlib.transforms as mtransforms
 
+from math import floor
+from statistics import mean 
+from scipy.signal import find_peaks
+
+def cycle(folder,dt,N_ts,save_path):
+    pressure = np.loadtxt(folder + '/PHistRCR.dat',skiprows=1)
+    T = round(dt*N_ts,3)
+    time = np.linspace(0,T,N_ts+1)
+    # Find peaks, keep only the maximums
+    peaks, _ = find_peaks(pressure[:,-1])
+    idx = np.where(pressure[peaks,-1] >= np.mean(pressure))
+    peaks_loc = peaks[idx]
+    P_peaks = pressure[peaks,-1][pressure[peaks,-1] >= np.mean(pressure)]
+    # Rmove the multiple picks found in the same location
+    idx_del  = []
+    t_pk_loc = time[peaks_loc].tolist()
+    peak_tdiff = [t_pk_loc[n]-t_pk_loc[n-1] for n in range(1,len(t_pk_loc))]
+    for i in range(0,len(peak_tdiff)):
+        if peak_tdiff[i]<= dt*10:
+            idx_del.append(i)
+            t_pk_loc[i] = []
+            peak_tdiff[i] = []    
+    t_pk_loc[:] = [x for x in t_pk_loc if x]
+    peak_tdiff[:] = [x for x in peak_tdiff if x]
+    P_peaks = np.delete(P_peaks,idx_del)
+
+    fig, ax = plt.subplots()
+    ax.plot(time,pressure[:,-1]/1333.22,'b')
+    ax.plot(t_pk_loc, P_peaks/1333.22, "or",label='peak location')
+    ax.set(xlabel='time [s]', ylabel='Pressure [mmHg]',
+       title='Pressure at last outlet')
+    ax.spines['right'].set_visible(False)
+    ax.spines['top'].set_visible(False)
+    ax.legend(loc=0)
+    # plt.show()
+    plt.savefig(save_path + '/Pressure_max_n_cycles.pdf')
+    return (floor(mean(peak_tdiff[1:])*1000)/1000,len(t_pk_loc))
 
 def error_plot(folder,t_step,r_criteria,save_path):
     # Load iterations and residual error
@@ -44,7 +75,7 @@ def error_plot(folder,t_step,r_criteria,save_path):
         title='Last nonlinear residual error for each time step')
     ax.spines['right'].set_visible(False)
     ax.spines['top'].set_visible(False)
-    plt.show()
+    # plt.show()
     plt.savefig(save_path + '/Last_nonlin_res_error.pdf')
     # Semilog scale
     fig, ax = plt.subplots()
@@ -54,14 +85,14 @@ def error_plot(folder,t_step,r_criteria,save_path):
         title='Log - Last nonlinear residual error for each time step')
     ax.spines['right'].set_visible(False)
     ax.spines['top'].set_visible(False)
-    plt.show()
+    # plt.show()
     plt.savefig(save_path + '/Log_Last_nonlin_res_error.pdf')
     fig.savefig(save_path + '/Log_Last_nonlin_res_error.jpg',dpi=150)
 
 
 def periodicity(project,folder,dt,T_cyc,n_cyc,save_path):
-    pressure = np.loadtxt(folder+'/PHistRCR.dat',skiprows=2,)
-    time = np.linspace(0,T_cyc*n_cyc,round(T_cyc/dt*n_cyc))
+    pressure = np.loadtxt(folder+'/PHistRCR.dat',skiprows=1,)
+    time = np.linspace(0,T_cyc*n_cyc,round(T_cyc/dt*n_cyc)+1)
     peak_P = []
     peak_P_pos = []
     Nc = round(T_cyc/dt)
@@ -74,24 +105,26 @@ def periodicity(project,folder,dt,T_cyc,n_cyc,save_path):
     
     fig, ax = plt.subplots()
     ax.plot(time,pressure[:,-1]/1333.22,'b')
-    ax.plot(time[peak_P_pos], peak_P,'ro',label='Cylce pike')
+    ax.plot(time[peak_P_pos], peak_P,'ro',label='Cycle pike')
     ax.set(xlabel='time [s]', ylabel='Pressure [mmHg]',
         title='Pressure at last outlet')
     ax.spines['right'].set_visible(False)
     ax.spines['top'].set_visible(False)
     ax.legend(loc=0)
-    plt.show()
+    # plt.show()
     plt.savefig(save_path + '/periodicity.pdf')
     fig.savefig(save_path + '/periodicity.jpg',dpi=150)
     
     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]))  
-        # txt = ['The numerical simulation \'{0}\' has achieve periodicity!'.format(project), 'Systolic Blood Pressure (SBP):', 'Second-last cycle = {0:.2f} mmHg'.format(peak_P[-2]), 'Last cycle = {0:.2f} mmHg'.format(peak_P[-1]), 'Delta_mmHg = {0:.2f} mmHg'.format(peak_Pdiff[-1])]
-        txt = ['The numerical simulation \'{0}\' has achieve periodicity!'.format(project), 'Systolic Blood Pressure (SBP):','Second-last cycle = {0:.2f} mmHg - Last cycle = {1:.2f} mmHg - Delta_mmHg = {2:.2f} mmHg'.format(peak_P[-2],peak_P[-1],peak_Pdiff[-1])]
+        txt = ['The numerical simulation \'{0}\' has achieved periodicity.'.format(project), 'Systolic Blood Pressure (SBP):', 'Second-last cycle = {0:.2f} mmHg'.format(peak_P[-2]), 'Last cycle = {0:.2f} mmHg'.format(peak_P[-1]), 'Delta_mmHg = {0:.2f} mmHg'.format(peak_Pdiff[-1])]
+    else:
+        print('The numerical simulation \'{0}\' has not 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]))  
+        txt = ['The numerical simulation \'{0}\' has achieved periodicity.'.format(project), 'Systolic Blood Pressure (SBP):', 'Second-last cycle = {0:.2f} mmHg'.format(peak_P[-2]), 'Last cycle = {0:.2f} mmHg'.format(peak_P[-1]), 'Delta_mmHg = {0:.2f} mmHg'.format(peak_Pdiff[-1])]
     return txt
     
 def pressure(folder,N_ts,T_cyc,dt,n_cyc,save_path):
-    pressure = np.loadtxt(folder+'/PHistRCR.dat',skiprows=2,)
+    pressure = np.loadtxt(folder+'/PHistRCR.dat',skiprows=1,)
     Nc = round(T_cyc/dt)
     time = np.linspace(0,T_cyc,Nc)
     fig, ax = plt.subplots()
@@ -109,28 +142,28 @@ def pressure(folder,N_ts,T_cyc,dt,n_cyc,save_path):
     ax.spines['right'].set_visible(False)
     ax.spines['top'].set_visible(False)
     ax.legend(loc=0)
-    plt.show()
+    # plt.show()
     plt.savefig(save_path + '/pressure.pdf')
     fig.savefig(save_path + '/pressure.jpg',dpi=150)
     return (DBP,MBP,SBP,PP)
 
 
 def flow(folder,N_ts,T_cyc,dt,n_cyc,save_path):
-    flow = np.loadtxt(folder+'/QHistRCR.dat',skiprows=2,)
+    flow = np.loadtxt(folder+'/QHistRCR.dat',skiprows=1,)
     Nc = round(T_cyc/dt)
     time = np.linspace(0,T_cyc,Nc)
     fig, ax = plt.subplots()
     Q = np.empty(flow.shape[1])
     for i in range(0,flow.shape[1]):
-        ax.plot(time,flow[N_ts-Nc:N_ts,i],label='ROI-'+str(i+2))
-        Q[i] = (mean(flow[N_ts-Nc:N_ts,i]))
+        ax.plot(time,flow[N_ts-Nc+1:N_ts+1,i],label='ROI-'+str(i+2))
+        Q[i] = (mean(flow[N_ts-Nc+1:N_ts+1,i]))
    
     ax.set(xlabel='time [s]', ylabel='Flow [mL/s]',
         title='Flow at each outlet')
     ax.spines['right'].set_visible(False)
     ax.spines['top'].set_visible(False)
     ax.legend(loc=0)
-    plt.show()
+    # plt.show()
     plt.savefig(save_path + '/flow.pdf')
     fig.savefig(save_path + '/flow.jpg',dpi=150)
     return Q
@@ -138,7 +171,10 @@ def flow(folder,N_ts,T_cyc,dt,n_cyc,save_path):
 def inlet_flow_waveform(project_folder,t_btw_rst,N_ts,dt,T_cyc,n_cyc,save_path):
     x = np.loadtxt(project_folder+'/ROI-1.flow')
     t = x[:,0]
-    Q = -x[:,1]
+    if (x[4,1] < 0):
+        Q = -x[:,1]
+    else:
+        Q = x[:,1]
     Nt_pts = np.linspace(t_btw_rst,N_ts,int(N_ts/t_btw_rst))
     t_pts = Nt_pts*dt
     # Put all the time values on a single cardiac cylce            
@@ -168,7 +204,7 @@ def inlet_flow_waveform(project_folder,t_btw_rst,N_ts,dt,T_cyc,n_cyc,save_path):
     for x, y, t in zip(t_pts[(-n_cyc-1):], Q_pts[(-n_cyc-1):], time):
         plt.text(x, y, t, transform=trans_offset, fontsize=12)    
     ax.legend(loc=0)
-    plt.show()
+    # plt.show()
     plt.savefig(save_path + '/inlet_waveform.pdf')
     fig.savefig(save_path + '/inlet_waveform.jpg',dpi=150)
     print('{0} time steps saved, available to visualize in ParaView.'.format(np.unique(np.round(t_pts,3)).shape[0]))
@@ -176,4 +212,92 @@ def inlet_flow_waveform(project_folder,t_btw_rst,N_ts,dt,T_cyc,n_cyc,save_path):
     return txt
       
 
+def rcr(project_folder):
+    rcr_lines = []
+    with open (project_folder + '/rcrt.dat', "rt") as myfile:
+        for myline in myfile:
+            rcr_lines.append(myline)
+    n_out = int((len(rcr_lines)-1)/6)
+    Rc_C_Rd = np.empty([n_out,3])
+    for i in range(0,n_out):
+        Rc_C_Rd[i][0] = float(rcr_lines[2+i*6][:-1])
+        Rc_C_Rd[i][1] = float(rcr_lines[3+i*6][:-1])
+        Rc_C_Rd[i][2] = float(rcr_lines[4+i*6][:-1])
+    return Rc_C_Rd
         
+def barPlot(project_folder,DBP,MBP,SBP,PP,Q_avg,save_path):
+    aimed_all = [1]     #No li agrada fer append en una llista buida, no recordo perquè i segur que es pot fer millor
+    f = open(project_folder+'/aimed.txt', 'r') #obrir el fitxer
+    for line in f:
+        aimed_all.append(line.split()) #fas una llista amb cada valor de la línia carregan-te els espais buits
+        
+    del aimed_all[0] #esborrar el primer valor que has posat perquè no estigues buida la llista
+        
+    def plot_bar(CFD,aimed,name,save_path,unit,decimals):
+        labels = []
+        for i in range(0,len(CFD)):
+            labels.append('ROI-'+str(i+2))
+        x = np.arange(len(labels))  # the label locations
+        width = 0.38  # the width of the bars
+        fig, ax = plt.subplots()
+        rects1 = ax.bar(x - width/2, CFD, width, label='CFD',color='k')
+        rects2 = ax.bar(x + width/2, aimed, width, label='Aimed',color='white', edgecolor='black', hatch='/')
+        ax.set_ylabel(name + unit)
+        ax.set_title(name)
+        ax.set_xticks(x)
+        ax.set_xticklabels(labels)
+        ax.legend()
+        high = max(max(aimed,CFD))
+        ax.set_ylim(top = 1.7*high)
+        def autolabel(rects,decimals):
+            """Attach a text label above each bar in *rects*, displaying its height."""
+            for rect in rects:
+                height = rect.get_height()
+                ax.annotate(decimals.format(height),
+                            xy=(rect.get_x() + rect.get_width() / 2, height),
+                            xytext=(0,3),  # 3 points vertical offset
+                            textcoords="offset points",
+                            ha='center', va='bottom',
+                            fontsize = 7.5)
+    
+        autolabel(rects1,decimals)
+        autolabel(rects2,decimals)
+        # Create labels
+        err = []
+        zip_object = zip(CFD, aimed)
+        for x1_i, x2_i in zip_object:
+            err.append(round((x1_i-x2_i)/x2_i*100,2))
+ 
+        # Text on the top of each barplot
+        for i in range(len(x)):
+            plt.text(x = x[i]- width/2, y = max(CFD[i],aimed[i])+0.15*max(max(CFD,aimed)), s = str(err[i])+'%', fontsize = 8, color='r')
+
+        fig.tight_layout()
+        # plt.show()
+        plt.savefig(save_path + '/'+name+'_bar.pdf')
+        fig.savefig(save_path + '/'+name+'_bar.jpg',dpi=150)        
+    
+    # DBP
+    CFD = [float(i) for i in DBP]
+    aimed = [float(i) for i in aimed_all[0]]
+    plot_bar(CFD,aimed,'DBP',save_path,' [mmHg]','{0:.1f}')
+    
+    # MBP
+    CFD = [float(i) for i in MBP]
+    aimed = [float(i) for i in aimed_all[1]]
+    plot_bar(CFD,aimed,'MBP',save_path,' [mmHg]','{0:.1f}')
+    
+    # SBP
+    CFD = [float(i) for i in SBP]
+    aimed = [float(i) for i in aimed_all[2]]
+    plot_bar(CFD,aimed,'SBP',save_path,' [mmHg]','{0:.1f}')
+    
+    # PP
+    CFD = [float(i) for i in PP]
+    aimed = [float(i) for i in aimed_all[3]]
+    plot_bar(CFD,aimed,'PP',save_path,' [mmHg]','{0:.1f}')
+    
+    # Q_avg
+    CFD = [float(i) for i in Q_avg]
+    aimed = [float(i) for i in aimed_all[4]]
+    plot_bar(CFD,aimed,'Q_avg',save_path,' [mL/s]','{0:.2f}')

generatePDF.py

@@ -7,7 +7,7 @@ Created on Fri Jul 17 05:34:59 2020
 
 from fpdf import FPDF
 
-def generatePDF(save_path,project,DBP,MBP,SBP,PP,Q_avg,txt1,txt2):
+def generatePDF(save_path,project,DBP,MBP,SBP,PP,Q_avg,txt1,txt2,Rc_C_Rd,T_cyc,n_cyc,N_ts,dt):
     class PDF(FPDF):
         def header(self):
             # Arial bold 15
@@ -78,24 +78,31 @@ def generatePDF(save_path,project,DBP,MBP,SBP,PP,Q_avg,txt1,txt2):
             self.add_page()
             self.section_title(num, title)
 
-    title = 'Project name: '+ project
+    title = 'Simulation Job name: '+ project
     pdf = PDF()
     pdf.set_title(title)
     pdf.set_author('Aloma Blanch Granada')
 
-    pdf.print_section(1, 'Cehcking convergency and periodicity')
-    pdf.image(save_path +'/Log_Last_nonlin_res_error.jpg', x = None, y = None, w = 140, h = 100, type = '', link = '')
-    pdf.image(save_path +'/periodicity.jpg', x = None, y = None, w = 140, h = 100, type = '', link = '')
-    pdf.cell(0, 10,txt1[0], 0, 1)
-    pdf.cell(0, 10,txt1[1], 0, 1)
-    pdf.cell(0, 10,txt1[2], 0, 1)
+    pdf.print_section(1, 'Checking convergency and periodicity')
+    pdf.image(save_path +'/Log_Last_nonlin_res_error.jpg', x = 40, y = None, w = 140, h = 100, type = '', link = '')
+    pdf.image(save_path +'/periodicity.jpg', x = 40, y = None, w = 140, h = 100, type = '', link = '')
+    pdf.set_font('Arial', '', 10)
+    pdf.cell(0, 10,' ', 0, 1)
+    pdf.cell(0, 5,txt1[0], 0, 1)
+    pdf.cell(0, 5,txt1[1], 0, 1)
+    pdf.cell(0, 5,txt1[2], 0, 1)
+    pdf.cell(0, 5,txt1[3], 0, 1)
+    pdf.cell(0, 5,txt1[4], 0, 1)
 
-    pdf.print_section(2, 'Cehcking Pressures at each outlet')
-    pdf.image(save_path +'/pressure.jpg', x = None, y = None, w = 140, h = 100, type = '', link = '')
+    pdf.print_section(2, 'Checking Pressures at each outlet')
+    pdf.image(save_path +'/pressure.jpg', x = 40, y = None, w = 140, h = 100, type = '', link = '')
 
     width_cell=[20,30,30,30,30,30,30];
     # pfd.SetFillColor(193,229,252); # Background color of header
+    pdf.set_font('Times', '', 10)
     # Header starts 
+    pdf.cell(0, 25,' ', 0, 1)
+    pdf.set_x(35)
     pdf.cell(width_cell[0],10,'ROI',1,0,'C') # First header column 
     pdf.cell(width_cell[1],10,'DBP [mmHg]',1,0,'C') # Second header column
     pdf.cell(width_cell[2],10,'MBP [mmHg]',1,0,'C') # Third header column 
@@ -103,6 +110,7 @@ def generatePDF(save_path,project,DBP,MBP,SBP,PP,Q_avg,txt1,txt2):
     pdf.cell(width_cell[4],10,'PP [mmHg]',1,1,'C') # Fourth header column
     # Rows
     for i in range(0,len(SBP)):
+        pdf.set_x(35)
         pdf.cell(width_cell[0],10,'ROI-'+ str(i+2),1,0,'C') # First column of row 1 
         pdf.cell(width_cell[1],10,str(round(DBP[i],2)),1,0,'C') # Second column of row 1 
         pdf.cell(width_cell[2],10,str(round(MBP[i],2)),1,0,'C') # Third column of row 1 
@@ -110,21 +118,64 @@ def generatePDF(save_path,project,DBP,MBP,SBP,PP,Q_avg,txt1,txt2):
         pdf.cell(width_cell[4],10,str(round(PP[i],2)),1,1,'C') # Fourth column of row 1 
 
 
-    pdf.print_section(3, 'Cehcking Flow Rate at each outlet')
-    pdf.image(save_path +'/flow.jpg', x = None, y = None, w = 140, h = 100, type = '', link = '')
+    pdf.print_section(3, 'Checking Flow Rate at each outlet')
+    pdf.image(save_path +'/flow.jpg', x = 40, y = None, w = 140, h = 100, type = '', link = '')
 
-    width_cell=[20,60];
+    width_cell=[20,45];
     # pfd.SetFillColor(193,229,252); # Background color of header
+    pdf.set_font('Times', '', 10)
     # Header starts 
+    pdf.cell(0, 25,' ', 0, 1)
+    pdf.set_x(80)
     pdf.cell(width_cell[0],10,'ROI',1,0,'C') # First header column 
     pdf.cell(width_cell[1],10,'Average Flow rate [mL/s]',1,1,'C') # Second header column
     # Rows
     for i in range(0,len(Q_avg)):
+        pdf.set_x(80)
         pdf.cell(width_cell[0],10,'ROI-'+ str(i+2),1,0,'C')
-        pdf.cell(width_cell[1],10,str(round(Q_avg[i],2)),1,1,'C')
+        pdf.cell(width_cell[1],10,str(round(Q_avg[i],3)),1,1,'C')
 
+    pdf.print_section(4, 'Checking CFD values vs aimed')
+    pdf.image(save_path +'/DBP_bar.jpg', x = 40, y = 50, w = 140, h = 100, type = '', link = '')
+    pdf.image(save_path +'/MBP_bar.jpg', x = 40, y = 170, w = 140, h = 100, type = '', link = '')
+    pdf.print_section(4, 'Checking CFD values vs aimed')
+    pdf.image(save_path +'/SBP_bar.jpg', x = 40, y = 50, w = 140, h = 100, type = '', link = '')
+    pdf.image(save_path +'/PP_bar.jpg', x = 40, y = 170, w = 140, h = 100, type = '', link = '')
+    pdf.print_section(4, 'Checking CFD values vs aimed')
+    pdf.image(save_path +'/Q_avg_bar.jpg', x = 40, y = 50, w = 140, h = 100, type = '', link = '')
+    pdf.cell(0, 10,'', 0, 1)
+    
+    pdf.print_section(5, 'Checking Inlet Flow Waveform and time steps saved')
+    pdf.image(save_path +'/inlet_waveform.jpg', x = 40, y = 50, w = 140, h = 100, type = '', link = '')
+    pdf.set_font('Arial', '', 10)
+    pdf.cell(0, 5,txt2, 0, 1)
+    
+    width_cell=[30,20,110];
+    pdf.print_section(6, 'Outlet Boundary Condition Applied')
+    pdf.set_font('Arial', '', 10)
+    pdf.cell(0, 10,'3-Element Windkessel model outlet boundary condition applied in SimVascular', 0, 1)
+    pdf.cell(0, 5,' ', 0, 1)
+    pdf.set_x(50)
+    pdf.set_font('Times', '', 10)
+    pdf.cell(width_cell[2],10,'3-Element Windkessel',1,1,'C') # First header column
+    pdf.set_x(50)
+    pdf.cell(width_cell[1],10,'ROI',1,0,'C') # Second header column
+    pdf.cell(width_cell[0],10,'Rc [g/cm^4 s]',1,0,'C') # Second header column
+    pdf.cell(width_cell[0],10,'C [cm^4 s^2/g]',1,0,'C') # Second header column
+    pdf.cell(width_cell[0],10,'Rd [g/cm^4 s]',1,1,'C') # Second header 
+    # Rows
+    for i in range(0,Rc_C_Rd.shape[0]):
+        pdf.set_x(50)
+        pdf.cell(width_cell[1],10,'ROI-'+ str(i+2),1,0,'C') # First column of row 1 
+        pdf.cell(width_cell[0],10,str(Rc_C_Rd[i][0]),1,0,'C') # Second column of row 1 
+        pdf.cell(width_cell[0],10,str(Rc_C_Rd[i][1]),1,0,'C') # Third column of row 1 
+        pdf.cell(width_cell[0],10,str(Rc_C_Rd[i][2]),1,1,'C') # Third column of row 1
+
+    pdf.print_section(7, 'Solver Parameters')
+    pdf.set_font('Arial', '', 10)
+    pdf.cell(0, 10,'Cardiac cycle period - ' + str(T_cyc) + 's', 0, 1)
+    pdf.cell(0, 10,'Number of cycles - ' + str(n_cyc), 0, 1)
+    pdf.cell(0, 10,'Number of time steps - ' + str(N_ts), 0, 1)
+    pdf.cell(0, 10,'Time step size - ' + str(dt) + 's', 0, 1)
     
-    pdf.print_section(4, 'Cehcking Inlet Flow Waveform and time steps saved')
-    pdf.image(save_path +'/inlet_waveform.jpg', x = None, y = None, w = 140, h = 100, type = '', link = '')
-    pdf.cell(0, 10,txt2, 0, 1)
     pdf.output(save_path + '/' + project + '-report.pdf', 'F')

main.py

@@ -5,21 +5,22 @@ Created on Thu Jul 16 15:08:27 2020
 @author: Aloma Blanch
 """
 
-import numpy as np
-import matplotlib.pyplot as plt
-import matplotlib.ticker as mtick
+# 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 tkinter.filedialog import askopenfilename, asksaveasfile, askdirectory
+from tkinter.filedialog import askdirectory
+# import pandas as pd
+# import tkinter as tk
 import os.path
-from scipy.signal import find_peaks_cwt
-from scipy import signal
-import statistics
-from fpdf import FPDF
+# from scipy.signal import find_peaks_cwt
+# from scipy import signal
+# import statistics
+# from fpdf import FPDF
 
 
-from functions import error_plot, periodicity, pressure, flow, inlet_flow_waveform
+from functions import error_plot, periodicity, pressure, flow, inlet_flow_waveform, rcr, cycle, barPlot
 from generatePDF import generatePDF  
 
 
@@ -32,11 +33,6 @@ save_path = folder+'/'+project+'-report'
 os.mkdir(save_path)
 save_pdf = save_path + '/' + project + '-report.pdf'
 
-# Input parameters
-T_cyc = 0.477
-n_cyc = 5
-
-
 # Read important parameters from - solver.inp file
 mylines = []                             # Declare an empty list named mylines.
 with open (project_folder + '/solver.inp', 'rt') as myfile: # Open lorem.txt for reading text data.
@@ -52,6 +48,12 @@ rc = float(mylines[26][18:-1])
 # Imesteps between Restarts - idx 6
 t_btw_rst = int(mylines[6][37:-1])
 
+# Extracting Outlet Boundary Conditions from - rcrt.dat fle
+Rc_C_Rd = rcr(project_folder)
+
+# Extracting number of cycles and period of cardiac cycle
+(T_cyc,n_cyc) = cycle(folder,dt,save_path)
+T_cyc=0.477    
 # Cehcking convergency and periodicity
 error_plot(folder,dt,rc,save_path)
 txt1 = periodicity(project,folder,dt,T_cyc,n_cyc,save_path)
@@ -65,6 +67,8 @@ txt1 = periodicity(project,folder,dt,T_cyc,n_cyc,save_path)
 # Inlet Flow Rate + and t saved
 txt2 = inlet_flow_waveform(project_folder,t_btw_rst,N_ts,dt,T_cyc,n_cyc,save_path)
 
+# Extract bar plots CFD vs. aimed
+barPlot(project_folder,DBP,MBP,SBP,PP,Q_avg,save_path)
 
 # Create PDF report
-generatePDF(save_path,project,DBP,MBP,SBP,PP,Q_avg,txt1,txt2)
+generatePDF(save_path,project,DBP,MBP,SBP,PP,Q_avg,txt1,txt2,Rc_C_Rd,T_cyc,n_cyc,N_ts,dt)

report_main.py

@@ -0,0 +1,99 @@
+ #!/usr/bin/env python
+"""
+Created on Thu Jul 16 15:08:27 2020
+
+@author: Aloma Blanch
+"""
+
+# 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
+from tkinter.filedialog import askdirectory
+# import pandas as pd
+# import tkinter as tk
+import os.path
+# from scipy.signal import find_peaks_cwt
+# from scipy import signal
+# import statistics
+# from fpdf import FPDF
+
+
+from functions import error_plot, periodicity, pressure, flow, inlet_flow_waveform, rcr, cycle, barPlot
+from generatePDF import generatePDF  
+
+
+# Selct dir 
+Tk().withdraw()
+folder = askdirectory()
+project_folder = os.path.dirname(folder)
+project = os.path.basename(project_folder)
+save_path = folder+'/'+project+'-report'
+os.mkdir(save_path)
+save_pdf = save_path + '/' + project + '-report.pdf'
+
+# Read important parameters from - solver.inp file
+mylines = []                             # Declare an empty list named mylines.
+with open (project_folder + '/solver.inp', 'rt') as myfile: # Open lorem.txt for reading text data.
+    for myline in myfile:                # For each line, stored as myline,
+        mylines.append(myline)           # add its contents to mylines.
+# Number of Timesteps - idx 3
+# Idea: remove the text to extract the number, the text part will be the same no matter the simulation
+N_ts = int(mylines[3][20:-1])
+# Time Step Size - idx 4
+dt = float(mylines[4][16:-1])
+# Residual criteria - idx 4
+rc = float(mylines[26][18:-1])
+# Imesteps between Restarts - idx 6
+t_btw_rst = int(mylines[6][37:-1])
+
+# Extracting Outlet Boundary Conditions from - rcrt.dat fle
+Rc_C_Rd = rcr(project_folder)
+
+# Extracting number of cycles and period of cardiac cycle
+(T_cyc,n_cyc) = cycle(folder,dt,N_ts,save_path)
+T_cyc=0.4769
+n_cyc=4   
+# Cehcking convergency and periodicity
+error_plot(folder,dt,rc,save_path)
+txt1 = periodicity(project,folder,dt,T_cyc,n_cyc,save_path)
+
+# Pressure - Outlets
+(DBP,MBP,SBP,PP) = pressure(folder,N_ts,T_cyc,dt,n_cyc,save_path)
+ 
+# Flow Rate - Outlets
+(Q_avg) = flow(folder,N_ts,T_cyc,dt,n_cyc,save_path)
+
+# Inlet Flow Rate + and t saved
+txt2 = inlet_flow_waveform(project_folder,t_btw_rst,N_ts,dt,T_cyc,n_cyc,save_path)
+
+# Extract bar plots CFD vs. aimed
+barPlot(project_folder,DBP,MBP,SBP,PP,Q_avg,save_path)
+
+# # Flow comparison specific for patient 120.
+# # Only ROI-5,6,8 because that is the PC-MRA data that we have.
+# def flow_comparison(folder,N_ts,T_cyc,dt,n_cyc,save_path):
+#     flow_sim = np.loadtxt(folder+'/QHistRCR.dat',skiprows=2,)
+#     flow_doc = np.loadtxt(os.path.dirname(os.path.dirname(os.path.dirname(os.path.dirname(os.path.dirname(folder)))))+'/Data/flow_waveforms.txt',skiprows=2,)
+    
+#     Nc = round(T_cyc/dt)
+#     time_sim = np.linspace(0,T_cyc,Nc)
+#     time_doc = np.linspace(0,T_cyc,flow_doc.shape[0])
+#     Q = np.empty(flow_sim.shape[1])
+#     for i in range(0,flow_sim.shape[1]):
+#         fig, ax = plt.subplots()
+#         ax.plot(time_sim,flow_sim[N_ts-Nc:N_ts,i],label='sim ROI-'+str(i+2))
+#         ax.plot(time_doc,flow_doc[:,1+i],label='doc ROI-'+str(i+2))
+   
+#     ax.set(xlabel='time [s]', ylabel='Flow [mL/s]',
+#         title='Flow at each outlet')
+#     ax.spines['right'].set_visible(False)
+#     ax.spines['top'].set_visible(False)
+#     ax.legend(loc=0)
+#     # plt.show()
+#     plt.savefig(save_path + '/flow_comparison.pdf')
+#     fig.savefig(save_path + '/flow_comparison.jpg',dpi=150)
+
+# Create PDF report
+generatePDF(save_path,project,DBP,MBP,SBP,PP,Q_avg,txt1,txt2,Rc_C_Rd,T_cyc,n_cyc,N_ts,dt)