Implement CMA-MAE archive thresholds

[itsdawei]

Description

In this PR, we modify archive_base to implement archive thresholds from CMA-MAE.

In addition to that, we also revive old implementation of add_single, which was removed in (#221) in favor of calling add with an array of a single element. We found it necessary to revive the old implementation so users of pyribs can easily "hack" the add_single function to use their own implementation. Otherwise, the users will either have to come up with their own implementation of add_single, or understand (and "hack") the much more complicated add function.

In #221 (comment), it was mentioned that add_batch without numba gives roughly 2x speed increase over add_single with numba. Since we did not revive the use of numba in add_single, add_batch now has a roughly 2.5x - 3x speed increase over add_single. For example, cma_me_imp takes ~100s on my computer with add_batch and ~300s with add_single.

TODO

• Revive old implementation of add_single in Implement batch addition in archives #221.
• Implement archive threshold single.
• Implement archive threshold batch.
• Verify CMA-MAE for add_single.
• Verify CMA-MAE for add.
• Verify CMA-MAEGA.

Testing

• Add testing for threshold updating.
• Fix the tests such that they call add_single and add with single solution

Can you also fix the tests such that they call add_single and add with single solution? We want to make sure we test both implementations, since we currently assume add_single calls add but that is no longer the case.

• Test for negative threshold_min.
• Re-run sphere.py test.
• add_single and add test with different learning rates.

Status

• I have read the guidelines in CONTRIBUTING.md
• I have formatted my code using yapf
• I have tested my code by running pytest
• I have linted my code with pylint
• I have added a one-line description of my change to the changelog in HISTORY.md
• This PR is ready to go

[itsdawei]

We ran python sphere.py <cma_me> 20 using add_mode = "single" for each of the implemented variant of CMA-ME. Our results are shown below:

Algorithm	New Implementation
cma_me_imp	92.09%, 17561542 QD
cma_me_imp_mu	93.67%, 17674516 QD
cma_me_rd	96.78%, 14260719 QD
cma_me_rd_mu	97.15%, 13483429 QD
cma_me_opt	10.04%, 2311070 QD
cma_me_mixed	97.22%, 18172492 QD

Since this is consistent with the results in #220 (comment), we verify that our new implementation of archive_base with learning_rate = 1.0 is equivalent to the CMA-ME algorithm, which is consistent with the theorems proven in Fontaine 2022.

[itsdawei]

We compare our implementation of CMA-MAE to that of the paper.

Experiments are ran with the following hyper-parameters: solution_dim = 100, archive_dim = 100, iterations=10000.

---

Our implementation (using add_mode="single")

Learning rate (alpha)	QD-Score	Coverage
0.000	5.5277	5.74%
0.001	62.7666	79.62%
0.010	64.9918	83.30%
0.100	60.5445	76.29%
1.000	37.7675	44.47%

Hence, we conclude that our implementation of CMA-MAE thresholds are consistent with that of the paper.

---

We further propose a method of batch adding to take advantage of the batching speed-up.

Our implementation (using add_mode="batch")

Learning rate (alpha)	QD-Score	Coverage
0.000	6.1481	6.44%
0.001	61.596	77.96%
0.010	63.3804	80.96%
0.100	61.7953	78.30%
1.000	39.1403	46.21%

---

We also run the algorithm CMA-MAEGA, achieving a QD score of 75.378 and coverage of 100%. This is consistent with the results in the paper.

[itsdawei]

Bug Expected behavior in original code:

if not already_occupied:
    threshold_values[new_index] = threshold_floor

# Update the threshold for this cell by the learning rate.
old_threshold = threshold_values[new_index]

if (not already_occupied or old_thresh < new_objective_value or 
    objective_values[new_index] < new_objective_value):
  
    # ...
  
    # Update the threshold if the threshold will increase.
    if old_thresh < new_objective_value:
        # ...
    
    # ...
 
# ...
value = objective_value - old_threshold

A problem with the orignial code is that threshold_floor defaults to 0. If the objective function have a negative codomain, e.g. $f(x) = -x^2-1$. Then, the algorithm breaks since no solutions will have objective greater than zero, i.e. no solutions will be added to the archive. For this reason, in this PR, we default threshold_floor to -np.inf. This way, when the learning_rate and threshold_floor are not set, the archive will not block any negative objective values.

How it affects us?

However, the code above will break as value will resolve to np.inf when using default setting of the archive. Hence, we should use a different old_threshold to compute the value when the cell is not already occupied.

[btjanaka]

To determine performance loss, I tested CMA-ME with improvement ranking both on master and on this branch (python examples/sphere.py cma_me_imp).

Master: 76.94 sec
This branch: 78.24 sec

Thus, we did not lose any significant performance in terms of wallclock time with this PR.

It's a bit harder to tell if CMA-MAE is doing alright in terms of time -- it's going to take twice as long just because it is inserting into two archives, and there are also threshold computations which occupy some time. Furthermore, CMA-MAE runs with different hyperparameters than CMA-ME does.

Nevertheless, I'll note that running python examples/sphere.py cma_mae gives a runtime of 319.91 sec.