From 07d2c4921dde67da38b003a1778cc76167271871 Mon Sep 17 00:00:00 2001 From: dimitrispapaxristos22-stack Date: Mon, 2 Feb 2026 20:10:34 +0200 Subject: [PATCH] Add files via upload MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Προσθήκη HHO και Alpine --- usermethod.cpp | 407 ++++++++++++++++++++++++++++++++++++++++++++++++ usermethod.h | 85 ++++++++++ userproblem.cpp | 47 ++++++ userproblem.h | 16 ++ 4 files changed, 555 insertions(+) create mode 100644 usermethod.cpp create mode 100644 usermethod.h create mode 100644 userproblem.cpp create mode 100644 userproblem.h diff --git a/usermethod.cpp b/usermethod.cpp new file mode 100644 index 00000000..bae69705 --- /dev/null +++ b/usermethod.cpp @@ -0,0 +1,407 @@ +#include "usermethod.h" +#include +#include +#include + +static inline double absd(double x) { return (x < 0.0) ? -x : x; } +static inline double clampd(double x, double lo, double hi) { + if (x < lo) return lo; + if (x > hi) return hi; + return x; +} + +UserMethod::UserMethod(Problem *p) +{ + myProblem = p; + + // Defaults (must be overridable via setters / external params) + popsize = 30; + max_iters = 500; + seed = 1; + sampling_mode = "uniform"; // "uniform" or "gaussian" + + max_evals = 0; // if 0 -> set in init() as popsize * max_iters + tol = 1e-8; + tol_window = 10; + max_stall = 50; + + p_local = 0.2; + local_every_R = 20; + top_k = 3; + local_steps = 50; + local_lr = 0.05; + + iter = 0; + stall = 0; + tol_hits = 0; + + besty = 1e100; + besty_old = 1e100; +} + +double UserMethod::rand01() const +{ + double r = (double)rand() / (double)RAND_MAX; + if (r < 1e-12) r = 1e-12; + if (r > 1.0 - 1e-12) r = 1.0 - 1e-12; + return r; +} + +double UserMethod::randUniform(double a, double b) const +{ + return a + (b - a) * rand01(); +} + +double UserMethod::randNormal(double mean, double stddev) const +{ + // Box–Muller transform + const double u1 = rand01(); + const double u2 = rand01(); + const double pi = 3.14159265358979323846; + const double z0 = std::sqrt(-2.0 * std::log(u1)) * std::cos(2.0 * pi * u2); + return mean + stddev * z0; +} + +void UserMethod::enforceBounds(Data &x) const +{ + for (int j = 0; j < (int)x.size(); ++j) + x[j] = clampd(x[j], L[j], R[j]); +} + +int UserMethod::bestIndex() const +{ + int b = 0; + for (int i = 1; i < (int)fit.size(); ++i) + if (fit[i] < fit[b]) b = i; + return b; +} + +Data UserMethod::levyFlight(int d, double beta) const +{ + // Mantegna's algorithm for Levy flight steps + const double pi = 3.14159265358979323846; + + const double num = std::tgamma(1.0 + beta) * std::sin(pi * beta / 2.0); + const double den = std::tgamma((1.0 + beta) / 2.0) * beta * std::pow(2.0, (beta - 1.0) / 2.0); + const double sigma_u = std::pow(num / den, 1.0 / beta); + + Data step; + step.resize(d); + + for (int j = 0; j < d; ++j) + { + const double u = randNormal(0.0, sigma_u); + const double v = randNormal(0.0, 1.0); + step[j] = u / std::pow(absd(v) + 1e-12, 1.0 / beta); + } + return step; +} + +void UserMethod::localOptimize(Data &x, double &y) +{ + // Simple projected gradient descent (hybrid local optimization) + for (int s = 0; s < local_steps; ++s) + { + Data g = myProblem->gradient(x); + + Data xn = x; + for (int j = 0; j < (int)xn.size(); ++j) + xn[j] = xn[j] - local_lr * g[j]; + + enforceBounds(xn); + + double yn = myProblem->statFunmin(xn); + if (yn + tol < y) + { + x = xn; + y = yn; + } + else + { + break; + } + } +} + +// ---------------- setters ---------------- +void UserMethod::setPopsize(int n) { if (n > 1) popsize = n; } +void UserMethod::setMaxIters(int it) { if (it > 0) max_iters = it; } +void UserMethod::setSeed(int s) { seed = s; } +void UserMethod::setSamplingMode(const std::string &m) { sampling_mode = m; } + +void UserMethod::setMaxEvals(int e) { if (e > 0) max_evals = e; } +void UserMethod::setTol(double t) { if (t > 0.0) tol = t; } +void UserMethod::setTolWindow(int g) { if (g > 0) tol_window = g; } +void UserMethod::setMaxStall(int s) { if (s > 0) max_stall = s; } + +void UserMethod::setPLocal(double p) { if (p >= 0.0 && p <= 1.0) p_local = p; } +void UserMethod::setLocalEveryR(int r) { if (r >= 0) local_every_R = r; } +void UserMethod::setTopK(int k) { if (k > 0) top_k = k; } +void UserMethod::setLocalSteps(int s) { if (s > 0) local_steps = s; } +void UserMethod::setLocalLr(double lr) { if (lr > 0.0) local_lr = lr; } + +// ---------------- lifecycle ---------------- +void UserMethod::init() +{ + srand(seed); + + const int d = myProblem->getDimension(); + L = myProblem->getLeftMargin(); + R = myProblem->getRightMargin(); + + if (max_evals <= 0) max_evals = popsize * max_iters; + + hawks.assign(popsize, Data(d, 0.0)); + fit.assign(popsize, 1e100); + + besty = 1e100; + besty_old = 1e100; + bestx.assign(d, 0.0); + + iter = 0; + stall = 0; + tol_hits = 0; + + // Initial sampling (two modes) + for (int i = 0; i < popsize; ++i) + { + Data x; + if (sampling_mode == "gaussian") + { + x.resize(d); + for (int j = 0; j < d; ++j) + { + const double mu = 0.5 * (L[j] + R[j]); + const double sd = (R[j] - L[j]) / 6.0; + x[j] = randNormal(mu, sd); + } + enforceBounds(x); + } + else + { + // uniform + x = myProblem->getSample(); + } + + const double y = myProblem->statFunmin(x); + hawks[i] = x; + fit[i] = y; + + if (y < besty) + { + besty = y; + bestx = x; + } + } +} + +void UserMethod::step() +{ + const int d = myProblem->getDimension(); + const Data rabbit = bestx; + + // mean position Xm + Data Xm(d, 0.0); + for (int i = 0; i < popsize; ++i) + for (int j = 0; j < d; ++j) + Xm[j] += hawks[i][j]; + for (int j = 0; j < d; ++j) + Xm[j] /= (double)popsize; + + const double T = (double)max_iters; + const double t = (double)iter; + + for (int i = 0; i < popsize; ++i) + { + const double E0 = 2.0 * rand01() - 1.0; // [-1,1] + const double E = 2.0 * E0 * (1.0 - t / T); // escaping energy + const double q = rand01(); + const double r = rand01(); + + Data Xnew = hawks[i]; + + if (absd(E) >= 1.0) + { + // Exploration phase + const int rand_idx = (int)(rand01() * popsize) % popsize; + const Data Xrand = hawks[rand_idx]; + + if (q >= 0.5) + { + const double r1 = rand01(); + const double r2 = rand01(); + for (int j = 0; j < d; ++j) + { + const double term = absd(Xrand[j] - 2.0 * r2 * hawks[i][j]); + Xnew[j] = Xrand[j] - r1 * term; + } + } + else + { + const double r3 = rand01(); + const double r4 = rand01(); + for (int j = 0; j < d; ++j) + { + const double randPos = L[j] + r4 * (R[j] - L[j]); + Xnew[j] = (rabbit[j] - Xm[j]) - r3 * randPos; + } + } + + enforceBounds(Xnew); + const double fnew = myProblem->statFunmin(Xnew); + + if (fnew < fit[i]) + { + hawks[i] = Xnew; + fit[i] = fnew; + } + continue; + } + + // Exploitation phase + const double J = 2.0 * (1.0 - rand01()); // jump strength in (0,2) + + if (r >= 0.5 && absd(E) >= 0.5) + { + // Soft besiege + for (int j = 0; j < d; ++j) + { + const double dX = rabbit[j] - hawks[i][j]; + Xnew[j] = dX - E * absd(J * rabbit[j] - hawks[i][j]); + } + enforceBounds(Xnew); + const double fnew = myProblem->statFunmin(Xnew); + if (fnew < fit[i]) { hawks[i] = Xnew; fit[i] = fnew; } + } + else if (r >= 0.5 && absd(E) < 0.5) + { + // Hard besiege + for (int j = 0; j < d; ++j) + { + const double dX = rabbit[j] - hawks[i][j]; + Xnew[j] = rabbit[j] - E * absd(dX); + } + enforceBounds(Xnew); + const double fnew = myProblem->statFunmin(Xnew); + if (fnew < fit[i]) { hawks[i] = Xnew; fit[i] = fnew; } + } + else if (r < 0.5 && absd(E) >= 0.5) + { + // Soft besiege + rapid dives + Data Y(d, 0.0); + for (int j = 0; j < d; ++j) + Y[j] = rabbit[j] - E * absd(J * rabbit[j] - hawks[i][j]); + enforceBounds(Y); + + Data Lf = levyFlight(d, 1.5); + Data Z(d, 0.0); + for (int j = 0; j < d; ++j) + Z[j] = Y[j] + rand01() * Lf[j] * (Y[j] - hawks[i][j]); + enforceBounds(Z); + + const double fY = myProblem->statFunmin(Y); + const double fZ = myProblem->statFunmin(Z); + + Data bestCand = hawks[i]; + double fbest = fit[i]; + if (fY < fbest) { bestCand = Y; fbest = fY; } + if (fZ < fbest) { bestCand = Z; fbest = fZ; } + + hawks[i] = bestCand; + fit[i] = fbest; + } + else + { + // Hard besiege + rapid dives + Data Y(d, 0.0); + for (int j = 0; j < d; ++j) + Y[j] = rabbit[j] - E * absd(J * rabbit[j] - Xm[j]); + enforceBounds(Y); + + Data Lf = levyFlight(d, 1.5); + Data Z(d, 0.0); + for (int j = 0; j < d; ++j) + Z[j] = Y[j] + rand01() * Lf[j] * (Y[j] - hawks[i][j]); + enforceBounds(Z); + + const double fY = myProblem->statFunmin(Y); + const double fZ = myProblem->statFunmin(Z); + + Data bestCand = hawks[i]; + double fbest = fit[i]; + if (fY < fbest) { bestCand = Y; fbest = fY; } + if (fZ < fbest) { bestCand = Z; fbest = fZ; } + + hawks[i] = bestCand; + fit[i] = fbest; + } + } + + // update best + besty_old = besty; + const int b = bestIndex(); + if (fit[b] < besty) + { + besty = fit[b]; + bestx = hawks[b]; + } + + // tol window + if (absd(besty - besty_old) <= tol) tol_hits++; + else tol_hits = 0; + + // stagnation + if (besty + tol < besty_old) stall = 0; + else stall++; + + // Hybrid local optimization every R iters, on top-k, with probability p_local + if (local_every_R > 0 && (iter % local_every_R == 0) && rand01() < p_local) + { + std::vector idx(popsize); + for (int i = 0; i < popsize; ++i) idx[i] = i; + + std::partial_sort(idx.begin(), + idx.begin() + std::min(top_k, popsize), + idx.end(), + [&](int a, int b){ return fit[a] < fit[b]; }); + + const int kk = std::min(top_k, popsize); + for (int k = 0; k < kk; ++k) + { + const int ii = idx[k]; + localOptimize(hawks[ii], fit[ii]); + if (fit[ii] < besty) { besty = fit[ii]; bestx = hawks[ii]; } + } + } + + iter++; +} + +bool UserMethod::terminated() +{ + if (iter >= max_iters) return true; + if (max_evals > 0 && myProblem->getFunctionCalls() >= max_evals) return true; + if (tol_hits >= tol_window) return true; + if (stall >= max_stall) return true; + return false; +} + +void UserMethod::done() +{ + // apply local optimization to the best point + localOptimize(bestx, besty); +} + +void UserMethod::solve() +{ + init(); + while (!terminated()) + step(); + done(); +} + +Data UserMethod::getBestx() const { return bestx; } +double UserMethod::getBesty() const { return besty; } + +UserMethod::~UserMethod() {} + diff --git a/usermethod.h b/usermethod.h new file mode 100644 index 00000000..5ea19199 --- /dev/null +++ b/usermethod.h @@ -0,0 +1,85 @@ +#ifndef USERMETHOD_H +#define USERMETHOD_H + +#include "problem.h" +#include +#include + +class UserMethod +{ +private: + Problem *myProblem; + + Matrix hawks; + std::vector fit; + + Data bestx; + double besty; + double besty_old; + + Data L, R; + + int popsize; + int max_iters; + int seed; + std::string sampling_mode; // "uniform" or "gaussian" + + int max_evals; + double tol; + int tol_window; + int max_stall; + + double p_local; + int local_every_R; + int top_k; + int local_steps; + double local_lr; + + int iter; + int stall; + int tol_hits; + + double rand01() const; + double randUniform(double a, double b) const; + double randNormal(double mean, double stddev) const; + + void enforceBounds(Data &x) const; + int bestIndex() const; + + Data levyFlight(int d, double beta) const; + void localOptimize(Data &x, double &y); + +public: + explicit UserMethod(Problem *p); + + void setPopsize(int n); + void setMaxIters(int it); + void setSeed(int s); + void setSamplingMode(const std::string &m); + + void setMaxEvals(int e); + void setTol(double t); + void setTolWindow(int g); + void setMaxStall(int s); + + void setPLocal(double p); + void setLocalEveryR(int r); + void setTopK(int k); + void setLocalSteps(int s); + void setLocalLr(double lr); + + void init(); + void step(); + bool terminated(); + void done(); + + void solve(); + + Data getBestx() const; + double getBesty() const; + + ~UserMethod(); +}; + +#endif + diff --git a/userproblem.cpp b/userproblem.cpp new file mode 100644 index 00000000..72ab7215 --- /dev/null +++ b/userproblem.cpp @@ -0,0 +1,47 @@ +#include "userproblem.h" +#include + +static inline double absd(double x) { return (x < 0.0) ? -x : x; } + +UserProblem::UserProblem(int n) : Problem(n) +{ + Data l(n, -10.0); + Data r(n, 10.0); + setLeftMargin(l); + setRightMargin(r); +} + +double UserProblem::funmin(Data &x) +{ + double s = 0.0; + for (int i = 0; i < (int)x.size(); ++i) + { + double xi = x[i]; + double gi = xi * std::sin(xi) + 0.1 * xi; + s += absd(gi); + } + return s; +} + +Data UserProblem::gradient(Data &x) +{ + Data g; + g.resize(x.size()); + const double eps = 1e-12; + + for (int i = 0; i < (int)x.size(); ++i) + { + double xi = x[i]; + double gi = xi * std::sin(xi) + 0.1 * xi; + + if (absd(gi) < eps) { g[i] = 0.0; continue; } + + double sign = (gi > 0.0) ? 1.0 : -1.0; + double dgi = std::sin(xi) + xi * std::cos(xi) + 0.1; + g[i] = sign * dgi; + } + return g; +} + +UserProblem::~UserProblem() {} + diff --git a/userproblem.h b/userproblem.h new file mode 100644 index 00000000..4a9180dd --- /dev/null +++ b/userproblem.h @@ -0,0 +1,16 @@ +#ifndef USERPROBLEM_H +#define USERPROBLEM_H + +#include "problem.h" + +class UserProblem : public Problem +{ +public: + explicit UserProblem(int n); + double funmin(Data &x) override; + Data gradient(Data &x) override; + ~UserProblem(); +}; + +#endif +