Enable numerical differentiation of vector-valued functions

so that multiple constraints at the same lattice site can share the
outcome of mass matrix diagonalization.

src/array_deriv.cpp

@@ -0,0 +1,114 @@
+// ====================================================================
+// This file is part of FlexibleSUSY.
+//
+// FlexibleSUSY is free software: you can redistribute it and/or modify
+// it under the terms of the GNU General Public License as published
+// by the Free Software Foundation, either version 3 of the License,
+// or (at your option) any later version.
+//
+// FlexibleSUSY is distributed in the hope that it will be useful, but
+// WITHOUT ANY WARRANTY; without even the implied warranty of
+// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
+// General Public License for more details.
+//
+// You should have received a copy of the GNU General Public License
+// along with FlexibleSUSY.  If not, see
+// <http://www.gnu.org/licenses/>.
+// ====================================================================
+
+/**
+ * @file  src/array_deriv.cpp
+ *
+ * @brief The following implementation is an extension of
+ * gsl_deriv_central() in GNU Scientific Library 1.16 to array-valued
+ * functions
+ *
+ */
+
+#include <gsl/gsl_math.h>
+#include "array_deriv.hpp"
+
+namespace flexiblesusy {
+
+using namespace std;
+using namespace Eigen;
+
+namespace {
+
+void
+central_deriv(function<ArrayXd(double x)> f,
+	      double x, double h,
+	      ArrayXd& result, ArrayXd& abserr_round, ArrayXd& abserr_trunc)
+{
+    /* Compute the derivative using the 5-point rule (x-h, x-h/2, x,
+       x+h/2, x+h). Note that the central point is not used.
+
+       Compute the error using the difference between the 5-point and
+       the 3-point rule (x-h,x,x+h). Again the central point is not
+       used. */
+
+    ArrayXd fm1 = f(x - h);
+    ArrayXd fp1 = f(x + h);
+
+    ArrayXd fmh = f(x - h / 2);
+    ArrayXd fph = f(x + h / 2);
+
+    ArrayXd r3 = 0.5 * (fp1 - fm1);
+    ArrayXd r5 = (4.0 / 3.0) * (fph - fmh) - (1.0 / 3.0) * r3;
+
+    ArrayXd e3 = (fp1.abs() + fm1.abs()) * GSL_DBL_EPSILON;
+    ArrayXd e5 = (fph.abs() + fmh.abs()) * (2.0 * GSL_DBL_EPSILON) + e3;
+
+    /* The next term is due to finite precision in x+h = O (eps * x) */
+
+    ArrayXd dy = r3.abs().cwiseMax(r5.abs()) *
+		 (fabs(x) / h / fabs(h) * GSL_DBL_EPSILON);
+
+    /* The truncation error in the r5 approximation itself is O(h^4).
+       However, for safety, we estimate the error from r5-r3, which is
+       O(h^2).  By scaling h we will minimise this estimated error, not
+       the actual truncation error in r5. */
+
+    result = r5 / h;
+    abserr_trunc = ((r5 - r3) / h).abs(); /* Estimated truncation error O(h^2) */
+    abserr_round = (e5 / h).abs() + dy;   /* Rounding error (cancellations) */
+}
+
+} // namespace
+
+void
+array_deriv_central(function<ArrayXd(double x)> f,
+		    double x, double h,
+		    ArrayXd& r_0, ArrayXd& error)
+{
+    ArrayXd round, trunc;
+    central_deriv(f, x, h, r_0, round, trunc);
+    error = round + trunc;
+    double round_sum = round.sum(), trunc_sum = trunc.sum();
+
+    if (round_sum < trunc_sum && round_sum > 0 && trunc_sum > 0)
+    {
+	ArrayXd r_opt, round_opt, trunc_opt, error_opt;
+
+	/* Compute an optimised stepsize to minimize the total error,
+	   using the scaling of the truncation error (O(h^2)) and
+	   rounding error (O(1/h)). */
+
+	double h_opt = h * pow(round_sum / (2.0 * trunc_sum), 1.0 / 3.0);
+	central_deriv(f, x, h_opt, r_opt, round_opt, trunc_opt);
+	error_opt = round_opt + trunc_opt;
+
+	/* Check that the new error is smaller, and that the new derivative
+	   is consistent with the error bounds of the original estimate. */
+
+	double error_sum = error.sum();
+	if (error_opt.sum() < error_sum &&
+	    (r_opt - r_0).abs().sum() < 4.0 * error_sum)
+        {
+	    r_0 = r_opt;
+	    error = error_opt;
+        }
+    }
+}
+
+} // namespace flexiblesusy

src/array_deriv.hpp

@@ -0,0 +1,46 @@
+// ====================================================================
+// This file is part of FlexibleSUSY.
+//
+// FlexibleSUSY is free software: you can redistribute it and/or modify
+// it under the terms of the GNU General Public License as published
+// by the Free Software Foundation, either version 3 of the License,
+// or (at your option) any later version.
+//
+// FlexibleSUSY is distributed in the hope that it will be useful, but
+// WITHOUT ANY WARRANTY; without even the implied warranty of
+// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
+// General Public License for more details.
+//
+// You should have received a copy of the GNU General Public License
+// along with FlexibleSUSY.  If not, see
+// <http://www.gnu.org/licenses/>.
+// ====================================================================
+
+#ifndef ARRAY_DERIV_H
+#define ARRAY_DERIV_H
+
+#include <functional>
+#include <Eigen/Dense>
+
+namespace flexiblesusy {
+
+/**
+ * Takes the numerical derivative of an array-valued function using an
+ * adaptive central difference algorithm.  The error in each component
+ * of the derivative is in general larger than what would result from
+ * applying gsl_deriv_central() on the same component of f.
+ *
+ * @param f	 function to differentiate
+ * @param x	 argument to derivative of f
+ * @param h	 initial step size
+ * @param result derivative of f
+ * @param abserr array of error estimates
+ *
+ */
+void array_deriv_central(std::function<Eigen::ArrayXd(double x)> f,
+			 double x, double h,
+			 Eigen::ArrayXd& result, Eigen::ArrayXd& abserr);
+
+}
+
+#endif

src/module.mk

@@ -5,6 +5,7 @@ LIBFLEXI_MK  := \
 		$(DIR)/module.mk
 
 LIBFLEXI_SRC := \
+		$(DIR)/array_deriv.cpp \
 		$(DIR)/betafunction.cpp \
 		$(DIR)/build_info.cpp \
 		$(DIR)/command_line_options.cpp \
@@ -28,6 +29,7 @@ LIBFLEXI_SRC := \
 		$(DIR)/wrappers.cpp
 
 LIBFLEXI_HDR := \
+		$(DIR)/array_deriv.hpp \
 		$(DIR)/betafunction.hpp \
 		$(DIR)/build_info.hpp \
 		$(DIR)/cextensions.hpp \

test/module.mk

@@ -14,6 +14,7 @@ LIBTEST_DEP := \
 LIBTEST     := $(DIR)/lib$(MODNAME)$(LIBEXT)
 
 TEST_SRC := \
+		$(DIR)/test_array_deriv.cpp \
 		$(DIR)/test_logger.cpp \
 		$(DIR)/test_betafunction.cpp \
 		$(DIR)/test_goldstones.cpp \
@@ -232,6 +233,9 @@ distclean::     distclean-$(MODNAME)
 
 $(DIR)/test_lowMSSM.sh.log: $(RUN_MSSM_EXE) $(RUN_lowMSSM_EXE)
 
+$(DIR)/test_array_deriv.x: $(DIR)/test_array_deriv.o $(LIBFLEXI)
+		$(CXX) -o $@ $(call abspathx,$^) $(BOOSTTESTLIBS) $(GSLLIBS)
+
 $(DIR)/test_logger.x: $(DIR)/test_logger.o $(LIBFLEXI) $(LIBLEGACY) $(filter-out -%,$(LOOPFUNCLIBS))
 		$(CXX) -o $@ $(call abspathx,$^) $(filter -%,$(LOOPFUNCLIBS)) $(BOOSTTESTLIBS) $(FLIBS)
 

test/test_array_deriv.cpp

@@ -0,0 +1,108 @@
+// ====================================================================
+// This file is part of FlexibleSUSY.
+//
+// FlexibleSUSY is free software: you can redistribute it and/or modify
+// it under the terms of the GNU General Public License as published
+// by the Free Software Foundation, either version 3 of the License,
+// or (at your option) any later version.
+//
+// FlexibleSUSY is distributed in the hope that it will be useful, but
+// WITHOUT ANY WARRANTY; without even the implied warranty of
+// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
+// General Public License for more details.
+//
+// You should have received a copy of the GNU General Public License
+// along with FlexibleSUSY.  If not, see
+// <http://www.gnu.org/licenses/>.
+// ====================================================================
+
+#include <cmath>
+#include <gsl/gsl_math.h>
+#include <gsl/gsl_deriv.h>
+#include "array_deriv.hpp"
+
+#define BOOST_TEST_DYN_LINK
+#define BOOST_TEST_MODULE test_array_deriv
+
+#include <boost/test/unit_test.hpp>
+#include <boost/test/test_case_template.hpp>
+#include <boost/test/floating_point_comparison.hpp>
+#include <boost/mpl/list.hpp>
+
+using namespace std;
+using namespace Eigen;
+using namespace flexiblesusy;
+
+ArrayXd foo(double x)
+{
+    return (ArrayXd(2) << 7 / sqrt(x), log(x)).finished();
+}
+
+ArrayXd bar(double x)
+{
+    return (ArrayXd(4) << sin(x), exp(x), atan(x), 11 / x).finished();
+}
+
+struct Param_gsl {
+    function<ArrayXd(double x)> f;
+    size_t i;
+};
+
+double f_wrap(double x, void *params)
+{
+    Param_gsl *p = (Param_gsl *)params;
+    return p->f(x)[p->i];
+}
+
+void gsl_array_deriv_central(function<ArrayXd(double x)> f,
+			     double x, double h,
+			     ArrayXd& result, ArrayXd& abserr)
+{
+    Param_gsl p;
+    p.f = f;
+
+    gsl_function f_gsl;
+    f_gsl.function = f_wrap;
+    f_gsl.params = &p;
+
+    size_t n = f(x).size();
+    result.resize(n);
+    abserr.resize(n);
+
+    for (size_t i = 0; i < n; i++) {
+	p.i = i;
+	gsl_deriv_central(&f_gsl, x, h, &result[i], &abserr[i]);
+    }
+}
+
+template<ArrayXd f_(double), int sign_>
+struct Test_deriv {
+    enum { sign = sign_ };
+    ArrayXd (*fp())(double) { return f_; };
+};
+
+typedef boost::mpl::list<
+    Test_deriv<foo,  1>,
+    Test_deriv<bar,  1>,
+    Test_deriv<bar, -1>
+> deriv_tests;
+
+BOOST_AUTO_TEST_CASE_TEMPLATE(test_array_deriv_central, T, deriv_tests)
+{
+    int sign = T::sign;
+
+    for (size_t n = 10000; n; n--) {
+	double x = sign * (1e2 * rand() / RAND_MAX + 1);
+
+	ArrayXd derivs_arr, abserrs_arr;
+	array_deriv_central
+	(T().fp(), x, GSL_SQRT_DBL_EPSILON, derivs_arr, abserrs_arr);
+
+	ArrayXd derivs_gsl, abserrs_gsl;
+	gsl_array_deriv_central
+	(T().fp(), x, GSL_SQRT_DBL_EPSILON, derivs_gsl, abserrs_gsl);
+
+	BOOST_CHECK(((derivs_arr - derivs_gsl).abs() < abserrs_arr).all());
+	BOOST_WARN (((derivs_arr - derivs_gsl).abs() < abserrs_gsl).all());
+    }
+}