try Gauss-Newton with barriers

test.jl

@@ -1,43 +1,38 @@
 import GLMakie as M
+import GaussNewton as GN
 import Metaheuristics as MH
 import DelimitedFiles: readdlm
+import LinearAlgebra: norm
 
 include("PLUV.jl")
 include("Utils.jl")
 
-
 # filtered_idx = d[:, 2] .> 0.0
 # @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")
 
-best_5k = begin
+best_5k_full = begin
   f_χ2 = open("POPC-test-5k/Results_collection-X2.dat", "r")
 
   d_χ2 = readdlm(f_χ2)
   close(f_χ2)
 
   idx_min_χ2 = argmin(d_χ2)[1]
-  @show idx_min_χ2
 
   f_params = open("POPC-test-5k/Results_collection.dat", "r")
   for i in 1:(idx_min_χ2-1)
     readline(f_params)
   end
   best_params = map(x -> parse(Float64, x), split(readline(f_params), ' '))
-  @show best_params
-
-  PLUV.reduce_to_free_parameters(meta, best_params)
 end
 
-mean_5k = begin
+mean_5k_full = begin
   f_params = open("POPC-test-5k/Results_collection.dat", "r")
   params_pop = readdlm(f_params)
   N_trials, _ = size(params_pop)
 
   mean_params = vec(sum(params_pop; dims=1) / N_trials)
-
-  PLUV.reduce_to_free_parameters(meta, mean_params)
 end
 
 # best_5k = PLUV.reduce_to_free_parameters(meta, [100000.0000000, 7.5004085, 23.5074449, 9.8664991, 0.0000000, 1.0000000, 0.2900000, 0.2203793,
@@ -48,6 +43,7 @@ end
 #   0.4158804, 0.3278631, 0.2296156, 0.2475607, 0.6664143, 0.1191859, 0.2618609, 0.9000000, 0.3050000, 0.2963982, 0.2770345, 0.4762528,
 #   0.8100000, 1.9706651, 37.0000000, 0.0299179, 0.0002167])
 
+best_5k_params = PLUV.reduce_to_free_parameters(meta, best_params)
 
 q, I_data, err = Utils.select_columns(qie_data, meta["q_min"], meta["q_max"], meta["binning"], meta["binning"])
 q_all, I_all, err_all = Utils.select_columns(qie_data, meta["q_min"], meta["q_max"], 1, meta["binning"])
@@ -58,39 +54,68 @@ I_all_lp = Utils.lowpass_filter(I_all; σ=2)
 w = Utils.compute_logscale_weights(q)
 w_all = Utils.compute_logscale_weights(q_all)
 
-I_init = PLUV.intensity(params_init, q)
+# I_init = PLUV.intensity(params_init, q)
 
 intensity_reduced, P_reduced, lb_reduced, ub_reduced = PLUV.reduce_to_free_parameters(meta, params_init, lower_bounds, upper_bounds, q)
 
+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
+
 bounds = MH.boxconstraints(lb=lb_reduced, ub=ub_reduced)
 
 function obj(P)
-  I_model = intensity_reduced(P)
+  I_model, neg_H20 = intensity_reduced(P)
-  return Utils.chi2(I_data, I_model, err)
+  χ2 = Utils.chi2(I_data, I_model, err)
+  factor = (neg_H20 == 0 ? 1.0 : 5.0 * (neg_H20 + 1))
+  return factor * χ2
 end
 
+barriered_obj = Utils.add_log_barriers(obj, simple_bounds; padding_factor=padding_factor)
+
 information = MH.Information(f_optimum=0.0)
-options = MH.Options(f_calls_limit=1_000_000, f_tol=1e-5);
+information = MH.Information()
+options = MH.Options(f_calls_limit=10_000, f_tol=1e-5);
 algorithm = MH.ECA(information=information, options=options)
+#algorithm = MH.PSO()
 
-I_best_5k = intensity_reduced(best_5k)
+I_best_5k, _ = PLUV.intensity(best_5k_full, q)
-I_mean_5k = intensity_reduced(mean_5k)
 
+I_mean_5k, _ = PLUV.intensity(mean_5k_full, q_all)
+
+# Gauss-Newton
+if true
+  _, result = GN.optimize(barriered_obj, simple_init)
+  if GN.has_converged(result)
+    P_best = result.minimizer
+    I_best, _ = intensity_reduced(P_best)
+  else
+    @error "Gauss-Newton did not converge"
+  end
+end
+
+# Metaheuristics
 if false
   result = MH.optimize(obj, bounds, algorithm)
+  @show MH.minimum(result)
 
   P_best = MH.minimizer(result)
-  I_best = intensity_reduced(P_best)
+  I_best, _ = intensity_reduced(P_best)
 end
 
 if true
   fig = M.Figure()
-  ax = M.Axis(fig[1, 1]; xscale=log10, yscale=log10)
+  #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")
   M.lines!(ax, q, I_best_5k, label="TSA best (5k)")
-  M.lines!(ax, q, I_mean_5k, label="TSA mean (5k)")
+  M.lines!(ax, q_all, I_mean_5k, label="TSA mean (5k)")
 
   M.axislegend()