use nlopt to find xopt

This commit updates both the chi2 and dm-search scripts to run nlopt at
the end of the MCMC to find the best fit point.

The reason for this is that for the dm-search script, we rely on the
best fit point being correct when doing the discovery threshold
analysis. To really get the best fit point with the MCMC we would need
to run a ridiculous number of steps, so it's better to run a fewer
amount of steps to get close and then run the minimization from there.

2 files changed, 30 insertions, 2 deletions

diff --git a/utils/chi2 b/utils/chi2
index 12b602a..0326e1a 100755
--- a/utils/chi2
+++ b/utils/chi2
@@ -35,6 +35,7 @@ from sddm.dc import estimate_errors, EPSILON, truncnorm_scaled
 import emcee
 from sddm import printoptions
 from sddm.utils import fast_cdf, correct_energy_bias
+import nlopt
 
 # Likelihood Fit Parameters
 # 0 - Atmospheric Neutrino Flux Scale
@@ -390,7 +391,20 @@ def do_fit(data,muon,data_mc,atmo_scale_factor,muon_scale_factor,bins,steps,prin
 
     samples = sampler.get_chain(flat=True,thin=thin)
 
-    return sampler.get_chain(flat=True)[sampler.get_log_prob(flat=True).argmax()], samples
+    # Now, we use nlopt to find the best set of parameters. We start at the
+    # best starting point from the MCMC and then run the SBPLX routine.
+    x0 = sampler.get_chain(flat=True)[sampler.get_log_prob(flat=True).argmax()]
+    opt = nlopt.opt(nlopt.LN_SBPLX, len(x0))
+    opt.set_min_objective(nll)
+    low = np.array(PRIORS_LOW)
+    high = np.array(PRIORS_HIGH)
+    opt.set_lower_bounds(low)
+    opt.set_upper_bounds(high)
+    opt.set_ftol_abs(1e-10)
+    opt.set_initial_step([0.01]*len(x0))
+    xopt = opt.optimize(x0)
+
+    return xopt, samples
 
 if __name__ == '__main__':
     import argparse
diff --git a/utils/dm-search b/utils/dm-search
index 2a008d5..4e8772e 100755
--- a/utils/dm-search
+++ b/utils/dm-search
@@ -36,6 +36,7 @@ from sddm.utils import fast_cdf, correct_energy_bias
 from scipy.integrate import quad
 from sddm.dm import *
 from sddm import SNOMAN_MASS
+import nlopt
 
 # Likelihood Fit Parameters
 # 0 - Atmospheric Neutrino Flux Scale
@@ -312,7 +313,20 @@ def do_fit(dm_particle_id,dm_mass,dm_energy,data,muon,data_mc,atmo_scale_factor,
 
     samples = sampler.get_chain(flat=True,thin=thin)
 
-    return sampler.get_chain(flat=True)[sampler.get_log_prob(flat=True).argmax()], samples
+    # Now, we use nlopt to find the best set of parameters. We start at the
+    # best starting point from the MCMC and then run the SBPLX routine.
+    x0 = sampler.get_chain(flat=True)[sampler.get_log_prob(flat=True).argmax()]
+    opt = nlopt.opt(nlopt.LN_SBPLX, len(x0))
+    opt.set_min_objective(nll)
+    low = np.array(PRIORS_LOW)
+    high = np.array(PRIORS_HIGH)
+    opt.set_lower_bounds(low)
+    opt.set_upper_bounds(high)
+    opt.set_ftol_abs(1e-10)
+    opt.set_initial_step([0.01]*len(x0))
+    xopt = opt.optimize(x0)
+
+    return xopt, samples
 
 def sample_priors():
     """