Mercurial > repos > pointsource_pde / changeset

--- a/.hgignore	Sun Oct 04 13:02:26 2026 -0500
+++ b/.hgignore	Sun Oct 04 16:30:24 2026 -0500
@@ -13,6 +13,10 @@
 .venv/
 basix-*/
 dolfinx-*/
-res
-res_*/
 results/
+results_*/
+tmp/
+__old*
+experiments/variant*
+REMAINING
+REMAINING.lock

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/stats/base.py	Sun Oct 04 16:30:24 2026 -0500
@@ -0,0 +1,60 @@
+import laser_and_mirrors_aux
+
+# Give name to the problem
+name = "laser_and_mirrors_aux"
+
+# Fine-tune algorithms
+alg_finetune = {
+    "tau0": 0.02,
+    "sigmap0": 0.0001,
+    "tautheta": 0.07,
+    "gradv_lipest_mult": 1,
+    "transport_tolerance": 1e5,
+    "inner_method": "Exact",
+    "merge": True,
+    "fitness_merging": True,
+    "merge_radius": 0.1,
+    "merge_interp": False,
+    "scaling_heuristic": None,
+    "extra_weight_optimisation_steps": None,
+}
+
+reg_params = {
+    "alpha": 2e-3,
+    "alpha_k": 60,
+    "alpha_c": 0.06,
+    "c_noise": 0.02,
+    "k_noise": 0.02,
+    "seed": 2393,
+    "true_init": False,
+    "true_target": False,
+    "noise": 0.01,
+    "sim_oversample": 2,
+    "sim_oversample_time": 2,
+    "scale_cstep": 300,
+    "do_non_sliding": True,
+}
+
+algorithm_overrides = {
+    "RadonFB": alg_finetune,
+    "RadonSlidingFB": alg_finetune,
+    "RadonPDPS": alg_finetune,
+    "RadonSlidingPDPS": alg_finetune,
+}
+
+
+def setup(prefix):
+    return laser_and_mirrors_aux.setup(
+        prefix,
+        alpha=reg_params["alpha"],
+        alpha_k=reg_params["alpha_k"],
+        alpha_c=reg_params["alpha_c"],
+        c_noise=reg_params["c_noise"],
+        k_noise=reg_params["k_noise"],
+        seed=reg_params["seed"],
+        true_init=reg_params["true_init"],
+        true_target=reg_params["true_target"],
+        sim_oversample=reg_params["sim_oversample"],
+        sim_oversample_time=reg_params["sim_oversample_time"],
+        scale_cstep=reg_params["scale_cstep"],
+    )

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/stats/extract_results.sh	Sun Oct 04 16:30:24 2026 -0500
@@ -0,0 +1,86 @@
+#!/bin/sh
+iter=${MAXITER:-5000}
+files=(radon_sliding_fb_log.txt radon_fb_log.txt radon_sliding_fb_config.json radon_fb_config.json radon_sliding_fb_stats.json radon_fb_stats.json)
+tmp=tmp
+
+this_script=$(realpath "$0")
+dir=$(dirname "$this_script")
+top="$dir/../"
+
+trap 'kill -- -$$' INT
+
+function generate {
+    i="$2"
+    job=$(basename "$i")
+    dest=$top/$tmp/"$i"
+    if echo "$i" | grep do_non_sliding=False  > /dev/null; then
+        do_non_sliding=false
+    else
+        do_non_sliding=true
+    fi
+
+    if ! test "$1" = partial; then
+        mkdir -p "$dest"
+        src1=$i/laser_and_mirrors_aux
+        src2=$i/laser_and_mirrors_aux2
+
+        $top/plot.py ${src1}/radon_sliding_fb --simple --iter ${iter} --save-only
+        $top/extract_aux.py ${src1}/radon_sliding_fb
+        cp ${src1}/radon_sliding_fb/res_${iter}.pdf ${dest}/aux_res_sliding_${iter}.pdf
+        cp ${src1}/radon_sliding_fb/true_${iter}.pdf ${dest}/aux_true_${iter}.pdf
+        cp ${src1}/radon_sliding_fb_log.txt ${dest}/aux_radon_sliding_fb_log.txt
+        cp ${src1}/radon_sliding_fb/res_aux.csv ${dest}/aux_radon_sliding_fb_aux.csv
+
+        if "$do_non_sliding"; then
+            $top/plot.py ${src1}/radon_fb --simple --iter ${iter} --save-only
+            $top/extract_aux.py ${src1}/radon_fb
+            cp ${src1}/radon_fb/res_${iter}.pdf ${dest}/aux_res_${iter}.pdf
+            cp ${src1}/radon_fb_log.txt ${dest}/aux_radon_fb_log.txt
+            cp ${src1}/radon_fb/res_aux.csv ${dest}/aux_radon_fb_aux.csv
+        fi
+
+        cp ${src1}/init_aux.csv ${dest}/aux_init_aux.csv
+        for file in $files; do
+            cp ${src1}/$file ${dest}/aux_$file
+        done
+
+        $top/plot.py ${src2}/radon_sliding_fb --simple --iter ${iter} --save-only
+        $top/extract_aux.py ${src2}/radon_sliding_fb
+        cp ${src2}/radon_sliding_fb/res_${iter}.pdf ${dest}/aux2_res_sliding_${iter}.pdf
+        cp ${src2}/radon_sliding_fb/true_${iter}.pdf ${dest}/aux2_true_${iter}.pdf
+        cp ${src2}/radon_sliding_fb_log.txt ${dest}/aux2_radon_sliding_fb_log.txt
+        cp ${src2}/radon_sliding_fb/res_aux.csv ${dest}/aux2_radon_sliding_fb_aux.csv
+
+        if "$do_non_sliding"; then
+            $top/plot.py ${src2}/radon_fb --simple --iter ${iter} --save-only
+            $top/extract_aux.py ${src2}/radon_fb
+            cp ${src2}/radon_fb/res_${iter}.pdf ${dest}/aux2_res_${iter}.pdf
+            cp ${src2}/radon_fb_log.txt ${dest}/aux2_radon_fb_log.txt
+            cp ${src2}/radon_fb/res_aux.csv ${dest}/aux2_radon_fb_aux.csv
+        fi
+
+        cp ${src2}/init_aux.csv ${dest}/aux2_init_aux.csv
+        for file in $files; do
+            cp ${src2}/$file ${dest}/aux2_$file
+        done
+    fi
+    (cd $dest; pdflatex -jobname="$job" '\def\maxiter{'$iter'}\input{../plotnew}')
+}
+
+for i in "$@"; do
+    echo "Generating ${i}…"
+    generate "$GEN_MANY_PARTIAL" "$i" &
+done
+
+wait
+
+files=()
+for i in "$@"; do
+    fn="$tmp/$i/$i.pdf"
+    if test -f "$fn"; then
+        files+=("$fn")
+    fi
+done
+#open "${files[@]}"
+echo "${files[@]}"
+pdfjam -o fused.pdf "${files[@]}"

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/stats/generate.sh	Sun Oct 04 16:30:24 2026 -0500
@@ -0,0 +1,45 @@
+#!/bin/sh
+
+function rep {
+    b="$1.py"
+    shift
+    p="$1"
+    shift
+    expr=""
+    name="$p"
+    comma=""
+    var="$1"
+    shift
+    while test "$var" != "--"; do
+        val="$1"
+        shift
+        expr="$expr s/\"$var\" *:.*/\"$var\": $val,/;"
+        name="$name$comma$var=$val"
+        comma=","
+        var="$1"
+        shift
+    done
+    var="$1"
+    shift
+    for val in "$@"; do
+        sed "$expr; s/\"$var\" *:.*/\"$var\": $val,/" $b > "$name$comma$var=$val].py"
+    done
+}
+
+this_script=$(realpath "$0")
+dir=$(dirname "$this_script")
+
+prefix="$dir/../experiments/variant["
+base="$dir/base"
+
+for seed in $(seq 1 8); do
+    # Most parameters halved or doubled.
+    # Step length parameters halved or divided by four, due to being chosen as high as works.
+    rep $base $prefix -- seed $seed
+    rep $base $prefix seed $seed -- tau0 0.01 0.005
+    rep $base $prefix seed $seed -- sigmap0 0.00005 0.0002
+    rep $base $prefix seed $seed -- alpha 1e-3 4e-3
+    rep $base $prefix seed $seed -- alpha_k 30 120
+    rep $base $prefix seed $seed -- alpha_c 0.03 0.12
+    rep $base $prefix seed $seed do_non_sliding False -- tautheta 0.0175 0.035
+done

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/stats/run.sh	Sun Oct 04 16:30:24 2026 -0500
@@ -0,0 +1,35 @@
+#!/bin/sh
+
+iters=5000
+
+set -e
+
+for f in "$@"; do
+    if ! test -f "$f"; then
+        echo "$f does not exist"
+        exit 1
+    fi
+done
+
+i=1
+n=$#
+for f in "$@"; do
+    a=$(basename "$f" .py)
+    #prefix="results_$(hg id -i)_$(hg --repo=../pointsource_algs id -i)_$(hg --repo=../measures id -i)"
+    prefix="results"
+    printf '\033]2;%s\033\\' "${prefix}_$a"
+    tmp=$(mktemp)
+    sed "s/laser_and_mirrors/laser_and_mirrors2/; s/mirrors2_aux/mirrors_aux2/;" "$f" > "$tmp"
+    if echo "$a" | grep do_non_sliding=False  > /dev/null; then
+        nice cargo run --release  --   -o "${prefix}_$a" -a radon_sliding_fb --max-iter $iters "$f" "$tmp"
+    else
+        nice cargo run --release  --   -o "${prefix}_$a" -a radon_sliding_fb -a radon_fb --max-iter $iters "$f" "$tmp"
+    fi
+    #echo "Results in ${prefix}_$a"
+    i=$(expr $i + 1)
+done
+
+for f in "$@"; do
+    a=$(basename "$f" .py)
+    echo "Results in ${prefix}_$a"
+done

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/stats/summarise.jl	Sun Oct 04 16:30:24 2026 -0500
@@ -0,0 +1,186 @@
+#!/usr/bin/env julia
+
+using CSV
+using DataFrames
+using Statistics
+using Distributions
+
+# --------------------------------------------------------------------------
+# CONFIGURATION
+#
+# Map from CSV file name (as it appears inside each parameter directory) to
+# the list of data columns in that file that should be summarized. Files may
+# have different columns from one another; the `iter` column is assumed to
+# be present in all of them and is NOT listed here.
+# --------------------------------------------------------------------------
+const LOG_COLS = ["value", "relative_value", "discrepancy", "aux_discrepancy", "joint_discrepancy"];
+const AUX_COLS = ["k", "c1", "c2", "k_dist", "c_dist"]
+
+const FILES_COLUMNS = Dict(
+    "aux_radon_fb_log.txt" => LOG_COLS,
+    "aux_radon_sliding_fb_log.txt" => LOG_COLS,
+    "aux2_radon_fb_log.txt" => LOG_COLS,
+    "aux2_radon_sliding_fb_log.txt" => LOG_COLS,
+    "aux_radon_fb_aux.csv" => AUX_COLS,
+    "aux_radon_sliding_fb_aux.csv" => AUX_COLS,
+    "aux2_radon_fb_aux.csv" => AUX_COLS,
+    "aux2_radon_sliding_fb_aux.csv" => AUX_COLS,
+)
+
+const ITER_COL = :iter
+const CONF_LEVEL = 0.95
+
+# --------------------------------------------------------------------------
+# Compute the half-width of a 95% (or CONF_LEVEL) confidence interval for
+# the mean of `values`, using the t-distribution. Falls back gracefully for
+# degenerate sample sizes.
+# --------------------------------------------------------------------------
+function conf_margin(values::AbstractVector{<:Real})
+    n = length(values)
+    if n < 2
+        return 0.0
+    end
+    s = std(values)
+    if s == 0
+        return 0.0
+    end
+    tdist = TDist(n - 1)
+    tcrit = quantile(tdist, 1 - (1 - CONF_LEVEL) / 2)
+    return tcrit * s / sqrt(n)
+end
+
+# --------------------------------------------------------------------------
+# Read one configured file from one parameter directory, if it exists.
+# Returns a DataFrame with at least `iter` and the configured columns, or
+# `nothing` if the file is missing or unreadable.
+# --------------------------------------------------------------------------
+function read_dir_file(dir::AbstractString, filename::AbstractString, cols::Vector{String})
+    path = joinpath(dir, filename)
+    if !isfile(path)
+        @warn "Missing file, skipping" path
+        return nothing
+    end
+
+    local df
+    try
+        df = CSV.read(path, DataFrame)
+    catch e
+        @warn "Could not read CSV, skipping" path exception = e
+        return nothing
+    end
+
+    if !hasproperty(df, ITER_COL)
+        @warn "File has no '$(ITER_COL)' column, skipping" path
+        return nothing
+    end
+
+    present_cols = String[]
+    for c in cols
+        if hasproperty(df, Symbol(c))
+            push!(present_cols, c)
+        else
+            @warn "Column '$c' not found in file, skipping that column for this directory" path
+        end
+    end
+    if isempty(present_cols)
+        return nothing
+    end
+
+    keep = vcat([ITER_COL], Symbol.(present_cols))
+    return df[:, keep]
+end
+
+# --------------------------------------------------------------------------
+# Build the summary DataFrame for one configured file, across all `dirs`.
+# --------------------------------------------------------------------------
+function summarize_file(dirs::Vector{String}, filename::AbstractString, cols::Vector{String})
+    parts = DataFrame[]
+    for dir in dirs
+        d = read_dir_file(dir, filename, cols)
+        d === nothing || push!(parts, d)
+    end
+
+    if isempty(parts)
+        @warn "No data found for file across any directory, skipping summary" filename
+        return nothing
+    end
+
+    combined = vcat(parts...; cols=:union)
+    grouped = groupby(combined, ITER_COL)
+
+    iters = sort(collect(keys(grouped)))
+    summary = DataFrame(iter=[k[ITER_COL] for k in iters])
+
+    for col in cols
+        colsym = Symbol(col)
+        mins = Vector{Union{Missing,Float64}}(undef, length(iters))
+        maxs = Vector{Union{Missing,Float64}}(undef, length(iters))
+        means = Vector{Union{Missing,Float64}}(undef, length(iters))
+        confmin = Vector{Union{Missing,Float64}}(undef, length(iters))
+        confmax = Vector{Union{Missing,Float64}}(undef, length(iters))
+
+        for (i, key) in enumerate(iters)
+            sub = grouped[key]
+            if !hasproperty(sub, colsym)
+                mins[i] = maxs[i] = means[i] = confmin[i] = confmax[i] = missing
+                continue
+            end
+            vals = collect(skipmissing(sub[!, colsym]))
+            if isempty(vals)
+                mins[i] = maxs[i] = means[i] = confmin[i] = confmax[i] = missing
+                continue
+            end
+            m = mean(vals)
+            margin = conf_margin(vals)
+            mins[i] = minimum(vals)
+            maxs[i] = maximum(vals)
+            means[i] = m
+            confmin[i] = m - margin
+            confmax[i] = m + margin
+        end
+
+        summary[!, Symbol(col, "_min")] = mins
+        summary[!, Symbol(col, "_max")] = maxs
+        summary[!, Symbol(col, "_confmax")] = confmax
+        summary[!, Symbol(col, "_confmin")] = confmin
+        summary[!, Symbol(col, "_mean")] = means
+    end
+
+    return summary
+end
+
+# --------------------------------------------------------------------------
+# Main
+# --------------------------------------------------------------------------
+function main()
+    if isempty(ARGS)
+        println(stderr, "Usage: julia summarize.jl DIR1 [DIR2 DIR3 ...]")
+        exit(1)
+    end
+
+    dirs = String[]
+    for d in ARGS
+        if isdir(d)
+            push!(dirs, d)
+        else
+            @warn "Not a directory, skipping" d
+        end
+    end
+
+    if isempty(dirs)
+        println(stderr, "No valid directories given.")
+        exit(1)
+    end
+
+    for (filename, cols) in FILES_COLUMNS
+        summary = summarize_file(dirs, filename, cols)
+        summary === nothing && continue
+
+        stem = first(splitext(filename))
+        outpath = joinpath(pwd(), "$(stem)_summary.csv")
+        CSV.write(outpath, summary)
+        println("Wrote $outpath ($(nrow(summary)) rows)")
+    end
+end
+
+main()