1"""GPOpt class."""
   2
   3# Authors: T. Moudiki
   4#
   5# License: BSD 3 Clause Clear
   6
   7import copy
   8import nnetsauce as ns
   9import numpy as np
  10import pandas as pd
  11import pickle
  12import shelve
  13from collections import namedtuple
  14from functools import partial
  15
  16try:
  17    from sklearn.utils.discovery import all_estimators
  18except:
  19    from sklearn.utils import all_estimators
  20from sklearn.base import RegressorMixin
  21from sklearn.gaussian_process import GaussianProcessRegressor
  22from sklearn.ensemble import RandomForestRegressor
  23from sklearn.gaussian_process.kernels import Matern
  24from sklearn.model_selection import GridSearchCV, RandomizedSearchCV
  25
  26import scipy.stats as st
  27from joblib import Parallel, delayed
  28from time import time
  29from tqdm import tqdm
  30from .config import REGRESSORS, REMOVED_REGRESSORS
  31from .utils import generate_sobol2
  32from .utils import Progbar
  33
  34
  35class GPOpt:
  36    """Class GPOpt.
  37
  38    # Arguments:
  39
  40        lower_bound: a numpy array;
  41            lower bound for researched minimum
  42
  43        upper_bound: a numpy array;
  44            upper bound for researched minimum
  45
  46        objective_func: a function;
  47            the objective function to be minimized
  48
  49        params_names: a list;
  50            names of the parameters of the objective function (optional)
  51
  52        surrogate_obj: an object;
  53            An ML model for estimating the uncertainty around the objective function
  54            Must be nnetsauce.CustomRegressor or nnetsauce.PredictionInterval
  55
  56        x_init:
  57            initial setting of points where `objective_func` is evaluated (optional)
  58
  59        y_init:
  60            initial setting values at points where `objective_func` is evaluated (optional)
  61
  62        n_init: an integer;
  63            number of points in the initial setting, when `x_init` and `y_init` are not provided
  64
  65        n_choices: an integer;
  66            number of points for the calculation of expected improvement
  67
  68        n_iter: an integer;
  69            number of iterations of the minimization algorithm
  70
  71        alpha: a float;
  72            Value added to the diagonal of the kernel matrix during fitting (for Matern 5/2 kernel)
  73
  74        n_restarts_optimizer: an integer;
  75            The number of restarts of the optimizer for finding the kernel’s parameters which maximize the log-marginal likelihood.
  76
  77        seed: an integer;
  78            reproducibility seed
  79
  80        save: a string;
  81            Specifies where to save the optimizer in its current state
  82
  83        n_jobs: an integer;
  84            number of jobs for parallel computing on initial setting or method `lazyoptimize` (can be -1)
  85
  86        acquisition: a string;
  87            acquisition function: "ei" (expected improvement) or "ucb" (upper confidence bound)
  88
  89        method: an str;
  90            "bayesian" (default) for Gaussian posteriors, "mc" for Monte Carlo posteriors,
  91            "splitconformal" (only for acquisition = "ucb") for conformalized surrogates
  92
  93        min_value: a float;
  94            minimum value of the objective function (default is None). For example,
  95            if objective function is accuracy, will be 1, and the algorithm will stop
  96
  97        per_second: a boolean;
  98            __experimental__, default is False (leave to default for now)
  99
 100        log_scale: a boolean;
 101            __experimental__, default is False (leave to default for now)
 102
 103    see also [Bayesian Optimization with GPopt](https://thierrymoudiki.github.io/blog/2021/04/16/python/misc/gpopt)
 104        or [Hyperparameters tuning with GPopt](https://thierrymoudiki.github.io/blog/2021/06/11/python/misc/hyperparam-tuning-gpopt) [Agnostic BayesOpt](https://thierrymoudiki.github.io/blog/2024/12/09/python/bayesconfoptim)
 105
 106    """
 107
 108    def __init__(
 109        self,
 110        lower_bound,
 111        upper_bound,
 112        objective_func=None,
 113        params_names=None,
 114        surrogate_obj=None,
 115        x_init=None,
 116        y_init=None,
 117        n_init=10,
 118        n_choices=100000,
 119        n_iter=190,
 120        alpha=1e-6,
 121        n_restarts_optimizer=25,
 122        seed=123,
 123        save=None,
 124        n_jobs=1,
 125        acquisition="ei",
 126        method="bayesian",
 127        min_value=None,
 128        per_second=False,  # /!\ very experimental
 129        log_scale=False,  # /!\ experimental
 130    ):
 131        n_dims = len(lower_bound)
 132
 133        assert n_dims == len(
 134            upper_bound
 135        ), "'upper_bound' and 'lower_bound' must have the same dimensions"
 136
 137        self.objective_func = objective_func
 138        self.params_names = params_names
 139        self.lower_bound = lower_bound
 140        self.upper_bound = upper_bound
 141        self.y_init = y_init
 142        self.log_scale = log_scale
 143        self.n_dims = n_dims
 144        self.n_init = n_init
 145        self.n_choices = n_choices
 146        self.n_iter = n_iter
 147        self.alpha = alpha
 148        self.n_restarts_optimizer = n_restarts_optimizer
 149        self.seed = seed
 150        self.save = save
 151        self.n_jobs = n_jobs  # for parallel case on initial design
 152        self.per_second = per_second
 153        self.x_min = None
 154        self.y_min = None
 155        self.y_mean = None
 156        self.y_std = None
 157        self.y_lower = None
 158        self.y_upper = None
 159        self.best_surrogate = None
 160        self.acquisition = acquisition
 161        self.min_value = min_value
 162        self.acq = np.array([])
 163        self.max_acq = []
 164        if surrogate_obj is None:
 165            self.surrogate_obj = GaussianProcessRegressor(
 166                kernel=Matern(nu=2.5),
 167                alpha=self.alpha,
 168                normalize_y=True,
 169                n_restarts_optimizer=self.n_restarts_optimizer,
 170                random_state=self.seed,
 171            )
 172        else:
 173            self.surrogate_obj = surrogate_obj
 174        assert method in (
 175            "bayesian",
 176            "mc",
 177            "splitconformal",
 178        ), "method must be in ('bayesian', 'mc', 'splitconformal')"
 179        self.method = method
 180        self.posterior_ = None
 181
 182        # Sobol seqs for initial design and choices
 183        sobol_seq_init = np.transpose(
 184            generate_sobol2(
 185                n_dims=self.n_dims,
 186                n_points=self.n_init,
 187                skip=2,
 188            )
 189        )
 190        sobol_seq_choices = np.transpose(
 191            generate_sobol2(
 192                n_dims=self.n_dims,
 193                n_points=self.n_choices,
 194                skip=self.n_init + 2,
 195            )
 196        )
 197
 198        # Sobol seqs for initial design and choices with bounds
 199        if self.log_scale == False:
 200            bounds_range = upper_bound - lower_bound
 201            self.x_init = (
 202                bounds_range * sobol_seq_init + lower_bound
 203                if x_init is None
 204                else x_init
 205            )
 206            self.x_choices = bounds_range * sobol_seq_choices + lower_bound
 207
 208        else:  # (!) experimental
 209            assert (
 210                lower_bound > 0
 211            ).all(), "all elements of `lower_bound` must be > 0"
 212            assert (
 213                upper_bound > 0
 214            ).all(), "all elements of `upper_bound` must be > 0"
 215
 216            log_lower_bound = np.log(lower_bound)
 217            log_upper_bound = np.log(upper_bound)
 218            log_bounds_range = log_upper_bound - log_lower_bound
 219            self.x_init = (
 220                np.minimum(
 221                    np.exp(log_bounds_range * sobol_seq_init + log_lower_bound),
 222                    1.7976931348623157e308,
 223                )
 224                if x_init is None
 225                else x_init
 226            )
 227            self.x_choices = np.minimum(
 228                np.exp(log_bounds_range * sobol_seq_choices + log_lower_bound),
 229                1.7976931348623157e308,
 230            )
 231
 232        # shelve for saving (not for loading)
 233        if self.save is not None:
 234            self.sh = shelve.open(filename=save, flag="c", writeback=True)
 235
 236        if self.per_second:
 237            self.timings = []
 238            self.rf_obj = RandomForestRegressor(
 239                n_estimators=250, random_state=self.seed
 240            )
 241
 242        self.n_jobs = n_jobs
 243
 244    # from sklearn.base
 245    def get_params(self):
 246        """Get object attributes.
 247
 248        Returns
 249        -------
 250        params : mapping of string to any
 251            Parameter names mapped to their values.
 252        """
 253        out = dict()
 254        param_names = dir(self)
 255        for key in param_names:
 256            if key.startswith("_") is False:
 257                out[key] = getattr(self, key, None)
 258
 259        return out
 260
 261    # for parallel case on initial design
 262    def eval_objective(self, arg):
 263        try:
 264            return self.objective_func(self.x_init[arg, :])
 265        except:
 266            return 1e06
 267
 268    # load data from stored shelve
 269    def load(self, path):
 270        """load data from stored shelve.
 271
 272        # Arguments
 273
 274        path : a string; path to stored shelve.
 275
 276        See also: [Bayesian Optimization with GPopt Part 2 (save and resume)](https://thierrymoudiki.github.io/blog/2021/04/30/python/misc/gpopt)
 277        """
 278
 279        self.sh = shelve.open(filename=path)
 280        for key, value in self.sh.items():
 281            setattr(self, key, value)
 282
 283    # update shelve in optimization loop
 284    def update_shelve(self):
 285        for key, value in self.get_params().items():
 286            if (callable(value) is False) & (key != "sh"):
 287                self.sh[key] = value
 288        self.sh.sync()
 289
 290    # closing shelve (can't be modified after)
 291    def close_shelve(self):
 292        """Close shelve.
 293
 294        # Arguments
 295
 296        No argument.
 297
 298        See also: [Bayesian Optimization with GPopt Part 2 (save and resume)](https://thierrymoudiki.github.io/blog/2021/04/30/python/misc/gpopt)
 299        """
 300
 301        self.sh.close()
 302
 303    # fit predict
 304    def surrogate_fit_predict(
 305        self,
 306        X_train,
 307        y_train,
 308        X_test,
 309        return_std=False,
 310        return_pi=False,
 311        param_search_init_design=False,
 312        param_distributions=None,
 313        **kwargs,
 314    ):
 315        if len(X_train.shape) == 1:
 316            X_train = X_train.reshape((-1, 1))
 317            X_test = X_test.reshape((-1, 1))
 318
 319        if (
 320            X_train.shape[0] <= self.n_init and param_search_init_design == True
 321        ):  # on initial design
 322            try:
 323                rs_obj = RandomizedSearchCV(
 324                    self.surrogate_obj,
 325                    param_distributions=param_distributions,
 326                    random_state=42,
 327                    cv=3,
 328                    **kwargs,
 329                )
 330                rs_obj.fit(X_train, y_train)
 331                self.surrogate_obj = rs_obj.best_estimator_
 332            except Exception as e:
 333                print(str(e))
 334
 335        # Get mean and standard deviation (+ lower and upper for not GPs)
 336        assert (
 337            return_std == True and return_pi == True
 338        ) == False, "must have either return_std == True or return_pi == True"
 339
 340        if return_std == True:
 341            self.posterior_ = "gaussian"
 342            self.surrogate_obj.fit(X_train, y_train)
 343            return self.surrogate_obj.predict(X_test, return_std=True)
 344
 345        elif (
 346            return_pi == True
 347        ):  # here, self.surrogate_obj must have `replications` not None
 348            if self.surrogate_obj.replications is not None:
 349                self.posterior_ = "mc"
 350                self.surrogate_obj.fit(X_train, y_train)
 351                try:  # it's a nnetsauce.CustomRegressor
 352                    res = self.surrogate_obj.predict(
 353                        X_test, return_pi=True, method="splitconformal"
 354                    )
 355                except Exception:  # it's a nnetsauce.PredictionInterval
 356                    try:
 357                        res = self.surrogate_obj.predict(X_test, return_pi=True)
 358                    except Exception as e:
 359                        print(
 360                            str(e)
 361                            + "Sureproof way is to encapsulate your surrogate in nnetsauce.PredictionInterval model"
 362                        )
 363
 364                self.y_sims = res.sims
 365                self.y_mean, self.y_std = (
 366                    np.mean(self.y_sims, axis=1),
 367                    np.std(self.y_sims, axis=1),
 368                )
 369                return self.y_mean, self.y_std, self.y_sims
 370
 371            else:  # self.surrogate_obj is conformalized (uses nnetsauce.PredictionInterval)
 372                assert (
 373                    self.acquisition == "ucb"
 374                ), "'acquisition' must be 'ucb' for conformalized surrogates"
 375                self.posterior_ = None
 376                self.surrogate_obj.fit(X_train, y_train)
 377                try:
 378                    res = self.surrogate_obj.predict(
 379                        X_test, return_pi=True, method="splitconformal"
 380                    )
 381                except Exception:
 382                    res = self.surrogate_obj.predict(X_test, return_pi=True)
 383                self.y_mean = res.mean
 384                self.y_lower = res.lower
 385                self.y_upper = res.upper
 386                return self.y_mean, self.y_lower, self.y_upper
 387
 388        else:
 389            raise NotImplementedError
 390
 391    # fit predict timings
 392    def timings_fit_predict(self, X_train, y_train, X_test):
 393        if len(X_train.shape) == 1:
 394            X_train = X_train.reshape((-1, 1))
 395            X_test = X_test.reshape((-1, 1))
 396
 397        # Get mean preds for timings
 398        return self.rf_obj.fit(X_train, y_train).predict(X_test)
 399
 400    # find next parameter by using expected improvement (ei)
 401    def next_parameter_by_acq(self, i, acq="ei"):
 402        if acq == "ei":
 403            if self.posterior_ == "gaussian":
 404                gamma_hat = (self.y_min - self.y_mean) / self.y_std
 405                self.acq = -self.y_std * (
 406                    gamma_hat * st.norm.cdf(gamma_hat) + st.norm.pdf(gamma_hat)
 407                )
 408            elif self.posterior_ == "mc":
 409                self.acq = -np.mean(
 410                    np.maximum(self.y_min - self.y_sims, 0), axis=1
 411                )
 412
 413        if acq == "ucb":
 414            if self.posterior_ == "gaussian":
 415                self.acq = self.y_mean - 1.96 * self.y_std
 416                self.ucb = self.y_mean + 1.96 * self.y_std
 417
 418            elif self.posterior_ is None:  # split conformal(ized) estimator
 419                self.acq = self.y_lower
 420                self.ucb = self.y_upper
 421
 422        # find max index -----
 423
 424        if self.per_second is False:
 425            # find index for max. ei
 426            max_index = self.acq.argmin()
 427
 428        else:  # self.per_second is True
 429            # predict timings on self.x_choices
 430            # train on X = self.parameters and y = self.timings
 431            # (must have same shape[0])
 432            timing_preds = self.timings_fit_predict(
 433                X_train=np.asarray(self.parameters),
 434                y_train=np.asarray(self.timings),
 435                X_test=self.x_choices,
 436            )
 437
 438            # find index for max. ei (and min. timings)
 439            max_index = (-self.acq / timing_preds).argmax()
 440
 441        self.max_acq.append(np.abs(self.acq[max_index]))
 442
 443        # Select next choice
 444        next_param = self.x_choices[max_index, :]
 445
 446        if next_param in np.asarray(self.parameters):
 447            if self.log_scale == False:
 448                np.random.seed(self.seed * i + 1000)
 449                next_param = (
 450                    self.upper_bound - self.lower_bound
 451                ) * np.random.rand(self.n_dims) + self.lower_bound
 452
 453            else:  # /!\ very... experimental
 454                np.random.seed(self.seed)
 455                log_upper_bound = np.log(self.upper_bound)
 456                log_lower_bound = np.log(self.lower_bound)
 457                log_bounds_range = log_upper_bound - log_lower_bound
 458
 459                next_param = np.minimum(
 460                    np.exp(
 461                        log_bounds_range * np.random.rand(self.n_dims)
 462                        + log_lower_bound
 463                    ),
 464                    1.7976931348623157e308,
 465                )
 466
 467        return next_param
 468
 469    # optimize the objective
 470    def optimize(
 471        self,
 472        verbose=1,
 473        n_more_iter=None,
 474        abs_tol=None,  # suggested 1e-4, for n_iter = 200
 475        ucb_tol=None,
 476        min_budget=50,  # minimum budget for early stopping
 477        func_args=None,
 478        param_search_init_design=False,
 479        param_distributions=None,
 480    ):
 481        """Launch optimization loop.
 482
 483        # Arguments:
 484
 485            verbose: an integer;
 486                verbose = 0: nothing is printed,
 487                verbose = 1: a progress bar is printed (longer than 0),
 488                verbose = 2: information about each iteration is printed (longer than 1)
 489
 490            n_more_iter: an integer;
 491                additional number of iterations for the optimizer (which has been run once)
 492
 493            abs_tol: a float;
 494                tolerance for convergence of the optimizer (early stopping based on acquisition function)
 495
 496            ucb_tol: a float;
 497                tolerance for convergence of the optimizer (early stopping based on length of prediction intervals)
 498                for UCB criterion
 499
 500            min_budget: an integer (default is 50);
 501                minimum number of iterations before early stopping controlled by `abs_tol`
 502
 503            func_args: a list;
 504                additional parameters for the objective function (if necessary)
 505
 506            param_search_init_design: a boolean;
 507                whether random search tuning must occur on the initial design or not
 508
 509            param_distributions: dict or list of dicts;
 510                Dictionary with parameters names (str) as keys and distributions or lists of
 511                parameters to try. Distributions must provide a rvs method for sampling
 512                (such as those from scipy.stats.distributions). If a list is given, it
 513                is sampled uniformly. If a list of dicts is given, first a dict is sampled
 514                uniformly, and then a parameter is sampled using that dict as above.
 515
 516        see also [Bayesian Optimization with GPopt](https://thierrymoudiki.github.io/blog/2021/04/16/python/misc/gpopt)
 517        and [Hyperparameters tuning with GPopt](https://thierrymoudiki.github.io/blog/2021/06/11/python/misc/hyperparam-tuning-gpopt)
 518
 519        """
 520
 521        # verbose = 0: nothing is printed
 522        # verbose = 1: a progress bar is printed (longer than 0)
 523        # verbose = 2: information about each iteration is printed (longer than 1)
 524        if func_args is None:
 525            func_args = []
 526
 527        if (
 528            n_more_iter is None
 529        ):  # initial optimization, before more iters are requested
 530            n_iter = self.n_iter
 531            # stopping iter for early stopping (default is total budget)
 532            iter_stop = n_iter  # potentially # got to check this
 533
 534            # initial design  ----------
 535
 536            if (verbose == 1) | (verbose == 2):
 537                print(f"\n Creating initial design... \n")
 538
 539            if verbose == 1:
 540                progbar = Progbar(target=self.n_init)
 541
 542            self.parameters = self.x_init.tolist()
 543            self.scores = []
 544
 545            if self.save is not None:
 546                self.update_shelve()
 547
 548            if self.y_init is None:  # calculate scores on initial design
 549                assert (
 550                    self.objective_func is not None
 551                ), "self.y_init is None: must have 'objective_func' not None"
 552
 553                if self.n_jobs == 1:
 554                    for i in range(self.n_init):
 555                        x_next = self.x_init[i, :]
 556
 557                        try:
 558                            if self.per_second is True:
 559                                start = time()
 560                                score = self.objective_func(x_next, *func_args)
 561                                if (np.isfinite(score) == False) or (
 562                                    np.isnan(score) == True
 563                                ):
 564                                    continue
 565                                self.timings.append(np.log(time() - start))
 566
 567                            else:  # self.per_second is False
 568                                score = self.objective_func(x_next, *func_args)
 569                                if (np.isfinite(score) == False) or (
 570                                    np.isnan(score) == True
 571                                ):
 572                                    continue
 573
 574                            self.scores.append(score)
 575
 576                            if self.save is not None:
 577                                self.update_shelve()
 578
 579                        except:
 580                            continue
 581
 582                        if verbose == 1:
 583                            progbar.update(i)  # update progress bar
 584
 585                        if verbose == 2:
 586                            print(f"point: {x_next}; score: {score}")
 587                    # end loop # calculate scores on initial design
 588
 589                    if verbose == 1:
 590                        progbar.update(self.n_init)
 591
 592                else:  # self.n_jobs != 1
 593                    assert (
 594                        self.per_second is False
 595                    ), "timings not calculated here"
 596
 597                    scores = Parallel(n_jobs=self.n_jobs, prefer="threads")(
 598                        delayed(self.objective_func)(self.x_init[i, :])
 599                        for i in range(self.n_init)
 600                    )
 601
 602                    self.scores = scores
 603
 604                    if self.save is not None:
 605                        self.update_shelve()
 606
 607            else:  # if self.y_init is not None:
 608                assert self.x_init.shape[0] == len(
 609                    self.y_init
 610                ), "must have: self.x_init.shape[0] == len(self.y_init)"
 611
 612                self.scores = pickle.loads(
 613                    pickle.dumps(self.y_init.tolist(), -1)
 614                )
 615
 616            # current best score on initial design
 617            min_index = (np.asarray(self.scores)).argmin()
 618            self.y_min = self.scores[min_index]
 619            self.x_min = self.x_init[min_index, :]
 620
 621            # current gp mean and std on initial design
 622            # /!\ if GP
 623            if param_search_init_design == False:
 624                if self.method == "bayesian":
 625                    self.posterior_ = "gaussian"
 626                    try:
 627                        y_mean, y_std = self.surrogate_fit_predict(
 628                            np.asarray(self.parameters),
 629                            np.asarray(self.scores),
 630                            self.x_choices,
 631                            return_std=True,
 632                            return_pi=False,
 633                        )
 634                    except ValueError:  # do not remove this
 635                        preds_with_std = self.surrogate_fit_predict(
 636                            np.asarray(self.parameters),
 637                            np.asarray(self.scores),
 638                            self.x_choices,
 639                            return_std=True,
 640                            return_pi=False,
 641                        )
 642                        y_mean, y_std = preds_with_std[0], preds_with_std[1]
 643                    self.y_mean = y_mean
 644                    self.y_std = np.maximum(2.220446049250313e-16, y_std)
 645
 646                elif self.method == "mc":
 647                    self.posterior_ = "mc"
 648                    assert self.surrogate_obj.__class__.__name__.startswith(
 649                        "CustomRegressor"
 650                    ) or self.surrogate_obj.__class__.__name__.startswith(
 651                        "PredictionInterval"
 652                    ), "for `method = 'mc'`, the surrogate must be a nnetsauce.CustomRegressor() or nnetsauce.PredictionInterval()"
 653                    assert (
 654                        self.surrogate_obj.replications is not None
 655                    ), "for `method = 'mc'`, the surrogate must be a nnetsauce.CustomRegressor() with a number of 'replications' provided"
 656                    preds_with_std = self.surrogate_fit_predict(
 657                        np.asarray(self.parameters),
 658                        np.asarray(self.scores),
 659                        self.x_choices,
 660                        return_std=False,
 661                        return_pi=True,
 662                    )
 663                    y_mean, y_std = preds_with_std[0], preds_with_std[1]
 664                    self.y_mean = y_mean
 665                    self.y_std = np.maximum(2.220446049250313e-16, y_std)
 666
 667                elif self.method == "splitconformal":
 668                    self.posterior_ = None
 669                    # assert self.surrogate_obj.__class__.__name__.startswith(
 670                    #    "PredictionInterval"
 671                    # ), "for `method = 'splitconformal'`, the surrogate must be a nnetsauce.PredictionInterval()"
 672                    preds_with_pi = self.surrogate_fit_predict(
 673                        np.asarray(self.parameters),
 674                        np.asarray(self.scores),
 675                        self.x_choices,
 676                        return_std=False,
 677                        return_pi=True,
 678                    )
 679                    y_lower = preds_with_pi[1]
 680                    self.lower = y_lower
 681
 682            else:
 683                assert (
 684                    param_distributions is not None
 685                ), "When 'param_search_init_design == False', 'param_distributions' must be provided"
 686
 687                if self.method == "bayesian":
 688                    self.posterior_ = "gaussian"
 689                    try:
 690                        y_mean, y_std = self.surrogate_fit_predict(
 691                            np.asarray(self.parameters),
 692                            np.asarray(self.scores),
 693                            self.x_choices,
 694                            return_std=True,
 695                            return_pi=False,
 696                            param_search_init_design=True,
 697                            param_distributions=param_distributions,
 698                        )
 699                    except ValueError:  # do not remove this
 700                        preds_with_std = self.surrogate_fit_predict(
 701                            np.asarray(self.parameters),
 702                            np.asarray(self.scores),
 703                            self.x_choices,
 704                            return_std=True,
 705                            return_pi=False,
 706                            param_search_init_design=True,
 707                            param_distributions=param_distributions,
 708                        )
 709                        y_mean, y_std = preds_with_std[0], preds_with_std[1]
 710                    self.y_mean = y_mean
 711                    self.y_std = np.maximum(2.220446049250313e-16, y_std)
 712
 713                elif self.method == "mc":
 714                    self.posterior_ = "mc"
 715                    assert self.surrogate_obj.__class__.__name__.startswith(
 716                        "CustomRegressor"
 717                    ) or self.surrogate_obj.__class__.__name__.startswith(
 718                        "PredictionInterval"
 719                    ), "for `method = 'mc'`, the surrogate must be a nnetsauce.CustomRegressor() or nnetsauce.PredictionInterval()"
 720                    assert (
 721                        self.surrogate_obj.replications is not None
 722                    ), "for `method = 'mc'`, the surrogate must be a nnetsauce.CustomRegressor() with a number of 'replications' provided"
 723                    preds_with_std = self.surrogate_fit_predict(
 724                        np.asarray(self.parameters),
 725                        np.asarray(self.scores),
 726                        self.x_choices,
 727                        return_std=False,
 728                        return_pi=True,
 729                        param_search_init_design=True,
 730                        param_distributions=param_distributions,
 731                    )
 732                    y_mean, y_std = preds_with_std[0], preds_with_std[1]
 733                    self.y_mean = y_mean
 734                    self.y_std = np.maximum(2.220446049250313e-16, y_std)
 735
 736                elif self.method == "splitconformal":
 737                    self.posterior_ = None
 738                    # assert self.surrogate_obj.__class__.__name__.startswith(
 739                    #    "PredictionInterval"
 740                    # ), "for `method = 'splitconformal'`, the surrogate must be a nnetsauce.PredictionInterval()"
 741                    preds_with_pi = self.surrogate_fit_predict(
 742                        np.asarray(self.parameters),
 743                        np.asarray(self.scores),
 744                        self.x_choices,
 745                        return_std=False,
 746                        return_pi=True,
 747                        param_search_init_design=True,
 748                        param_distributions=param_distributions,
 749                    )
 750                    y_lower = preds_with_pi[1]
 751                    self.lower = y_lower
 752
 753            # saving after initial design computation
 754            if self.save is not None:
 755                self.update_shelve()
 756
 757        else:  # if n_more_iter is not None
 758            assert self.n_iter > 5, "you must have n_iter > 5"
 759            n_iter = n_more_iter
 760            iter_stop = len(self.max_acq) + n_more_iter  # potentially
 761
 762        if (verbose == 1) | (verbose == 2):
 763            print(f"\n ...Done. \n")
 764            try:
 765                print(np.hstack((self.x_init, self.y_init.reshape(-1, 1))))
 766            except:
 767                pass
 768
 769        # end init design ----------
 770
 771        # if n_more_iter is None: # initial optimization, before more iters are requested
 772
 773        if (verbose == 1) | (verbose == 2):
 774            print(f"\n Optimization loop... \n")
 775
 776        # early stopping?
 777        if abs_tol is not None:
 778            assert (
 779                min_budget > 20
 780            ), "With 'abs_tol' provided, you must have 'min_budget' > 20"
 781            self.abs_tol = abs_tol
 782
 783        if ucb_tol is not None:
 784            assert (
 785                min_budget > 20
 786            ), "With 'ucb_tol' provided, you must have 'min_budget' > 20"
 787            assert (
 788                self.acquisition == "ucb"
 789            ), "With 'ucb_tol' provided, you must have 'acquisition' == 'ucb'"
 790            self.ucb_tol = ucb_tol
 791
 792        if verbose == 1:
 793            progbar = Progbar(target=n_iter)
 794
 795        # main loop ----------
 796
 797        for i in range(n_iter):
 798            # find next set of parameters (vector), maximizing acquisition function
 799            next_param = self.next_parameter_by_acq(i=i, acq=self.acquisition)
 800
 801            try:
 802                if self.per_second is True:
 803                    start = time()
 804
 805                    if self.objective_func is not None:
 806                        score_next_param = self.objective_func(
 807                            next_param, *func_args
 808                        )
 809
 810                        if (np.isfinite(score_next_param) == False) or (
 811                            np.isnan(score_next_param) == True
 812                        ):
 813                            continue
 814
 815                    else:
 816                        assert (self.x_init is not None) and (
 817                            self.y_init is not None
 818                        ), "self.objective_func is not None: must have (self.x_init is not None) and (self.y_init is not None)"
 819
 820                        print(f"\n next param: {next_param} \n")
 821                        score_next_param = float(
 822                            input("get new score: \n")
 823                        )  # or an API response
 824
 825                        if (np.isfinite(score_next_param) == False) or (
 826                            np.isnan(score_next_param) == True
 827                        ):
 828                            continue
 829
 830                    self.timings.append(np.log(time() - start))
 831
 832                else:  # self.per_second is False:
 833                    if self.objective_func is not None:
 834                        score_next_param = self.objective_func(
 835                            next_param, *func_args
 836                        )
 837
 838                        if (np.isfinite(score_next_param) == False) or (
 839                            np.isnan(score_next_param) == True
 840                        ):
 841                            continue
 842
 843                    else:
 844                        assert (self.x_init is not None) and (
 845                            self.y_init is not None
 846                        ), "self.objective_func is not None: must have (self.x_init is not None) and (self.y_init is not None)"
 847
 848                        print(f"\n next param: {next_param} \n")
 849                        score_next_param = float(
 850                            input("get new score: \n")
 851                        )  # or an API response
 852
 853                        if (np.isfinite(score_next_param) == False) or (
 854                            np.isnan(score_next_param) == True
 855                        ):
 856                            continue
 857
 858            except:
 859                continue
 860
 861            self.parameters.append(next_param.tolist())
 862
 863            self.scores.append(score_next_param)
 864
 865            if self.save is not None:
 866                self.update_shelve()
 867
 868            if verbose == 2:
 869                print(f"iteration {i + 1} -----")
 870                print(f"current minimum:  {self.x_min}")
 871                print(f"current minimum score:  {self.y_min}")
 872                print(f"next parameter: {next_param}")
 873                print(f"score for next parameter: {score_next_param} \n")
 874
 875            if score_next_param < self.y_min:
 876                self.x_min = next_param
 877                self.y_min = score_next_param
 878                if self.save is not None:
 879                    self.update_shelve()
 880                if self.y_min == self.min_value:
 881                    break
 882
 883            if self.posterior_ == "gaussian" and self.method == "bayesian":
 884                try:
 885                    self.y_mean, self.y_std = self.surrogate_fit_predict(
 886                        np.asarray(self.parameters),
 887                        np.asarray(self.scores),
 888                        self.x_choices,
 889                        return_std=True,
 890                        return_pi=False,
 891                    )
 892                except:
 893                    (
 894                        self.y_mean,
 895                        self.y_std,
 896                        lower,
 897                        upper,
 898                    ) = self.surrogate_fit_predict(
 899                        np.asarray(self.parameters),
 900                        np.asarray(self.scores),
 901                        self.x_choices,
 902                        return_std=True,
 903                        return_pi=False,
 904                    )
 905
 906            elif self.posterior_ in (None, "mc") and self.method in (
 907                "mc",
 908                "splitconformal",
 909            ):
 910                (
 911                    self.y_mean,
 912                    self.y_lower,
 913                    self.y_upper,
 914                ) = self.surrogate_fit_predict(
 915                    np.asarray(self.parameters),
 916                    np.asarray(self.scores),
 917                    self.x_choices,
 918                    return_std=False,
 919                    return_pi=True,
 920                )
 921
 922            else:
 923                return NotImplementedError
 924
 925            if self.save is not None:
 926                self.update_shelve()
 927
 928            if verbose == 1:
 929                progbar.update(i + 1)  # update progress bar
 930
 931            # early stopping
 932
 933            if abs_tol is not None:
 934                # if self.max_acq.size > (self.n_init + self.n_iter * min_budget_pct):
 935                if len(self.max_acq) > min_budget:
 936                    diff_max_acq = np.abs(np.diff(np.asarray(self.max_acq)))
 937
 938                    if diff_max_acq[-1] <= abs_tol:
 939                        iter_stop = len(self.max_acq)  # index i starts at 0
 940
 941                        break
 942
 943            if ucb_tol is not None:
 944                if len(self.max_acq) > min_budget:
 945                    # print(f"self.ucb: {self.ucb}")
 946                    # print(f"self.acq: {self.acq}")
 947                    # print(f"mean(self.ucb/self.acq): {np.mean(self.ucb/self.acq)/100}")
 948
 949                    if (
 950                        np.abs(np.mean(self.ucb / self.acq) / 100) <= ucb_tol
 951                    ):  # self.ucb is the upper confidence bound for UCB criterion
 952                        iter_stop = len(self.max_acq)
 953
 954                        break
 955
 956        # end main loop ----------
 957
 958        if (verbose == 1) & (i < (n_iter - 1)):
 959            progbar.update(n_iter)
 960
 961        self.n_iter = iter_stop
 962        if self.save is not None:
 963            self.update_shelve()
 964
 965        DescribeResult = namedtuple(
 966            "DescribeResult", ("best_params", "best_score")
 967        )
 968
 969        if self.params_names is None:
 970            return DescribeResult(self.x_min, self.y_min)
 971
 972        else:
 973            return DescribeResult(
 974                dict(zip(self.params_names, self.x_min)), self.y_min
 975            )
 976
 977    # optimize the objective
 978    def lazyoptimize(
 979        self,
 980        verbose=1,
 981        abs_tol=None,  # suggested 1e-4, for n_iter = 200
 982        ucb_tol=None,
 983        min_budget=50,  # minimum budget for early stopping
 984        func_args=None,
 985        estimators="all",
 986        type_pi="kde",  # for now, 'kde', 'bootstrap', 'splitconformal'
 987        type_exec="independent",  # "queue" or "independent" (default)
 988    ):
 989        """Launch optimization loop.
 990
 991        # Arguments:
 992
 993            verbose: an integer;
 994                verbose = 0: nothing is printed,
 995                verbose = 1: a progress bar is printed (longer than 0),
 996                verbose = 2: information about each iteration is printed (longer than 1)
 997
 998            n_more_iter: an integer;
 999                additional number of iterations for the optimizer (which has been run once)
1000
1001            abs_tol: a float;
1002                tolerance for convergence of the optimizer (early stopping based on expected improvement)
1003
1004            ucb_tol: a float;
1005                tolerance for convergence of the optimizer (early stopping based on length of prediction intervals)
1006
1007            min_budget: an integer (default is 50);
1008                minimum number of iterations before early stopping controlled by `abs_tol`
1009
1010            func_args: a list;
1011                additional parameters for the objective function (if necessary)
1012
1013            estimators: an str or a list of strs (estimators names)
1014                if "all", then 30 models are fitted. Otherwise, only those provided in the list
1015                are adjusted; for example ["RandomForestRegressor", "Ridge"]
1016
1017            type_pi: an str;
1018                "kde" (default) or, "splitconformal"; type of prediction intervals for the surrogate
1019                model
1020
1021            type_exec: an str;
1022                "independent" (default) is when surrogate models are adjusted independently on
1023                the same design set and the best model is chosen eventually; "queue" is when
1024                surrogate models are adjusted one after the other, on a design set with
1025                increasing size;
1026
1027        """
1028
1029        # Base case: Gaussian Process
1030        gp_opt_obj = GPOpt(
1031            objective_func=self.objective_func,
1032            lower_bound=self.lower_bound,
1033            upper_bound=self.upper_bound,
1034            n_init=self.n_init,
1035            n_iter=self.n_iter,
1036            alpha=self.alpha,
1037            n_restarts_optimizer=self.n_restarts_optimizer,
1038            seed=self.seed,
1039            n_jobs=self.n_jobs,
1040            acquisition="ei",
1041            surrogate_obj=GaussianProcessRegressor(
1042                kernel=Matern(nu=2.5),
1043                alpha=self.alpha,
1044                normalize_y=True,
1045                n_restarts_optimizer=self.n_restarts_optimizer,
1046                random_state=self.seed,
1047            ),
1048        )
1049        if verbose == 2:
1050            print(
1051                f"\n adjusting surrogate model # {0} (GaussianProcessRegressor(Matern(5/2)))... \n"
1052            )
1053        res_base = gp_opt_obj.optimize(
1054            verbose=verbose,
1055            abs_tol=abs_tol,  # suggested 1e-4, for n_iter = 200
1056            min_budget=min_budget,  # minimum budget for early stopping
1057            func_args=func_args,
1058        )
1059
1060        if estimators == "all":
1061            if type_pi == "kde":
1062                self.regressors = REGRESSORS
1063
1064            else:
1065                self.regressors = [
1066                    (
1067                        est[0],
1068                        ns.PredictionInterval(
1069                            est[1](), type_pi="splitconformal"
1070                        ),
1071                    )
1072                    for est in all_estimators()
1073                    if (
1074                        issubclass(est[1], RegressorMixin)
1075                        and (est[0] not in REMOVED_REGRESSORS)
1076                    )
1077                ]
1078
1079        else:
1080            if type_pi == "kde":
1081                self.regressors = [
1082                    (
1083                        "CustomRegressor(" + est[0] + ")",
1084                        ns.CustomRegressor(
1085                            est[1](), replications=150, type_pi=type_pi
1086                        ),
1087                    )
1088                    for est in all_estimators()
1089                    if (
1090                        issubclass(est[1], RegressorMixin)
1091                        and (est[0] not in REMOVED_REGRESSORS)
1092                        and (est[0] in estimators)
1093                    )
1094                ]
1095
1096            elif type_pi == "splitconformal":
1097                self.regressors = [
1098                    (
1099                        est[0],
1100                        ns.PredictionInterval(
1101                            est[1](), type_pi="splitconformal"
1102                        ),
1103                    )
1104                    for est in all_estimators()
1105                    if (
1106                        issubclass(est[1], RegressorMixin)
1107                        and (est[0] not in REMOVED_REGRESSORS)
1108                        and (est[0] in estimators)
1109                    )
1110                ]
1111
1112        df_res = pd.DataFrame(
1113            np.empty((len(self.regressors) + 1, 2)), columns=["Model", "Score"]
1114        )
1115        df_res.iloc[0, 0] = "GaussianProcessRegressor"
1116        df_res.iloc[0, 1] = res_base.best_score
1117
1118        self.surrogate_fit_predict = partial(
1119            self.surrogate_fit_predict, return_pi=True
1120        )
1121
1122        if (
1123            type_exec == "queue"
1124        ):  # when models are adjusted one after the other on a design set with increasing size
1125            self.x_min = None
1126
1127            self.y_min = np.inf
1128
1129            score_next_param = np.inf
1130
1131            DescribeResult = namedtuple(
1132                "DescribeResult", ("best_params", "best_score", "scores")
1133            )
1134
1135            if verbose == 2:
1136                print(
1137                    f"\n adjusting surrogate model # {1} ({self.regressors[0][0]})... \n"
1138                )
1139
1140            df_res.iloc[0, 0] = self.regressors[0][0]
1141
1142            gp_opt_obj_prev = GPOpt(
1143                objective_func=self.objective_func,
1144                lower_bound=self.lower_bound,
1145                upper_bound=self.upper_bound,
1146                n_init=self.n_init,
1147                n_iter=self.n_iter,
1148                alpha=self.alpha,
1149                n_restarts_optimizer=self.n_restarts_optimizer,
1150                seed=self.seed,
1151                n_jobs=self.n_jobs,
1152                acquisition=self.acquisition,
1153                method=self.method,
1154                min_value=self.min_value,
1155                surrogate_obj=copy.deepcopy(self.regressors[0][1]),
1156            )
1157
1158            gp_opt_obj_prev.optimize(
1159                verbose=verbose,
1160                abs_tol=abs_tol,  # suggested 1e-4, for n_iter = 200
1161                ucb_tol=ucb_tol,
1162                min_budget=min_budget,  # minimum budget for early stopping
1163                func_args=func_args,
1164            )
1165
1166            score_next_param = gp_opt_obj_prev.y_min
1167
1168            df_res.iloc[0, 1] = score_next_param
1169
1170            if self.n_jobs is None:  # sequential optimization
1171                for i in range(len(self.regressors)):
1172                    try:
1173                        if verbose == 2:
1174                            print(
1175                                f"\n adjusting surrogate model # {i + 2} ({self.regressors[i][0]})... \n"
1176                            )
1177
1178                        df_res["Model"][i + 1] = self.regressors[i][0]
1179
1180                        gp_opt_obj = GPOpt(
1181                            objective_func=self.objective_func,
1182                            lower_bound=self.lower_bound,
1183                            upper_bound=self.upper_bound,
1184                            n_init=self.n_init,
1185                            n_iter=self.n_iter,
1186                            alpha=self.alpha,
1187                            n_restarts_optimizer=self.n_restarts_optimizer,
1188                            seed=self.seed,
1189                            n_jobs=self.n_jobs,
1190                            acquisition=self.acquisition,
1191                            method=self.method,
1192                            min_value=self.min_value,
1193                            surrogate_obj=copy.deepcopy(self.regressors[i][1]),
1194                            x_init=np.asarray(gp_opt_obj_prev.parameters),
1195                            y_init=np.asarray(gp_opt_obj_prev.scores),
1196                        )
1197
1198                        gp_opt_obj.optimize(
1199                            verbose=verbose,
1200                            abs_tol=abs_tol,  # suggested 1e-4, for n_iter = 200
1201                            ucb_tol=ucb_tol,
1202                            min_budget=min_budget,  # minimum budget for early stopping
1203                            func_args=func_args,
1204                        )
1205
1206                        score_next_param = gp_opt_obj.y_min
1207
1208                        df_res.iloc[i, 1] = score_next_param
1209
1210                        if score_next_param < self.y_min:
1211                            self.x_min = gp_opt_obj.x_min
1212                            self.y_min = score_next_param
1213                            if self.y_min == self.min_value:
1214                                break
1215
1216                        if verbose == 2:
1217                            print(f"Global iteration #{i + 1} -----")
1218                            print(f"current minimum:  {self.x_min}")
1219                            print(f"current minimum score:  {self.y_min}")
1220                            print(
1221                                f"score for next parameter: {score_next_param} \n"
1222                            )
1223
1224                        gp_opt_obj_prev = copy.deepcopy(gp_opt_obj)
1225
1226                    except ValueError:
1227                        continue
1228
1229            elif self.n_jobs >= 2 or self.n_jobs == -1:  # parallel optimization
1230                pass
1231            else:
1232                raise ValueError(
1233                    "n_jobs must be either None or >= 2 or equal to -1"
1234                )
1235            return DescribeResult(
1236                self.x_min, self.y_min, df_res.sort_values(by="Score")
1237            )
1238
1239        elif (
1240            type_exec == "independent"
1241        ):  # when models are adjusted independently on the same design set and the best model is chosen eventually
1242            self.x_min = None
1243
1244            self.y_min = np.inf
1245
1246            score_next_param = np.inf
1247
1248            DescribeResult = namedtuple(
1249                "DescribeResult",
1250                ("best_params", "best_score", "best_surrogate", "scores"),
1251            )
1252
1253            if verbose == 2:
1254                print(
1255                    f"\n adjusting surrogate model # {1} ({self.regressors[0][0]})... \n"
1256                )
1257
1258            if self.n_jobs is None:  # sequential optimization
1259                for i in range(len(self.regressors)):
1260                    # try:
1261
1262                    if verbose == 2:
1263                        print(
1264                            f"\n adjusting surrogate model # {i + 1} ({self.regressors[i][0]})... \n"
1265                        )
1266
1267                    df_res["Model"][i + 1] = self.regressors[i][0]
1268
1269                    gp_opt_obj = GPOpt(
1270                        objective_func=self.objective_func,
1271                        lower_bound=self.lower_bound,
1272                        upper_bound=self.upper_bound,
1273                        n_init=self.n_init,
1274                        n_iter=self.n_iter,
1275                        alpha=self.alpha,
1276                        n_restarts_optimizer=self.n_restarts_optimizer,
1277                        seed=self.seed,
1278                        n_jobs=self.n_jobs,
1279                        acquisition=self.acquisition,
1280                        method=self.method,
1281                        min_value=self.min_value,
1282                        surrogate_obj=copy.deepcopy(self.regressors[i][1]),
1283                    )
1284
1285                    gp_opt_obj.optimize(
1286                        verbose=verbose,
1287                        abs_tol=abs_tol,  # suggested 1e-4, for n_iter = 200
1288                        ucb_tol=ucb_tol,
1289                        min_budget=min_budget,  # minimum budget for early stopping
1290                        func_args=func_args,
1291                    )
1292
1293                    score_next_param = gp_opt_obj.y_min
1294
1295                    df_res.iloc[i, 1] = score_next_param
1296
1297                    if score_next_param < self.y_min:
1298                        self.x_min = gp_opt_obj.x_min
1299                        self.y_min = score_next_param
1300                        self.best_surrogate = copy.deepcopy(
1301                            gp_opt_obj.surrogate_obj
1302                        )
1303                        if self.y_min == self.min_value:
1304                            break
1305
1306                    if verbose == 2:
1307                        print(f"Global iteration #{i + 1} -----")
1308                        print(f"current minimum:  {self.x_min}")
1309                        print(f"current minimum score:  {self.y_min}")
1310                        print(
1311                            f"score for next parameter: {score_next_param} \n"
1312                        )
1313
1314                    # except ValueError:
1315
1316                    #    continue
1317
1318            elif self.n_jobs >= 2 or self.n_jobs == -1:  # parallel optimization
1319
1320                def foo(i):
1321                    df_res["Model"][i + 1] = self.regressors[i][0]
1322
1323                    gp_opt_obj = GPOpt(
1324                        objective_func=self.objective_func,
1325                        lower_bound=self.lower_bound,
1326                        upper_bound=self.upper_bound,
1327                        n_init=self.n_init,
1328                        n_iter=self.n_iter,
1329                        alpha=self.alpha,
1330                        n_restarts_optimizer=self.n_restarts_optimizer,
1331                        seed=self.seed,
1332                        n_jobs=self.n_jobs,
1333                        acquisition=self.acquisition,
1334                        method=self.method,
1335                        min_value=self.min_value,
1336                        surrogate_obj=copy.deepcopy(self.regressors[i][1]),
1337                    )
1338
1339                    try:
1340                        gp_opt_obj.optimize(
1341                            verbose=0,  # important
1342                            abs_tol=abs_tol,  # suggested 1e-4, for n_iter = 200
1343                            ucb_tol=ucb_tol,
1344                            min_budget=min_budget,  # minimum budget for early stopping
1345                            func_args=func_args,
1346                        )
1347
1348                        return gp_opt_obj
1349                    except ValueError:
1350                        return None
1351
1352                tmp_results = Parallel(n_jobs=self.n_jobs)(
1353                    delayed(foo)(i) for i in tqdm(range(len(self.regressors)))
1354                )
1355
1356                for i in tqdm(range(len(tmp_results))):
1357                    if tmp_results[i] is not None:
1358                        gp_opt_obj = copy.deepcopy(tmp_results[i])
1359
1360                        score_next_param = gp_opt_obj.y_min
1361
1362                        df_res.iloc[i, 1] = score_next_param
1363
1364                        if score_next_param < self.y_min:
1365                            self.x_min = gp_opt_obj.x_min
1366                            self.y_min = score_next_param
1367                            self.best_surrogate = copy.deepcopy(
1368                                gp_opt_obj.surrogate_obj
1369                            )
1370                            if self.y_min == self.min_value:
1371                                break
1372
1373                        if verbose == 2:
1374                            print(f"Global iteration #{i + 1} -----")
1375                            print(f"current minimum:  {self.x_min}")
1376                            print(f"current minimum score:  {self.y_min}")
1377                            print(
1378                                f"score for next parameter: {score_next_param} \n"
1379                            )
1380
1381                    else:
1382                        continue
1383
1384            else:
1385                raise ValueError(
1386                    "n_jobs must be either None or >= 2 or equal to -1"
1387                )
1388            return DescribeResult(
1389                self.x_min,
1390                self.y_min,
1391                self.best_surrogate,
1392                df_res.sort_values(by="Score"),
1393            )
1394
1395        else:
1396            NotImplementedError(
1397                "type_exec must be either 'queue' or 'independent'"
1398            )