Color Harmony Algorithm (CHA)

Contents

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

Introduction

The Color Harmony Algorithm (CHA), proposed by researchers Zaeimi and Goddosian in 2020, draws on a source of inspiration that seems too far removed from the computational nature of the problem. At first glance, the connection between Munsell's rules for color combinations in painting and numerical optimization of functions seems far-fetched. But if you look closely, the metaphor turns out to be less of an embellishment and more of a structural framework for the entire algorithm — and it is precisely this structure that makes it an interesting subject for analysis.

In standard metaheuristics, a population is a “crowd” in which all solutions are equal and interact according to the same rules. CHA introduces an explicit hierarchy: the population is organized into a hue circle divided into 100 sectors across 10 color groups, with one leader agent in each group; rules governing harmonious combinations apply between groups. Candidates are constructed not as mutations of individual solutions, but as color combinations — linear mixtures based on strictly defined rules for selecting parents.

From the perspective of a trader or researcher selecting a metaheuristic for their task, CHA has two appealing features. The first is a deliberately structured population: the algorithm does not collapse to a single point due to random clustering, because the geometry of the hue circle maintains population diversity by design. The second is explicit switching between the concentration (exploitation) and dispersion (exploration) phases based on a clear criterion rather than hidden hyperparameters.

However, these advantages come at a cost. By the end of this article, readers will have a complete picture: the algorithm’s framework, a working implementation in MQL5 using the standard C_AO framework, objective results from a standardized benchmark, and an understanding of the scenarios in which CHA makes sense to use and those in which it is better to choose something else.

Implementation of the Algorithm

If you were an artist looking for the perfect combination of colors for your painting, you would start with a few random colors on your palette and gradually mix them together until you achieved the desired shade. Sometimes you might experiment with completely new colors so you do not get stuck in the same color scheme. The color harmony algorithm works in much the same way. Each potential solution to a problem is represented as a “color” in a color space, and harmonious colors in art point to promising avenues for finding better solutions.

The algorithm’s geometry is based on the Munsell hue circle, divided into 100 sectors. These sectors are grouped into 10 “color groups,” each consisting of 10 sectors (red, orange, yellow, and so on). Each group has one special sector — the purest sector — which contains the group’s best solution, known as an agent. The purest sectors are numbered 5, 15, 25, ..., 95 — these are the centers of the groups. The remaining 90 solutions in the population are distributed across the remaining 90 sectors.

The algorithm consists of two main phases, which alternate as needed.

The concentration phase — searching for better solutions by combining those already found. The algorithm takes pairs of solutions selected according to the rules of harmony and forms perturbed linear combinations of them. This is similar to how an artist mixes harmonious colors in the hope of creating a promising new shade. For example, if we have a good red solution and a good blue solution, we try a “purple” solution between them, and if it is better than either parent, it earns a place in the population.

The dispersion phase — introducing diversity. When the population becomes too clustered (all solutions are too similar), the algorithm replaces several of the solutions closest to the center with new ones taken from the “color memory” — a repository of previously found promising but unused solutions. It is like how an artist might add an unexpected color to their palette to break away from a well-worn color scheme.

Phase switching is controlled by the measured population diversity Dpop and the current threshold Dth. If Dpop < Dth, the next iteration uses the dispersion phase, after which the threshold Dth is reduced by multiplying it by the coefficient damp (usually 0.5 or 0.96). Over time, the threshold becomes increasingly strict, and the algorithm focuses on pinpointing the optimum with ever-greater precision.

Harmony Templates

In Munsell’s theory, harmonious colors are not arbitrary pairs, but combinations that fit specific geometric templates on the hue circle; the authors of the article selected five such harmony templates:

• Type X: two arcs, each covering 26% of a circle, spaced 180° apart
• V-type: a single 26% arc
• Y-type: a small 5% arc and a large 26% arc, 180° apart
• T-type: one large 50% arc
• L-type: a small 5% arc and a medium 22% arc, 90° apart

When a template is selected for a specific agent, its gray arcs define which sectors of the hue circle are considered harmonious for that agent. The second parent for creating a new color (combination) can only be selected from solutions that fall within these sectors.

Rules for Forming Combinations

In each concentration phase, each of the 10 agents generates Ncomb combinations (4 by default). The type of each combination is determined by the PA (Power of Agent) parameter, which increases from 0 to Ncomb as the algorithm progresses. The higher the PA, the more intensive the exploitation, and the more combinations include agents as parents.

There are three types of combinations:

• Type (i) — an anchor agent is combined with a non-agent from the harmonious sectors. This is a local exploitation mechanism: the group leader pulls the best of its rank-and-file members toward itself.
• Type (ii) — the global best agent is combined with the anchor agent. This serves as a bridge between the groups: the lagging group can “latch onto” the leader’s achievements through its agent.
• Type (iii) — two non-agents from harmonious sectors. A pure exploration mechanism that is not tied to any leaders.

The combination is calculated using the following formula:

X_new = r1·X_a + r2·X_b + perturbation,

where r1 + r2 = 1, r1 is drawn uniformly from [0.1, 0.9] (we exclude degenerate cases where the candidate almost coincides with one of the parents), and the perturbation is a Gaussian perturbation whose amplitude gradually decreases over the course of the algorithm from 10% to 1% of the coordinate range.

It is worth pausing here. The perturbation is not described in the original article. This is a critically important extension; without it, the algorithm simply does not work on a complex benchmark. The explanation is simple: a pure linear combination of two parents always lies on the line segment between them — that is, within the convex hull of the parents. If the global optimum of the objective function lies outside the convex hull of the initial population, no mixing operations can reach it — the combinations converge toward the center but do not leave the initial volume.

For the smooth unimodal functions from the original article, this is not a problem: the optimum is usually inside the initial cloud. For functions with sharp peaks outside the initial cloud (Forest and Megacity in the benchmark), this completely blocks the algorithm. Perturbation is a minimal extension that allows candidates to move beyond the convex hull.

Dispersion phase

When Dpop &lt; Dth, the dispersion phase begins. The algorithm determines the center of the population (the coordinate-wise mean), finds the Ns solutions closest to the center (excluding agents — they are never replaced), and replaces them with new colors using the following formula: X_new = rcm·X_CM + (1 - rcm)·X_old,

where X_CM is a point from the color memory, rcm ∈ [rcm0, 1] (by default, rcm0 = 0.7, meaning that the new color is primarily determined by the memory but retains some connection to the old position).

Color Memory (CM) is a repository of Ksets sets, each containing Ns points. It stores combinations that were generated but not accepted into the population (because they were worse than the solutions occupying their slots during the concentration phase). The idea is that candidates rejected in one phase may prove useful in the dispersion phase: they already provide some diversity beyond the current center of mass.

Figure 1. The algorithm's hue circle.

Hue circle: 10 large, filled petal-like groups labeled with color names (Red, Orange, ..., Magenta) arranged around it. The thin lines outside the circle indicate the boundaries of the 100 sectors. The purest sectors are shown as agent points: circles with red outlines and white fill. Inside the circle are small white dots (standard solutions), and in the very center is a yellow circle (the global best).

Figure 2. Five mini-circles

Five Harmony Templates. Five mini-circles, each with dark “gray arcs,” as in the article. The anchor agent (red dot) is at the top, making it easy to compare the templates as your eyes move between them. The arc dimensions are exactly as specified in Appendix A: X — two at 26%, V — one at 26%, Y — 5% + 26%, T — one at 50%, L — 5% + 22%. Below is an explanation of the dimensions and angles.

Figure 3. Phases of the Algorithm's Operation

Two Phases. On the left, "Concentration": two parents (a red agent and a gray non-agent), a dotted line between them, with the linear combination point on the line, surrounded by a blue perturbation cloud containing candidates. On the right, "Dispersion": the black square represents the center of the population; the gray circles around it represent the points closest to the center (to be replaced); the purple diamonds represent points from CM memory; and the purple arrows indicate a jump. The large red circles represent agents at a distance; they are not affected. Below is a legend for the 5 types of symbols.

Let's write out the general pseudocode for the algorithm:

1. Initialization
   1.1. Set parameters (popSize=50, Ncomb=4, Ns=10, Ksets=4,
        damp=0.5, Dthf=0.01, rcm0=0.7)
   1.2. Set diversity thresholds:
        Dthi = 0.5 × max(rangeMax - rangeMin)
        Dth_cur = Dthi
        N_Dth = floor(log(Dthf/Dthi) / log(damp))
   1.3. Generate initial population uniformly in [rangeMin, rangeMax]
   1.4. Evaluate FF, place top-10 on the purest sectors as agents,
        the rest — uniformly by groups

2. Main loop
   2.1. If Dpop < Dth_cur:
        launch dispersion phase
        Dth_cur ← Dth_cur × damp
        iter_dp ← iter_dp + 1
        iter_cp ← 0
   2.2. Otherwise:
        launch concentration phase
        iter_cp ← iter_cp + 1
   2.3. Update the best solution fB, cB
   2.4. If the stop criteria is not reached — go to step 2.1

3. Concentration phase
   3.1. Calculate PA by (8)-(10) Zaeimi and Goddosian equations
   3.2. For each of 10 agents:
        - iterate through 5 harmony patterns
        - for each collect an H set of harmonious solutions
        - generate Ncomb combinations (i/ii/iii types by PA)
        - select the template set with maximum diversity
   3.3. Place candidates to slots with suitable parent sectors
   3.4. Evaluate FF
   3.5. Greedy selection: each slot compares the candidate with its
        snapshot, the best one wins
   3.6. Discarded — to CM_temp
   3.7. Re-assign group agents

4. Dispersion phase
   4.1. Find population center
   4.2. Find the Ns non-agents closest to the population center
   4.3. SCM source:
        if CM is filled — take the most diverse set
        otherwise — random points from CM_temp or random ones in space
   4.4. Generate new colors using the Munsell-based formula
   4.5. Set SCM to the main CM (replacing the least diverse set)

Let's move on to the implementation. S_CHA_Vector. This is a technical wrapper around a one-dimensional dynamic array v[]. The Init method takes the problem dimensionality and allocates the internal array v accordingly. The emulation of two-dimensional indexing occurs entirely at the point of access: X[i].v[c] is used instead of X[i][c]. The structure is used in three places within the class — in the population snapshot (snap_c[]), in the color memory buffer (CMtempC[]), and in the points fields of the S_CHA_CMSet structure.

//+------------------------------------------------------------------+
struct S_CHA_Vector
 {
  double  v          [];

  void               Init(int dim)
   {
    ArrayResize(v, dim);
   }
 };
//+------------------------------------------------------------------+

S_CHA_NewColor. An entry for a single candidate in the concentration phase. In addition to coordinates and fitness, the structure stores the indices of the hue circle sectors occupied by the parents used in the combination (sec_a and sec_b). These sectors are needed when placing a candidate in a population slot: according to the second-type rule from the article by Zaemi and Goddosian, a new color must be placed in the sector of one of its components. The Init method initializes the coordinate array; fitness is set to -DBL_MAX (the minimum value for a maximization problem), and the sector indices are set to -1 (unassigned).

//+------------------------------------------------------------------+
//| Entry for one candidate in the concentration phase:              |
//| coordinates + parent sector indices on the hue circle.           |
//+------------------------------------------------------------------+
struct S_CHA_NewColor
 {
  double  c          [];      // candidate coordinates
  double             f;       // fitness (set externally, OUTSIDE Moving)
  int                sec_a;   // first parent sector (0..99)
  int                sec_b;   // second parent sector (0..99)

  void               Init(int dim)
   {
    ArrayResize(c, dim);
    f = -DBL_MAX;
    sec_a = -1;
    sec_b = -1;
   }
 };

S_CHA_CMSet. A single set of poor but diverse colors in the algorithm's main memory. The points field is an array of Ns S_CHA_Vector wrappers, each containing a vector whose length equals the problem dimensionality. The diversity field stores the precomputed diversity of the set relative to the center of the search space. This allows the set with the lowest diversity to be quickly identified for replacement during the dispersion phase, without recalculating the metric each time. The `Init` method takes the number of points in the set and the problem dimension, allocates memory, and sets the diversity to zero.

//+------------------------------------------------------------------+
//| A single set of "poor but diverse" colors in main memory.        |
//+------------------------------------------------------------------+
struct S_CHA_CMSet
 {
  S_CHA_Vector points [];   // [Ns] wrappers, each of length dim
  double             diversity;   // Diversity of the set relative to the search center

  void               Init(int n_pts, int dim)
   {
    ArrayResize(points, n_pts);
    for(int i = 0; i < n_pts; i++)
      points [i].Init(dim);
    diversity = 0.0;
   }
 };

The C_AO_CHA class inherits from the standard C_AO class and implements the Init/Moving/Revision interface. The algorithm itself is divided into six basic subroutines, each responsible for a specific block of logic: placement on the hue circle after the first fitness evaluation, calculation of Power of Agent, generation and application of the concentration and dispersion phases, and updating of group leaders. These subroutines are private and are called from the Moving and Revision dispatcher methods.

Internal state. The iter_cp and iter_dp counters track the number of concentration-phase iterations since the last dispersion phase and the total number of dispersion phases completed, respectively. iter_cp is reset to zero every time a dispersion phase is executed — this is important for the correct calculation of Power of Agent. The parameters Dthi (initial threshold) and Dth_cur (current threshold) determine the switching point between phases; N_Dth is an estimate of the expected number of dispersion phases over the entire run, which is required for the Power of Agent formula. The PA field stores the current value of this parameter, which is recalculated before each concentration phase. The `nextIsDispersion` flag determines which phase will be run in the next call to `Moving`.

The geometry of the hue circle is stored in two arrays. For each individual in the population, sectorOf[] stores its sector, from 0 to 99; its group is determined by sectorOf[i] / 10. For each group, agentIdx[10] stores the index of its agent in the a[] array, or −1 if no agent has been assigned to the group.

Concentration-phase buffers. Concentration-phase candidates are stored separately from the population, in an array `candidates[]` of size 10·Ncomb. At the time of generation, candidates are not yet assigned to population slots — only after the candidates have been constructed and their parent sectors recorded does the algorithm decide which slots to place them in for FF evaluation. The `candidateUsed[]` flag marks the candidates accepted into the population based on the results of greedy selection; those that were not accepted are sent to the temporary color memory `CMtempC[]`.

A population snapshot (snap_c, snap_f, snap_sectorOf) is taken at the beginning of ConcentrationPhase_Generate and is used for two purposes. First, it preserves the old slot coordinates in case a candidate temporarily placed in a slot for fitness evaluation turns out to be worse than the slot's original contents — in which case the slot is restored from the snapshot. Second, the snapshot of fitness values and sectors is used in all hue-circle rules to prevent confusion between the “state before Moving” and the “state with a temporarily placed candidate.”

The association between a population slot and a candidate's index is stored in `slotCandidate[]`: if a candidate occupies a slot, its index in `candidates[]` is stored there; otherwise, it is -1. A parallel Boolean array `isCandidate[]` is needed to quickly check whether a slot is occupied by a candidate without comparing it to −1.

Color memory. The main memory, CM[], has a capacity of Ksets sets, each containing Ns points. The CMcount field shows how many sets have actually been stored — as long as CMcount < Ksets, new sets are added; once the memory is full, the least diverse set is overwritten by a new one if the new set has greater diversity. The temporary memory CMtempC[] is a flat circular buffer with a capacity of Ksets·Ns that stores rejected candidates from the concentration phase. The CMtempCnt counter is incremented with each write; when it reaches twice the buffer capacity, it is reset to CMtempCap to prevent it from growing indefinitely. During the dispersion phase, if CM is not yet full, CMtempC serves as the source of diverse points.

//+------------------------------------------------------------------+
//|                       Color Harmony Algorithm                    |
//+------------------------------------------------------------------+
class C_AO_CHA : public C_AO
 {
public:
                    ~C_AO_CHA() {}
                     C_AO_CHA()
   {
    ao_name = "CHA";
    ao_desc = "Color Harmony Algorithm";
    ao_link = "https://www.mql5.com/en/articles/18132";

    popSize = 50;
    Ncomb   = 4;       // combinations per template (per agent)
    Ns      = 10;      // replacements during the dispersion phase
    Ksets   = 4;       // number of sets in CM
    damp    = 0.5;     // diversity threshold compression factor
    Dthf    = 0.01;    // final diversity threshold
    rcm0    = 0.7;     // lower bound on the CM contribution fraction during the dispersion phase

    ArrayResize(params, 7);
    params [0].name = "popSize"; params [0].val = popSize;
    params [1].name = "Ncomb";   params [1].val = Ncomb;
    params [2].name = "Ns";      params [2].val = Ns;
    params [3].name = "Ksets";   params [3].val = Ksets;
    params [4].name = "damp";    params [4].val = damp;
    params [5].name = "Dthf";    params [5].val = Dthf;
    params [6].name = "rcm0";    params [6].val = rcm0;
   }

  void               SetParams()
   {
    popSize = (int)params [0].val;
    Ncomb   = (int)params [1].val;
    Ns      = (int)params [2].val;
    Ksets   = (int)params [3].val;
    damp    = params      [4].val;
    Dthf    = params      [5].val;
    rcm0    = params      [6].val;

    //--- safeguards 
    if(popSize < 20)
      popSize = 20;          // at least 10 agents + 10 non-agents
    if(Ncomb   < 1)
      Ncomb   = 1;
    if(Ncomb * 10 > popSize - 10)
      Ncomb   = (popSize - 10) / 10;
    if(Ns      < 1)
      Ns      = 1;
    if(Ns      > popSize - 10)
      Ns      = popSize - 10;
    if(Ksets   < 1)
      Ksets   = 1;
    if(damp    <= 0 || damp >= 1)
      damp    = 0.5;
    if(Dthf    <= 0)
      Dthf    = 0.01;
    if(rcm0    < 0 || rcm0 >= 1)
      rcm0    = 0.7;
   }

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

  void               Moving();
  void               Revision();

  //--- visible parameters 
  int                Ncomb;  // combinations per template (per agent)
  int                Ns;     // replacements during the dispersion phase
  int                Ksets;  // number of sets in CM
  double             damp;   // diversity threshold compression factor
  double             Dthf;   // final diversity threshold
  double             rcm0;   // lower bound on the CM contribution fraction during the dispersion phase

private: //---
  //--- global algorithm state 
  int                iter_cp;       // consecutive concentration phase iterations
  int                iter_dp;       // number of iterations in the dispersion phase (= number of Dth updates)
  int                N_Dth;         // floor(log(Dthf/Dthi)/log(damp))
  double             Dthi;          // initial diversity threshold
  double             Dth_cur;       // current threshold
  int                PA;            // agent power at the current concentration iteration
  bool               nextIsDispersion; // flag: the next Moving call is the dispersion phase

  //--- hue circle 
  // sectorOf [i] ∈ [0..99] — the hue circle sector for the i-th individual in a[].
  // Sectors 5, 15, 25, ..., 95 are the “purest sectors”; agents are placed in them.
  // A group is defined as sectorOf[i] / 10 (0..9).
  int    sectorOf    [];          // [popSize]

  // agentIdx [g] — the index in a[] of the agent in group g (0..9). −1 if not specified.
  int    agentIdx    [10];

  //--- concentration phase candidate buffer 
  // 10·Ncomb candidates from all 10 agents. They are stored SEPARATELY from a[],
  // because at the generation stage, we do not yet know which slot
  // to place each candidate in — this is decided in Revision according to the hue circle rules.
  // During Moving(), candidates are "temporarily" placed in slots
  // occupied by non-agents for FF evaluation; in Revision, the candidates are assigned according to the rules.
  S_CHA_NewColor candidates [];     // [10·Ncomb]
  bool           candidateUsed [];  // [10·Ncomb], true = accepted onto the hue circle

  //--- Snapshot of the current population during the concentration phase 
  // Needed because Moving overwrites the coordinates of non-agents with
  // candidates. If a candidate is rejected, the slot is restored
  // from the snapshot.
  S_CHA_Vector snap_c        [];   // [popSize] wrappers, each of length coords
  double       snap_f        [];   // [popSize]
  int          snap_sectorOf [];   // [popSize]

  //--- the mapping "slot in a[] → candidate index" during the phase
  // -1 if the slot is not occupied by a candidate (an agent or a free non-agent).
  int     slotCandidate [];         // [popSize]
  bool    isCandidate   [];         // [popSize]

  //--- color memory 
  S_CHA_CMSet  CM     [];       // [&lt;= Ksets]
  int                CMcount;         // Actual CM size

  // Temporary memory — a flat buffer of candidates not placed on the hue circle
  S_CHA_Vector CMtempC [];      // [&lt;= Ksets*Ns] wrappers, each of length coords
  int                CMtempCnt;
  int                CMtempCap;

  //--- internal subroutines 
  bool               InitialPopulationWithDiversity();
  void               PlaceOnHueCircle();                     // after the initial FF evaluation

  void               ConcentrationPhase_Generate();          // in Moving()
  void               ConcentrationPhase_UpdateCircle();      // in Revision()

  void               DispersionPhase_Generate();             // in Moving()
  void               DispersionPhase_UpdateCircle();         // in Revision()

  void               ComputePA();                            // PA_t according to equations (8)–(10)
  void               UpdateAgentsByGroup();                  // Best in the group → agent

  //--- Utilities 
  double             DiversityPopulation();
  void               PopulationCenter(double &out []);

  // Templates: returns an array of sectors for each agent
  // (0..99) that fall within the gray arcs of the template when it is rotated to align with this agent.
  // templateType: 0=X, 1=V, 2=Y, 3=T, 4=L
  void               GreySectorsForTemplate(int templateType,
      int anchorSector,
      int &out []);

  int                GroupOf(int sector) { return sector / 10; }
  int                PurestSector(int group)  { return group * 10 + 5; }
 };
//+------------------------------------------------------------------+

Init is the main initialization procedure. After the standard initialization of the base class, the maximum coordinate range is calculated, and the initial diversity threshold Dthi is set to 0.5 · max(rangeMax − rangeMin). Next, the expected number of dispersion phases for the entire run is estimated using the formula N_Dth = floor(log(Dthf/Dthi) / log(damp)); this number is used in the Power of Agent formulas. All counters are reset to zero, and the `nextIsDispersion` flag is set to `false`. Next, memory is allocated for all buffers: the hue circle, concentration phase candidates, the population snapshot, slot-candidate associations, and the main and temporary color memory. Finally, `InitialPopulationWithDiversity` is called to generate the initial population.

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

//--- initial diversity threshold Dthi = 0.5 · max(rangeMax − rangeMin) 
  double maxRange = 0.0;
  for(int c = 0; c < coords; c++)
   {
    double r = rangeMax [c] - rangeMin [c];
    if(r > maxRange)
      maxRange = r;
   }
  Dthi    = 0.5 * maxRange;
  Dth_cur = Dthi;

//--- N_Dth = floor( log(Dthf/Dthi) / log(damp) )
  if(Dthi > 0.0 && Dthf > 0.0 && damp > 0.0 && damp != 1.0)
    N_Dth = (int)MathFloor(MathLog(Dthf / Dthi) / MathLog(damp));
  else
    N_Dth = 1;
  if(N_Dth < 1)
    N_Dth = 1;

//--- counters and flags 
  iter_cp = 0;
  iter_dp = 0;
  PA      = 0;
  nextIsDispersion = false;

//--- hue circle
  ArrayResize(sectorOf, popSize);
  ArrayInitialize(sectorOf, -1);
  for(int g = 0; g < 10; g++)
    agentIdx [g] = -1;

//--- candidate buffer: 10 agents · Ncomb combinations 
  int Kcand = 10 * Ncomb;
  ArrayResize(candidates, Kcand);
  for(int i = 0; i < Kcand; i++)
    candidates [i].Init(coords);

  ArrayResize(candidateUsed, Kcand);
  ArrayInitialize(candidateUsed, false);

//--- snapshot and slot↔candidate associations 
  ArrayResize(snap_c, popSize);
  for(int i = 0; i < popSize; i++)
    snap_c [i].Init(coords);
  ArrayResize(snap_f,        popSize);
  ArrayResize(snap_sectorOf, popSize);

  ArrayResize(slotCandidate, popSize);
  ArrayResize(isCandidate,   popSize);
  ArrayInitialize(slotCandidate, -1);
  ArrayInitialize(isCandidate,   false);

//--- color memory 
  ArrayResize(CM, Ksets);
  for(int s = 0; s < Ksets; s++)
    CM [s].Init(Ns, coords);
  CMcount   = 0;

  CMtempCap = Ksets * Ns;
  ArrayResize(CMtempC, CMtempCap);
  for(int i = 0; i < CMtempCap; i++)
    CMtempC [i].Init(coords);
  CMtempCnt = 0;

//--- initial population with Dth_i check 
  if(!InitialPopulationWithDiversity())
    return false;

  return true;
 }
//+------------------------------------------------------------------+

InitialPopulationWithDiversity generates an initial population with a population diversity check. On each attempt (up to 50), a population is sampled uniformly within the ranges [rangeMin, rangeMax], and its diversity is then calculated. If it reaches or exceeds Dthi, the attempts are stopped. In practice, when the distribution is uniform over large ranges, a single attempt is sufficient. If population diversity has not reached the threshold after all 50 attempts, the last generated population is retained — this is a scenario for very narrow or specific search ranges where it is geometrically impossible to achieve high diversity.

//+------------------------------------------------------------------------------------+
//| Generate an initial population with population diversity ≥ Dthi.                   |
//| Stores the coordinates in a[i].cP[] (they will be transferred to a[i].c via Moving)|
//+------------------------------------------------------------------------------------+
bool C_AO_CHA::InitialPopulationWithDiversity()
 {
// Regeneration loop; in practice, it converges in 1–2 attempts with a
// uniform distribution over [rangeMin, rangeMax].
  const int maxAttempts = 50;
  for(int attempt = 0; attempt < maxAttempts; attempt++)
   {
    for(int i = 0; i < popSize; i++)
      for(int c = 0; c < coords; c++)
        a [i].cP [c] = u.RNDfromCI(rangeMin [c], rangeMax [c]);

    // Calculating diversity using cP[]
    double D = 0.0;
    for(int c = 0; c < coords; c++)
     {
      double mean = 0.0;
      for(int i = 0; i < popSize; i++)
        mean += a [i].cP [c];
      mean /= popSize;
      double dj = 0.0;
      for(int i = 0; i < popSize; i++)
        dj += MathAbs(a [i].cP [c] - mean);
      dj /= popSize;
      D += dj;
     }
    D /= coords;

    if(D >= Dthi)
      return true;
   }

// If convergence is not achieved within maxAttempts, keep the last population.
  return true;
 }
//+------------------------------------------------------------------+

Moving — the main dispatcher. On the first call (when the `revision` flag is still `false`), the method simply transfers the initial population from `cP[]` to `c[]` via `SeInDiSp` for discretization — this is necessary so that the outer loop of the test bench can calculate the FF for the initial population. The placement on the hue circle and the assignment of agents are done in the first Revision, not here. In all subsequent calls, the method checks the `nextIsDispersion` flag and calls either `DispersionPhase_Generate` or a combination of `ComputePA` and `ConcentrationPhase_Generate`. The method itself does not do any substantive work — it just handles routing.

//+------------------------------------------------------------------+
//|                            Moving                                |
//+------------------------------------------------------------------+
//  Dispatcher logic:
//   • First call (revision == false): copy cP → c so that the outer
//     loop can calculate the FF for the initial population. Placement on the hue circle
//     and assignment of agents take place in the first Revision().
//   • Subsequent calls: either the concentration phase (generation of
//     10·Ncomb_eff candidates in a[].c), or the dispersion phase
//     (generating Ns candidates to replace the individuals closest to the center).
//+------------------------------------------------------------------+
void C_AO_CHA::Moving()
 {
  if(!revision)
   {
    // First pass—we simply copy the initial population into c[]
    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]);
    return;
   }

  if(nextIsDispersion)
   {
    DispersionPhase_Generate();
   }
  else
   {
    ComputePA();
    ConcentrationPhase_Generate();
   }
 }
//+------------------------------------------------------------------+

Revision — the second main dispatcher, called immediately after the FF evaluation in each slot. First, a standard loop updates each individual's personal best and the global best for the entire population. Next, if this is the first call to `Revision`, `PlaceOnHueCircle` is called to place the population — sorted by fitness — on the hue circle and assign 10 agents; after that, the `revision` flag is set to `true`, and the method terminates.

On subsequent calls, the hue circle update rules for the phase just evaluated are applied: either ConcentrationPhase_UpdateCircle or DispersionPhase_UpdateCircle. After the rules are applied, the counters and the diversity threshold are updated, and the method determines which phase to run next: if the measured population diversity has fallen below the current threshold Dth_cur, the nextIsDispersion flag is set.

//+------------------------------------------------------------------+
//|                           Revision                               |
//+------------------------------------------------------------------+
//  Dispatcher logic:
//   • If revision == false: this is the first evaluation of the population—update
//     fB/cB, copy f→fP in each slot, place the individuals on the
//     hue circle, assign agents. revision = true.
//   • Otherwise: apply the hue circle update rules for the most recent phase
//     (concentration or dispersion), switch phases based on the
//     diversity threshold, update Dth_cur as needed.
//+------------------------------------------------------------------+
void C_AO_CHA::Revision()
 {
//--- Update personal bests and the global best 
  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);
     }
   }

  if(!revision)
   {
    // First Revision: Simply place the population on the hue circle
    PlaceOnHueCircle();
    revision = true;
    return;
   }

//--- update the circle using the results of the newly evaluated phase 
  if(nextIsDispersion)
   {
    DispersionPhase_UpdateCircle();
    iter_dp ++;
    Dth_cur *= damp;
    iter_cp = 0;                  // Reset — a new series of concentration phases begins
    nextIsDispersion = false;
   }
  else
   {
    ConcentrationPhase_UpdateCircle();
    iter_cp ++;
   }

//--- Decide which phase comes next 
  if(DiversityPopulation() < Dth_cur)
    nextIsDispersion = true;
 }
//+------------------------------------------------------------------+

DiversityPopulation — an implementation of equations (2)–(4) from the original article. For each coordinate, the population mean is calculated, followed by the mean absolute deviation from that mean. The final diversity value is the average of the deviations across all coordinates. The metric is robust to outliers and performs well across a wide range of problem dimensionalities.

//+-------------------------------------------------------------------+
//|              Utilities: Population Diversity and Population Center|
//+-------------------------------------------------------------------+
double C_AO_CHA::DiversityPopulation()
 {
  double D = 0.0;
  for(int c = 0; c < coords; c++)
   {
    double mean = 0.0;
    for(int i = 0; i < popSize; i++)
      mean += a [i].c [c];
    mean /= popSize;

    double dj = 0.0;
    for(int i = 0; i < popSize; i++)
      dj += MathAbs(a [i].c [c] - mean);
    dj /= popSize;

    D += dj;
   }
  return D / coords;
 }
//+------------------------------------------------------------------+

PopulationCenter — a simple calculation of the arithmetic mean for each coordinate. It is used in the dispersion phase to identify the slots closest to the center of the population — these are the ones to be replaced.

//+------------------------------------------------------------------+
void C_AO_CHA::PopulationCenter(double &out [])
 {
  ArrayResize(out, coords);
  for(int c = 0; c < coords; c++)
   {
    double s = 0.0;
    for(int i = 0; i < popSize; i++)
      s += a [i].c [c];
    out [c] = s / popSize;
   }
 }
//+------------------------------------------------------------------+

ComputePA — calculation of the Power of Agent parameter using equations (8)–(10) in the article. First, Ni = floor(iter_dp · Ncomb / N_Dth) is calculated, capped at Ncomb — this is the “exploitation pressure,” which increases over time. Then, step = (Ncomb + 1) / (Ncomb + Ni) is the step size with which iter_cp approaches saturation. The final value of PA = floor(iter_cp / step) is also capped at Ncomb. Meaning: The higher the PA, the more combinations are generated in the concentration phase with the agent as one of the parents — that is, the stronger the exploitation of group leaders.

//+--------------------------------------------------------------------------+
//|       Power of Agent: Equations (8)–(10) from the article (Zaeimi 2020)  |
//+--------------------------------------------------------------------------+
void C_AO_CHA::ComputePA()
 {
  int Ni;
  if(N_Dth > 0)
   {
    double denom = (double)N_Dth / (double)Ncomb;
    Ni = (int)MathFloor((double)iter_dp / denom);
    if(Ni > Ncomb)
      Ni = Ncomb;
   }
  else
    Ni = 0;
  if(Ni < 0)
    Ni = 0;

  double step = (double)(Ncomb + 1) / (double)(Ncomb + Ni);
  if(step <= 0)
    step = 1.0;

  PA = (int)MathFloor((double)iter_cp / step);
  if(PA > Ncomb)
    PA = Ncomb;
  if(PA < 0)
    PA = 0;
 }
//+------------------------------------------------------------------+

UpdateAgentsByGroup — is called after the hue circle rules have been applied in both phases. For each of the 10 groups, the individual with the best fitness score is selected from among those whose sectors belong to that group. If the newly found leader is not the same as the group's previous agent, a sector swap occurs — the previous agent moves to the former leader's sector, and the new agent takes over the group's purest sector. This implements the rule stated in the article: “if replaced colors have better fitness value than their related agents, their positions on the hue circle will be exchanged”.

//+---------------------------------------------------------------------------------------+
//|  Group agent = individual with the best fitness. If the new best individual is not    |
//|  in the purest sector, swap the sectors of the old and new agents.                    |
//+---------------------------------------------------------------------------------------+
void C_AO_CHA::UpdateAgentsByGroup()
 {
  for(int g = 0; g < 10; g++)
   {
    int    oldAg  = agentIdx [g];
    int    bestI  = oldAg;
    double bestF  = (oldAg >= 0) ? a [oldAg].f : -DBL_MAX;

    for(int i = 0; i < popSize; i++)
     {
      if(i == oldAg)
        continue;
      if(sectorOf [i] / 10 != g)
        continue;
      if(a [i].f > bestF)
       {
        bestF = a [i].f;
        bestI = i;
       }
     }

    if(bestI >= 0 && bestI != oldAg)
     {
      // Sector swap: new agent → the purest sector in the group
      if(oldAg >= 0)
       {
        int tmp = sectorOf [bestI];
        sectorOf [bestI] = sectorOf [oldAg];
        sectorOf [oldAg] = tmp;
       }
      else
       {
        sectorOf [bestI] = PurestSector(g);
       }
      agentIdx [g] = bestI;
     }
   }
 }
//+------------------------------------------------------------------+

PlaceOnHueCircle — the initial placement of the population on the hue circle following the initial FF evaluation. First, the individuals are sorted by fitness in descending order (bubble sort based on indices, without rearranging the a[] structures themselves). The top 10 are then assigned to the purest sectors, with the groups visited in random order — the group order is shuffled using the classic Fisher–Yates algorithm. The remaining popSize − 10 individuals are distributed among the groups as evenly as possible: first, the base number of individuals per group is calculated as (popSize − 10) / 10, then the remainder, (popSize − 10) % 10, is distributed among the first groups. Within each group, the available sectors (all except the purest sector) are shuffled and assigned to individuals in turn. If there are more than nine individuals in a group, sectors are reused.

//+------------------------------------------------------------------+
//|     Initial placement of the population on the hue circle after  |
//|     the initial FF evaluation.                                   |
//+------------------------------------------------------------------+
void C_AO_CHA::PlaceOnHueCircle()
 {
// Sorting indices by fitness (in descending order: best first)
  int idx [];
  ArrayResize(idx, popSize);
  for(int i = 0; i < popSize; i++)
    idx [i] = i;

  for(int i = 0; i < popSize - 1; i++)
    for(int j = i + 1; j < popSize; j++)
      if(a [idx [j]].f > a [idx [i]].f)
       {
        int t = idx [i];
        idx [i] = idx [j];
        idx [j] = t;
       }

// Top 10 → agents in the purest sectors, with groups in random order
  int groupOrder [10];
  for(int g = 0; g < 10; g++)
    groupOrder [g] = g;
  for(int g = 9; g > 0; g--)
   {
    int j = (int)u.RNDfromCI(0.0, (double)(g + 1));
    if(j > g)
      j = g;
    int t = groupOrder [g];
    groupOrder [g] = groupOrder [j];
    groupOrder [j] = t;
   }

  for(int k = 0; k < 10; k++)
   {
    int slot = idx [k];
    int g    = groupOrder [k];
    sectorOf [slot] = PurestSector(g);
    agentIdx [g]    = slot;
   }

// The remaining popSize-10 are evenly distributed among the groups
  int rem        = popSize - 10;
  int perGroup   = rem / 10;
  int extra      = rem % 10;

  int sectorList [];
  ArrayResize(sectorList, rem);
  int pos = 0;

  for(int g = 0; g < 10; g++)
   {
    int n = perGroup + (g < extra ? 1 : 0);

    // Free sectors in the group (excluding the purest sector)
    int freeSecs [9];
    int fsc = 0;
    for(int s = 0; s < 10; s++)
     {
      int sec = g * 10 + s;
      if(sec == PurestSector(g))
        continue;
      freeSecs [fsc++] = sec;
     }

    // Shuffle
    for(int s = fsc - 1; s > 0; s--)
     {
      int j = (int)u.RNDfromCI(0.0, (double)(s + 1));
      if(j > s)
        j = s;
      int t = freeSecs [s];
      freeSecs [s] = freeSecs [j];
      freeSecs [j] = t;
     }

    // Assign n sectors (with repeats if n &gt; 9)
    for(int s = 0; s < n; s++)
      sectorList [pos++] = freeSecs [s % fsc];
   }

// Shuffle sectorList and assign the sectors
  for(int i = ArraySize(sectorList) - 1; i > 0; i--)
   {
    int j = (int)u.RNDfromCI(0.0, (double)(i + 1));
    if(j > i)
      j = i;
    int t = sectorList [i];
    sectorList [i] = sectorList [j];
    sectorList [j] = t;
   }

  for(int k = 10; k < popSize; k++)
    sectorOf [idx [k]] = sectorList [k - 10];
 }
//+------------------------------------------------------------------+

GreySectorsForTemplate. Returns an array of sectors that fall within the “gray arcs” of one of the five harmony templates when the template is rotated to align with the anchor sector anchor. The arc sizes correspond exactly to Appendix A of the original article: X-type — two arcs, each covering 26% of the circle, 180° apart; V-type—one 26% arc; Y-type — a small 5% arc and a large 26% arc, 180° apart; T-type — one 50% arc; L-type — a small 5% arc and a medium 22% arc, 90° apart. All sector offsets are taken modulo 100 to handle wraparound at the circle's boundary. This method is used in the concentration phase to select candidate parents for a combination.

//+------------------------------------------------------------------+
//|  Gray sectors for the template rotated to align with anchor.     |
//|  Arc sizes: X/V — 26%, T — 50%, Y — 5%+26%, L — 5%+22%.          |
//|  Angles between arcs: X/Y — 180°, L — 90°.                       |
//+------------------------------------------------------------------+
void C_AO_CHA::GreySectorsForTemplate(int templateType, int anchor, int &out [])
 {
  ArrayResize(out, 0);

// helper lambda emulation: adding an interval (offset ± half) relative to anchor
// range [center - low, center + high], inclusive, modulo 100
// To avoid creating an extra function, we write it inline.

  if(templateType == 0)
   {
    // Type X: two arcs of 26% each (26 sectors each), with an angle of 180° between their bisectors
    for(int d = -13; d <= 12; d++)
     {
      int s1 = ((anchor + d) % 100 + 100) % 100;
      int s2 = ((anchor + 50 + d) % 100 + 100) % 100;
      int n = ArraySize(out);
      ArrayResize(out, n + 2);
      out [n]     = s1;
      out [n + 1] = s2;
     }
   }
  else
    if(templateType == 1)
     {
      // Type V: one 26% arc
      for(int d = -13; d <= 12; d++)
       {
        int s = ((anchor + d) % 100 + 100) % 100;
        int n = ArraySize(out);
        ArrayResize(out, n + 1);
        out [n] = s;
       }
     }
    else
      if(templateType == 2)
       {
        // Type Y: a small 5% arc centered at anchor + a large 26% arc centered at (anchor + 50)
        for(int d = -2; d <= 2; d++)
         {
          int s = ((anchor + d) % 100 + 100) % 100;
          int n = ArraySize(out);
          ArrayResize(out, n + 1);
          out [n] = s;
         }
        for(int d = -13; d <= 12; d++)
         {
          int s = ((anchor + 50 + d) % 100 + 100) % 100;
          int n = ArraySize(out);
          ArrayResize(out, n + 1);
          out [n] = s;
         }
       }
      else
        if(templateType == 3)
         {
          // Type T: one 50% arc
          for(int d = -25; d <= 24; d++)
           {
            int s = ((anchor + d) % 100 + 100) % 100;
            int n = ArraySize(out);
            ArrayResize(out, n + 1);
            out [n] = s;
           }
         }
        else
          if(templateType == 4)
           {
            // Type L: a small 5% arc centered at anchor + a large 22% arc centered at (anchor + 25)
            for(int d = -2; d <= 2; d++)
             {
              int s = ((anchor + d) % 100 + 100) % 100;
              int n = ArraySize(out);
              ArrayResize(out, n + 1);
              out [n] = s;
             }
            for(int d = -11; d <= 10; d++)
             {
              int s = ((anchor + 25 + d) % 100 + 100) % 100;
              int n = ArraySize(out);
              ArrayResize(out, n + 1);
              out [n] = s;
             }
           }
 }
//+------------------------------------------------------------------+

ConcentrationPhase_Generate — the main procedure for generating candidates. The first step is to create a snapshot of the entire population. The second step is to reset the slot↔candidate mappings. The third and main step is a loop over the 10 agents.

For each agent, all 5 harmony templates are examined. For each template, a set H is formed — the indices of the slots whose sectors fall within the gray arcs. This set is divided into Hagents (agents in gray arcs) and Hnon (non-agents in gray arcs). If Hnon is empty, the template is skipped — there is no material to construct combinations from.

Before generating combinations, the adaptive perturbation amplitude (perturbScale) is calculated. It is tied to the current diversity threshold, Dth_cur: as long as the threshold is high (the population is sparse), the amplitude is about 10% of the range; as the threshold decreases toward Dthf, the amplitude drops to 1%. This is not a static constant from the original formulation, but a critical extension of the algorithm — without it, all combinations lie strictly within the convex hull of their parents and cannot leave it.

Next, Ncomb combinations are generated. If the combination index m &lt; PA, it is constructed according to type (i) or (ii): with a probability of 1/3, the combination “global best agent × anchor agent” is selected (Type ii — a bridge between leaders); otherwise, the “anchor agent × random non-agent from the gray arcs” combination is chosen (Type i — local exploitation). If m ≥ PA, a combination of type (iii) is constructed — two distinct non-agents from the gray arcs.

The combination itself is a linear combination of the two parents' coordinates with coefficients r1 and r2 = 1 − r1, where r1 is drawn uniformly from [0.1, 0.9] (to avoid degenerate cases where the candidate almost coincides with one of the parents). A Gaussian perturbation with amplitude `perturbScale` is added to the result. The final value is clamped to the range [rangeMin, rangeMax].

After generating all Ncomb combinations for the current template, their aggregate diversity relative to the anchor agent is calculated. The template with the greatest diversity wins — its combinations are saved in bestSet. If no template returns any results (an exceptional case), a fallback is used with a perturbation of the anchor agent's coordinates within a radius of 10% of the range.

The best set is stored in the shared candidates[] array. After the loop over agents is complete, candidates are placed in non-agent slots for FF evaluation. The placement follows the rule stated in the article: a new color is assigned to the sector of one of its parents. First, the system searches for a non-agent slot with a suitable sector; if none is found, it falls back to any available non-agent slot. Agent slots are left unchanged during the concentration phase — their coordinates are restored from the snapshot so that their fitness isn't recalculated unnecessarily during the next FF evaluation.

//+------------------------------------------------------------------+
//|                    CONCENTRATION PHASE — generation              |
//+------------------------------------------------------------------+
void C_AO_CHA::ConcentrationPhase_Generate()
 {
//--- 1. Save a population snapshot 
  for(int i = 0; i < popSize; i++)
   {
    for(int c = 0; c < coords; c++)
      snap_c [i].v [c] = a [i].c [c];
    snap_f        [i] = a [i].f;
    snap_sectorOf [i] = sectorOf [i];
   }

//--- 2. Reset slot↔candidate associations and usage flags 
  ArrayInitialize(slotCandidate, -1);
  ArrayInitialize(isCandidate,   false);
  ArrayInitialize(candidateUsed, false);

//--- 3. For each agent: 5 templates → choose the best set 
  int Kcand   = 10 * Ncomb;
  int candPos = 0;

// a temporary buffer for the set of Ncomb combinations
  S_CHA_NewColor combo [];
  ArrayResize(combo, Ncomb);
  for(int i = 0; i < Ncomb; i++)
    combo [i].Init(coords);

  S_CHA_NewColor bestSet [];
  ArrayResize(bestSet, Ncomb);
  for(int i = 0; i < Ncomb; i++)
    bestSet [i].Init(coords);

  for(int g = 0; g < 10; g++)
   {
    int agSlot = agentIdx [g];
    if(agSlot < 0)
      continue;
    int anchorSec = snap_sectorOf [agSlot];

    double bestDiv = -1.0;

    for(int t = 0; t < 5; t++)
     {
      int greySectors [];
      GreySectorsForTemplate(t, anchorSec, greySectors);

      // Collect H — the indices of the slots whose sectors are within the gray arcs
      int H [];
      ArrayResize(H, 0);
      for(int i = 0; i < popSize; i++)
       {
        for(int s = 0; s < ArraySize(greySectors); s++)
         {
          if(snap_sectorOf [i] == greySectors [s])
           {
            int n = ArraySize(H);
            ArrayResize(H, n + 1);
            H [n] = i;
            break;
           }
         }
       }

      // Divide H into agents and non-agents
      int Hagents [], Hnon [];
      ArrayResize(Hagents, 0);
      ArrayResize(Hnon,    0);
      for(int i = 0; i < ArraySize(H); i++)
       {
        bool ag = false;
        for(int gg = 0; gg < 10; gg++)
          if(agentIdx [gg] == H [i])
           {
            ag = true;
            break;
           }
        int n;
        if(ag)
         {
          n = ArraySize(Hagents);
          ArrayResize(Hagents, n + 1);
          Hagents [n] = H [i];
         }
        else
         {
          n = ArraySize(Hnon);
          ArrayResize(Hnon, n + 1);
          Hnon [n] = H [i];
         }
       }

      if(ArraySize(Hnon) == 0)
        continue;  // There is no "material" for combinations

      // Find the global best agent — required for Type (ii)
      int    bestAg  = agSlot;
      double bestAgF = snap_f [agSlot];
      for(int gg = 0; gg < 10; gg++)
       {
        int s = agentIdx [gg];
        if(s < 0)
          continue;
        if(snap_f [s] > bestAgF)
         {
          bestAgF = snap_f [s];
          bestAg  = s;
         }
       }

      // Adaptive perturbation amplitude: tied to the current
      // population diversity threshold Dth_cur, rather than the iter_dp counter. Logic:
      // When the population is dispersed (Dth_cur is high), apply a strong perturbation;
      // when it is clustered (Dth_cur is low), focus precisely.
      // 10% when Dth_cur = Dthi (start), 1% when Dth_cur = Dthf (finish).
      double ratio = (Dthi > 0.0)
                     ? MathMax(0.0, MathMin(1.0, (Dth_cur - Dthf) / (Dthi - Dthf)))
                     : 0.0;
      double perturbScale = 0.01 + 0.09 * ratio;

      // Generate Ncomb combinations
      for(int m = 0; m < Ncomb; m++)
       {
        int    ia = -1, ib = -1;

        if(m < PA)
         {
          // With a probability of 1/3 — Type (ii): global best × anchor agent.
          // This provides strong exploitation, forming a bridge between the best group leaders.
          // Otherwise — Type (i): anchor agent × non-agent (local exploitation).
          if(bestAg != agSlot && u.RNDfromCI(0.0, 1.0) < 0.33)
           {
            ia = bestAg;
            ib = agSlot;
           }
          else
           {
            ia = agSlot;
            int kk = (int)u.RNDfromCI(0.0, (double)ArraySize(Hnon));
            if(kk >= ArraySize(Hnon))
              kk = ArraySize(Hnon) - 1;
            ib = Hnon [kk];
           }
         }
        else
         {
          // Type (iii): two non-agents, exploration
          if(ArraySize(Hnon) >= 2)
           {
            int k1 = (int)u.RNDfromCI(0.0, (double)ArraySize(Hnon));
            if(k1 >= ArraySize(Hnon))
              k1 = ArraySize(Hnon) - 1;
            int k2;
            int safety = 0;
            do
             {
              k2 = (int)u.RNDfromCI(0.0, (double)ArraySize(Hnon));
              if(k2 >= ArraySize(Hnon))
                k2 = ArraySize(Hnon) - 1;
              safety++;
              if(safety > 20)
               {
                k2 = (k1 + 1) % ArraySize(Hnon);
                break;
               }
             }
            while(k2 == k1);
            ia = Hnon [k1];
            ib = Hnon [k2];
           }
          else
           {
            ia = (ArraySize(Hagents) > 0) ? Hagents [0] : agSlot;
            ib = Hnon [0];
           }
         }

        // r1 in [0.1, 0.9] — we exclude degenerate combinations
        double r1 = u.RNDfromCI(0.1, 0.9);
        double r2 = 1.0 - r1;

        for(int c = 0; c < coords; c++)
         {
          double v = r1 * snap_c [ia].v [c] + r2 * snap_c [ib].v [c];

          // A perturbation with an adaptive amplitude based on the cooling schedule
          double range = rangeMax [c] - rangeMin [c];
          v += u.RNDfromCI(-1.0, 1.0) * range * perturbScale;

          if(v < rangeMin [c])
            v = rangeMin [c];
          if(v > rangeMax [c])
            v = rangeMax [c];
          combo [m].c [c] = v;
         }
        combo [m].sec_a = snap_sectorOf [ia];
        combo [m].sec_b = snap_sectorOf [ib];
       }

      // Diversity of combinations relative to the agent (center = agent's coordinates)
      double divAgg = 0.0;
      for(int c = 0; c < coords; c++)
       {
        double dj = 0.0;
        for(int m = 0; m < Ncomb; m++)
          dj += MathAbs(combo [m].c [c] - snap_c [agSlot].v [c]);
        dj /= Ncomb;
        divAgg += dj;
       }
      divAgg /= coords;

      if(divAgg > bestDiv)
       {
        bestDiv = divAgg;
        for(int m = 0; m < Ncomb; m++)
         {
          for(int c = 0; c < coords; c++)
            bestSet [m].c [c] = combo [m].c [c];
          bestSet [m].sec_a = combo [m].sec_a;
          bestSet [m].sec_b = combo [m].sec_b;
         }
       }
     } // End of the loop over templates

    // Fallback: no template found suitable material → perturbations around the agent
    if(bestDiv < 0)
     {
      for(int m = 0; m < Ncomb; m++)
       {
        for(int c = 0; c < coords; c++)
         {
          double range = rangeMax [c] - rangeMin [c];
          double v = snap_c [agSlot].v [c] + u.RNDfromCI(-1.0, 1.0) * range * 0.1;
          if(v < rangeMin [c])
            v = rangeMin [c];
          if(v > rangeMax [c])
            v = rangeMax [c];
          bestSet [m].c [c] = v;
         }
        bestSet [m].sec_a = anchorSec;
        bestSet [m].sec_b = anchorSec;
       }
     }

    // Add bestSet to the shared candidates list
    for(int m = 0; m < Ncomb; m++)
     {
      for(int c = 0; c < coords; c++)
        candidates [candPos].c [c] = bestSet [m].c [c];
      candidates [candPos].sec_a = bestSet [m].sec_a;
      candidates [candPos].sec_b = bestSet [m].sec_b;
      candidates [candPos].f     = -DBL_MAX;
      candPos++;
     }
   } // end of the loop over agents

//--- 4. Assign candidates to non-agent slots for FF evaluation 
// According to rule (ii) in the article, a new color is placed in the sector of one of its
// parents. First, we look for a slot with a matching sector; if it is occupied or
// no such slot exists, fall back to any free non-agent slot.

  bool isAgentSlot [];
  ArrayResize(isAgentSlot, popSize);
  ArrayInitialize(isAgentSlot, false);
  for(int g = 0; g < 10; g++)
    if(agentIdx [g] >= 0)
      isAgentSlot [agentIdx [g]] = true;

  bool slotTaken [];
  ArrayResize(slotTaken, popSize);
  ArrayInitialize(slotTaken, false);

  for(int m = 0; m < Kcand; m++)
   {
    int sec_a = candidates [m].sec_a;
    int sec_b = candidates [m].sec_b;
    int chosenSlot = -1;

    // 1) attempt: a non-agent slot with a matching sector
    for(int i = 0; i < popSize; i++)
     {
      if(isAgentSlot [i] || slotTaken [i])
        continue;
      if(snap_sectorOf [i] == sec_a || snap_sectorOf [i] == sec_b)
       {
        chosenSlot = i;
        break;
       }
     }

    // 2) fallback: any available non-agent slot
    if(chosenSlot < 0)
     {
      for(int i = 0; i < popSize; i++)
       {
        if(isAgentSlot [i] || slotTaken [i])
          continue;
        chosenSlot = i;
        break;
       }
     }

    if(chosenSlot < 0)
      break; // Just in case—all slots are taken

    slotTaken     [chosenSlot] = true;
    slotCandidate [chosenSlot] = m;
    isCandidate   [chosenSlot] = true;

    for(int c = 0; c < coords; c++)
      a [chosenSlot].c [c] = u.SeInDiSp(candidates [m].c [c],
                                        rangeMin  [c],
                                        rangeMax  [c],
                                        rangeStep [c]);
   }

//--- 5. Agent slots — no changes (a[i].c = snap_c[i]) 
  for(int g = 0; g < 10; g++)
   {
    int s = agentIdx [g];
    if(s < 0)
      continue;
    for(int c = 0; c < coords; c++)
      a [s].c [c] = u.SeInDiSp(snap_c [s].v [c],
                               rangeMin  [c],
                               rangeMax  [c],
                               rangeStep [c]);
   }
 }
//+------------------------------------------------------------------+

ConcentrationPhase_UpdateCircle — applying the FF evaluation results. First, the candidate fitness values are copied from a[i].f back to candidates[m].f via the slotCandidate[] links. Next, a greedy selection is performed: for each slot into which a candidate has been placed, the candidate's fitness is compared to the fitness in the snapshot. If a candidate is better, it remains in the slot, and the `candidateUsed[m]` flag is set. If the candidate is worse, the slot is restored from the snapshot. This is a simplification of the original hue circle rule, in which the cross-filtering in steps 4–12 allowed only a small fraction of candidates into the population; greedy selection gives all candidates an equal chance and significantly speeds up the algorithm's progress.

After selection, the agent slots are fully restored from the snapshot — this is necessary in case SeInDiSp in Moving slightly shifted the coordinates during discretization. All rejected candidates are stored in the circular buffer CMtempC[]. At the end, `UpdateAgentsByGroup` is called to potentially change the group leaders.

//+------------------------------------------------------------------+
//|         CONCENTRATION PHASE — hue circle update (greedy + sector)|
//+------------------------------------------------------------------+
void C_AO_CHA::ConcentrationPhase_UpdateCircle()
 {
  int Kcand = ArraySize(candidates);

//--- 1. Read the candidates' fitness values from their temporary slots in a[] 
  for(int i = 0; i < popSize; i++)
   {
    int m = slotCandidate [i];
    if(m >= 0)
      candidates [m].f = a [i].f;
   }

//--- 2. Greedy selection: for each slot with a candidate, compare
// its fitness with the fitness in snapshot. The better candidate is retained; the worse one is
// rolled back. This gives all 10·Ncomb candidates an equal chance of
// entry into the population (as opposed to strict crossover rules
// in the original article, which discard most candidates
// and leave the population virtually unchanged).
  for(int i = 0; i < popSize; i++)
   {
    int m = slotCandidate [i];
    if(m < 0)
      continue;

    if(candidates [m].f > snap_f [i])
     {
      // accept — the coordinates are already in a[i].c, and the fitness is already in a[i].f
      candidateUsed [m] = true;
      // sectorOf[i] remains unchanged (rule (ii) of the article)
     }
    else
     {
      // reject — roll back the slot
      for(int c = 0; c < coords; c++)
        a [i].c [c] = snap_c [i].v [c];
      a [i].f = snap_f [i];
     }
   }

//--- 3. Agent slots — restore from snapshot in case
// SeInDiSp in Moving slightly shifted the coordinates during discretization.
  for(int g = 0; g < 10; g++)
   {
    int s = agentIdx [g];
    if(s < 0)
      continue;
    for(int c = 0; c < coords; c++)
      a [s].c [c] = snap_c [s].v [c];
    a [s].f = snap_f [s];
   }

//--- 4. Rejected candidates → CMtemp 
  for(int m = 0; m < Kcand; m++)
   {
    if(candidateUsed [m])
      continue;

    int dst = CMtempCnt % CMtempCap;
    for(int c = 0; c < coords; c++)
      CMtempC [dst].v [c] = candidates [m].c [c];
    CMtempCnt++;
    if(CMtempCnt >= 2 * CMtempCap)
      CMtempCnt = CMtempCap;
   }

//--- 5. Update group agents (exchange sectors when 
//       a new leader emerges in the group)
  UpdateAgentsByGroup();
 }
//+------------------------------------------------------------------+
        

DispersionPhase_Generate — implements the Munsell formula for population dispersion. First, the center of the population is calculated; then, the Ns non-agents closest to the center are identified — these are the ones to be replaced. The source of material for replacement is selected as follows: if the main color memory (CM) is already full (CMcount ≥ Ksets), the most diverse set from it is taken; otherwise, if there are any accumulated candidates in CMtempC, a random sample is taken from them; otherwise, a fallback to uniform generation of points across the entire search space is used.

The selected source SCM is evaluated for diversity relative to the center of the search space. If the main memory is not full, the new set is simply added to it; if it is full, the system searches for an existing set with the lowest diversity — and replaces it with the new one if the new set has greater diversity.

The final step is to generate replacements using the formula X_new = rcm·X_SCM + (1 − rcm)·X_old, where the coefficient rcm is drawn uniformly from [rcm0, 1]. The rcm0 parameter (default 0.7) ensures that the new color is primarily determined by the point in memory, while maintaining some connection to the old position. Replacements are written directly into the slots a[sxIdx[k]].c[] using SeInDiSp for discretization.

//+------------------------------------------------------------------+
//|                    Dispersion phase — generation                 |
//+------------------------------------------------------------------+
void C_AO_CHA::DispersionPhase_Generate()
 {
//--- 1. Center of the population 
  double center [];
  PopulationCenter(center);

//--- 2. SX = the Ns non-agents closest to the center 
  double dist [];
  ArrayResize(dist, popSize);
  for(int i = 0; i < popSize; i++)
   {
    bool ag = false;
    for(int g = 0; g < 10; g++)
      if(agentIdx [g] == i)
       {
        ag = true;
        break;
       }
    if(ag)
     {
      dist [i] = DBL_MAX;
      continue;
     }

    double d = 0.0;
    for(int c = 0; c < coords; c++)
      d += MathAbs(a [i].c [c] - center [c]);
    dist [i] = d;
   }

  int sxIdx [];
  ArrayResize(sxIdx, Ns);
  for(int k = 0; k < Ns; k++)
   {
    double mn = DBL_MAX;
    int    mi = -1;
    for(int i = 0; i < popSize; i++)
      if(dist [i] < mn)
       {
        mn = dist [i];
        mi = i;
       }
    sxIdx [k] = mi;
    if(mi >= 0)
      dist [mi] = DBL_MAX;
   }

//--- 3. SCM source: CM (if full) or CMtemp 
  double centerSpace [];
  ArrayResize(centerSpace, coords);
  for(int c = 0; c < coords; c++)
    centerSpace [c] = (rangeMin [c] + rangeMax [c]) * 0.5;

  S_CHA_Vector scmPoints [];
  ArrayResize(scmPoints, Ns);
  for(int i = 0; i < Ns; i++)
    scmPoints [i].Init(coords);

  bool scmReady = false;

  if(CMcount >= Ksets)
   {
    // Choose the most diverse set from CM
    int    bi = 0;
    double bd = CM [0].diversity;
    for(int s = 1; s < CMcount; s++)
      if(CM [s].diversity > bd)
       {
        bd = CM [s].diversity;
        bi = s;
       }

    for(int i = 0; i < Ns; i++)
      for(int c = 0; c < coords; c++)
        scmPoints [i].v [c] = CM [bi].points [i].v [c];
    scmReady = true;
   }
  else
    if(CMtempCnt > 0)
     {
      int avail = MathMin(CMtempCnt, CMtempCap);
      if(avail >= Ns)
       {
        // Random selection of Ns points
        for(int k = 0; k < Ns; k++)
         {
          int idx = (int)u.RNDfromCI(0.0, (double)avail);
          if(idx >= avail)
            idx = avail - 1;
          for(int c = 0; c < coords; c++)
            scmPoints [k].v [c] = CMtempC [idx].v [c];
         }
       }
      else
       {
        for(int k = 0; k < Ns; k++)
         {
          int src = k % avail;
          for(int c = 0; c < coords; c++)
            scmPoints [k].v [c] = CMtempC [src].v [c];
         }
       }
      scmReady = true;
     }

  if(!scmReady)
   {
    // Fallback: random points throughout the entire search space
    for(int i = 0; i < Ns; i++)
      for(int c = 0; c < coords; c++)
        scmPoints [i].v [c] = u.RNDfromCI(rangeMin [c], rangeMax [c]);
   }

//--- 4. Save SCM in CM (step 2 of the dispersion phase) 
// Diversity of SCM relative to the center of the search space
  double scmDiv = 0.0;
  for(int c = 0; c < coords; c++)
   {
    double dj = 0.0;
    for(int i = 0; i < Ns; i++)
      dj += MathAbs(scmPoints [i].v [c] - centerSpace [c]);
    dj /= Ns;
    scmDiv += dj;
   }
  scmDiv /= coords;

  if(CMcount < Ksets)
   {
    int idx = CMcount;
    for(int i = 0; i < Ns; i++)
      for(int c = 0; c < coords; c++)
        CM [idx].points [i].v [c] = scmPoints [i].v [c];
    CM [idx].diversity = scmDiv;
    CMcount++;
   }
  else
   {
    int    mi = 0;
    double md = CM [0].diversity;
    for(int s = 1; s < CMcount; s++)
      if(CM [s].diversity < md)
       {
        md = CM [s].diversity;
        mi = s;
       }

    if(scmDiv > md)
     {
      for(int i = 0; i < Ns; i++)
        for(int c = 0; c < coords; c++)
          CM [mi].points [i].v [c] = scmPoints [i].v [c];
      CM [mi].diversity = scmDiv;
     }
   }

//--- 5. Generate new colors and write them to the SX slots 
  for(int k = 0; k < Ns; k++)
   {
    int slot = sxIdx [k];
    if(slot < 0)
      continue;
    double rcm = u.RNDfromCI(rcm0, 1.0);
    for(int c = 0; c < coords; c++)
     {
      double v = rcm * scmPoints [k].v [c] + (1.0 - rcm) * a [slot].c [c];
      a [slot].c [c] = u.SeInDiSp(v, rangeMin [c], rangeMax [c], rangeStep [c]);
     }
   }
 }
//+------------------------------------------------------------------+

DispersionPhase_UpdateCircle is the shortest method in the algorithm: it simply calls UpdateAgentsByGroup. According to the article, after the dispersion phase, “if replaced colors have better fitness value than their related agents, their positions on the hue circle will be exchanged” — and that is exactly what UpdateAgentsByGroup does: for each group, it checks whether a new leader has emerged and, if necessary, swaps the sectors of the new leader and the previous agent.

//+---------------------------------------------------------------------------+
//|       DISPERSION PHASE — updating the hue circle after the FF evaluation  |
//+---------------------------------------------------------------------------+
void C_AO_CHA::DispersionPhase_UpdateCircle()
 {
// According to the article: “if replaced colors have better fitness value than their
// related agents, their positions on the hue circle will be exchanged.”
// UpdateAgentsByGroup does exactly that (swapping sectors when the leader changes).
  UpdateAgentsByGroup();
 }
//+------------------------------------------------------------------+

Test Results

Testing was conducted using the standard C_AO benchmark: the Hilly, Forest, and Megacity functions in 10, 50, and 1,000 dimensions (5, 25, and 500 functions), with a budget of 10,000 FF evaluations and 10 repetitions. Default algorithm parameters. The CHA algorithm passes the anti-cheat test without any issues, scoring a comparable number of points at the lower dimensionalities of the main test.

CHA|Color Harmony Algorithm|50.0|4.0|10.0|4.0|0.5|0.01|0.7|
=============================
5 Hilly's; Func runs: 10000; result: 0.5861221565962496
25 Hilly's; Func runs: 10000; result: 0.46909923573046725
500 Hilly's; Func runs: 10000; result: 0.30514641914713303
=============================
5 Forest's; Func runs: 10000; result: 0.4702465987009819
25 Forest's; Func runs: 10000; result: 0.21050099785621573
500 Forest's; Func runs: 10000; result: 0.07964154815476658
=============================
5 Megacity's; Func runs: 10000; result: 0.42000000000000004
25 Megacity's; Func runs: 10000; result: 0.2437333333333333
500 Megacity's; Func runs: 10000; result: 0.12732000000000107
=============================
All score: 2.91181 (32.35%)

Visualization of the algorithm running on the test bench and on selected functions.

CHA on the Hilly test function

CHA on the Forest test function

CHA on the Megacity test function

CHA on the Ackley function

CHA on the Paraboloid function

The CHA algorithm is presented in the ranking table for informational purposes only.

cc	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	ANS	across neighbourhood search	1.00000	0.88228	0.40138	2.28366	1.00000	0.95281	0.28092	2.23373	0.94667	0.85733	0.22389	2.02789	6.545	72.72
2	AMOm	animal migration optimization M	0.91624	0.83603	0.46790	2.22017	0.98482	0.92010	0.36391	2.26883	0.91733	0.81707	0.25177	1.98617	6.475	71.94
3	CLA	code lock algorithm (joo)	0.95139	0.86199	0.37879	2.19217	0.99349	0.93500	0.26497	2.19346	0.93600	0.84267	0.24060	2.01927	6.405	71.17
4	(P+O)ES	(P+O) evolution strategies	0.86571	0.89539	0.39740	2.15850	0.97761	0.89820	0.26878	2.14459	0.92133	0.80240	0.23952	1.96325	6.266	69.62
5	SDSm	stochastic diffusion search M	0.95195	0.84944	0.36249	2.16388	0.98061	0.88457	0.22112	2.08630	0.92267	0.79013	0.21380	1.92660	6.177	68.63
6	AAm	archery algorithm M	0.84685	0.73320	0.42590	2.00595	0.96709	0.77837	0.27789	2.02335	0.86133	0.77707	0.28712	1.92552	5.955	66.17
7	SIA	simulated isotropic annealing (joo)	0.93543	0.86504	0.38483	2.18530	0.94069	0.80609	0.23835	1.98513	0.86400	0.66160	0.19536	1.72096	5.891	65.46
8	TETA	time evolution travel algorithm (joo)	0.91452	0.86369	0.25579	2.03400	0.99654	0.91291	0.14394	2.05339	0.85467	0.82213	0.10443	1.78123	5.869	65.21
9	ESG	evolution of social groups (joo)	0.98111	0.79857	0.31167	2.09135	0.98954	0.82270	0.15032	1.96256	0.92133	0.73440	0.15315	1.80888	5.863	65.14
10	CTA	comet tail algorithm (joo)	0.92435	0.86786	0.27838	2.07059	0.99039	0.84571	0.19448	2.03058	0.95467	0.69680	0.11008	1.76155	5.863	65.14
11	ECBO	enhanced colliding bodies optimization	0.94024	0.72363	0.32356	1.98743	0.99477	0.80291	0.13056	1.92824	0.87600	0.70160	0.17433	1.75193	5.668	62.98
12	DA	dialectical algorithm	0.93117	0.75400	0.26205	1.94722	0.98925	0.81375	0.08662	1.88962	0.92667	0.68107	0.11315	1.72089	5.558	61.76
13	BBO	biogeography-based optimization	0.95876	0.70609	0.35752	2.02237	0.92981	0.70660	0.16970	1.80611	0.87467	0.63013	0.20813	1.71293	5.541	61.57
14	BHAm	black hole algorithm M	0.79558	0.76207	0.34682	1.90447	0.99836	0.75798	0.13826	1.89460	0.85067	0.64427	0.17020	1.66514	5.464	60.71
15	HS	harmony search	0.91420	0.69049	0.29924	1.90393	0.97627	0.73373	0.14193	1.85193	0.91733	0.62720	0.15364	1.69817	5.454	60.60
16	RFO	royal flush optimization (joo)	0.80989	0.74481	0.34546	1.90016	0.95251	0.77926	0.15185	1.88362	0.80400	0.66427	0.19071	1.65898	5.443	60.48
17	BOAm	billiards optimization algorithm M	0.76177	0.72421	0.25275	1.73873	0.90890	0.81960	0.28853	2.01703	0.83733	0.74613	0.09763	1.68109	5.437	60.41
18	ASO	anarchy society optimization	0.73070	0.73713	0.31195	1.77978	0.99732	0.87700	0.17619	2.05051	0.72000	0.68773	0.18988	1.59761	5.428	60.31
19	EOm	Extremal Optimization_M	0.76527	0.75205	0.31908	1.83640	0.99999	0.76426	0.12437	1.88862	0.84133	0.64133	0.15247	1.63513	5.360	59.56
20	ACS	artificial cooperative search	0.75545	0.77162	0.31653	1.84360	1.00000	0.80488	0.10705	1.91193	0.76933	0.60800	0.14157	1.51890	5.274	58.60
21	SSG	saplings sowing and growing	0.75436	0.63206	0.35935	1.74577	0.91907	0.69694	0.19755	1.81356	0.81867	0.60533	0.21347	1.63747	5.197	57.74
22	AOSm	atomic orbital search M	0.76184	0.68435	0.31344	1.75963	0.90015	0.80044	0.11501	1.81560	0.82800	0.63280	0.15696	1.61776	5.193	57.70
23	TSEA	turtle shell evolution algorithm (joo)	0.95809	0.64852	0.29571	1.90232	0.99522	0.58104	0.10542	1.68168	0.92133	0.52160	0.14567	1.58860	5.173	57.48
24	DE	differential evolution	0.96398	0.62346	0.26089	1.84833	0.98482	0.77018	0.11459	1.86959	0.93067	0.36213	0.11000	1.40280	5.121	56.90
25	BIO	blood inheritance optimization (joo)	0.72580	0.66522	0.31228	1.70330	0.99995	0.68125	0.11540	1.79660	0.85467	0.59333	0.15364	1.60164	5.102	56.69
26	(PO)ES	(PO) evolution strategies	0.73972	0.58190	0.38896	1.71058	0.91199	0.59975	0.21262	1.72436	0.82400	0.56240	0.23432	1.62072	5.056	56.18
27	BO	bonobo optimizer	0.75555	0.64366	0.32657	1.72578	0.94332	0.70442	0.13999	1.78773	0.73467	0.61440	0.16728	1.51635	5.030	55.89
28	SRA	successful restaurateur algorithm (joo)	0.89010	0.63359	0.29115	1.81484	0.96634	0.55285	0.08914	1.60833	0.89333	0.52800	0.13911	1.56044	4.984	55.38
29	CRO	chemical reaction optimisation	0.91281	0.65681	0.29866	1.86828	0.90513	0.56020	0.10939	1.57472	0.82800	0.50133	0.14149	1.47082	4.914	54.60
30	BCOm	bacterial chemotaxis optimization M	0.82589	0.61733	0.31584	1.75906	0.95296	0.63718	0.11984	1.70998	0.76533	0.51653	0.15800	1.43986	4.909	54.54
31	DOA	dream optimization algorithm	0.78522	0.78121	0.36036	1.92679	0.61584	0.42117	0.12254	1.15955	0.86667	0.72587	0.21127	1.80381	4.890	54.33
32	ABO	african buffalo optimization	0.92295	0.62528	0.29885	1.84708	0.92992	0.57468	0.09372	1.59832	0.73333	0.51333	0.14324	1.38990	4.835	53.72
33	BSA	bird swarm algorithm	0.94432	0.67941	0.26401	1.88774	0.91649	0.65619	0.12054	1.69322	0.80933	0.33547	0.10652	1.25132	4.832	53.69
34	TSm	tabu search M	0.87806	0.61040	0.28993	1.77839	0.98116	0.52165	0.08544	1.58825	0.82667	0.49547	0.13552	1.45766	4.824	53.60
35	BSA	backtracking search algorithm	0.87128	0.53190	0.28675	1.68993	0.92408	0.51602	0.09153	1.53163	0.96000	0.47253	0.13760	1.57013	4.792	53.24
36	WOAm	whale optimization algorithm M	0.93893	0.59477	0.26695	1.80065	0.98036	0.53873	0.07112	1.59021	0.78667	0.47600	0.11892	1.38159	4.772	53.02
37	CSO	competitive swarm optimizer	0.85151	0.60786	0.29896	1.75833	0.84085	0.58491	0.11974	1.54550	0.80000	0.48560	0.14184	1.42744	4.731	52.57
38	FBA	fractal-based algorithm	0.69419	0.64267	0.28955	1.62641	0.99812	0.54905	0.08705	1.63422	0.76133	0.51253	0.13689	1.41075	4.671	51.90
39	ECOi	eco-inspired evolutionary algorithm	0.78817	0.54402	0.29360	1.62579	0.88996	0.46592	0.09747	1.45335	0.78533	0.45173	0.14295	1.38001	4.459	49.54
40	BSO	brain storm optimization	0.92207	0.57625	0.29732	1.79564	0.80764	0.42508	0.09448	1.32720	0.77200	0.36533	0.13065	1.26798	4.391	48.79
41	CAm	camel algorithm M	0.71534	0.56917	0.35985	1.64436	0.84094	0.47174	0.10850	1.42118	0.70400	0.41947	0.19563	1.31910	4.385	48.72
42	ACOm	ant colony optimization_M	0.71885	0.48410	0.30990	1.51285	0.75792	0.48639	0.11871	1.36302	0.83600	0.48667	0.16148	1.48415	4.360	48.44
43	BSO	brain storm optimization	0.91331	0.55813	0.29705	1.76849	0.85329	0.44038	0.09447	1.38814	0.72267	0.32880	0.13325	1.18472	4.341	48.23
44	AEFA	artificial electric field algorithm	0.88798	0.65769	0.25276	1.79843	0.95550	0.63602	0.03946	1.63098	0.60800	0.16000	0.09845	0.86645	4.296	47.73
45	SOA	simple optimization algorithm	0.91664	0.47040	0.27095	1.65799	0.88264	0.31546	0.06044	1.25854	0.82933	0.31733	0.11579	1.26245	4.179	46.43
	CHA	color_harmony_algorithm	0.58612	0.46909	0.30514	1.36035	0.47024	0.21050	0.07964	0.76038	0.42000	0.24373	0.12732	0.79105	2.912	32.35
	RW	random walk	0.49970	0.32333	0.25791	1.08094	0.30754	0.11470	0.04400	0.46624	0.36133	0.17013	0.10244	0.63390	2.181	24.23

Conclusions

CHA is an example of how an appealing metaphor and the mathematical nature of an algorithm can diverge, and this becomes apparent during the implementation process. The structure of the hue circle and the rules for harmonious color combinations seem meaningful, but the basic crossover operator (linear mixture) reduces the entire algorithm to movement within the convex hull of the initial population. Without an explicit mechanism to move offspring outside this convex hull, CHA cannot optimize functions whose optimum lies outside the initial population cloud. We added this mechanism — a Gaussian perturbation with a cooling schedule — and the algorithm started working. But it is important to understand that this is not part of the original formulation.

While fine-tuning the implementation to make it fully functional, we had to deviate from the exact description in the article in several places. In the original, the rule for updating the hue circle involves fifteen steps with cross-filtering of candidates between groups. In practice, these filters sort out most candidates in the initial phase and significantly slow down the population's movement. Step 6 (“the agent in the other parent’s group must be better than the agent in the current slot’s group”) has been removed entirely; step 8 (“the group with the minimum average fitness”) has been replaced with simply selecting the candidate with the maximum fitness from among the eligible ones.

Steps 4–12 of the original circle rules establish a global mapping between candidates and slots. This results in approximately 5 to 10 accepted candidates out of 40 — three-quarters of the FF evaluations are wasted. This is replaced with a greedy approach: each slot compares the candidate placed in it with its snapshot taken before Moving; the better of the two wins. All 40 candidates have a chance of being selected for the population. The main extension is an adaptively varying perturbation amplitude.

The algorithm’s strengths include a structured population (the geometry of the hue circle prevents premature clustering), a clear distinction between the exploitation and exploration phases based on an intuitive diversity criterion, and the ability to recall previously found solutions for use in restarting the process. These properties make CHA a suitable choice for problems where the optimum is known to lie within a relatively compact region of the search space — typical problems include model calibration, parameter tuning for systems already in operation, and local optimization around a known approximation.

Weaknesses include low efficiency on functions with a highly uneven fitness landscape. For high-dimensional global optimization problems with irregular objective functions, it is advisable to choose algorithms from the top of the ranking table.

CHA’s final score on our benchmark is 32.35%, which places the algorithm outside our ranking. This is a fair result for a method whose main value lies not in performance, but in the conceptual freshness of its approach and in reminding us to check the mathematical soundness behind every aesthetically appealing metaphor. The complete source code for the algorithm and the test script is available in the attachment to the article.

Figure 4. Color coding of algorithms by test

Figure 5. A 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 Algorithm

Pros:

1. It solves simple problems.

Cons:

1. A large number of external parameters.
2. Tendency to stagnate.
3. The algorithm is very difficult to configure and debug.

An archive containing the latest versions of the algorithm source code is attached to the article. The author of this article assumes no responsibility for the absolute accuracy of the descriptions of canonical algorithms; many of these algorithms 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

Files used in this 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	TestStand3D.mqh	Include file	3D visualization panel for the test bench
6	Utilities.mqh	Include file	Utility functions library
7	CalculationTestResults.mqh	Include file	Script for calculating results for the comparison table
8	Test_AO_All.mq5	Script	Unified test bench for all population-based optimization algorithms
9	Test_AO_AntiCheat	Script	Test for cheating in optimization algorithms
10	Simple use of population optimization algorithms.mq5	Script	Simple example of using population-based optimization algorithms without visualization
11	Test_AO_CHA.mq5	Script	Test bench for CHA