From f7917094b74fd3efe1b8532e394f3e4571bcfa2e Mon Sep 17 00:00:00 2001 From: Gaspard Jankowiak Date: Thu, 16 May 2024 12:30:37 +0200 Subject: [PATCH] fix initial conditions, check bounds --- PLUV.jl | 2 +- Utils.jl | 43 +++++++++++++++++++++++++++++++++++++------ test.jl | 27 +++++++++++++++++++-------- 3 files changed, 57 insertions(+), 15 deletions(-) diff --git a/PLUV.jl b/PLUV.jl index 2c49bb0..7763872 100644 --- a/PLUV.jl +++ b/PLUV.jl @@ -333,7 +333,7 @@ end function reduce_to_free_parameters(metadata::Dict, params_init::NamedTuple, lb, ub, q) default_params = Float64[params_init...] free_idx = get_free_parameters_idx(metadata) - P_reduced = view(default_params, free_idx) + P_reduced = default_params[free_idx] lb_reduced = lb[free_idx] ub_reduced = ub[free_idx] function reduced(P::AbstractArray{Float64}) diff --git a/Utils.jl b/Utils.jl index 2756c42..379b1c8 100644 --- a/Utils.jl +++ b/Utils.jl @@ -3,12 +3,23 @@ module Utils import TOML: parsefile import DelimitedFiles: readdlm import DSP +import Printf: @sprintf +import Crayons as C function load_init_params(filename::String) meta = parsefile(filename)["metadata"] - params = NamedTuple(Dict(Symbol(p["symbol"]) => p["initial_value"] for p in meta["parameters"])) + param_names = [p["symbol"] for p in meta["parameters"]] + param_initial = [p["initial_value"] for p in meta["parameters"]] + params = NamedTuple{Tuple(map(Symbol, param_names))}(param_initial) lower_bounds = [p["min"] for p in meta["parameters"]] upper_bounds = [p["max"] for p in meta["parameters"]] + println("Got $(length(param_names)) parameters:") + for (i, k) in enumerate(eachindex(param_names)) + lb_str = @sprintf "%.3g" lower_bounds[i] + ub_str = @sprintf "%.3g" upper_bounds[i] + p = @sprintf "%.3g" param_initial[i] + println("$(lpad(param_names[i], 7)): $(lpad(lb_str, 7)) <= $(lpad(p, 7)) <= $(lpad(ub_str, 7)) ") + end return (meta=meta, params=params, lower_bounds=lower_bounds, upper_bounds=upper_bounds) end @@ -85,9 +96,11 @@ function chi2(I_data, I_model, err) return sum((I_data .- I_model ./ err) .^ 2) end -function add_log_barriers(f::Function, constraints::Vector{Tuple{T,T}}; padding_factor=0.1) where {T<:Real} +function add_log_barriers(f::Function, constraints::Vector{Tuple{T,T}}; padding_factor=0.1, mode::Symbol=:outer) where {T<:Real} n = length(constraints) + @assert mode in [:outer, :inner] + lower_shifts = Float64[] upper_shifts = Float64[] @@ -124,21 +137,23 @@ function add_log_barriers(f::Function, constraints::Vector{Tuple{T,T}}; padding_ end end + barrier_sign = mode == :outer ? 1.0 : -1.0 + function std_lower_barrier(s::Float64) - return max(-log(s + 1.0), 0.0)^2 + return max(-log(s + barrier_sign), 0.0)^2 end function barrier(x::Vector{Float64}) lb = (x .- lower_shifts) ./ barrier_widths ub = (upper_shifts .- x) ./ barrier_widths - lb_violations = findall(<=(-1.0), lb) - ub_violations = findall(<=(-1.0), ub) + lb_violations = findall(<=(-barrier_sign), lb) + ub_violations = findall(<=(-barrier_sign), ub) if !isempty(lb_violations) err_msg = "" for i in lb_violations - err_msg *= " · x[$i] = $(x[i]) < $(constraints[i][1]), width = $(barrier_widths[i])" + err_msg *= "\n · x[$i] = $(x[i]) < $(constraints[i][1]), width = $(barrier_widths[i])" end @error err_msg throw("lower bound violation for indice(s) $(lb_violations)") @@ -160,4 +175,20 @@ function add_log_barriers(f::Function, constraints::Vector{Tuple{T,T}}; padding_ return x::Vector{Float64} -> f(x) .+ barrier(x) end +function print_check_bounds(param_names, param_values, lb, ub) + for (i, v) in enumerate(param_values) + lb_str = @sprintf "%.3g" lb[i] + ub_str = @sprintf "%.3g" ub[i] + p = @sprintf "%.3g" v + if lb[i] <= v <= ub[i] + color = :white + bold = false + else + color = :light_red + bold = true + end + println("$(lpad(param_names[i], 7)): ", C.Crayon(foreground=color, bold=bold), "$(lpad(lb_str, 7)) <= $(lpad(p, 7)) <= $(lpad(ub_str, 7)) ") + end +end + end # module Utils diff --git a/test.jl b/test.jl index b3e4885..3604034 100644 --- a/test.jl +++ b/test.jl @@ -11,6 +11,7 @@ include("Utils.jl") # @views q, I, err = d[filtered_idx, 1], d[filtered_idx, 2], d[filtered_idx, 3] meta, params_init, lower_bounds, upper_bounds, qie_data = Utils.load_config("Data/PAR_POPC-test.toml") +param_names = keys(params_init) best_5k_full = begin f_χ2 = open("POPC-test-5k/Results_collection-X2.dat", "r") @@ -61,10 +62,7 @@ intensity_reduced, P_reduced, lb_reduced, ub_reduced = PLUV.reduce_to_free_param simple_bounds = collect(zip(lb_reduced, ub_reduced)) simple_init = 0.5 * (lb_reduced .+ ub_reduced) -padding_factor = fill(0.1, length(simple_init)) - -padding_factor[17] = 10.0 -padding_factor[18] = 10.0 +padding_factor = fill(1e-1, length(simple_init)) bounds = MH.boxconstraints(lb=lb_reduced, ub=ub_reduced) @@ -82,7 +80,7 @@ function obj_residuals(P) return factor * residuals end -barriered_obj = Utils.add_log_barriers(obj_residuals, simple_bounds; padding_factor=padding_factor) +barriered_obj = Utils.add_log_barriers(obj_residuals, simple_bounds; padding_factor=padding_factor, mode=:outer) information = MH.Information(f_optimum=0.0) information = MH.Information() @@ -96,10 +94,23 @@ I_mean_5k, _ = PLUV.intensity(mean_5k_full, q_all) # Gauss-Newton if true - _, result = GN.optimize(barriered_obj, simple_init) + initial_guess = PLUV.reduce_to_free_parameters(meta, collect(Float64, values(params_init))) + initial_guess = best_5k_params + + _, result = GN.optimize(barriered_obj, initial_guess; show_trace=true, iscale=1, ZCP=1e-2, ZCPMIN=1e-2) + #_, result = GN.optimize(barriered_obj, initial_guess) + if GN.has_converged(result) + @info "Gauss-Newton converged" + @show result P_best = result.minimizer I_best, _ = intensity_reduced(P_best) + + @info "Simulated Annealing 5k" + Utils.print_check_bounds(param_names, best_5k_params, lb_reduced, ub_reduced) + + @info "Our result" + Utils.print_check_bounds(param_names, P_best, lb_reduced, ub_reduced) else @error "Gauss-Newton did not converge" end @@ -116,8 +127,8 @@ end if true fig = M.Figure() - #ax = M.Axis(fig[1, 1]; xscale=log10, yscale=log10) - ax = M.Axis(fig[1, 1]) + ax = M.Axis(fig[1, 1]; xscale=log10, yscale=log10) + #ax = M.Axis(fig[1, 1]) M.lines!(ax, q, I_best, label="MH best (julia)") M.scatter!(ax, q, I_data, label="data")