//+——————————————————————————————————————————————————————————————————+
//|                                                          C_AO_DO |
//|                                  Copyright 2007-2025, Andrey Dik |
//|                                https://www.mql5.com/ru/users/joo |
//+——————————————————————————————————————————————————————————————————+

#include "#C_AO.mqh"

//————————————————————————————————————————————————————————————————————
struct S_DO_Coord
{
    double v;
};
//————————————————————————————————————————————————————————————————————

//————————————————————————————————————————————————————————————————————
class C_AO_DO : public C_AO
{
  public:
  ~C_AO_DO () { }

  C_AO_DO ()
  {
    ao_name = "DO";
    ao_desc = "Dandelion Optimizer";
    ao_link = "https://www.mql5.com/ru/articles/20540";

    popSize = 50;

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

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

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

  void Moving   ();
  void Revision ();

  private: //—————————————————————————————————————————————————————————
  int    epochs;       // maximum iterations
  int    currentEpoch; // current iteration
  double sigma_u;      // precomputed Levy sigma for beta=1.5

  S_DO_Coord mean   [];  // mean position
  S_DO_Coord levy   [];  // levy flight steps
  S_DO_Coord center [];  // center of search range
  S_DO_Coord range  [];  // range width

  void   LevyFlight      ();
  double LognormalPDF    (double x, double mu, double sigma);
  void   BoundaryControl (int idx);
};
//————————————————————————————————————————————————————————————————————

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

  //------------------------------------------------------------------
  epochs       = epochsP;
  currentEpoch = 0;

  ArrayResize (mean,   coords);
  ArrayResize (levy,   coords);
  ArrayResize (center, coords);
  ArrayResize (range,  coords);

  // Предвычисление центра и ширины диапазона для каждой координаты
  for (int c = 0; c < coords; c++)
  {
    center [c].v = (rangeMax [c] + rangeMin [c]) * 0.5;
    range  [c].v =  rangeMax [c] - rangeMin [c];
  }

  // Предвычисление sigma_u для Levy flight с beta=1.5
  sigma_u = 0.6966;

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

//————————————————————————————————————————————————————————————————————
void C_AO_DO::Moving ()
{
  //------------------------------------------------------------------
  // Первая итерация: инициализация популяции
  if (!revision)
  {
    for (int i = 0; i < popSize; i++)
    {
      for (int c = 0; c < coords; c++)
      {
        a [i].c [c] = u.RNDfromCI (rangeMin [c], rangeMax [c]);
        a [i].c [c] = u.SeInDiSp (a [i].c [c], rangeMin [c], rangeMax [c], rangeStep [c]);
      }
    }

    revision = true;
    return;
  }

  //------------------------------------------------------------------
  currentEpoch++;
  double t = (double)currentEpoch;
  double T = (double)epochs;

  // Параметр alpha (eq. 8): alpha = rand * (t?/T? - 2t/T + 1)
  double alpha = u.RNDfromCI (0.0, 1.0) * ((t * t) / (T * T) - 2.0 * t / T + 1.0);

  // Параметры для k (eq. 11)
  double denom = T * T - 2.0 * T + 1.0;
  if (MathAbs (denom) < 1e-10) denom = 1e-10;
  double aa = 1.0 / denom;
  double bb = -2.0 * aa;
  double cc = 1.0 - aa - bb;
  double k  = 1.0 - u.RNDfromCI (0.0, 1.0) * (cc + aa * t * t + bb * t);

  //==================================================================
  // Rising stage (Фаза подъёма)
  //==================================================================
  if (u.RNDprobab () < 0.8)
  {
    // Вихревой подъём (eq. 5)
    for (int i = 0; i < popSize; i++)
    {
      double theta = (2.0 * u.RNDfromCI (0.0, 1.0) - 1.0) * M_PI;
      double row   = 1.0 / MathExp (theta);
      double vx    = row * MathCos (theta);
      double vy    = row * MathSin (theta);
      double vxvy  = vx * vy;

      for (int c = 0; c < coords; c++)
      {
        double lamb    = MathAbs (u.GaussDistribution (0, 1, -3, 3));
        double lognPDF = LognormalPDF (lamb, 0.0, 1.0);
        double NEW     = u.RNDfromCI (rangeMin [c], rangeMax [c]);

        a [i].c [c] = a [i].c [c] + alpha * vxvy * lognPDF * (NEW - a [i].c [c]);
      }

      BoundaryControl (i);
    }
  }
  else
  {
    // Линейный подъём (eq. 10) - масштабирование относительно центра диапазона
    for (int i = 0; i < popSize; i++)
    {
      for (int c = 0; c < coords; c++)
      {
        // Смещение относительно центра, масштабирование, возврат
        double offset = a [i].c [c] - center [c].v;
        a [i].c [c] = center [c].v + offset * k;
      }

      BoundaryControl (i);
    }
  }

  //==================================================================
  // Decline stage (Фаза снижения)
  //==================================================================
  // Вычисление среднего положения (eq. 14)
  for (int c = 0; c < coords; c++)
  {
    mean [c].v = 0.0;
  }

  for (int i = 0; i < popSize; i++)
  {
    for (int c = 0; c < coords; c++)
    {
      mean [c].v += a [i].c [c];
    }
  }

  for (int c = 0; c < coords; c++)
  {
    mean [c].v /= popSize;
  }

  // Обновление позиций (eq. 13)
  for (int i = 0; i < popSize; i++)
  {
    for (int c = 0; c < coords; c++)
    {
      // Ограничиваем beta типичным диапазоном нормального распределения
      double beta = u.GaussDistribution (0, 1, -3, 3);
      double betaAlpha = beta * alpha;

      double delta = -betaAlpha * (mean [c].v - betaAlpha * a [i].c [c]);

      // Ограничиваем смещение
      double maxDelta = range [c].v * 0.3;
      if (delta >  maxDelta) delta =  maxDelta;
      if (delta < -maxDelta) delta = -maxDelta;

      a [i].c [c] = a [i].c [c] + delta;
    }

    BoundaryControl (i);
  }

  //==================================================================
  // Landing stage (Фаза приземления)
  //==================================================================
  double ratio = 2.0 * t / T;
  if (ratio > 2.0) ratio = 2.0; // Ограничение ratio

  for (int i = 0; i < popSize; i++)
  {
    LevyFlight ();

    for (int c = 0; c < coords; c++)
    {
      double elite   = cB [c];
      double current = a [i].c [c];

      // eq. 15: x = Elite + levy * alpha * (Elite - x * ratio)
      double delta = levy [c].v * alpha * (elite - current * ratio);

      // Ограничиваем итоговое смещение половиной диапазона
      double maxDelta = range [c].v * 0.5;

      if (delta >  maxDelta) delta =  maxDelta;
      if (delta < -maxDelta) delta = -maxDelta;

      a [i].c [c] = elite + delta;
    }

    BoundaryControl (i);
  }
}
//————————————————————————————————————————————————————————————————————

//————————————————————————————————————————————————————————————————————
void C_AO_DO::LevyFlight ()
{
  // Levy flight с beta = 1.5
  for (int c = 0; c < coords; c++)
  {
    double uu = u.GaussDistribution (0, sigma_u, -3.0 * sigma_u, 3.0 * sigma_u);
    double vv = u.GaussDistribution (0, 1, -3, 3);

    if (MathAbs (vv) < 1e-10) vv = 1e-10;

    levy [c].v = uu / MathPow (MathAbs (vv), 0.6667);
  }
}
//————————————————————————————————————————————————————————————————————

//————————————————————————————————————————————————————————————————————
double C_AO_DO::LognormalPDF (double x, double mu, double sigma)
{
  if (x <= 0.0) return 0.0;

  double logx  = MathLog (x);
  double diff  = logx - mu;
  double coeff = 1.0 / (x * sigma * MathSqrt (2.0 * M_PI));
  double expon = MathExp (-diff * diff / (2.0 * sigma * sigma));

  return coeff * expon;
}
//————————————————————————————————————————————————————————————————————

//————————————————————————————————————————————————————————————————————
void C_AO_DO::BoundaryControl (int idx)
{
  for (int c = 0; c < coords; c++)
  {
    double val = a [idx].c [c];
    double min = rangeMin [c];
    double max = rangeMax [c];

    // Итеративное отражение от границ (максимум 10 итераций)
    int iter = 0;
    while ((val < min || val > max) && iter < 10)
    {
      if (val < min) val = min + (min - val);
      if (val > max) val = max - (val - max);
      iter++;
    }

    // Если отражение не помогло, случайная позиция
    if (val < min || val > max)
    {
      val = u.RNDfromCI (min, max);
    }

    a [idx].c [c] = u.SeInDiSp (val, min, max, rangeStep [c]);
  }
}
//————————————————————————————————————————————————————————————————————

//————————————————————————————————————————————————————————————————————
void C_AO_DO::Revision ()
{
  int    bestIdx = 0;
  double bestFit = a [0].f;

  for (int i = 1; i < popSize; i++)
  {
    if (a [i].f > bestFit)
    {
      bestFit = a [i].f;
      bestIdx = i;
    }
  }

  if (bestFit > fB)
  {
    fB = bestFit;
    ArrayCopy (cB, a [bestIdx].c, 0, 0, coords);
  }
}
//————————————————————————————————————————————————————————————————————