Adding CMA-ES (both inside MVRSM and separate) and Hospital simulator problem

Author: lbliek

expensiveoptimbenchmark/problems/DockerHospitalBenchmark.py

@@ -0,0 +1,62 @@
+from .base import BaseProblem
+
+import subprocess
+import platform
+import numpy as np
+import os
+
+class DockerHospitalBenchmarkProblem(BaseProblem):
+
+    def __init__(self, name, d, lbs, ubs, vartype, direction, errval):
+        self.name = name
+        self.d = d
+        self.lb = np.asarray(lbs)
+        self.ub = np.asarray(ubs)
+        self.vt = np.asarray(vartype)
+        self.direction = 1 if direction == "min" else -1
+        self.errval = errval
+
+        if os.path.exists("./evaluate.sh"):
+            self.evalCommand = ["./evaluate.sh"]
+        elif os.path.exists("/evaluate.sh"):
+            self.evalCommand = ["/evaluate.sh"]
+        else:
+            # Note: experiment should be ran with sudo under linux to allow this command to work,
+            # or under an user in the docker group.
+            self.evalCommand = ["docker", "run", "--rm", "mrebolle/r-geccoc:Track1", "-c"]
+    
+    def evaluate(self, xs):
+        concatenatedCandidate = ",".join(["%.8f" % x if xvt == 'cont' else "%i" % x for (x, xvt) in zip(xs, self.vt)])
+        parsedCandidate = f"Rscript objfun.R \"{concatenatedCandidate}\""
+        cmd = f"{self.evalCommand + [parsedCandidate]}"
+        print(f"Running '{cmd}'")
+        # res = subprocess.check_output(cmd, shell=True)
+        res = subprocess.check_output(self.evalCommand + [parsedCandidate])
+        reslast = res.strip().split(b"\n")[-1]
+        print(f"Result: {res}. Objective: {reslast}")
+        try:
+            return self.direction * float(reslast)
+        except:
+            return self.errval
+
+    def lbs(self):
+        return self.lb
+
+    def ubs(self):
+        return self.ub
+
+    def vartype(self):
+        return self.vt
+
+    def dims(self):
+        return self.d
+
+    def __str__(self):
+        return f"Hospital()"
+
+Hospital = DockerHospitalBenchmarkProblem("hospital", 29, 
+    np.array([ 6,  7, 3, 3, 3, 5, 3, 3, 25, 17, 2, 1, 0.25, 0.05, 0.07, 0.005, 0.07, 1e-04, 0.08, 0.25, 0.08, 0.5, 1e-6, 2, 1e-6, 1e-6, 1, 2,  0.5 ]), 
+    np.array([14, 13, 7, 9, 7, 9, 5, 7, 35, 25, 5, 7, 2   , 0.15, 0.11, 0.02 , 0.13, 0.002, 0.12, 0.35, 0.12, 0.9, 0.01, 4, 1.1, 0.0625, 2, 5, 0.75]), 
+    ['cont'] * 29, "min", 10000.0)
+
+dockerhospitalsimbenches = [Hospital]

expensiveoptimbenchmark/run_experiment.py

@@ -292,6 +292,16 @@ def generate_synthetic(args):
     for dockerproblem in dockersimbenches
 })
 
+from problems.DockerHospitalBenchmark import dockerhospitalsimbenches
+problems.update({
+    dockerproblem.name.lower(): {
+        'args': set(),
+        'defaults': {},
+        'constructor': generate_construct_synthetic(dockerproblem)
+    }
+    for dockerproblem in dockerhospitalsimbenches
+})
+
 def nop(*x, **y):
     pass
 
@@ -325,6 +335,7 @@ def execute_IDONE(params, problem, max_eval, log):
     assert rand_evals >= 1, "IDONE requires at least one initial random evaluation."
     return optimize_IDONE(problem, max_eval, rand_evals=rand_evals, model=type_model, binarize_categorical=binarize_categorical, binarize_int=binarize_int, sampling=sampling, enable_scaling=enable_scaling, log=log, exploration_prob=expl_prob, idone_log=idone_log)
 
+## MVRSM
 def execute_MVRSM(params, problem, max_eval, log):
     from solvers.MVRSM.wMVRSM import optimize_MVRSM
     if params['--model'] not in ['basic', 'advanced']:
@@ -333,21 +344,27 @@ def execute_MVRSM(params, problem, max_eval, log):
         raise ValueError("--binarize-categorical should be a boolean.")
     if params['--scaling'] not in ['true', 't', 'yes', 'y', 'false', 'f', 'no', 'n']:
         raise ValueError("--scaling should be a boolean.")
+    if params['--optimizer'] not in ['adam', 'L-BFGS-B', 'CMA-ES']:
+        raise ValueError("Valid optimizers are `adam` and 'L-BFGS-B'")
 
     type_model = params['--model']
     binarize_categorical = params['--binarize-categorical'] in ['true','t', 'yes', 'y']
     enable_scaling = params['--scaling'] in ['true','t', 'yes', 'y']
     rand_evals = int(params['--rand-evals']) - 1
+    optimizer = params.get('--optimizer')
+    bound_h = params['--boundary-handler']
     assert rand_evals >= 0, "MVRSM requires at least one initial random evaluation."
 
-    return optimize_MVRSM(problem, max_eval, rand_evals=rand_evals, model=type_model, binarize_categorical=binarize_categorical, enable_scaling=enable_scaling, log=log)
+    return optimize_MVRSM(problem, max_eval, rand_evals=rand_evals, model=type_model, binarize_categorical=binarize_categorical, enable_scaling=enable_scaling, optimizer=optimizer, bound_h=bound_h, log=log)
 
-# SA
+## SA
 def execute_SA(params, problem, max_eval, log):
     from solvers.SA.wSA import optimize_SA
 
     return optimize_SA(problem, max_eval, log=log)
 
+
+## DONE
 def check_DONEjl():
     from solvers.DONEjl.wDONEjl import minimize_DONEjl
     import numpy as np
@@ -376,7 +393,7 @@ def execute_DONEjl(params, problem, max_eval, log):
 
     return optimize_DONEjl(problem, rand_evals, max_eval, hyperparams, log=log)
 
-# Hyperopt TPE
+## Hyperopt TPE
 def execute_hyperopt(params, problem, max_eval, log):
     from solvers.hyperopt.whyperopt import optimize_hyperopt_tpe
     rand_evals = int(params['--rand-evals'])
@@ -396,7 +413,7 @@ def execute_hyperopt_rnd(params, problem, max_eval, log):
 
     return optimize_hyperopt_rnd(problem, max_eval, cparams=conversion_params, log=log)
 
-# pyGPGO
+## pyGPGO
 def execute_pygpgo(params, problem, max_eval, log):
     from solvers.pyGPGO.wpyGPGO import optimize_pyGPGO
     from pyGPGO.covfunc import matern32
@@ -410,14 +427,14 @@ def execute_pygpgo(params, problem, max_eval, log):
     acq = Acquisition(mode='ExpectedImprovement')
     return optimize_pyGPGO(problem, max_eval, gp, acq, random_init_evals=rand_evals, log=log)
 
-# bayesian-optimization
+## bayesian-optimization
 def execute_bayesianoptimization(params, problem, max_eval, log):
     from solvers.bayesianoptimization.wbayesianoptimization import optimize_bayesian_optimization
     rand_evals = int(params['--rand-evals'])
     # TODO: Allow picking different configurations?
     return optimize_bayesian_optimization(problem, max_eval, random_init_evals=rand_evals, log=log)
 
-# smac
+## smac
 def execute_smac(params, problem, max_eval, log):
     from solvers.smac.wsmac import optimize_smac
     rand_evals = int(params['--rand-evals'])
@@ -428,12 +445,20 @@ def check_smac():
     from solvers.smac.wsmac import optimize_smac
     pass
 
-# CoCaBO
+## CoCaBO
 def execute_cocabo(params, problem, max_eval, log):
     from solvers.CoCaBO.wCoCaBo import optimize_CoCaBO
     rand_evals = int(params['--rand-evals'])
     return optimize_CoCaBO(problem, max_eval, init_points=rand_evals, log=log)
 
+## CMA
+def execute_cma(params, problem, max_eval, log):
+    from solvers.CMA.wCMA import optimize_CMA
+    bound_h = params['--boundary-handler']
+
+    return optimize_CMA(problem, max_eval, bound_h=bound_h, binarize_categorical=False, log=log)
+
+## Basinhopping
 def execute_scipy_basinhopping(params, problem, max_eval, log):
     from solvers.scipy.wscipy import optimize_basinhopping
     T = float(params['--T'])
@@ -442,6 +467,7 @@ def execute_scipy_basinhopping(params, problem, max_eval, log):
 
     return optimize_basinhopping(problem, max_eval, T=T, stepsize=stepsize, localmethod=method, log=log)
 
+## Local search
 def execute_scipy_local(params, problem, max_eval, log):
     from solvers.scipy.wscipy import optimize_scipy_local
     method = params['--method']
@@ -466,12 +492,14 @@ def execute_scipy_local(params, problem, max_eval, log):
         'check': nop
     },
     'mvrsm': {
-        'args': {'--model', '--binarize-categorical', '--rand-evals', '--scaling'},
+        'args': {'--model', '--binarize-categorical', '--rand-evals', '--scaling', '--optimizer', '--boundary-handler'},
         'defaults': {
             '--model': 'advanced',
             '--binarize-categorical': 'false',
             '--rand-evals': '5',
-            '--scaling': 'true'
+            '--scaling': 'true',
+            '--optimizer': 'L-BFGS-B',
+            '--boundary-handler': 'transform'
         },
         'executor': execute_MVRSM,
         'check': nop
@@ -544,6 +572,14 @@ def execute_scipy_local(params, problem, max_eval, log):
         'executor': execute_cocabo,
         'check': nop
     },
+    'cma': {
+        'args': {'--boundary-handler'},
+        'defaults': {
+            '--boundary-handler': 'penalty'
+        },
+        'executor': execute_cma,
+        'check': nop
+    },
     'scipy-basin': {
         'args': {'--T', '--stepsize', '--method'},
         'defaults': {
@@ -645,6 +681,7 @@ def execute_scipy_local(params, problem, max_eval, log):
     print(f" --binarize-categorical=<t|true|f|false> \t Whether to binarize categorical variables. (default: false)")
     print(f" --scaling=<t|true|f|false> \t Whether scaling is applied. (default: true)")
     print(f" --rand-evals=<int> \t Number of random evaluations. (default: 1)")
+    print(f" --optimizer=<adam|LBFGSB> Optimizer for the surrogate model (default: LBFGSB)")
     print()
     # HyperOpt
     print(f" hyperopt")

expensiveoptimbenchmark/solvers/CMA/wCMA.py

@@ -0,0 +1,74 @@
+import numpy as np
+
+from cma import BoundPenalty, BoundTransform, fmin, CMAOptions
+from ..utils import Monitor, Binarizer
+
+def optimize_CMA(problem, max_evals, bound_h='transform', binarize_categorical=False, log=None):
+    d = problem.dims()
+
+    vartypes = problem.vartype()
+
+    if binarize_categorical:
+        raise Exception("While sort-of implemented, this implementation has issues. Please do not enable binarization.")
+
+    # Bounds
+    lb = problem.lbs()
+    ub = problem.ubs()
+
+    # Generate initial point, round the integers.
+
+    if binarize_categorical:
+        b = Binarizer(problem.vartype()[perm] == 'cat', lb, ub)
+        x0 = b.binarize(x0)
+        # Get new bounds.
+        lb = b.lbs()
+        ub = b.ubs()
+        # All new variables introduced by the binarized are 'integer'.
+        # print(f"Binarization introduced {b.dout - b.din} new variables.")
+        vartypes = vartypes[self.origin_mapping]
+    
+    cmascale = (ub - lb)
+    x0 = (0.5*(ub-lb) + lb) / cmascale
+    not_continuous_variables = vartypes != "cont"
+    x0[not_continuous_variables] = np.round(x0[not_continuous_variables])
+    sigma0 = 1/4
+
+    mon = Monitor(f"CMAES{'/cbinarized' if binarize_categorical else ''}/{bound_h}", problem, log=log)
+    def f(x):
+        x = x * cmascale
+        if binarize_categorical:
+            xred = b.unbinarize(x)
+        else:
+            xred = x
+        x = x[vartypes != "cont"]
+        # print(xred)
+        mon.commit_start_eval()
+        r = problem.evaluate(xred)
+        mon.commit_end_eval(xred, r)
+        return r
+    
+    # lb, ub, max_evals
+    opts = CMAOptions()
+    opts['bounds'] = [list(lb / cmascale), list(ub / cmascale)]
+    opts['BoundaryHandler'] = BoundPenalty if bound_h == "penalty" else BoundTransform
+    opts['maxfevals'] = max_evals
+    opts['verbose'] = 0 # Supress printing!
+    opts['integer_variables'] = [i for i,t in enumerate(vartypes) if t != 'cont']
+    mon.start()
+    res = fmin(f, x0 / cmascale, sigma0, options=opts)
+    mon.end()
+    
+    # """
+    # - `res[0]` (`xopt`) -- best evaluated solution  
+    # - `res[1]` (`fopt`) -- respective function value  
+    # - `res[2]` (`evalsopt`) -- respective number of function evaluations  
+    # - `res[3]` (`evals`) -- number of overall conducted objective function evaluations  
+    # - `res[4]` (`iterations`) -- number of overall conducted iterations  
+    # - `res[5]` (`xmean`) -- mean of the final sample distribution  
+    # - `res[6]` (`stds`) -- effective stds of the final sample distribution  
+    # - `res[-3]` (`stop`) -- termination condition(s) in a dictionary  
+    # - `res[-2]` (`cmaes`) -- class `CMAEvolutionStrategy` instance  
+    # - `res[-1]` (`logger`) -- class `CMADataLogger` instance
+    # """
+
+    return mon.best_x, mon.best_fitness, mon