Sampling

.gitignore

@@ -108,4 +108,11 @@ deploy.sh
 x.py
 run.py
 install
-.vs/
+.vs/
+
+scripts/distributions/*pdf
+scripts/distributions/*png
+scripts/distributions/*pkl
+scripts/distributions/figures
+
+data_old

include/bounds.hpp

@@ -11,11 +11,13 @@ namespace parameters
 	struct Parameters;
 }
 
-
 namespace bounds
 {
 	using Mask = Eigen::Array<bool, Eigen::Dynamic, 1>;
 
+	Mask is_out_of_bounds(const Vector &xi, const Vector &lb, const Vector &ub);
+	bool any_out_of_bounds(const Vector &xi, const Vector &lb, const Vector &ub);
+
 	struct BoundCorrection
 	{
 		virtual ~BoundCorrection() = default;
@@ -26,68 +28,73 @@ namespace bounds
 		BoundCorrection(const Vector &lb, const Vector &ub) : lb(lb), ub(ub), db(ub - lb),
 															  diameter((ub - lb).norm()) {}
 
-		void correct(const Eigen::Index i, parameters::Parameters& p);
+		void correct(const Eigen::Index i, parameters::Parameters &p);
 
-		virtual Vector correct_x(const Vector& xi, const Mask& oob) = 0;
+		virtual Vector correct_x(const Vector &xi, const Mask &oob) = 0;
 
-		[[nodiscard]] Mask is_out_of_bounds(const Vector& xi) const;
+		[[nodiscard]] Mask is_out_of_bounds(const Vector &xi) const;
 
-		[[nodiscard]] Vector delta_out_of_bounds(const Vector& xi, const Mask& oob) const;
+		[[nodiscard]] Vector delta_out_of_bounds(const Vector &xi, const Mask &oob) const;
 
 		[[nodiscard]] bool any_out_of_bounds() const
 		{
 			return n_out_of_bounds > 0;
-		}		
+		}
 	};
 
-	struct NoCorrection final : BoundCorrection
+	struct NoCorrection : BoundCorrection
 	{
 		using BoundCorrection::BoundCorrection;
 
-		Vector correct_x(const Vector& xi, const Mask& oob) override
+		Vector correct_x(const Vector &xi, const Mask &oob) override
 		{
 			return xi;
 		}
 	};
 
+	struct Resample final : NoCorrection
+	{
+		using NoCorrection::NoCorrection;
+	};
+
 	struct COTN final : BoundCorrection
 	{
 		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.)) {}
 
-		Vector correct_x(const Vector& xi, const Mask& oob) override;
+		Vector correct_x(const Vector &xi, const Mask &oob) override;
 	};
 
 	struct Mirror final : BoundCorrection
 	{
 		using BoundCorrection::BoundCorrection;
 
-		Vector correct_x(const Vector& xi, const Mask& oob) override;
+		Vector correct_x(const Vector &xi, const Mask &oob) override;
 	};
 
 	struct UniformResample final : BoundCorrection
 	{
-		sampling::Random<std::uniform_real_distribution<>> sampler;
+		sampling::Uniform sampler;
 
 		UniformResample(Eigen::Ref<const Vector> lb, Eigen::Ref<const Vector> ub) : BoundCorrection(lb, ub), sampler(static_cast<size_t>(lb.size())) {}
 
-		Vector correct_x(const Vector& xi, const Mask& oob) override;
+		Vector correct_x(const Vector &xi, const Mask &oob) override;
 	};
 
 	struct Saturate final : BoundCorrection
 	{
 		using BoundCorrection::BoundCorrection;
 
-		Vector correct_x(const Vector& xi, const Mask& oob) override;
+		Vector correct_x(const Vector &xi, const Mask &oob) override;
 	};
 
 	struct Toroidal final : BoundCorrection
 	{
 		using BoundCorrection::BoundCorrection;
 
-		Vector correct_x(const Vector& xi, const Mask& oob) override;
-	};	
+		Vector correct_x(const Vector &xi, const Mask &oob) override;
+	};
 
 	inline std::shared_ptr<BoundCorrection> get(const parameters::CorrectionMethod &m, const Vector &lb, const Vector &ub)
 	{
@@ -104,7 +111,8 @@ namespace bounds
 			return std::make_shared<Saturate>(lb, ub);
 		case CorrectionMethod::TOROIDAL:
 			return std::make_shared<Toroidal>(lb, ub);
-
+		case CorrectionMethod::RESAMPLE:
+			return std::make_shared<Resample>(lb, ub);
 		default:
 		case CorrectionMethod::NONE:
 			return std::make_shared<NoCorrection>(lb, ub);

include/common.hpp

@@ -273,7 +273,7 @@ namespace rng
 	struct uniform
 	{
 		/**
-		 * @brief Generate a random uniform number in the closed interval [-1, 1]
+		 * @brief Generate a random uniform number in the closed interval [0, 1]
 		 *
 		 * @tparam G the type of the generator
 		 * @param gen the generator instance
@@ -282,7 +282,7 @@ namespace rng
 		template <typename G>
 		T operator()(G &gen)
 		{
-			return static_cast<T>(2.0 * gen() - gen.min()) / gen.max() - gen.min() - 1;
+			return static_cast<T>(gen() - gen.min()) / gen.max() - gen.min();
 		}
 	};
 
@@ -320,8 +320,8 @@ namespace rng
 
 			if (generate)
 			{
-				T u1 = std::abs(rng(gen));
-				T u2 = std::abs(rng(gen));
+				T u1 = rng(gen);
+				T u2 = rng(gen);
 				const T root_log_u1 = std::sqrt(-2.0 * std::log(u1));
 				const T two_pi_u2 = 2.0 * M_PI * u2;
 				r1 = (sigma * (root_log_u1 * std::sin(two_pi_u2))) + mu;

include/es.hpp

@@ -0,0 +1,104 @@
+#pragma once
+
+#include "sampling.hpp"
+#include "stats.hpp"
+#include "bounds.hpp"
+
+namespace es
+{
+    struct OnePlusOneES
+    {
+        OnePlusOneES(
+            const size_t d,
+            const Vector &x0,
+            const double f0,
+            const double sigma0,
+            const size_t budget,
+            const double target,
+            const parameters::Modules &modules)
+            : d(d), sigma(sigma0), decay(1.0 / std::sqrt(static_cast<double>(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)),
+              corrector(bounds::get(modules.bound_correction, Vector::Ones(d) * -5.0, Vector::Ones(d) * 5.0))
+        {
+        }
+
+        Vector sample();
+        void step(FunctionType &objective);
+        void operator()(FunctionType &objective);
+
+        size_t d;
+        double sigma;
+        double decay;
+        Vector x;
+        double f;
+        size_t t;
+        size_t budget;
+        double target;
+        bool rejection_sampling;
+
+        std::shared_ptr<sampling::Sampler> sampler;
+        std::shared_ptr<bounds::BoundCorrection> corrector;
+    };
+
+    struct MuCommaLambdaES
+    {
+
+        MuCommaLambdaES(
+            const size_t d,
+            const Vector &x0,
+            const double sigma0, 
+            const size_t budget,
+            const double 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)), 
+              mu_inv(1.0 / mu),
+              m(x0), sigma(sigma0 * Vector::Ones(d)),
+              f(Vector::Constant(lambda, std::numeric_limits<double>::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())),
+              t(0), e(0), budget(budget), target(target),
+              sampler(sampling::get(d, modules, lambda)),
+              sigma_sampler(std::make_shared<sampling::Gaussian>(d)),
+              rejection_sampling(modules.bound_correction == parameters::CorrectionMethod::RESAMPLE),
+              corrector(bounds::get(modules.bound_correction, Vector::Ones(d) * -5.0, Vector::Ones(d) * 5.0))
+        {
+            // tau = 1.0 / sampler->expected_length();
+            // tau_i = 1.0 / std::sqrt(sampler->expected_length());
+        }
+
+        Vector sample(const Vector si);
+        void step(FunctionType &objective);
+        void operator()(FunctionType &objective);
+
+        size_t d;
+        size_t lambda;
+        size_t mu;
+        double tau;
+        double tau_i;
+        double mu_inv;
+
+        Vector m;
+        Vector sigma;
+        Vector f;
+        Matrix X;
+        Matrix S;
+
+        double f_min;
+        Vector x_min;
+        size_t t;
+        size_t e;
+        size_t budget;
+        double target;
+
+        std::shared_ptr<sampling::Sampler> sampler;
+        std::shared_ptr<sampling::Sampler> sigma_sampler;
+
+        bool rejection_sampling;
+        std::shared_ptr<bounds::BoundCorrection> corrector;
+    };
+}

include/matrix_adaptation.hpp

@@ -12,14 +12,13 @@ namespace matrix_adaptation
 	{
 		Vector m, m_old, dm, ps;
 		double dd;
-		double chiN;
+		double expected_length_z;
 		Matrix inv_C;
 
-		Adaptation(const size_t dim, const Vector& x0, const Vector& ps) :
+		Adaptation(const size_t dim, const Vector& x0, const Vector& ps, const double expected_length_z) :
 			m(x0), m_old(dim), dm(Vector::Zero(dim)),
 			ps(ps), dd(static_cast<double>(dim)),
-			chiN(sqrt(dd)* (1.0 - 1.0 / (4.0 * dd) + 1.0
-				/ (21.0 * pow(dd, 2.0)))),
+			expected_length_z(expected_length_z),
 			inv_C(Matrix::Identity(dim, dim))
 		{
 		}
@@ -43,7 +42,7 @@ namespace matrix_adaptation
 
 	struct None final : Adaptation
 	{
-		None(const size_t dim, const Vector& x0) : Adaptation(dim, x0, Vector::Ones(dim))
+		None(const size_t dim, const Vector& x0, const double expected_length_z) : Adaptation(dim, x0, Vector::Ones(dim), expected_length_z)
 		{
 		}
 
@@ -72,7 +71,7 @@ namespace matrix_adaptation
 
 		bool hs = true;
 
-		CovarianceAdaptation(const size_t dim, const Vector& x0) : Adaptation(dim, x0, Vector::Zero(dim)),
+		CovarianceAdaptation(const size_t dim, const Vector& x0, const double expected_length_z) : Adaptation(dim, x0, Vector::Zero(dim), expected_length_z),
 			pc(Vector::Zero(dim)),
 			d(Vector::Ones(dim)),
 			B(Matrix::Identity(dim, dim)),
@@ -112,7 +111,7 @@ namespace matrix_adaptation
 		Matrix M;
 		Matrix M_inv;
 
-		MatrixAdaptation(const size_t dim, const Vector& x0) : Adaptation(dim, x0, Vector::Ones(dim)),
+		MatrixAdaptation(const size_t dim, const Vector& x0, const double expected_length_z) : Adaptation(dim, x0, Vector::Ones(dim), expected_length_z),
 			M(Matrix::Identity(dim, dim)),
 			M_inv(Matrix::Identity(dim, dim))
 		{
@@ -132,20 +131,20 @@ namespace matrix_adaptation
 		Vector invert_y(const Vector&) override;
 	};
 
-	inline std::shared_ptr<Adaptation> get(const parameters::Modules& m, const size_t dim, const Vector& x0)
+	inline std::shared_ptr<Adaptation> get(const parameters::Modules& m, const size_t dim, const Vector& x0, const double expected_z)
 	{
 		using namespace parameters;
 		switch (m.matrix_adaptation)
 		{
 		case MatrixAdaptationType::MATRIX:
-			return std::make_shared<MatrixAdaptation>(dim, x0);
+			return std::make_shared<MatrixAdaptation>(dim, x0, expected_z);
 		case MatrixAdaptationType::NONE:
-			return std::make_shared<None>(dim, x0);
+			return std::make_shared<None>(dim, x0, expected_z);
 		case MatrixAdaptationType::SEPERABLE:
-			return std::make_shared<SeperableAdaptation>(dim, x0);
+			return std::make_shared<SeperableAdaptation>(dim, x0, expected_z);
 		default:
 		case MatrixAdaptationType::COVARIANCE:
-			return std::make_shared<CovarianceAdaptation>(dim, x0);
+			return std::make_shared<CovarianceAdaptation>(dim, x0, expected_z);
 		}
 	}
 } // namespace parameters