file operations for circle test

apps/toy/generate_chains.jl

+include("utils/all.jl")
+
+μ, σ = .5, .1
+π_ref = MixtureModel([
+    TruncatedNormal(-μ, σ, -1, 0),
+    TruncatedNormal( μ, σ,  0, 1),
+])
+
+single_peak = TruncatedNormal(μ, σ, -1, 1)
+dist_single_peak = (; 
+    ks = exact_distance_ks(single_peak, π_ref),
+    cm = exact_distance_cm(single_peak, π_ref)
+)
+
+z = halton_points(2000, 3)
+f = normalmix_proxy(z, μ, σ)
+l = MultilevelSampledLogDensity(f, [100, 2000])
+
+w = CyclicWalk(-1,1,1.0)
+
+
+begin 
+    println("Sampling chains...\n")
+
+    println("MetropolisHastings")
+    s_mh = MetropolisHastings(w)
+    @time c_mh = sample(f, s_mh, MCMCSerial(), 100000, 100)
+    display(c_mh)
+    @time write_chain("$(@__DIR__)/data/mh.csv", c_mh)
+
+    println("ChristenFox")
+    s_cf = ChristenFox(w)
+    @time c_cf = sample(l, s_cf, MCMCSerial(), 100000, 100)
+    display(c_cf)
+    @time write_chain("$(@__DIR__)/data/cf.csv", c_cf)
+    
+    saved = mean(get_info(c_cf).evaluations) / mean(get_info(c_mh).evaluations)
+    println("\n => Saved costs ", round(saved * 100, digits=2), " % \n")
+end
+
+
+begin
+    @time write_metrics("$(@__DIR__)/data/mh.csv", π_ref) 
+    @time write_metrics("$(@__DIR__)/data/cf.csv", π_ref) 
+end
+
+
+#= 
+fig = Figure(); ax = Axis(fig[1,1])
+plot_pdf!(ax, π_ref, color=:black, linestyle=:dash, label=L"π_{∞}")
+plot_pdf!(ax, f)
+#hist!(ax, c.states[:,1], normalization=:pdf, label="MH")
+Legend(fig[1,2], ax)
+display(fig)
+=#
\ No newline at end of file

apps/toy/utils/all.jl

 using Revise
 using MultilevelChainSampler
 
+using ProgressBars
+using Distributions
+using DataFrames, CSV
+using Random
+using HaltonSequences, Primes
+using CairoMakie
+using AbstractMCMC: MCMCSerial
+
+include("files.jl")
 include("proxies.jl")
 include("analysis.jl")
 include("visualize.jl")

apps/toy/utils/analysis.jl

-using Distributions
 
 # Komolgorov - Smirnov 
 # ( adapted from HypothesisTests.jl )
@@ -13,12 +12,18 @@ end
 
 function exact_distance_ks(f::Function, d::Distribution, h = 1e-4)
     x = [-1:h:1 ... ]
-    F = cdf.(d, x)
-    y = f.(x)
-    p = exp.(y); p = p / (2 * mean(p))
-    Δx = x[2:end].-x[1:end-1]
-    P = cumsum(p .* h)
-    return maximum(abs.(P .- F))
+    f = f.(x); f = exp.(f); f = f / (2 * mean(f))
+    F = cumsum(f .* h)
+    P = cdf.(d, x)
+    δ = maximum(abs.(F .- P))
+    return δ
+end 
+
+function exact_distance_ks(f::Distribution, d::Distribution, h = 1e-4)
+    x = [-1:h:1 ... ]
+    F = cdf.(f, x)
+    P = cdf.(d, x)
+    return maximum(abs.(F .- P))
 end 
 
 
@@ -28,25 +33,35 @@ function distance_cm(x::Vector{<:Real}, d::UnivariateDistribution)
     x = sort(x)
     cdfs = cdf.(d, x)
     ω² = mean( ((1:n) / n .+ .5 - cdfs ).^2 ) + 1/(12*n^2)
-    return ω²
+    return sqrt(ω²)
 end
 
 function exact_distance_cm(f::Function, d::Distribution, h = 1e-4)
     x = [-1:h:1 ... ]
-    y = f.(x)
-    F = cdf.(d, x)
-    f = pdf.(d, x)
-    p = exp.(y); p = p / (2 * mean(p))
-    P = cumsum(p .* h)
-    return 2 * mean( (P .- F).^2 .* p )   
+    f = f.(x); f = exp.(f); f = f / (2 * mean(f))
+    F = cumsum(f .* h)
+    P = cdf.(d, x)
+    p = pdf.(d, x)
+    ω² =  2 * mean( (F .- P).^2 .* p ) 
+    return sqrt(ω²)  
 end 
 
-# Quantify convergence rate
+function exact_distance_cm(f::Distribution, d::Distribution, h = 1e-4)
+    x = [-1:h:1 ... ]
+    F = cdf.(f, x)
+    P = cdf.(d, x)
+    p = pdf.(d, x)
+    ω² =  2 * mean( (F .- P).^2 .* p ) 
+    return sqrt(ω²)  
+end 
+
+
+# Quantify convergence rate  y ≈ C⋅x^α
 function asymptotic_rate(x::Vector, y::Vector)
     log_x = log.(x)
     log_y = log.(y) 
     A = hcat( ones( length(x)), log_x )
 
     (a,b) = (A'A) \ (A' * log_y)
-    return b
+    return exp(a),b
 end

apps/toy/utils/files.jl

+
+function write_chain(csv_path::String, data::MultilevelChainSampler.SimpleChains)   
+    states = data.states
+    info   = get_info(data)
+
+    # write to .csv files
+    name, ext = splitext(csv_path)
+    ext == ".csv" || error("Only CSV files allowed.")
+    CSV.write("$name.states.csv",  DataFrame(states, :auto))
+    CSV.write("$name.info.csv",    DataFrame(info))
+    # ( ... ignore energies )
+
+    # check if files were created
+    ispath("$name.states.csv") && ispath("$name.info.csv")
+end
+
+function load_data(csv_path::String, tabs = [ "states", "info", "metrics" ])
+
+    name, ext = splitext(csv_path)
+    ext == ".csv" || error("Only CSV files allowed.")
+
+    vals = Dict( (v=>nothing for v in tabs)... )
+    for k in tabs
+        f = "$name.$k.csv"
+        if ispath(f)
+            vals[k] = CSV.read(f, DataFrame)
+        else
+            print("File '$f' not found")
+        end
+    end
+    return vals
+end
+
+function calculate_metrics(csv_path::String, d::Distribution, Ns=nothing)
+    states = load_data(csv_path, ["states"])
+
+    # Logarithmic time scale by default
+    if isnothing(Ns) 
+        Nmax = size(states,1)
+        Ns = round.(Int, exp.([0:.01:1 ... ]  .* log(Nmax)) )
+    end
+
+    # Calculate the metrics
+    n_chains = size(states, 2)
+    ks = zeros(length(Ns), n_chains)
+    cm = zeros(length(Ns), n_chains)
+    for n=ProgressBar(1:n_chains)
+        for (i,N) in enumerate(Ns)
+            ks[i,n] = distance_ks(states[1:N, n], d)
+            cm[i,n] = distance_ks(states[1:N, n], d)
+        end
+    end
+    
+    # Write to file
+    columns = [ map(i->"$l_$i", 1:n_chains) for l in ["ks", "cm"] ]
+    df = DataFrame( hcat(ks, cm), columns )
+    df.N = Ns 
+    name, ext = splitext(csv_path)
+    CSV.write("$name.metrics.csv", df)
+end
\ No newline at end of file

apps/toy/utils/proxies.jl

-using Random
-using HaltonSequences, Primes
 
 function rand_points(n::Int, d::Int=2)
     x = rand(n, d)
@@ -30,19 +28,20 @@ function as_tuples(z::Matrix)
     return [ tuple(z[i,:] ... ) for i=1:n]
 end
 
-function gauss_proxy(z::Matrix; μ::Number=0.0, σ::Number=1.0)
-    z2 = sum(z .^ 2, dims=2)[:, 1]
-    f = x -> - ( (sign(x) * sqrt( 4/pi * mean( z2 .< x^2 ) ) - μ ) / σ )^2 / 2
-    return LogDensity(f)
+function normal_proxy(z::Matrix, μ::Number=0.0, σ::Number=1.0)
+    f = function (x, z)
+        X2 = 4/pi * ( z[1]^2 + z[2]^2  < x^2 )
+        X1 = 2 * (z[3]^2 < x^2) * sign(x)
+        return - 1 / (2 * σ^2) * X2  + μ / (σ^2) * X1 
+    end
+    return SampledLogDensity( as_tuples(z), f)
 end
 
-function gauss_proxy(z::Matrix, nsamples::Vector{Int}; μ::Number=0.0, σ::Number=1.0)
-    fs = [ gauss_proxy(z[1:n, :]; μ, σ) for n=nsamples ]  
-    return MultilevelLogDensity(fs)
+function normalmix_proxy(z::Matrix, μ::Number=0.0, σ::Number=1.0)
+    f = function (x, z)
+        X2 = 4/pi * ( z[1]^2 + z[2]^2  < x^2 )
+        X1 = (z[3]^2 < x^2) # leave out the sign
+        return - 1 / (2 * σ^2) * X2  + μ / (σ^2) * X1 
+    end
+    return SampledLogDensity( as_tuples(z), f)
 end
-
-
-function gaussmix_proxy(z::Matrix, μ::Vector=[-.5, .5], σ::Vector=[.1, .1])
-    fs = [ sample_gauss_proxy(z, m, s) for (m,s) in zip(μ, σ) ]
-    return LogDensity( x -> log( sum(exp.([logdensity(f, x) for f in fs]))) )
-end
\ No newline at end of file

apps/toy/utils/visualize.jl

-using CairoMakie
 
 function plot_pdf!(ax, d::Distribution, x=[-1:.1e-4:1 ... ]; kwargs...)
     y = pdf.(d, x)
@@ -12,6 +11,9 @@ function plot_pdf!(ax, f::Function, x=[-1:.1e-4:1 ... ]; kwargs... )
     return lines!(ax, x, y; kwargs...)  
 end
 
+function plot_pdf!(ax, f::MultilevelChainSampler.AbstractLogDensity, x=[-1:.1e-4:1 ...]; kwargs...)
+    return plot_pdf!(ax, x->logdensity(f, x), x; kwargs...)
+end
 
 
 function errorlines!(ax, x, m, s; clamp=-Inf, rate=false, args...)