icaros-usc / pyribs

A bare-bones Python library for quality diversity optimization.
https://pyribs.org
MIT License
205 stars 31 forks source link

[BUG or Rather Improvement?] #392

Open lunjohnzhang opened 9 months ago

lunjohnzhang commented 9 months ago

Description

I am trying to use EvolutionStrategyEmitter with CMAEvolutionStrategy to generate high dimensional solutions (i.e. > 3000) where all dimensions are within the bounds [1e-3, None] (i.e. larger than a small positive number). However, since the current implementation of CMAEvolutionStrategy simply keeps resampling until all solutions are within the bounds, it ends up iterating forever with high dimensional solutions (in theory it'll eventually sample solutions within the bounds but it's taking too long).

Would it be possible to change the implementation such that it only resamples for a fixed number of iterations and clip the solution within the bounds after that?

Steps to Reproduce

The following code reproduce the above issue:

import fire
import numpy as np

from threadpoolctl import threadpool_limits
from ribs.emitters import EvolutionStrategyEmitter
from ribs.emitters.opt import CMAEvolutionStrategy
from ribs.archives import GridArchive
from ribs.schedulers import Scheduler
from ribs._utils import readonly

class DebugCMAEvolutionStrategy(CMAEvolutionStrategy):
    # Limit OpenBLAS to single thread. This is typically faster than
    # multithreading because our data is too small.
    @threadpool_limits.wrap(limits=1, user_api="blas")
    def ask(self, lower_bounds, upper_bounds, batch_size=None):
        """Samples new solutions from the Gaussian distribution.

        Args:
            lower_bounds (float or np.ndarray): scalar or (solution_dim,) array
                indicating lower bounds of the solution space. Scalars specify
                the same bound for the entire space, while arrays specify a
                bound for each dimension. Pass -np.inf in the array or scalar to
                indicated unbounded space.
            upper_bounds (float or np.ndarray): Same as above, but for upper
                bounds (and pass np.inf instead of -np.inf).
            batch_size (int): batch size of the sample. Defaults to
                ``self.batch_size``.
        """
        if batch_size is None:
            batch_size = self.batch_size

        self._solutions = np.empty((batch_size, self.solution_dim),
                                   dtype=self.dtype)
        self.cov.update_eigensystem(self.current_eval, self.lazy_gap_evals)
        transform_mat = self.cov.eigenbasis * np.sqrt(self.cov.eigenvalues)

        # Resampling method for bound constraints -> sample new solutions until
        # all solutions are within bounds.
        remaining_indices = np.arange(batch_size)
        while len(remaining_indices) > 0:
            unscaled_params = self._rng.normal(
                0.0,
                self.sigma,
                (len(remaining_indices), self.solution_dim),
            ).astype(self.dtype)
            new_solutions, out_of_bounds = self._transform_and_check_sol(
                unscaled_params, transform_mat, self.mean, lower_bounds,
                upper_bounds)
            self._solutions[remaining_indices] = new_solutions
            print("new sol: ", new_solutions)

            # Find indices in remaining_indices that are still out of bounds
            # (out_of_bounds indicates whether each value in each solution is
            # out of bounds).
            remaining_indices = remaining_indices[np.any(out_of_bounds,
                                                         axis=1)]
            print("Remaining indices:", remaining_indices)
        return readonly(self._solutions)

def infinite_sampling(sol_size, seed):
    archive = GridArchive(
        solution_dim=sol_size,
        dims=[100, 50],
        ranges=[[10, 30], [0, 10]],
        seed=seed,
        dtype=np.float32,
    )

    bounds = [(1e-3, None) for _ in range(sol_size)]

    initial_solution = np.ones(sol_size)
    emitters = [
        EvolutionStrategyEmitter(
            archive,
            es=DebugCMAEvolutionStrategy,
            x0=initial_solution,
            sigma0=1,
            seed=seed,
            bounds=bounds,
            batch_size=10,
        )
    ]

    scheduler = Scheduler(
        archive,
        emitters,
    )

    sols = scheduler.ask()
    print(sols)

    return emitters, scheduler

if __name__ == "__main__":
    fire.Fire(infinite_sampling)

To reproduce the issue, run with python <script_name> 3000 42.

btjanaka commented 9 months ago

Hi @lunjohnzhang, thank you for the suggestions! I can confirm that I am able to reproduce the behavior. This seems somewhere between a bug report and a feature request.

The new behavior seems reasonable, but I'm not sure what the API should be for it. Perhaps there could be a bounds_iters parameter to CMAEvolutionStrategy.__init__ that is used in ask as part of the behavior you describe? bounds_iters could default to None to indicate the behavior should not be activated. If someone wants the behavior, they can pass bounds_iters in the es_kwargs of EvolutionStrategyEmitter.

Would you be willing to write a PR for this? No worries if not; I can always circle back to it later.

lunjohnzhang commented 9 months ago

Thanks for the reply! The solution with bounds_iters makes sense to me and I will be happy to write a PR for it. I guess you would want this feature in other ES as well?

btjanaka commented 9 months ago

Yes, that would be great. Thanks!