Better floating point precision

include/bounds.hpp

@@ -22,7 +22,7 @@ namespace bounds
 	{
 		virtual ~BoundCorrection() = default;
 		Vector lb, ub, db;
-		double diameter;
+		Float diameter;
 		size_t n_out_of_bounds = 0;
 
 		BoundCorrection(const Vector &lb, const Vector &ub) : lb(lb), ub(ub), db(ub - lb),
@@ -61,7 +61,7 @@ namespace bounds
 	{
 		sampling::Gaussian sampler;
 
-		COTN(Eigen::Ref<const Vector> lb, Eigen::Ref<const Vector> ub) : BoundCorrection(lb, ub), sampler(static_cast<size_t>(lb.size()), rng::normal<double>(0, 1.0 / 3.)) {}
+		COTN(Eigen::Ref<const Vector> lb, Eigen::Ref<const Vector> ub) : BoundCorrection(lb, ub), sampler(static_cast<size_t>(lb.size()), rng::normal<Float>(0, 1.0 / 3.)) {}
 
 		Vector correct_x(const Vector &xi, const Mask &oob) override;
 	};

include/common.hpp

@@ -22,45 +22,49 @@
 #include <Eigen/Eigenvalues>
 #include <Eigen/QR>
 
-using Matrix = Eigen::MatrixXd;
-using Vector = Eigen::VectorXd;
-using Array = Eigen::ArrayXd;
+using Float = long double;
+using Matrix = Eigen::Matrix<Float, -1, -1>;
+using Vector = Eigen::Matrix<Float, -1, 1>;
+using Array = Eigen::Array<Float, -1, 1>;
 using size_to = std::optional<size_t>;
 
 template <typename T>
 std::ostream &operator<<(std::ostream &os, const std::vector<T> &x);
 
-using FunctionType = std::function<double(const Vector &)>;
+using FunctionType = std::function<Float(const Vector &)>;
 
 namespace constants
 {
-	extern double tolup_sigma;
-	extern double tol_condition_cov;
-	extern double tol_min_sigma;
-	extern double stagnation_quantile;
-	extern double sigma_threshold;
+	extern Float tolup_sigma;
+	extern Float tol_condition_cov;
+	extern Float tol_min_sigma;
+	extern Float stagnation_quantile;
+	extern Float sigma_threshold;
 	extern size_t cache_max_doubles;
 	extern size_t cache_min_samples;
 	extern bool cache_samples;
+	extern Float lb_sigma;
+	extern Float ub_sigma; 
+	extern bool clip_sigma; 
 }
 
 /**
  * @brief Cdf of a standard normal distribution.
  *
  * see: ndtr_ndtri.cpp
  * @param x lower tail of the probabilty
- * @return double quantile corresponding to the lower tail probability q
+ * @return Float quantile corresponding to the lower tail probability q
  */
-double cdf(const double x);
+Float cdf(const Float x);
 
 /**
  * @brief Percent point function (inverse of cdf) of a standard normal distribution.
  *
  * see: ndtri.cpp
  * @param x lower tail of the probabilty
- * @return double quantile corresponding to the lower tail probability q
+ * @return Float quantile corresponding to the lower tail probability q
  */
-double ppf(const double x);
+Float ppf(const Float x);
 
 /**
  * @brief Generate a sobol sequence using 8 byte integer numbers.
@@ -70,24 +74,24 @@ double ppf(const double x);
  * @param seed The current seed of the sobol sequence
  * @param quasi the vector of random numbers in which to place the output
  */
-void i8_sobol(int dim_num, long long int *seed, double quasi[]);
+void i8_sobol(int dim_num, long long int *seed, Float quasi[]);
 
 struct Solution
 {
 	//! Coordinates
 	Vector x;
 	//! Function value
-	double y;
+	Float y;
 	//! Generation
 	size_t t;
 	//! Evaluations
 	size_t e;
 
-	Solution(const Vector &x, const double y, const size_t t = 0, const size_t e = 0) : x(x), y(y), t(t), e(e)
+	Solution(const Vector &x, const Float y, const size_t t = 0, const size_t e = 0) : x(x), y(y), t(t), e(e)
 	{
 	}
 
-	Solution() : Solution({}, std::numeric_limits<double>::infinity()) {}
+	Solution() : Solution({}, std::numeric_limits<Float>::infinity()) {}
 
 	[[nodiscard]] size_t n() const
 	{
@@ -171,9 +175,9 @@ namespace utils
 	 *
 	 * @param running_times the vector of measured running times
 	 * @param budget the maximum budget allocated to each run
-	 * @return std::pair<double, size_t> (ERT, number of successfull runs)
+	 * @return std::pair<Float, size_t> (ERT, number of successfull runs)
 	 */
-	std::pair<double, size_t> compute_ert(const std::vector<size_t> &running_times, size_t budget);
+	std::pair<Float, size_t> compute_ert(const std::vector<size_t> &running_times, size_t budget);
 
 	/**
 	 * \brief calculate the nearest power of two
@@ -184,7 +188,7 @@ namespace utils
 	template<typename T>
 	T nearest_power_of_2(const T value)
 	{
-		const double val = static_cast<double>(value);
+		const Float val = static_cast<Float>(value);
 		return static_cast<T>(pow(2.0, std::floor(std::log2(val))));
 	}
 
@@ -251,14 +255,14 @@ namespace rng
 		size_t dim;
 		size_t n_samples;
 
-		std::vector<double> cache;
+		std::vector<Float> cache;
 		Shuffler shuffler;
 
 		CachedShuffleSequence(const size_t d);
 
-		void fill(const std::vector<double>& c);
+		void fill(const std::vector<Float>& c);
 
-		void transform(const std::function<double(double)>& f);
+		void transform(const std::function<Float(Float)>& f);
 
 		Vector get_index(const size_t idx);
 
@@ -269,7 +273,7 @@ namespace rng
 	 * @brief distribution which in combination with mt19997 produces the same
 	 * random numbers for gcc and msvc
 	 */
-	template <typename T = double>
+	template <typename T = Float>
 	struct uniform
 	{
 		/**
@@ -290,7 +294,7 @@ namespace rng
 	 * @brief Box-Muller random normal number generator. Ensures similar numbers generated
 	 * on different operating systems.
 	 */
-	template <typename T = double>
+	template <typename T = Float>
 	struct normal
 	{
 		T mu;
@@ -314,7 +318,7 @@ namespace rng
 		template <typename G>
 		T operator()(G &gen)
 		{
-			static uniform<double> rng;
+			static uniform<Float> rng;
 			static T r1, r2;
 			static bool generate = true;
 
@@ -338,6 +342,7 @@ namespace rng
 
 namespace functions
 {
-	double sphere(const Vector &x);
-	double rastrigin(const Vector &x);
+	Float sphere(const Vector &x);
+	Float rastrigin(const Vector &x);
+	Float ellipse(const Vector& x);
 }

include/es.hpp

@@ -11,12 +11,12 @@ namespace es
         OnePlusOneES(
             const size_t d,
             const Vector &x0,
-            const double f0,
-            const double sigma0,
+            const Float f0,
+            const Float sigma0,
             const size_t budget,
-            const double target,
+            const Float target,
             const parameters::Modules &modules)
-            : d(d), sigma(sigma0), decay(1.0 / std::sqrt(static_cast<double>(d) + 1)),
+            : d(d), sigma(sigma0), decay(1.0 / std::sqrt(static_cast<Float>(d) + 1)),
               x(x0), f(f0), t(1), budget(budget), target(target),
               rejection_sampling(modules.bound_correction == parameters::CorrectionMethod::RESAMPLE),
               sampler(sampling::get(d, modules, 1)),
@@ -29,13 +29,13 @@ namespace es
         void operator()(FunctionType &objective);
 
         size_t d;
-        double sigma;
-        double decay;
+        Float sigma;
+        Float decay;
         Vector x;
-        double f;
+        Float f;
         size_t t;
         size_t budget;
-        double target;
+        Float target;
         bool rejection_sampling;
 
         std::shared_ptr<sampling::Sampler> sampler;
@@ -48,19 +48,19 @@ namespace es
         MuCommaLambdaES(
             const size_t d,
             const Vector &x0,
-            const double sigma0, 
+            const Float sigma0, 
             const size_t budget,
-            const double target,
+            const Float target,
             const parameters::Modules &modules)
             : d(d), lambda(d * 5), mu(std::floor(lambda / 4)),
-              tau(1.0 / std::sqrt(static_cast<double>(d))), // v1 ->
-              tau_i(1.0 / pow(static_cast<double>(d), .25)), 
+              tau(1.0 / std::sqrt(static_cast<Float>(d))), // v1 ->
+              tau_i(1.0 / pow(static_cast<Float>(d), .25)), 
               mu_inv(1.0 / mu),
               m(x0), sigma(sigma0 * Vector::Ones(d)),
-              f(Vector::Constant(lambda, std::numeric_limits<double>::infinity())),
+              f(Vector::Constant(lambda, std::numeric_limits<Float>::infinity())),
               X(d, lambda), S(d, lambda),
-              f_min(std::numeric_limits<double>::infinity()),
-              x_min(Vector::Constant(d, std::numeric_limits<double>::signaling_NaN())),
+              f_min(std::numeric_limits<Float>::infinity()),
+              x_min(Vector::Constant(d, std::numeric_limits<Float>::signaling_NaN())),
               t(0), e(0), budget(budget), target(target),
               sampler(sampling::get(d, modules, lambda)),
               sigma_sampler(std::make_shared<sampling::Gaussian>(d)),
@@ -78,22 +78,22 @@ namespace es
         size_t d;
         size_t lambda;
         size_t mu;
-        double tau;
-        double tau_i;
-        double mu_inv;
+        Float tau;
+        Float tau_i;
+        Float mu_inv;
 
         Vector m;
         Vector sigma;
         Vector f;
         Matrix X;
         Matrix S;
 
-        double f_min;
+        Float f_min;
         Vector x_min;
         size_t t;
         size_t e;
         size_t budget;
-        double target;
+        Float target;
 
         std::shared_ptr<sampling::Sampler> sampler;
         std::shared_ptr<sampling::Sampler> sigma_sampler;