diff -r e091766f556d -r aca9c90f151c src/Stats.jl
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/src/Stats.jl	Thu Apr 25 11:16:24 2024 -0500
@@ -0,0 +1,78 @@
+__precompile__()
+
+module Stats
+
+
+########################
+# Load external modules
+########################
+
+using CSV, DataFrames
+using Statistics
+
+export calculate_statistics
+
+function calculate_statistics()
+    ImName = ("lighthouse200x300", "shepplogan256x256", "brainphantom256x256")
+    AlgName = ("dualscaling", "greedy","noprediction","primalonly","proximal","rotation","zerodual")
+    mystart = 41 # Corresponds to the 500th iterate
+    
+    # Define an array to store results
+    results = DataFrame(experiment = String[], α = Float64[], algorithm = String[], psnr_mean1 = Float64[], psnr_mean500 = Float64[], psnr_ci = String[], ssim_mean1 = Float64[], ssim_mean500 = Float64[], ssim_ci = String[])
+
+    for imname in ImName
+        for algname in AlgName
+            directory_path = "./img/"
+            files = readdir(directory_path)
+            filtered_files = filter(file -> startswith(file, "$(imname)_pdps_known_$(algname)") && endswith(file, ".txt"), files)
+
+            for file in filtered_files
+                filename = directory_path * file
+                data = CSV.File(filename, delim='\t', header=2) |> DataFrame
+
+                # Extract α from filename
+                α, _, _ = extract_parameters(filename)
+
+
+                # Extract SSIM and PSNR columns starting from 1st iteration
+                ssim_values1 = Float64.(data[:, :ssim])
+                psnr_values1 = Float64.(data[:, :psnr])
+
+                # Extract SSIM and PSNR columns starting from 500th iteration
+                ssim_values500 = Float64.(data[mystart:end, :ssim])
+                psnr_values500 = Float64.(data[mystart:end, :psnr])
+
+                # Calculate mean and confidence intervals
+                ssim_mean1 = round(mean(ssim_values1), digits=4)
+                psnr_mean1 = round(mean(psnr_values1), digits=4)
+
+                ssim_mean500 = round(mean(ssim_values500), digits=4)
+                psnr_mean500 = round(mean(psnr_values500), digits=4)
+                ssim_std500 = round(std(ssim_values500), digits=4)
+                psnr_std500 = round(std(psnr_values500), digits=4)
+                n = length(ssim_values500)
+
+                ssim_ci_lower = round(ssim_mean500 - 1.96 * ssim_std500 / sqrt(n), digits=4)
+                ssim_ci_upper = round(ssim_mean500 + 1.96 * ssim_std500 / sqrt(n), digits=4)
+                psnr_ci_lower = round(psnr_mean500 - 1.96 * psnr_std500 / sqrt(n), digits=4)
+                psnr_ci_upper = round(psnr_mean500 + 1.96 * psnr_std500 / sqrt(n), digits=4)
+
+                ssim_ci = "$(ssim_ci_lower) - $(ssim_ci_upper)"
+                psnr_ci = "$(psnr_ci_lower) - $(psnr_ci_upper)"
+                experiment = "$(imname)"
+                algorithm = "$(algname)"
+
+                # Append results to DataFrame
+                push!(results, (experiment, α, algorithm, psnr_mean1, psnr_mean500, psnr_ci, ssim_mean1, ssim_mean500, ssim_ci))
+            end
+        end
+    end
+    sort!(results, [:experiment, :α])
+    csv_path = "./img/summarystats.csv"
+    if isfile(csv_path)
+    	rm(csv_path)
+    end
+        CSV.write(csv_path, results)
+end
+
+end