fix initial conditions, check bounds

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})

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

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")