//+——————————————————————————————————————————————————————————————————+
//|                                                         C_AO_HHO |
//|                                  Copyright 2007-2025, Andrey Dik |
//|                                https://www.mql5.com/ru/users/joo |
//+——————————————————————————————————————————————————————————————————+
#include "#C_AO.mqh"


//────────────────────────────────────────────────────────────────────
class C_AO_HHO : public C_AO
{
  public: //------------------------------------------------------------
  ~C_AO_HHO () { }
  C_AO_HHO ()
  {
    ao_name = "HHO";
    ao_desc = "Harris Hawks Optimization";
    ao_link = "https://www.mql5.com/ru/articles/21898";

    popSize = 50;
    beta    = 1.5;

    ArrayResize (params, 2);
    params [0].name = "popSize"; params [0].val = popSize;
    params [1].name = "beta";    params [1].val = beta;
  }

  void SetParams ()
  {
    popSize = (int)params [0].val;
    beta    =      params [1].val;
  }

  bool Init (const double &rangeMinP  [],
             const double &rangeMaxP  [],
             const double &rangeStepP [],
             const int     epochsP = 0);

  void Moving   ();
  void Revision ();

  //------------------------------------------------------------------
  double beta; // показатель хвоста Леви (β = 1.5 в оригинале)

  private: //---------------------------------------------------------
  int    t;      // текущая итерация (счёт с 1)
  int    T;      // всего итераций
  double sigma;  // σ для формулы Леви (Eq. 9), вычисляется один раз в Init

  double Xm   []; // [coords]  среднее положение стаи (Eq. 2)
  double Ytmp []; // [coords]  кандидат Y  (стратегии 3, 4)

  double LGamma   (double x); // логарифм Г(x), аппроксимация Ланцоша
  double LevyStep ();         // [L2] скалярный шаг Леви: u~N(0,σ), v~N(0,1)

  double Clamp (double v, double lo, double hi)
  {
    return v < lo ? lo : (v > hi ? hi : v);
  }
};
//————————————————————————————————————————————————————————————————————

//————————————————————————————————————————————————————————————————————
bool C_AO_HHO::Init (const double &rangeMinP  [],
                     const double &rangeMaxP  [],
                     const double &rangeStepP [],
                     const int     epochsP = 0)
{
  if (!StandardInit (rangeMinP, rangeMaxP, rangeStepP)) return false;

  //------------------------------------------------------------------
  T = (epochsP > 0) ? epochsP : 1;
  t = 0;

  ArrayResize (Xm,   coords);
  ArrayResize (Ytmp, coords);

  //--- Eq. (9): σ = [Г(1+β)·sin(πβ/2) / (Г((1+β)/2)·β·2^((β-1)/2))]^(1/β)
  double num = exp (LGamma (1.0 + beta)) * sin (M_PI * beta * 0.5);
  double den = exp (LGamma ((1.0 + beta) * 0.5)) * beta * pow (2.0, (beta - 1.0) * 0.5);
  sigma = (den > 1e-15) ? pow (fabs (num / den), 1.0 / beta) : 1.0;

  //--- начальные позиции: равномерно случайно в [LB, UB]
  for (int i = 0; i < popSize; i++)
  {
    for (int c = 0; c < coords; c++)
    {
      a [i].cP [c] = u.RNDfromCI (rangeMin [c], rangeMax [c]);
    }
    a [i].fP = -DBL_MAX; // сброс истории улучшений для [L3]
  }

  return true;
}
//————————————————————————————————————————————————————————————————————

//————————————————————————————————————————————————————————————————————
void C_AO_HHO::Moving ()
{
  for (int i = 0; i < popSize; i++)
  {
    for (int c = 0; c < coords; c++)
    {
      a [i].c [c] = u.SeInDiSp (a [i].cP [c], rangeMin [c], rangeMax [c], rangeStep [c]);
    }
  }
}
//————————————————————————————————————————————————————————————————————

//————————————————————————————————————————————————————————————————————
// Revision(): обновляет X_prey, X_m и вычисляет cP для следующего шага.
//
// Структура соответствует оригиналу:
//   |E| >= 1 → исследование:
//       q < 0.5  → X_rand formula (Eq. 1a)      [L1]
//       q >= 0.5 → X_prey − X_m  (Eq. 1b)
//   |E| <  1 → эксплуатация:
//       r>=0.5, |E|>=0.5 → стратегия 1: мягкая осада          Eq. (4-5)
//       r>=0.5, |E|< 0.5 → стратегия 2: жёсткая осада         Eq. (6)
//       r< 0.5, |E|>=0.5 → стратегия 3: мягкая осада + Lévy  Eq. (7-10)
//       r< 0.5, |E|< 0.5 → стратегия 4: жёсткая осада + Lévy Eq. (11-13)
//
// [L3] В оригинале стратегии 3/4 делают два вызова fobj (для Y и Z).
//      Прокси в нашем фреймворке: если агент улучшился в прошлой итерации
//      (a[i].f > a[i].fP) — доверяем направленному шагу Y;
//      иначе добавляем Lévy → Z (как запасной вариант).
void C_AO_HHO::Revision ()
{
  t++;

  //─── 1. Обновление личных и глобального лучших ──────────────────
  for (int i = 0; i < popSize; i++)
  {
    if (a [i].f > a [i].fB)
    {
      a [i].fB = a [i].f;
      ArrayCopy (a [i].cB, a [i].c, 0, 0, coords);
    }
    if (a [i].f > fB)
    {
      fB = a [i].f;
      ArrayCopy (cB, a [i].c, 0, 0, coords);
    }
  }

  //─── 2. Среднее положение стаи Xm (Eq. 2) ───────────────────────
  ArrayInitialize (Xm, 0.0);
  for (int i = 0; i < popSize; i++)
  {
    for (int c = 0; c < coords; c++)
    {
      Xm [c] += a [i].c [c];
    }
  }
  for (int c = 0; c < coords; c++)
  {
    Xm [c] /= (double)popSize;
  }

  //─── 3. E1 — общий для всей итерации  ─────────────
  double E1 = 2.0 * (1.0 - (double)t / (double)T);

  //─── 4. Обновление позиции каждого ястреба ──────────────────────
  for (int i = 0; i < popSize; i++)
  {
    //--- [L3] Фиксируем признак улучшения ДО обновления fP
    bool improvedLast = (a [i].f > a [i].fP);
    a [i].fP = a [i].f; // сохраняем для следующей итерации

    //--- Eq. (3): E = E1·E0, E0 ~ U(−1,1)  [E0 не нужен отдельно]
    double E    = E1 * (2.0 * u.RNDfromCI (0.0, 1.0) - 1.0);
    double absE = fabs (E);

    if (absE >= 1.0)
    {
      //══ ИССЛЕДОВАНИЕ ══════════════════════════════════════════════
      // q генерируется только здесь
      double q = u.RNDfromCI (0.0, 1.0);

      if (q < 0.5)
      {
        //─ [L1] Eq. (1a): ориентация на случайного члена стаи ──────
        //  Xi(t+1) = X_rand − r1·|X_rand − 2·r2·Xi|
        int    ri = u.RNDminusOne (popSize);
        double r1 = u.RNDfromCI (0.0, 1.0);
        double r2 = u.RNDfromCI (0.0, 1.0);
        for (int c = 0; c < coords; c++)
        {
          a [i].cP [c] = Clamp (a [ri].c [c] - r1 * fabs (a [ri].c [c] - 2.0 * r2 * a [i].c [c]), rangeMin [c], rangeMax [c]);
        }
      }
      else
      {
        //─ Eq. (1b): ориентация на разность X_prey − X_m ───────────
        //  Xi(t+1) = (X_prey − X_m) − r3·((UB−LB)·r4 + LB)
        double r3 = u.RNDfromCI (0.0, 1.0);
        double r4 = u.RNDfromCI (0.0, 1.0);
        for (int c = 0; c < coords; c++)
        {
          a [i].cP [c] = Clamp ((cB [c] - Xm [c]) - r3 * (rangeMin [c] + r4 * (rangeMax [c] - rangeMin [c])), rangeMin [c], rangeMax [c]);
        }
      }
    }
    else
    {
      //══ ЭКСПЛУАТАЦИЯ ══════════════════════════════════════════════
      // r генерируется только здесь
      double r = u.RNDfromCI (0.0, 1.0);

      if (r >= 0.5 && absE >= 0.5)
      {
        //─ Стратегия 1: мягкая осада ─────────────────────────────────
        //  ΔX = X_prey − Xi                                   Eq. (5)
        //  Xi(t+1) = ΔX − E·|J·X_prey − Xi|                  Eq. (4)
        //  J генерируется здесь (Jump_strength)
        double J = 2.0 * (1.0 - u.RNDfromCI (0.0, 1.0));
        for (int c = 0; c < coords; c++)
        {
          a [i].cP [c] = Clamp ((cB [c] - a [i].c [c]) - E * fabs (J * cB [c] - a [i].c [c]), rangeMin [c], rangeMax [c]);
        }
      }
      else
      {
        if (r >= 0.5 && absE < 0.5)
        {
          //─ Стратегия 2: жёсткая осада ────────────────────────────────
          //  Xi(t+1) = X_prey − E·|X_prey − Xi|                Eq. (6)
          //  J не используется в этой стратегии
          for (int c = 0; c < coords; c++)
          {
            a [i].cP [c] = Clamp (cB [c] - E * fabs (cB [c] - a [i].c [c]), rangeMin [c], rangeMax [c]);
          }
        }
        else
        {
          if (r < 0.5 && absE >= 0.5)
          {
            //─ Стратегия 3: мягкая осада + пикирование Леви ──────────────
            //  Y = X_prey − E·|J·X_prey − Xi|                    Eq. (7)
            //  Z = Y + rand·LF(D)                                 Eq. (8)
            //  Xi = Y если F(Y)<F(Xi), иначе Z если F(Z)<F(Xi),
            //          иначе Xi остаётся (stay-put — фреймворк не поддерживает).
            //  [L3] прокси: improvedLast → Y; иначе → Z.
            //  J генерируется здесь
            double J = 2.0 * (1.0 - u.RNDfromCI (0.0, 1.0));
            for (int c = 0; c < coords; c++)
            {
              Ytmp [c] = Clamp (cB [c] - E * fabs (J * cB [c] - a [i].c [c]), rangeMin [c], rangeMax [c]);
            }

            if (improvedLast)
            {
              for (int c = 0; c < coords; c++) a [i].cP [c] = Ytmp [c];
            }
            else
            {
              // Z = Y + rand·LF(D), rand(1,dim).*Levy(dim)  [L2]
              for (int c = 0; c < coords; c++)
              {
                a [i].cP [c] = Clamp (Ytmp [c] + u.RNDfromCI (0.0, 1.0) * LevyStep (), rangeMin [c], rangeMax [c]);
              }
            }
          }
          else   // r < 0.5 && absE < 0.5
          {
            //─ Стратегия 4: жёсткая осада + пикирование Леви ─────────────
            //  Y' = X_prey − E·|J·X_prey − X_m|                  Eq. (12)
            //  Z' = Y' + rand·LF(D)                               Eq. (13)
            //  Xi = Y' если F(Y')<F(Xi), иначе Z' если F(Z')<F(Xi),
            //          иначе Xi остаётся (stay-put — фреймворк не поддерживает).
            //  [L3] прокси: improvedLast → Y'; иначе → Z'.
            //  J генерируется здесь
            double J = 2.0 * (1.0 - u.RNDfromCI (0.0, 1.0));
            for (int c = 0; c < coords; c++)
            {
              Ytmp [c] = Clamp (cB [c] - E * fabs (J * cB [c] - Xm [c]), rangeMin [c], rangeMax [c]);
            }

            if (improvedLast)
            {
              for (int c = 0; c < coords; c++)
              {
                a [i].cP [c] = Ytmp [c];
              }
            }
            else
            {
              // Z' = Y' + rand·LF(D)  [L2]
              for (int c = 0; c < coords; c++)
              {
                a [i].cP [c] = Clamp (Ytmp [c] + u.RNDfromCI (0.0, 1.0) * LevyStep (), rangeMin [c], rangeMax [c]);
              }
            }
          }
        }
      }
    }
  } // for i
}
//————————————————————————————————————————————————————————————————————

//————————————————————————————————————————————————————————————————————
// [L2] Скалярный шаг Леви:
//   u = randn()*σ,  v = randn()
//   step = u / |v|^(1/β)
// Нормальные сэмплы генерируются методом Бокса-Мюллера.
// Levy() возвращает вектор — здесь скаляр; вызывающий код
// вызывает функцию per-coordinate, rand(1,dim).*Levy(dim).
double C_AO_HHO::LevyStep ()
{
  // Box-Muller: два независимых N(0,1)
  double r1 = u.RNDfromCI (1e-15, 1.0);
  double r2 = u.RNDfromCI (0.0,   1.0);
  double n1 = sqrt (-2.0 * log (r1)) * MathCos (2.0 * M_PI * r2);
  double n2 = sqrt (-2.0 * log (r1)) * sin (2.0 * M_PI * r2);

  double uu = n1 * sigma; // u ~ N(0, sigma^2)
  double vv = n2;         // v ~ N(0, 1)
  double av = fabs (vv);
  if (av < 1e-15) av = 1e-15;
  return fabs (uu / pow (av, 1.0 / beta));
}
//————————————————————————————————————————————————————————————————————

//————————————————————————————————————————————————————————————————————
// Логарифм гамма-функции — аппроксимация Ланцоша (g=7, n=8).
// Точность ~15 значащих цифр для Re(x) > 0.
double C_AO_HHO::LGamma (double x)
{
  static const double cL [8] =
  {
    676.5203681218851,    -1259.1392167224028,
    771.32342877765313,   -176.61502916214059,
    12.507343278686905,    -0.13857109526572012,
    9.9843695780195716e-6,  1.5056327351493116e-7
  };
  if (x < 0.5) return log (M_PI / sin (M_PI * x)) - LGamma (1.0 - x);
  x -= 1.0;
  double ag = 0.99999999999980993;
  for (int k = 0; k < 8; k++) ag += cL [k] / (x + k + 1.0);
  double tt = x + 7.5;
  return 0.5 * log (2.0 * M_PI) + (x + 0.5) * log (tt) - tt + log (ag);
}
//————————————————————————————————————————————————————————————————————