From 9014cc08fe883f0b09423720d7273f56b2595241 Mon Sep 17 00:00:00 2001
From: Orso Meneghini <orso82@gmail.com>
Date: Mon, 27 Nov 2023 11:50:08 -0800
Subject: [PATCH] Playing around with Surrogates.jl

---
 Project.toml                               |  2 ++
 src/actors/transport/flux_matcher_actor.jl | 28 ++++++++++++++++++++++
 2 files changed, 30 insertions(+)

diff --git a/Project.toml b/Project.toml
index 8bd8059b9..255f01843 100644
--- a/Project.toml
+++ b/Project.toml
@@ -56,6 +56,8 @@ Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
 SimulationParameters = "32ab26d3-25d6-405d-a295-367385e16093"
 SnoopPrecompile = "66db9d55-30c0-4569-8b51-7e840670fc0c"
 SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b"
+Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
+Surrogates = "6fc51010-71bc-11e9-0e15-a3fcc6593c49"
 TAUENN = "4f5ffc88-c87f-11eb-142e-4d2ba0394a67"
 TEQUILA = "a60c9cbd-e72f-4185-96b6-b8fc312c4d37"
 TGLFNN = "558c7b13-fd9f-4806-b461-592296cfa9d0"
diff --git a/src/actors/transport/flux_matcher_actor.jl b/src/actors/transport/flux_matcher_actor.jl
index c7b8990a8..868590214 100644
--- a/src/actors/transport/flux_matcher_actor.jl
+++ b/src/actors/transport/flux_matcher_actor.jl
@@ -1,5 +1,7 @@
 import NLsolve
 using LinearAlgebra
+import Surrogates
+import Statistics
 
 #= ================ =#
 #  ActorFluxMatcher  #
@@ -112,6 +114,32 @@ function _step(actor::ActorFluxMatcher)
         actor_logging(dd, old_logging)
     end
 
+    #### SURROGATES ####
+    old_logging = actor_logging(dd, false)
+    try
+        tmp = hcat(map(collect, z_scaled_history)...)
+
+        # Calculate the mean and standard deviation for each set of corresponding elements.
+        means = Statistics.mean(tmp; dims=2)  # Result is a 1-row matrix
+        std_devs = Statistics.std(tmp; dims=2)  # Result is a 1-row matrix
+
+        N = 5  # Number of sigmas (change this to your specific number of sigmas)
+
+        # Calculate the lower and upper bounds.
+        lower_bound = means .- N * std_devs
+        upper_bound = means .+ N * std_devs
+        lower_bound = scale_z_profiles(means .* 0.0 .- 4.0)
+        upper_bound = scale_z_profiles(means .* 0.0 .+ 0.0)
+        @show lower_bound
+        @show upper_bound
+
+        radial_basis = Surrogates.SecondOrderPolynomialSurrogate(z_scaled_history, err_history, lower_bound, upper_bound)
+        Surrogates.surrogate_optimize(z -> norm(flux_match_errors(actor, z)), Surrogates.SRBF(), lower_bound, upper_bound, radial_basis, Surrogates.SobolSample(); maxiters=50)
+    finally
+        actor_logging(dd, old_logging)
+    end
+    #### SURROGATES ####
+
     if par.do_plot
         display(res)
     end