Benchmarking

using Statistics: quantile"""Compare algorithm accuracy against SPA reference.Returns DataFrame with error statistics including percentiles."""function compare_accuracy(obs::Observer, times::Vector{DateTime}, algo)    # Get SPA reference positions    spa_pos = solar_position(obs, times, SPA())    # Get algorithm positions    algo_pos = solar_position(obs, times, algo)    # Calculate errors for all positions    elev_errors = abs.(algo_pos.elevation .- spa_pos.elevation)    azim_errors = abs.(algo_pos.azimuth .- spa_pos.azimuth)    # Handle azimuth wraparound (0° and 360° are the same)    azim_errors = min.(azim_errors, 360.0 .- azim_errors)    return (        elevation_mean = mean(elev_errors),        elevation_p2_5 = quantile(elev_errors, 0.025),        elevation_p97_5 = quantile(elev_errors, 0.975),        elevation_max = maximum(elev_errors),        azimuth_mean = mean(azim_errors),        azimuth_p2_5 = quantile(azim_errors, 0.025),        azimuth_p97_5 = quantile(azim_errors, 0.975),        azimuth_max = maximum(azim_errors),        n_samples = length(times),    )end

# Algorithms to compare (excluding SPA which is the reference)algorithms = [    ("PSA", PSA()),    ("NOAA", NOAA()),    ("Walraven", Walraven()),    ("USNO", USNO()),]# Collect accuracy dataaccuracy_results = DataFrame(    Algorithm = String[],    Location = String[],    Elevation_Mean_Error = Float64[],    Elevation_P2_5 = Float64[],    Elevation_P97_5 = Float64[],    Elevation_Max_Error = Float64[],    Azimuth_Mean_Error = Float64[],    Azimuth_P2_5 = Float64[],    Azimuth_P97_5 = Float64[],    Azimuth_Max_Error = Float64[],)for (algo_name, algo) in algorithms    for loc in locations        stats = compare_accuracy(loc.obs, times, algo)        push!(accuracy_results, (            Algorithm = algo_name,            Location = loc.name,            Elevation_Mean_Error = stats.elevation_mean,            Elevation_P2_5 = stats.elevation_p2_5,            Elevation_P97_5 = stats.elevation_p97_5,            Elevation_Max_Error = stats.elevation_max,            Azimuth_Mean_Error = stats.azimuth_mean,            Azimuth_P2_5 = stats.azimuth_p2_5,            Azimuth_P97_5 = stats.azimuth_p97_5,            Azimuth_Max_Error = stats.azimuth_max,        ))    endendaccuracy_results

# Aggregate results by algorithm (mean across all locations)algo_stats = combine(    groupby(accuracy_results, :Algorithm),    :Elevation_Mean_Error => mean => :Elev_Mean,    :Elevation_P2_5 => mean => :Elev_P2_5,    :Elevation_P97_5 => mean => :Elev_P97_5,    :Azimuth_Mean_Error => mean => :Azim_Mean,    :Azimuth_P2_5 => mean => :Azim_P2_5,    :Azimuth_P97_5 => mean => :Azim_P97_5,)# Sort by algorithm orderalgo_order = ["PSA", "NOAA", "Walraven", "USNO"]algo_stats = algo_stats[sortperm([findfirst(==(algo), algo_order) for algo in algo_stats.Algorithm]), :]xs = 1:nrow(algo_stats)fig = Figure(size = (900, 400), backgroundcolor = :transparent, fontsize = 12, textcolor = "#f5ab35")# Elevation error plot with error barsax1 = Axis(fig[1, 1],    title = "Elevation Error vs SPA (95% CI)",    xlabel = "Algorithm",    ylabel = "Error (degrees)",    xticks = (xs, algo_stats.Algorithm),    backgroundcolor = :transparent,)# Error bars showing 95% intervalerrorbars!(ax1, xs, algo_stats.Elev_Mean,    algo_stats.Elev_Mean .- algo_stats.Elev_P2_5,    algo_stats.Elev_P97_5 .- algo_stats.Elev_Mean,    color = :steelblue, linewidth = 2, whiskerwidth = 10)scatter!(ax1, xs, algo_stats.Elev_Mean, color = :steelblue, markersize = 12)# Azimuth error plot with error barsax2 = Axis(fig[1, 2],    title = "Azimuth Error vs SPA (95% CI)",    xlabel = "Algorithm",    ylabel = "Error (degrees)",    xticks = (xs, algo_stats.Algorithm),    backgroundcolor = :transparent,)errorbars!(ax2, xs, algo_stats.Azim_Mean,    algo_stats.Azim_Mean .- algo_stats.Azim_P2_5,    algo_stats.Azim_P97_5 .- algo_stats.Azim_Mean,    color = :coral, linewidth = 2, whiskerwidth = 10)scatter!(ax2, xs, algo_stats.Azim_Mean, color = :coral, markersize = 12)

fig_err = Figure(size = (900, 500), backgroundcolor = :transparent, fontsize = 12, textcolor = "#f5ab35")# Convert to day of year for x-axisday_of_year = [Dates.dayofyear(t) for t in hourly_times]ax1 = Axis(fig_err[1, 1],    title = "PSA Elevation Error vs SPA (2024, London)",    xlabel = "Day of Year",    ylabel = "Error (degrees)",    backgroundcolor = :transparent,)scatter!(ax1, day_of_year, elev_errors, markersize = 1.5, color = (:steelblue, 0.5))hlines!(ax1, [0.0], color = :gray, linestyle = :dash)ax2 = Axis(fig_err[2, 1],    title = "PSA Azimuth Error vs SPA (2024, London)",    xlabel = "Day of Year",    ylabel = "Error (degrees)",    backgroundcolor = :transparent,)scatter!(ax2, day_of_year, azim_errors, markersize = 1.5, color = (:coral, 0.5))hlines!(ax2, [0.0], color = :gray, linestyle = :dash)

fig2 = Figure(size = (900, 500), backgroundcolor = :transparent, fontsize = 11, textcolor = "#f5ab35")for (i, loc) in enumerate(locations)    row = (i - 1) ÷ 2 + 1    col = (i - 1) % 2 + 1    ax = Axis(fig2[row, col],        title = loc.name,        xlabel = "Algorithm",        ylabel = "Mean Elevation Error (°)",        xticks = (1:length(algorithms), [a[1] for a in algorithms]),        backgroundcolor = :transparent,    )    loc_data = filter(r -> r.Location == loc.name, accuracy_results)    barplot!(ax, 1:length(algorithms), loc_data.Elevation_Mean_Error, color = :teal)endLabel(fig2[0, :], "Elevation Error by Location", fontsize = 14, font = :bold)

# Single position benchmarksobs = Observer(51.5074, -0.1278, 11.0)  # Londondt = DateTime(2024, 6, 21, 12, 0, 0)single_benchmarks = DataFrame(    Algorithm = String[],    Time_ns = Float64[],    Allocations = Int[],)for (name, algo) in [("PSA", PSA()), ("NOAA", NOAA()), ("Walraven", Walraven()),                      ("USNO", USNO()), ("SPA", SPA())]    b = @benchmark solar_position($obs, $dt, $algo) samples=100 evals=10    push!(single_benchmarks, (        Algorithm = name,        Time_ns = median(b.times),        Allocations = b.allocs,    ))end# Add relative timingsingle_benchmarks.Time_μs = single_benchmarks.Time_ns ./ 1000spa_time = filter(row -> row.Algorithm == "SPA", single_benchmarks).Time_ns[1]single_benchmarks.Relative_to_SPA = single_benchmarks.Time_ns ./ spa_timesingle_benchmarks[:, [:Algorithm, :Time_μs, :Allocations, :Relative_to_SPA]]

# Vector benchmarks for different sizessizes = [100, 1_000, 10_000]vector_benchmarks = DataFrame(    Algorithm = String[],    N = Int[],    Time_ms = Float64[],    Throughput = Float64[],  # positions per second)for n in sizes    times_vec = collect(DateTime(2024, 1, 1):Hour(1):(DateTime(2024, 1, 1) + Hour(n-1)))    for (name, algo) in [("PSA", PSA()), ("NOAA", NOAA()), ("Walraven", Walraven()),                          ("USNO", USNO()), ("SPA", SPA())]        b = @benchmark solar_position($obs, $times_vec, $algo) samples=10 evals=1        time_ms = median(b.times) / 1e6        push!(vector_benchmarks, (            Algorithm = name,            N = n,            Time_ms = time_ms,            Throughput = n / (time_ms / 1000),        ))    endend# Pivot for displayvector_pivot = unstack(vector_benchmarks, :Algorithm, :N, :Time_ms)vector_pivot

fig3 = Figure(size = (900, 400), backgroundcolor = :transparent, fontsize = 12, textcolor = "#f5ab35")# Scaling plot (log-log)ax1 = Axis(fig3[1, 1],    title = "Computation Time vs Input Size",    xlabel = "Number of Timestamps",    ylabel = "Time (ms)",    xscale = log10,    yscale = log10,    backgroundcolor = :transparent,)colors = [:blue, :orange, :green, :purple, :red]algo_names = ["PSA", "NOAA", "Walraven", "USNO", "SPA"]for (i, algo) in enumerate(algo_names)    data = filter(r -> r.Algorithm == algo, vector_benchmarks)    lines!(ax1, data.N, data.Time_ms, label = algo, color = colors[i], linewidth = 2)    scatter!(ax1, data.N, data.Time_ms, color = colors[i], markersize = 8)endaxislegend(ax1, position = :rb, framevisible = false, labelsize = 10)# Throughput plotax2 = Axis(fig3[1, 2],    title = "Throughput at N=10,000",    xlabel = "Algorithm",    ylabel = "Positions per Second",    xticks = (1:length(algo_names), algo_names),    backgroundcolor = :transparent,)throughput_10k = filter(r -> r.N == 10_000, vector_benchmarks)barplot!(ax2, 1:length(algo_names), throughput_10k.Throughput ./ 1e6, color = colors)ax2.ylabel = "Million Positions / Second"

# Benchmark sizessizes = [100, 1_000, 10_000]# Results storagecomparison_results = DataFrame(    Algorithm = String[],    N = Int[],    Julia_ms = Float64[],    Python_ms = Float64[],    Speedup = Float64[],)for n in sizes    # Create time vectors    julia_times_vec = collect(DateTime(2024, 1, 1):Hour(1):(DateTime(2024, 1, 1) + Hour(n-1)))    py_times_idx = create_pandas_times(n)    for (algo_name, algo) in [("PSA", PSA()), ("NOAA", NOAA()), ("Walraven", Walraven()),                               ("USNO", USNO()), ("SPA", SPA())]        # Julia benchmark        julia_bench = @benchmark solar_position($obs, $julia_times_vec, $algo) samples=5 evals=1        julia_time_ms = median(julia_bench.times) / 1e6        # Python benchmark        py_func, py_kwargs = solposx_algorithms[algo_name]        # Benchmark Python function using BenchmarkTools        if isempty(py_kwargs)            py_bench = @benchmark $py_func($py_times_idx, $lat, $lon) samples=5 evals=1        else            py_bench = @benchmark $py_func($py_times_idx, $lat, $lon; $py_kwargs...) samples=5 evals=1        end        python_time_ms = median(py_bench.times) / 1e6        speedup = python_time_ms / julia_time_ms        push!(comparison_results, (            Algorithm = algo_name,            N = n,            Julia_ms = round(julia_time_ms, digits=3),            Python_ms = round(python_time_ms, digits=3),            Speedup = round(speedup, digits=1),        ))    endendcomparison_results

fig5 = Figure(size = (600, 750), backgroundcolor = :transparent, fontsize = 12, textcolor = "#f5ab35")# Group by algorithm for plottingalgo_names = ["PSA", "NOAA", "Walraven", "USNO", "SPA"]colors_julia = [:blue, :green, :purple, :orange, :red]ax1 = Axis(fig5[1, 1],    title = "Computation Time: Julia vs Python",    xlabel = "Number of Timestamps",    ylabel = "Time (ms)",    xscale = log10,    yscale = log10,    backgroundcolor = :transparent,)# Store line objects for legendsjulia_lines = []python_lines = []for (i, algo) in enumerate(algo_names)    data = filter(r -> r.Algorithm == algo, comparison_results)    # Julia times (solid lines)    l1 = lines!(ax1, data.N, data.Julia_ms,        color = colors_julia[i],        linewidth = 2)    scatter!(ax1, data.N, data.Julia_ms, color = colors_julia[i], markersize = 6)    push!(julia_lines, l1)    # Python times (dashed lines)    l2 = lines!(ax1, data.N, data.Python_ms,        color = colors_julia[i],        linewidth = 2,        linestyle = :dash)    scatter!(ax1, data.N, data.Python_ms, color = colors_julia[i], markersize = 6,        marker = :utriangle)    push!(python_lines, l2)end# Algorithm legend (colors)leg1 = Legend(fig5[1, 2][1, 1], julia_lines, algo_names, "Algorithm",    framevisible = true, backgroundcolor = :transparent, labelsize = 10, titlesize = 11)# Line style legend (solid vs dashed)style_lines = [LineElement(color = :gray, linewidth = 2, linestyle = :solid),               LineElement(color = :gray, linewidth = 2, linestyle = :dash)]leg2 = Legend(fig5[1, 2][2, 1], style_lines, ["Julia", "Python"], "Library",    framevisible = true, backgroundcolor = :transparent, labelsize = 10, titlesize = 11)# Speedup bar chartax2 = Axis(fig5[2, 1:2],    title = "Julia vs Python at N=1,000",    xlabel = "Algorithm",    ylabel = "Speedup Factor (×)",    xticks = (1:length(algo_names), algo_names),    backgroundcolor = :transparent,)speedup_1k = filter(r -> r.N == 1_000, comparison_results)barplot!(ax2, 1:length(algo_names), speedup_1k.Speedup, color = colors_julia)hlines!(ax2, [1.0], color = :gray, linestyle = :dash)

# Combine accuracy and performance data for Pareto analysis at N=1000# Get mean elevation error across all locations for each algorithmpareto_data = combine(    groupby(accuracy_results, :Algorithm),    :Elevation_Mean_Error => mean => :Mean_Error)# Add SPA as reference (error = 0)push!(pareto_data, (Algorithm = "SPA", Mean_Error = 0.0))# Get Julia timing at N=1000julia_times_1k = filter(r -> r.N == 1_000, comparison_results)pareto_julia = leftjoin(pareto_data,    select(julia_times_1k, :Algorithm, :Julia_ms),    on = :Algorithm)pareto_julia.Library = fill("Julia", nrow(pareto_julia))# Get Python timing at N=1000pareto_python = leftjoin(pareto_data,    select(julia_times_1k, :Algorithm, :Python_ms),    on = :Algorithm)rename!(pareto_python, :Python_ms => :Julia_ms)  # Reuse column name for conveniencepareto_python.Library = fill("Python", nrow(pareto_python))# Combine bothpareto_combined = vcat(pareto_julia, pareto_python)rename!(pareto_combined, :Julia_ms => :Time_ms)# Sort by error then timesort!(pareto_combined, [:Mean_Error, :Time_ms])pareto_combined

fig_pareto = Figure(size = (700, 550), backgroundcolor = :transparent, fontsize = 12, textcolor = "#f5ab35")ax = Axis(fig_pareto[1, 1],    title = "Algorithm Selection: Accuracy vs Performance (N=1,000)",    xlabel = "Computation Time (ms)",    ylabel = "Mean Elevation Error vs SPA (°)",    backgroundcolor = :transparent,)# Color scheme for algorithmscolors_pareto = Dict(    "PSA" => :steelblue,    "NOAA" => :orange,    "Walraven" => :green,    "USNO" => :purple,    "SPA" => :red)# Marker scheme for librarymarkers_lib = Dict(    "Julia" => :circle,    "Python" => :utriangle)# Plot pointsjulia_points = []python_points = []for row in eachrow(pareto_combined)    marker = markers_lib[row.Library]    color = colors_pareto[row.Algorithm]    s = scatter!(ax, [row.Time_ms], [row.Mean_Error],            color = color,            marker = marker,            markersize = 15)    if row.Library == "Julia"        push!(julia_points, s)    else        push!(python_points, s)    endend# Add algorithm labels for Julia implementationsfor row in eachrow(filter(r -> r.Library == "Julia", pareto_combined))    text!(ax, row.Time_ms, row.Mean_Error,          text = "  " * row.Algorithm,          align = (:left, :center),          fontsize = 10)end# Add algorithm labels for Python implementationsfor row in eachrow(filter(r -> r.Library == "Python", pareto_combined))    text!(ax, row.Time_ms, row.Mean_Error,          text = "  " * row.Algorithm,          align = (:left, :center),          fontsize = 9,          color = (:gray40, 0.8))end# Connect Julia points to show the trade-off curvejulia_data = sort(filter(r -> r.Library == "Julia", pareto_combined), :Time_ms)lines!(ax, julia_data.Time_ms, julia_data.Mean_Error,       color = (:steelblue, 0.3), linestyle = :dash, linewidth = 1.5, label = "Julia Frontier")# Connect Python pointspython_data = sort(filter(r -> r.Library == "Python", pareto_combined), :Time_ms)lines!(ax, python_data.Time_ms, python_data.Mean_Error,       color = (:coral, 0.3), linestyle = :dot, linewidth = 1.5, label = "Python Frontier")# Set y-axis limits for better discriminationylims!(ax, nothing, 0.005)# Add guidelines for accuracy levelshlines!(ax, [0.001, 0.002], color = (:gray, 0.2), linestyle = :dot, linewidth = 1)# Create legend for library typeslibrary_elements = [    MarkerElement(marker = :circle, color = :gray, markersize = 12),    MarkerElement(marker = :utriangle, color = :gray, markersize = 12)]Legend(fig_pareto[1, 2], library_elements, ["Julia", "Python"], "Library",       framevisible = true, backgroundcolor = :transparent, labelsize = 10)