src/Stats.jl

__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