Ebola Optimization Search Algorithm (EOSA)

Table of Contents

1. Introduction
2. Implementation of the Algorithm
3. Test Results
4. Conclusions

Introduction

In this article, we will examine another optimization algorithm — or rather, translate its idea into an implementation based on the natural phenomenon of viral spread.

The Ebola virus is one of the deadliest pathogens known to humankind. The 2014–2016 outbreak in West Africa demonstrated just how effectively the virus can spread through a population. It was precisely this grim efficiency that inspired researchers Oyelade and Ezugwu to develop a new optimization method.

EOSA is a new bio-inspired metaheuristic optimization algorithm based on a model of Ebola virus transmission in a population and published in 2021. Imagine a village where the first person to fall ill — "patient zero" — appeared. That person becomes a source of infection, a superspreader. The virus is transmitted in two ways:

• Short-distance transmission — through close contact between people, such as a handshake, a hug, or caring for a sick person. An infected person stays in their neighborhood but actively infects their neighbors. In optimization terms, this is exploitation: an intensive search in the neighborhood of the best solution.
• Long-distance transmission — through travel and migration. An infected person gets on a bus and travels to a neighboring city, carrying the virus into new, previously unaffected areas. This exploration is a global search in unexplored areas of the solution space.

Quarantine isolates some of the infected, preventing runaway spread. In the algorithm, this preserves the diversity of the population, preventing all agents from converging at a single point. We will see what came of this below.

Implementation of the Algorithm

The SEIR-HDVQ Mathematical Model

The authors developed a complex epidemiological model with eight compartments (subpopulations):

Symbol	State	Description
S	Susceptible	Susceptible — not yet infected
E	Exposed	Exposed individuals — in the incubation period
I	Infected	Infected individuals are active spreaders
H	Hospitalized	Hospitalized
R	Recovered	Recovered
D	Dead	Dead
V	Vaccinated	Vaccinated
Q	Quarantine	In quarantine

The model is described by a system of differential equations (Eq. 6–12), which define transitions between states:

(6) ∂S/∂t = π - (β₁I + β₃D + β₄R + β₂(PE)η)S - (τS + ΓI)

(7) ∂I/∂t = (β₁I + β₃D + β₄R + β₂(PE)λ)S − (Γ + γ)I − τS

(8) ∂H/∂t = αI − (γ + ω)H

(9) ∂R/∂t = γI − ΓR

(10) ∂V/∂t = γI − (μ + θ)V

(11) ∂D/∂t = (τS + ΓI) − δD

(12) ∂Q/∂t = (πI − (γR + ΓD)) − ξQ

where π is the recruitment rate, β₁–β₄ are the contact rates, η and λ are the decay rates, τ is the natural mortality rate, Γ is the disease-induced mortality rate, α is the hospitalization rate, γ is the recovery rate, ω is the treatment rate, μ is the vaccination rate, θ is the vaccine response, δ is the burial rate, and ξ is the quarantine rate.

Let's take a look at the main equations for updating agent positions in the original article:

Eq. 1 - Position update: mIᵢᵗ⁺¹ = mIᵢᵗ + ρ × M(I); where ρ is the scaling coefficient and M(I) is the movement speed.

Eq. 2 - Exploitation (short distance): M(I) = srate × rand(0,1) + M(Ind_best).

Eq. 3 - Exploration (long distance): M(s) = lrate × rand(0,1) + M(Ind_best).

Eq. 4 - Initialization: individualᵢ = Lᵢ + rand(0,1) × (Uᵢ + Lᵢ).

Problems with the Original Formulas

The analysis of the equations revealed critical discrepancies:

• Eq. 2 and Eq. 3 do not involve any directed motion. The formula M(I) = srate × rand(0,1) + M(Ind_best) simply adds a random number to the position of the best solution. The direction vector (target - current), which is required for the agent to move toward the target, is missing. Without this component, the agent does not move toward the best solution but instead ends up in a position unrelated to its current location.
• Eq. 4 contains an obvious typographical error. The formula L + rand × (U + L) yields an incorrect range. With bounds U=10 and L=5, the result falls within the interval [5, 20], which exceeds the upper bound. The standard initialization formula is: L + rand × (U - L).
• Eq. 2 and Eq. 3 are practically identical. The only difference is the srate vs. lrate coefficient. However, the text in the article refers to movement toward the best solution (exploitation) and movement via a random transmitter agent (exploration), which is not reflected in the formulas.

The main equations for updating agent positions are now presented in a new form:

• Eq. 2 — Exploitation (short distance): M(I) = srate × rand(0,1) × (IndBest - current). A direction vector (Ind_best - current) has been added to guide the agent toward the best solution. Without this component, the formula is meaningless for optimization.
• Eq. 3 — Exploration (long distance): M(s) = lrate × rand(0,1) × (transmitter - current) + srate × rand(0,1) × (IndBest - current). A movement component toward a random transmitter agent has been added (transmitter - current), which corresponds to the textual description of long-distance transmission via contact with other individuals.
• Eq. 4 — Initialization: individualᵢ = Lᵢ + rand(0,1) × (Uᵢ - Lᵢ). An obvious typo has been corrected: the plus sign has been replaced with a minus sign to ensure correct generation within the specified range [L, U].

Pseudocode for the algorithm:

1. Initialize S, E, I, H, R, V, Q ← ∅
2. Create the initial population S according to Eq. 4
3. Select the index case (patient zero) → I
4.
5. WHILE epoch ≤ max_epoch AND |I| > 0:
6. Q ← quarantine a portion of I according to Eq. 12
7. fracI ← I \ Q (active infected)
8.
9. FOR each agent in fracI:
10. IF the incubation period has elapsed:
11. Compute neighborhood
12. IF neighborhood < 0.5:
13. Exploitation according to Eq. 2
14. ELSE:
15. Exploration according to Eq. 3
16.
17. // State transitions:
18. H ← hospitalize a portion of I according to Eq. 8
19. R ← agents recovered from I according to Eq. 9
20. V ← vaccinate a portion of H according to Eq. 10
21. D ← deaths from I according to Eq. 11
22.
23. I ← I - R - D
24. S ← S + R // those who have recovered return
25. S ← S - D + new // the deceased are replaced by new agents
26.
27. Update the best solution

Problem: a gap between the description and optimization. When attempting a literal implementation of the algorithm, several critical problems emerged, we ran into serious problems:

• Incorrect motion equations. The original Eq. 2 and Eq. 3 do not involve directed motion toward a target. Without a vector (target - current), agents cannot systematically improve their positions.
• The population dies out. At the start, only one "index case" agent is infected (active). With a high mortality rate (δ = 0.5) and a low infection rate, the population rapidly declines to zero.
• Agent freezing. Hospitalized (H), vaccinated (V), and deceased (D) agents do not participate in the search. They are "switched off" from optimization, which drastically reduces efficiency.
• Weak generation of new solutions. The infection formula numNew = β₁ × |I| × srate yields zero new agents when |I| is small.
• Uncertainty around "neighborhood." The article does not provide an explicit formula for calculating the "neighborhood" parameter, on which the choice of exploitation or exploration strategy depends.

Final Implementation

The philosophy of adaptation: I had to preserve the spirit of the epidemiological metaphor, but rethink it for optimization. Instead of literally modeling the spread of a disease, we use key mechanisms as a source of inspiration. Imagine not the epidemic itself, but information about a good solution that spreads through the population like a virus:

• Agents that are close to the best solution "catch" its quality — they refine their positions locally;
• Distant agents undertake "journeys" — Lévy flights — in search of new, promising areas.

Let's introduce a simplified structure. Instead of eight states, all agents are active:

Parameter	Value	Description
popSize	50	Population size
srate	1.5	Exploitation intensity
lrate	1.0	Exploration intensity
quarantine	0.05	Probability that an agent skips the current iteration (quarantine rate)

Strategy Selection Mechanism

The "neighborhood" parameter is calculated randomly, taking the agent's fitness into account: neighborhood = rand() × (1 - fitness_ratio × 0.5).

Better agents enter exploitation mode more often, while poorer-performing agents explore more often. This makes intuitive sense: if you have struck gold, dig deeper; if your mine is empty, look for a new spot.

Exploitation: local search. When (neighborhood < 0.5), the agent is in "close contact" mode:

// Target point: a mixture of the global best and personal best
target = w × gBest + (1-w) × pBest

// Movement toward the target with local noise
M = srate × rand() × (target - current)
M += Gaussian_noise × range × 0.05

new_position = current + ρ × M

Let's draw an analogy: an infected person in their neighborhood. He visits his neighbors (gBest, pBest), but doesn't leave the neighborhood. A little Gaussian noise is like chance encounters on the street.

Exploration: global search. When (neighborhood ≥ 0.5), the agent "travels". Two options. Option A — Lévy flights (50% of cases):

levy = LevyFlight()  // Lévy distribution
step = lrate × levy × range × 0.1
attraction = 0.1 × rand() × (gBest - current)

new_position = current + step + attraction

Example: an infected person boards a plane and flies to a random city. The Lévy distribution is many "short flights plus rare transcontinental flights." This is exactly how real epidemics spread in the age of globalization.

Option B — transmission via an agent (50% of cases):

j = random_agent ≠ i

M = lrate × rand() × (agent[j] - current)
M += 0.3 × srate × rand() × (gBest - current)

new_position = current + ρ × M

Example: an infected person meets a traveler from another city and learns about the situation there.

Lévy flights are the mathematics of long-distance travel. The Lévy distribution is a heavy-tailed distribution characteristic of animal movements in search of food, human migration, and the spread of epidemics.

It is characterized by many small steps, interspersed with occasional giant leaps.

double LevyFlight()
{
    // Mantegna's algorithm, β = 1.5
    double sigma_u = 0.6966;  // precomputed constant
    
    double u = Gaussian() × sigma_u;
    double v = Gaussian();
    
    double step = u / |v|^0.6667;
    
    return clamp(step, -3, 3);
}

Quarantine: preserving diversity. With a 5% probability, the agent skips an iteration:

if (rand() < quarantine)
    continue;  // agent in isolation

Why? Without quarantine, all agents rush toward the best solution and may merge into a single point. Quarantine is a form of “inertia”: some of the population retain their positions, ensuring diversity for future iterations.

Personal memory (pBest). Each agent remembers its own best position:

if (current.fitness > pBest.fitness)
    pBest = current;

In exploitation mode, the agent moves toward a combination of gBest and pBest. This prevents a situation where all agents blindly follow a single leader. An analogy: you remember where you personally felt good, even if everyone says things are better somewhere else right now.

Why were these changes made? Occam and metaheuristics. The complexity of the epidemiological model does not translate into optimization quality. Eight states and sixteen parameters create a large number of “degrees of freedom” and chaotic dynamics, in which most agents are frozen in inactive states. This situation can be described using Occam's razor: of two explanations, choose the simplest one. Of the two algorithms, choose the one with fewer parameters that need tuning or meta-optimization.

Lévy flights are a scientifically sound alternative. The original study refers to "long-distance movements" but uses a simple linear combination. Lévy flights are a mathematically grounded model of long-distance movements.

Fitness bias — adaptive balance. Instead of a fixed threshold (neighborhood = 0.5), we use an adaptive mechanism in which good solutions are exploited, while poor solutions continue to be explored. This is only natural: a successful businessman expands on what works, while a beginner tries different things.
EOSA is an example of how a biological metaphor can inspire the creation of an optimization algorithm, but its literal implementation is not always optimal and can be interpreted in several ways. In the end, the key ideas were retained:

• two search modes (exploitation/exploration), analogous to short- and long-distance virus transmission;
• quarantine — as a mechanism for preserving diversity;
• the scale factor "ρ" — as the attenuation of the epidemic over time.

Proven techniques were added: "Lévy flights" for global exploration, personal memory pBest for stability, and "Fitness bias" for adaptive strategy balancing. The result is a compact algorithm with four parameters that preserves the epidemiological intuition of the original. Here are the formulas for the final implementation:

Initialization (Eq. 4): aᵢ = rangeMin + rand(0,1) × (rangeMax - rangeMin)

Scale factor: ρ = 1.0 - epoch/max_epoch × 0.5 (from 1.0 to 0.5)

Exploitation (neighborhood < 0.5):

target = w × gBest + (1-w) × pBest, w ∈ [0.3, 0.7]

M = srate × rand() × (target - aᵢ) + Gaussian() × range × 0.05 × ρ

aᵢ = aᵢ + ρ × M

Exploration — Lévy flights (50%):

step = lrate × LevyFlight() × range × 0.1

attraction = 0.1 × rand() × (gBest - aᵢ)

aᵢ = aᵢ + step + attraction × ρ

Exploration — via an agent (50%):

M = lrate × rand() × (aⱼ - aᵢ) + 0.3 × srate × rand() × (gBest - aᵢ)

aᵢ = aᵢ + ρ × M

Lévy flight (β = 1.5):

σᵤ = 0.6966

u = Gaussian() × σᵤ

v = Gaussian()

step = u / |v|^0.6667

return clamp(step, -3, 3)

Let's take a look at the diagram illustrating one implementation of the EOSA algorithm.

Figure 1. Schematic illustration of how the EOSA algorithm works

The top part of the illustration shows two operating modes:

EXPLOITATION (left): the agent moves toward the target (a combination of gBest and pBest) + Gaussian noise for local search; EXPLORATION (right): Lévy flights or transmission via a random agent j (50/50).

Let us prepare detailed pseudocode for the subsequent implementation of the algorithm in code.

Algorithm EOSA
Input: popSize, srate, lrate, ξ, epochs, rangeMin[], rangeMax[]
Output: gBest — best solution found

BEGIN
// Initialization
FOR each agent i DO
a[i] ← random position in [rangeMin, rangeMax]
pBest[i] ← a[i]
END FOR
gBest ← best of population

// Main loop
FOR epoch = 1 TO epochs DO
ρ ← 1.0 − epoch/epochs × 0.5 // scale factor: 1.0 → 0.5

FOR each agent i DO
IF rand() < ξ THEN CONTINUE // quarantine

neighborhood ← rand() × (1 − fitnessBias)

IF neighborhood < 0.5 THEN
//——————————————————————————————————————
// EXPLOITATION: local search
//——————————————————————————————————————
w ← rand(0.3, 0.7)
target ← w·gBest + (1−w)·pBest[i]
M ← srate·rand()·(target − a[i]) + GaussianNoise
a[i] ← a[i] + ρ·M

ELSE
//——————————————————————————————————————
// EXPLORATION: global search
//——————————————————————————————————————
IF rand() < 0.5 THEN
// Lévy flight
step ← lrate · LévyFlight() · range · 0.1
a[i] ← a[i] + step + 0.1·rand()·(gBest − a[i])
ELSE
// Transmission via a random agent j
j ← random agent, j ≠ i
M ← lrate·rand()·(a[j] − a[i]) + 0.3·srate·rand()·(gBest − a[i])
a[i] ← a[i] + ρ·M
END IF
END IF

a[i] ← clamp(a[i], rangeMin, rangeMax) // boundary control
END FOR

// Update best solutions
FOR each agent i DO
IF f(a[i]) > f(pBest[i]) THEN pBest[i] ← a[i]
IF f(a[i]) > f(gBest) THEN gBest ← a[i]
END FOR
END FOR

RETURN gBest
END
// Lévy Flight (Mantegna’s algorithm, β = 1.5)

Function LévyFlight()
u ← Gaussian() × 0.6966
v ← Gaussian()
RETURN clamp(u / |v|^0.6667, −3, 3)
END
// Default parameters: popSize=50, srate=1.5, lrate=1.0, ξ=0.05

Let us move on to describing the implementation of the EOSA algorithm.

The C_AO_EOSA class implements the EOSA (Ebola Optimization Search Algorithm) optimization algorithm, which is inspired by the mechanisms of Ebola virus transmission. The class inherits from the base class C_AO and contains four configurable parameters: popSize — population size (default: 50), srate — exploitation intensity (1.5), lrate — exploration intensity using Lévy flights (1.0), and quarantine — the probability of an agent being quarantined (0.05). The internal variables "epochs" and "currentEpoch" track the total number of iterations and the current epoch, respectively.

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

  C_AO_EOSA ()
  {
    ao_name = "EOSA";
    ao_desc = "Ebola Optimization Search Algorithm";
    ao_link = "https://www.mql5.com/en/articles/20932";

    popSize    = 50;
    srate      = 1.5;
    lrate      = 1.0;
    quarantine = 0.05;

    ArrayResize (params, 4);
    params [0].name = "popSize";    params [0].val = popSize;
    params [1].name = "srate";      params [1].val = srate;      // exploitation intensity
    params [2].name = "lrate";      params [2].val = lrate;      // exploration intensity (Lévy)
    params [3].name = "quarantine"; params [3].val = quarantine; // quarantine rate
  }

  void SetParams ()
  {
    popSize    = (int)params [0].val;
    srate      = params      [1].val;
    lrate      = params      [2].val;
    quarantine = params      [3].val;
  }

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

  void Moving   ();
  void Revision ();

  private: //—————————————————————————————————————————————————————————
  double srate;      // exploitation intensity
  double lrate;      // exploration intensity (Lévy)
  double quarantine; // quarantine probability

  int    epochs;
  int    currentEpoch;

  void   BoundaryControl (int idx);
  double LevyFlight      ();
  double RandGauss       ();
};
//————————————————————————————————————————————————————————————————————

The Init method initializes the algorithm. It calls the standard initialization of the base class, passing the search boundary arrays rangeMin and rangeMax and the discretization rangeStep. It then saves the total number of epochs and resets the current epoch counter to zero. The method returns 'true' if initialization is successful.

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

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

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

The Moving method is the core of the algorithm and implements the agent movement logic. On the first iteration, when the "revision" flag is 'false', the method initializes the population with random positions within the specified bounds and then returns. In subsequent iterations, the adaptive scale factor "ρ" is computed; it decreases linearly from 1.0 to 0.5 as the optimization progresses, ensuring a gradual transition from intensive exploration to fine-tuning the solutions. For each agent, the quarantine condition is checked first: with probability "quarantine," the agent skips the current iteration, which prevents premature convergence and preserves population diversity.

Next, the "neighborhood" parameter is calculated, which determines the agent's behavior strategy. This parameter is a random value adjusted to account for the quality of the agent's current solution: the best agents are more likely to enter exploitation mode, while the worst agents are more likely to enter exploration mode. If the "neighborhood" value is less than 0.5, the agent performs exploitation — a local search around the best known solutions. The target point is formed as a weighted combination of the global best solution gBest and the agent's personal best solution pBest, with a weight "w" randomly selected from the interval [0.3, 0.7]. Movement toward the target is performed with intensity "srate", with Gaussian noise added for local exploration of the neighborhood.

If the "neighborhood" value is at least 0.5, the agent performs exploration — a global search using one of two strategies selected with equal probability. The first strategy is a Lévy flight, which simulates the long-distance movements characteristic of epidemic spread through travel and migration. The Lévy step is scaled by the coefficient "lrate", and a weak attraction toward the global best solution is added. The second strategy is transmission via a random agent, in which the current agent moves toward another randomly selected agent in the population, with additional attraction toward the best solution.

//————————————————————————————————————————————————————————————————————
void C_AO_EOSA::Moving ()
{
  //------------------------------------------------------------------
  // FIRST ITERATION: Population initialization
  //------------------------------------------------------------------
  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++;

  // Adaptive scale: decreases over time
  double progress = (double)currentEpoch / (double)MathMax (epochs, 1);
  double rho = 1.0 - progress * 0.5;  // 1.0 → 0.5

  //------------------------------------------------------------------
  // Updating the positions of all agents
  //------------------------------------------------------------------
  for (int i = 0; i < popSize; i++)
  {
    //================================================================
    // QUARANTINE: an agent in isolation skips the iteration
    // Epidemiological interpretation: isolation prevents spread
    //================================================================
    if (u.RNDprobab () < quarantine)
    {
      continue;
    }

    //================================================================
    // Determine neighborhood (proximity to the "epicenter"—the best solution)
    // Use random selection with a fitness-based bias
    //================================================================
    double neighborhood = u.RNDprobab ();

    // Bias: the best agents perform exploitation more often
    double fitnessRatio = 0.5;
    if (fB > -DBL_MAX + 1e10 && a [i].f > -DBL_MAX + 1e10)
    {
      // Normalize fitness for the bias
      double range_f = MathMax (fB - a [i].f, 1e-10);
      fitnessRatio = MathExp (-range_f / MathMax (MathAbs (fB), 1e-10));
    }

    // Combine randomness with fitness bias
    neighborhood = neighborhood * (1.0 - fitnessRatio * 0.5);

    //================================================================
    // UPDATE STRATEGY SELECTION
    //================================================================
    if (neighborhood < 0.5)
    {
      //==============================================================
      // EXPLOITATION: Short-distance transmission (Eq. 2)
      // Epidemiology: close contact → local spread
      // Optimization: intensive search around the best solutions
      //==============================================================

      // Choose a target point: a blend of the global best and personal best
      double w = u.RNDfromCI (0.3, 0.7);

      for (int c = 0; c < coords; c++)
      {
        double target;

        // Target point: a weighted combination of gBest and pBest
        if (a [i].fB > -DBL_MAX + 1e10)
        {
          target = w * cB [c] + (1.0 - w) * a [i].cB [c];
        }
        else
        {
          target = cB [c];
        }

        // Interpretation of Eq. 2: M(I) = srate * rand * direction + noise
        double r1 = u.RNDfromCI (0.0, 1.0);
        double direction = target - a [i].c [c];

        // Movement toward the target + local randomness
        double M = srate * r1 * direction;

        // Add Gaussian noise for local exploration
        double range = rangeMax [c] - rangeMin [c];
        M += RandGauss () * range * 0.05 * rho;

        // Eq. 1: position update
        a [i].c [c] = a [i].c [c] + rho * M;
      }
    }
    else
    {
      //==============================================================
      // EXPLORATION: Long-distance transmission (Eq. 3)
      // Epidemiology: distant contact → global spread
      // Optimization: Lévy flights for global search
      //==============================================================

      // Select a random agent j ≠ i (transmitter)
      int j = u.RNDminusOne (popSize);
      if (j == i) j = (i + 1) % popSize;

      // Decide: Lévy flight or movement via the transmitter
      if (u.RNDprobab () < 0.5)
      {
        //------------------------------------------------------------
        // LÉVY FLIGHT: simulation of long-distance movements
        // In epidemiology: travel, migration
        //------------------------------------------------------------
        for (int c = 0; c < coords; c++)
        {
          double range = rangeMax [c] - rangeMin [c];
          double levy = LevyFlight ();  // β = 1.5

          // Lévy step from the current position
          double step = lrate * levy * range * 0.1;

          // Weak attraction toward the best solution
          double attraction = 0.1 * u.RNDfromCI (0.0, 1.0) * (cB [c] - a [i].c [c]);

          a [i].c [c] = a [i].c [c] + step + attraction * rho;
        }
      }
      else
      {
        //------------------------------------------------------------
        // TRANSMISSION VIA AN AGENT: movement toward another agent
        // In epidemiology: transmission through contact
        //------------------------------------------------------------
        for (int c = 0; c < coords; c++)
        {
          double r1 = u.RNDfromCI (0.0, 1.0);
          double r2 = u.RNDfromCI (0.0, 1.0);

          // Eq. 3: M(s) = lrate * rand * (a[j] - a[i]) + attraction_to_best
          double M = lrate * r1 * (a [j].c [c] - a [i].c [c]);

          // A weak attraction toward the best solution (not dominant)
          M += 0.3 * srate * r2 * (cB [c] - a [i].c [c]);

          a [i].c [c] = a [i].c [c] + rho * M;
        }
      }
    }

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

The LevyFlight method generates a random step according to the Lévy distribution using Mantegna’s algorithm. A fixed value of β = 1.5 is used, for which the constant σ_u = 0.6966 has been precomputed. The method generates two Gaussian random numbers and computes the step using the formula u/|v|^0.6667, clipping the result to the interval [-3, 3] to prevent extreme outliers.

//————————————————————————————————————————————————————————————————————
double C_AO_EOSA::LevyFlight ()
{
  // Mantegna’s algorithm for the Lévy distribution
  // beta = 1.5 (fixed)

  // Precomputed sigma_u for beta = 1.5:
  // sigma_u = [Γ(2.5) × sin(3π/4) / (Γ(1.25) × 1.5 × 2^0.25)]^(2/3)
  //         = [1.3293 × 0.7071 / (0.9064 × 1.5 × 1.1892)]^0.6667
  //         ≈ 0.6966
  double sigma_u = 0.6966;

  double u_val = RandGauss () * sigma_u;
  double v_val = RandGauss ();

  // Division-by-zero protection
  if (MathAbs (v_val) < 1e-10) v_val = 1e-10;

  // step = u / |v|^(1/beta) = u / |v|^0.6667 for beta = 1.5
  double step = u_val / MathPow (MathAbs (v_val), 0.6667);

  // Limit extreme values
  if (step > 3.0)  step = 3.0;
  if (step < -3.0) step = -3.0;

  return step;
}
//————————————————————————————————————————————————————————————————————

The RandGauss method generates a Gaussian random variable with a mean of zero and a variance of one using the Box-Muller method. Two uniformly distributed numbers are transformed into a normally distributed value through logarithmic and trigonometric transformations.

//————————————————————————————————————————————————————————————————————
double C_AO_EOSA::RandGauss ()
{
  // Box-Muller transform
  double u1 = u.RNDfromCI (1e-10, 1.0);
  double u2 = u.RNDfromCI (0.0, 1.0);

  return MathSqrt (-2.0 * MathLog (u1)) * MathCos (2.0 * M_PI * u2);
}
//————————————————————————————————————————————————————————————————————

The BoundaryControl method controls the search boundaries using the reflection method. Unlike simple clipping, reflection preserves the agent's momentum: when the agent goes beyond the lower boundary, its position is reflected upward; when it goes beyond the upper boundary, it is reflected downward. If, after several reflections, the position remains outside the valid range, the coordinate is randomly reinitialized.

//————————————————————————————————————————————————————————————————————
void C_AO_EOSA::BoundaryControl (int idx)
{
  for (int c = 0; c < coords; c++)
  {
    // Reflection at boundaries instead of clipping
    double val = a [idx].c [c];
    double min = rangeMin [c];
    double max = rangeMax [c];
    double range = max - min;

    // Reflection
    while (val < min || val > max)
    {
      if (val < min) val = min + (min - val);
      if (val > max) val = max - (val - max);

      // Infinite Loop Protection
      if (val < min - range || val > max + range)
      {
        val = min + u.RNDfromCI (0.0, 1.0) * range;
        break;
      }
    }

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

The Revision method updates the best solutions found. For each agent, the system checks whether its current fitness exceeds its personal best, and updates pBest if necessary. At the same time, the best agent in the population is tracked, and if its fitness exceeds the global best, gBest is updated. This approach ensures that information about promising areas of the search space is retained.

//————————————————————————————————————————————————————————————————————
void C_AO_EOSA::Revision ()
{
  //------------------------------------------------------------------
  // Update Personal and Global Bests
  //------------------------------------------------------------------
  int    bestIdx = 0;
  double bestFit = -DBL_MAX;

  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 > bestFit)
    {
      bestFit = a [i].f;
      bestIdx = i;
    }
  }

  //------------------------------------------------------------------
  // Update the global best (Eq. 5)
  //------------------------------------------------------------------
  if (bestFit > fB)
  {
    fB = bestFit;
    ArrayCopy (cB, a [bestIdx].c, 0, 0, coords);
  }
}
//————————————————————————————————————————————————————————————————————

Test Results

Now we can look at the results: they are not outstanding, but they show that the algorithm is capable of tackling the tasks. The EOSA algorithm was tested on a standard set of test functions: Hilly (smooth multimodal), Forest (sharp peaks), and Megacity (discrete stepwise). Testing was conducted for problem dimensions of 5, 25, and 500, with a limit of 10,000 objective function evaluations.

The overall result was about 38.5% of the maximum possible score, which is approximately 10% below the threshold for inclusion in the table of the best population-based optimization methods.

EOSA|Ebola Optimization Search Algorithm|50.0|3.0|2.0|0.01|
=============================
5 Hilly; Func runs: 10000; result: 0.710222611733512
25 Hilly; Func runs: 10000; result: 0.4632808052522758
500 Hilly; Func runs: 10000; result: 0.29182585953654855
=============================
5 Forest; Func runs: 10000; result: 0.6237714734740634
25 Forest; Func runs: 10000; result: 0.38209042951228955
500 Forest; Func runs: 10000; result: 0.20125437332932763
=============================
5 Megacity; Func runs: 10000; result: 0.4676923076923078
25 Megacity; Func runs: 10000; result: 0.2132307692307692
500 Megacity; Func runs: 10000; result: 0.11516923076923188
=============================
Overall score: 3.46854 (38.54%)

The visualization shows that the algorithm has more difficulty handling discrete functions such as Megacity within a finite number of iterations. A stepwise landscape with flat plateaus does not provide gradient information, and the mechanism for moving toward the best solution becomes ineffective when many points have the same fitness.

EOSA on the Hilly test function

EOSA on the Forest test function

EOSA on the Megacity test function

You can also view a visualization of the algorithm's operation on simpler standard functions; the parabolic function also has some issues at high problem dimensionality, and convergence is not optimal.

EOSA on the test function GoldsteinPrice

EOSA on the test function Paraboloid

In the ranking table of the best population-based optimization methods, the EOSA algorithm is included for reference.

No.	AO	Description	Hilly			Hilly Final	Forest			Forest Final	Megacity (discrete)			Megacity Final	Final Result	% of MAX
			10 p (5 F)	50 p (25 F)	1,000 p (500 F)		10 p (5 F)	50 p (25 F)	1,000 p (500 F)		10 p (5 F)	50 p (25 F)	1,000 p (500 F)
1	DOAdingom	dingo_optimization_algorithm_M	0.47968	0.45367	0.46369	1.39704	0.94145	0.87909	0.91454	2.73508	0.78615	0.86061	0.84805	2.49481	6.627	73.63
2	ANS	across neighbourhood search	0.94948	0.84776	0.43857	2.23581	1.00000	0.92334	0.39988	2.32323	0.70923	0.63477	0.23091	1.57491	6.134	68.15
3	CLA	code lock algorithm (joo)	0.95345	0.87107	0.37590	2.20042	0.98942	0.91709	0.31642	2.22294	0.79692	0.69385	0.19303	1.68380	6.107	67.86
4	AMOm	animal migration optimization M	0.90358	0.84317	0.46284	2.20959	0.99001	0.92436	0.46598	2.38034	0.56769	0.59132	0.23773	1.39675	5.987	66.52
5	(P+O)ES	(P+O) evolution strategies	0.92256	0.88101	0.40021	2.20379	0.97750	0.87490	0.31945	2.17185	0.67385	0.62985	0.18634	1.49003	5.866	65.17
6	CTA	comet tail algorithm (joo)	0.95346	0.86319	0.27770	2.09435	0.99794	0.85740	0.33949	2.19484	0.88769	0.56431	0.10512	1.55712	5.846	64.96
7	TETA	time evolution travel algorithm (joo)	0.91362	0.82349	0.31990	2.05701	0.97096	0.89532	0.29324	2.15952	0.73462	0.68569	0.16021	1.58052	5.797	64.41
8	SDSm	stochastic diffusion search M	0.93066	0.85445	0.39476	2.17988	0.99983	0.89244	0.19619	2.08846	0.72333	0.61100	0.10670	1.44103	5.709	63.44
9	BOAm	billiards optimization algorithm M	0.95757	0.82599	0.25235	2.03590	1.00000	0.90036	0.30502	2.20538	0.73538	0.52523	0.09563	1.35625	5.598	62.19
10	AAm	archery algorithm M	0.91744	0.70876	0.42160	2.04780	0.92527	0.75802	0.35328	2.03657	0.67385	0.55200	0.23738	1.46323	5.548	61.64
11	ESG	evolution of social groups (joo)	0.99906	0.79654	0.35056	2.14616	1.00000	0.82863	0.13102	1.95965	0.82333	0.55300	0.04725	1.42358	5.529	61.44
12	SIA	simulated isotropic annealing (joo)	0.95784	0.84264	0.41465	2.21513	0.98239	0.79586	0.20507	1.98332	0.68667	0.49300	0.09053	1.27020	5.469	60.76
13	EOm	extremal_optimization_M	0.76166	0.77242	0.31747	1.85155	0.99999	0.76751	0.23527	2.00277	0.74769	0.53969	0.14249	1.42987	5.284	58.71
14	BBO	biogeography-based optimization	0.94912	0.69456	0.35031	1.99399	0.93820	0.67365	0.25682	1.86867	0.74615	0.48277	0.17369	1.40261	5.265	58.50
15	ACS	artificial cooperative search	0.75547	0.74744	0.30407	1.80698	1.00000	0.88861	0.22413	2.11274	0.69077	0.48185	0.13322	1.30583	5.226	58.06
16	DA	dialectical algorithm	0.86183	0.70033	0.33724	1.89940	0.98163	0.72772	0.28718	1.99653	0.70308	0.45292	0.16367	1.31967	5.216	57.95
17	BHAm	black hole algorithm M	0.75236	0.76675	0.34583	1.86493	0.93593	0.80152	0.27177	2.00923	0.65077	0.51646	0.15472	1.32195	5.196	57.73
18	ASO	anarchy society optimization	0.84872	0.74646	0.31465	1.90983	0.96148	0.79150	0.23803	1.99101	0.57077	0.54062	0.16614	1.27752	5.178	57.54
19	RFO	royal flush optimization (JOO)	0.83361	0.73742	0.34629	1.91733	0.89424	0.73824	0.24098	1.87346	0.63154	0.50292	0.16421	1.29867	5.089	56.55
20	AOSm	atomic orbital search M	0.80232	0.70449	0.31021	1.81702	0.85660	0.69451	0.21996	1.77107	0.74615	0.52862	0.14358	1.41835	5.006	55.63
21	TSEA	turtle shell evolution algorithm (JOO)	0.96798	0.64480	0.29672	1.90949	0.99449	0.61981	0.22708	1.84139	0.69077	0.42646	0.13598	1.25322	5.004	55.60
22	BSA	backtracking_search_algorithm	0.97309	0.54534	0.29098	1.80941	0.99999	0.58543	0.21747	1.80289	0.84769	0.36953	0.12978	1.34700	4.959	55.10
23	DE	differential evolution	0.95044	0.61674	0.30308	1.87026	0.95317	0.78896	0.16652	1.90865	0.78667	0.36033	0.02953	1.17653	4.955	55.06
24	SRA	successful restaurateur algorithm (JOO)	0.96883	0.63455	0.29217	1.89555	0.94637	0.55506	0.19124	1.69267	0.74923	0.44031	0.12526	1.31480	4.903	54.48
25	BO	bonobo_optimizer	0.77565	0.63805	0.32908	1.74278	0.88088	0.76344	0.25573	1.90005	0.61077	0.49846	0.14246	1.25169	4.895	54.38
26	CRO	chemical reaction optimization	0.94629	0.66112	0.29853	1.90593	0.87906	0.58422	0.21146	1.67473	0.75846	0.42646	0.12686	1.31178	4.892	54.36
27	BIO	blood inheritance optimization (JOO)	0.81568	0.65336	0.30877	1.77781	0.89937	0.65319	0.21760	1.77016	0.67846	0.47631	0.13902	1.29378	4.842	53.80
28	DOA	dream_optimization_algorithm	0.85556	0.70085	0.37280	1.92921	0.73421	0.48905	0.24147	1.46473	0.77231	0.47354	0.18561	1.43146	4.825	53.62
29	BSA	bird swarm algorithm	0.89306	0.64900	0.26250	1.80455	0.92420	0.71121	0.24939	1.88479	0.69385	0.32615	0.10012	1.12012	4.809	53.44
30	DEA	dolphin_echolocation_algorithm	0.75995	0.67572	0.34171	1.77738	0.89582	0.64223	0.23941	1.77746	0.61538	0.44031	0.15115	1.20684	4.762	52.91
31	HS	harmony search	0.86509	0.68782	0.32527	1.87818	0.99999	0.68002	0.09590	1.77592	0.62000	0.42267	0.05458	1.09725	4.751	52.79
32	SSG	saplings sowing and growing	0.77839	0.64925	0.39543	1.82308	0.85973	0.62467	0.17429	1.65869	0.64667	0.44133	0.10598	1.19398	4.676	51.95
33	BCOm	bacterial chemotaxis optimization M	0.75953	0.62268	0.31483	1.69704	0.89378	0.61339	0.22542	1.73259	0.65385	0.42092	0.14435	1.21912	4.649	51.65
34	ABO	african buffalo optimization	0.83337	0.62247	0.29964	1.75548	0.92170	0.58618	0.19723	1.70511	0.61000	0.43154	0.13225	1.17378	4.634	51.49
35	(PO)ES	(PO) evolution strategies	0.79025	0.62647	0.42935	1.84606	0.87616	0.60943	0.19591	1.68151	0.59000	0.37933	0.11322	1.08255	4.610	51.22
36	FBA	fractal-based algorithm	0.79000	0.65134	0.28965	1.73099	0.87158	0.56823	0.18877	1.62858	0.61077	0.46062	0.12398	1.19537	4.555	50.61
37	TSm	tabu search M	0.87795	0.61431	0.29104	1.78330	0.92885	0.51844	0.19054	1.63783	0.61077	0.38215	0.12157	1.11449	4.536	50.40
38	BSO	brain storm optimization	0.93736	0.57616	0.29688	1.81041	0.93131	0.55866	0.23537	1.72534	0.55231	0.29077	0.11914	0.96222	4.498	49.98
39	WOAm	whale optimization algorithm M	0.84521	0.56298	0.26263	1.67081	0.93100	0.52278	0.16365	1.61743	0.66308	0.41138	0.11357	1.18803	4.476	49.74
40	AEFA	artificial electric field algorithm	0.87700	0.61753	0.25235	1.74688	0.92729	0.72698	0.18064	1.83490	0.66615	0.11631	0.09508	0.87754	4.459	49.55
41	AEO	artificial ecosystem-based optimization algorithm	0.91380	0.46713	0.26470	1.64563	0.90223	0.43705	0.21400	1.55327	0.66154	0.30800	0.28563	1.25517	4.454	49.49
42	CAm	camel algorithm M	0.78684	0.56042	0.35133	1.69859	0.82772	0.56041	0.24336	1.63149	0.64846	0.33092	0.13418	1.11356	4.444	49.37
43	ACOm	ant colony optimization M	0.88190	0.66127	0.30377	1.84693	0.85873	0.58680	0.15051	1.59604	0.59667	0.37333	0.02472	0.99472	4.438	49.31
44	CMAES	covariance_matrix_adaptation_evolution_strategy	0.76258	0.72089	0.00000	1.48347	0.82056	0.79616	0.00000	1.61672	0.75846	0.49077	0.00000	1.24923	4.349	48.33
45	DA_duelist	duelist_algorithm	0.92782	0.53778	0.27792	1.74352	0.86957	0.47536	0.18193	1.52686	0.62153	0.33569	0.11715	1.07437	4.345	48.28
	EOSA	ebola_optimization_search_algorithm	0.71022	0.46328	0.29182	1.46532	0.62377	0.38209	0.20125	1.20711	0.46769	0.21323	0.11516	0.79608	3.469	38.54
	RW	random walk	0.48754	0.32159	0.25781	1.06694	0.37554	0.21944	0.15877	0.75375	0.27969	0.14917	0.09847	0.52734	2.348	26.09

Conclusions

The original article proposes a complex epidemiological model with eight agent states and sixteen parameters. However, the mathematical framework of differential equations used to describe the spread of disease cannot be directly translated into effective search operators. An attempt at literal implementation leads to the extinction of the active population and the freezing of most agents in inactive states. My simplified implementation retained the key metaphors — short- and long-distance transmission and quarantine — but lost some of the potential variety of mechanisms. Perhaps the potential for a more effective search lay precisely in the complex transitions between states, but that potential could not be realized.

The very concept of the “spread of infection” as a metaphor for optimization has a fundamental limitation: in a real epidemic, the “goal” is to spread the virus as widely as possible, whereas in optimization the goal is to focus on the best solution. These two tasks are fundamentally opposite in nature. Despite its modest practical results, the EOSA algorithm is of significant value as educational material. It clearly illustrates the process of transforming a biological metaphor into a computational algorithm, along with all the associated challenges and trade-offs.

The example of EOSA clearly illustrates a typical problem with “metaphorical” algorithms: an appealing analogy does not guarantee efficiency. The complexity of the original model may end up being a hindrance rather than an advantage. The process of simplification and adaptation is an integral part of metaheuristic development.

EOSA is a viable optimization algorithm. Its performance is insufficient to qualify it as one of the best methods, especially for discrete and high-dimensional problems. The epidemiological metaphor on which the approach was based turned out to be less than ideal for optimization purposes, and the complex mathematical model in the original article required significant simplification to produce a working algorithm.

Nevertheless, EOSA is valuable as an educational case study that illustrates the entire process from a biological idea to a software implementation: analyzing the original concept, identifying problems with a literal interpretation, finding compromises, integrating proven techniques, and objectively evaluating the results. This experience is valuable for understanding the nature of metaheuristic algorithms and developing the skills needed to design them.

Figure 3. Color coding of the algorithms for the corresponding tests

Figure 3. Histogram of the algorithm test results (on a scale from 0 to 100; the higher the score, the better, where 100 is the maximum possible theoretical result; the archive contains a script for calculating the ranking table)

Pros and cons of the EOSA algorithm:

Pros:

1. It can handle some simple types of problems.

Cons:

1. It has the most difficulty handling discrete functions.

An archive containing the latest versions of the algorithm code is attached to the article. The author of this article does not guarantee absolute accuracy in the descriptions of the canonical algorithms; many of them have been modified to improve their search capabilities. The conclusions and judgments presented in the articles are based on the results of the experiments conducted.

• GitHub: https://github.com/JQSakaJoo/Population-optimization-algorithms-MQL5

Software used in the article

#	Name	Type	Description
1	#C_AO.mqh	Include file	Parent class for population-based optimization algorithms
2	#C_AO_enum.mqh	Include file	Enumeration of population-based optimization algorithms
3	TestFunctions.mqh	Include file	Test function library
4	TestStandFunctions.mqh	Include file	Test bench function library
5	Utilities.mqh	Include file	Utility function library
6	CalculationTestResults.mqh	Include file	Script for calculating results for the comparison table
7	Testing AOs.mq5	Script	Unified test bench for all population-based optimization algorithms
8	Simple Use of Population Optimization Algorithms.mq5	Script	A simple example of using population-based optimization algorithms without visualization
9	Test_AO_EOSA.mq5	Script	Test bench for EOSA