nnetsauce

 21class AdaBoostClassifier(Boosting, ClassifierMixin):
 22    """AdaBoost Classification (SAMME) model class derived from class Boosting
 23
 24    Parameters:
 25
 26        obj: object
 27            any object containing a method fit (obj.fit()) and a method predict
 28            (obj.predict())
 29
 30        n_estimators: int
 31            number of boosting iterations
 32
 33        learning_rate: float
 34            learning rate of the boosting procedure
 35
 36        n_hidden_features: int
 37            number of nodes in the hidden layer
 38
 39        reg_lambda: float
 40            regularization parameter for weights
 41
 42        reg_alpha: float
 43            controls compromize between l1 and l2 norm of weights
 44
 45        activation_name: str
 46            activation function: 'relu', 'tanh', 'sigmoid', 'prelu' or 'elu'
 47
 48        a: float
 49            hyperparameter for 'prelu' or 'elu' activation function
 50
 51        nodes_sim: str
 52            type of simulation for the nodes: 'sobol', 'hammersley', 'halton',
 53            'uniform'
 54
 55        bias: boolean
 56            indicates if the hidden layer contains a bias term (True) or not
 57            (False)
 58
 59        dropout: float
 60            regularization parameter; (random) percentage of nodes dropped out
 61            of the training
 62
 63        direct_link: boolean
 64            indicates if the original predictors are included (True) in model's
 65            fitting or not (False)
 66
 67        n_clusters: int
 68            number of clusters for 'kmeans' or 'gmm' clustering (could be 0:
 69                no clustering)
 70
 71        cluster_encode: bool
 72            defines how the variable containing clusters is treated (default is one-hot)
 73            if `False`, then labels are used, without one-hot encoding
 74
 75        type_clust: str
 76            type of clustering method: currently k-means ('kmeans') or Gaussian
 77            Mixture Model ('gmm')
 78
 79        type_scaling: a tuple of 3 strings
 80            scaling methods for inputs, hidden layer, and clustering respectively
 81            (and when relevant).
 82            Currently available: standardization ('std') or MinMax scaling ('minmax')
 83
 84        col_sample: float
 85            percentage of covariates randomly chosen for training
 86
 87        row_sample: float
 88            percentage of rows chosen for training, by stratified bootstrapping
 89
 90        seed: int
 91            reproducibility seed for nodes_sim=='uniform'
 92
 93        verbose: int
 94            0 for no output, 1 for a progress bar (default is 1)
 95
 96        method: str
 97            type of Adaboost method, 'SAMME' (discrete) or 'SAMME.R' (real)
 98
 99        backend: str
100            "cpu" or "gpu" or "tpu"
101
102    Attributes:
103
104        alpha_: list
105            AdaBoost coefficients alpha_m
106
107        base_learners_: dict
108            a dictionary containing the base learners
109
110    Examples:
111
112    See also [https://github.com/Techtonique/nnetsauce/blob/master/examples/adaboost_classification.py](https://github.com/Techtonique/nnetsauce/blob/master/examples/adaboost_classification.py)
113
114    ```python
115    import nnetsauce as ns
116    import numpy as np
117    from sklearn.datasets import load_breast_cancer
118    from sklearn.linear_model import LogisticRegression
119    from sklearn.model_selection import train_test_split
120    from sklearn import metrics
121    from time import time
122
123    breast_cancer = load_breast_cancer()
124    Z = breast_cancer.data
125    t = breast_cancer.target
126    np.random.seed(123)
127    X_train, X_test, y_train, y_test = train_test_split(Z, t, test_size=0.2)
128
129    # SAMME.R
130    clf = LogisticRegression(solver='liblinear', multi_class = 'ovr',
131                            random_state=123)
132    fit_obj = ns.AdaBoostClassifier(clf,
133                                    n_hidden_features=int(11.22338867),
134                                    direct_link=True,
135                                    n_estimators=250, learning_rate=0.01126343,
136                                    col_sample=0.72684326, row_sample=0.86429443,
137                                    dropout=0.63078613, n_clusters=2,
138                                    type_clust="gmm",
139                                    verbose=1, seed = 123,
140                                    method="SAMME.R")
141
142    start = time()
143    fit_obj.fit(X_train, y_train)
144    print(f"Elapsed {time() - start}")
145
146    start = time()
147    print(fit_obj.score(X_test, y_test))
148    print(f"Elapsed {time() - start}")
149
150    preds = fit_obj.predict(X_test)
151
152    print(metrics.classification_report(preds, y_test))
153
154    ```
155
156    """
157
158    # construct the object -----
159    _estimator_type = "classifier"
160
161    def __init__(
162        self,
163        obj,
164        n_estimators=10,
165        learning_rate=0.1,
166        n_hidden_features=1,
167        reg_lambda=0,
168        reg_alpha=0.5,
169        activation_name="relu",
170        a=0.01,
171        nodes_sim="sobol",
172        bias=True,
173        dropout=0,
174        direct_link=False,
175        n_clusters=2,
176        cluster_encode=True,
177        type_clust="kmeans",
178        type_scaling=("std", "std", "std"),
179        col_sample=1,
180        row_sample=1,
181        seed=123,
182        verbose=1,
183        method="SAMME",
184        backend="cpu",
185    ):
186        self.type_fit = "classification"
187        self.verbose = verbose
188        self.method = method
189        self.reg_lambda = reg_lambda
190        self.reg_alpha = reg_alpha
191
192        super().__init__(
193            obj=obj,
194            n_estimators=n_estimators,
195            learning_rate=learning_rate,
196            n_hidden_features=n_hidden_features,
197            activation_name=activation_name,
198            a=a,
199            nodes_sim=nodes_sim,
200            bias=bias,
201            dropout=dropout,
202            direct_link=direct_link,
203            n_clusters=n_clusters,
204            cluster_encode=cluster_encode,
205            type_clust=type_clust,
206            type_scaling=type_scaling,
207            col_sample=col_sample,
208            row_sample=row_sample,
209            seed=seed,
210            backend=backend,
211        )
212
213        self.alpha_ = []
214        self.base_learners_ = dict.fromkeys(range(n_estimators))
215
216    def fit(self, X, y, sample_weight=None, **kwargs):
217        """Fit Boosting model to training data (X, y).
218
219        Parameters:
220
221            X: {array-like}, shape = [n_samples, n_features]
222                Training vectors, where n_samples is the number
223                of samples and n_features is the number of features.
224
225            y: array-like, shape = [n_samples]
226                Target values.
227
228            **kwargs: additional parameters to be passed to
229                    self.cook_training_set or self.obj.fit
230
231        Returns:
232
233             self: object
234        """
235
236        assert mx.is_factor(y), "y must contain only integers"
237
238        assert self.method in (
239            "SAMME",
240            "SAMME.R",
241        ), "`method` must be either 'SAMME' or 'SAMME.R'"
242
243        assert (self.reg_lambda <= 1) & (
244            self.reg_lambda >= 0
245        ), "must have self.reg_lambda <= 1 &  self.reg_lambda >= 0"
246
247        assert (self.reg_alpha <= 1) & (
248            self.reg_alpha >= 0
249        ), "must have self.reg_alpha <= 1 &  self.reg_alpha >= 0"
250
251        # training
252        n, p = X.shape
253        self.n_classes = len(np.unique(y))
254        self.classes_ = np.unique(y)  # for compatibility with sklearn
255        self.n_classes_ = len(self.classes_)  # for compatibility with sklearn
256
257        if sample_weight is None:
258            w_m = np.repeat(1.0 / n, n)
259        else:
260            w_m = np.asarray(sample_weight)
261
262        base_learner = CustomClassifier(
263            self.obj,
264            n_hidden_features=self.n_hidden_features,
265            activation_name=self.activation_name,
266            a=self.a,
267            nodes_sim=self.nodes_sim,
268            bias=self.bias,
269            dropout=self.dropout,
270            direct_link=self.direct_link,
271            n_clusters=self.n_clusters,
272            type_clust=self.type_clust,
273            type_scaling=self.type_scaling,
274            col_sample=self.col_sample,
275            row_sample=self.row_sample,
276            seed=self.seed,
277        )
278
279        if self.verbose == 1:
280            pbar = Progbar(self.n_estimators)
281
282        if self.method == "SAMME":
283            err_m = 1e6
284            err_bound = 1 - 1 / self.n_classes
285            self.alpha_.append(1.0)
286            x_range_n = range(n)
287
288            for m in range(self.n_estimators):
289                preds = base_learner.fit(
290                    X, y, sample_weight=w_m.ravel(), **kwargs
291                ).predict(X)
292
293                self.base_learners_.update({m: deepcopy(base_learner)})
294
295                cond = [y[i] != preds[i] for i in x_range_n]
296
297                err_m = max(
298                    sum([elt[0] * elt[1] for elt in zip(cond, w_m)]),
299                    2.220446049250313e-16,
300                )  # sum(w_m) == 1
301
302                if self.reg_lambda > 0:
303                    err_m += self.reg_lambda * (
304                        (1 - self.reg_alpha) * 0.5 * sum([x**2 for x in w_m])
305                        + self.reg_alpha * sum([abs(x) for x in w_m])
306                    )
307
308                err_m = min(err_m, err_bound)
309
310                alpha_m = self.learning_rate * log(
311                    (self.n_classes - 1) * (1 - err_m) / err_m
312                )
313
314                self.alpha_.append(alpha_m)
315
316                w_m_temp = [exp(alpha_m * cond[i]) for i in x_range_n]
317
318                sum_w_m = sum(w_m_temp)
319
320                w_m = np.asarray([w_m_temp[i] / sum_w_m for i in x_range_n])
321
322                base_learner.set_params(seed=self.seed + (m + 1) * 1000)
323
324                if self.verbose == 1:
325                    pbar.update(m)
326
327            if self.verbose == 1:
328                pbar.update(self.n_estimators)
329
330            self.n_estimators = len(self.base_learners_)
331            self.classes_ = np.unique(y)
332
333            return self
334
335        if self.method == "SAMME.R":
336            Y = mo.one_hot_encode2(y, self.n_classes)
337
338            if sample_weight is None:
339                w_m = np.repeat(1.0 / n, n)  # (N, 1)
340
341            else:
342                w_m = np.asarray(sample_weight)
343
344            for m in range(self.n_estimators):
345                probs = base_learner.fit(
346                    X, y, sample_weight=w_m.ravel(), **kwargs
347                ).predict_proba(X)
348
349                np.clip(
350                    a=probs, a_min=2.220446049250313e-16, a_max=1.0, out=probs
351                )
352
353                self.base_learners_.update({m: deepcopy(base_learner)})
354
355                w_m *= np.exp(
356                    -1.0
357                    * self.learning_rate
358                    * (1.0 - 1.0 / self.n_classes)
359                    * xlogy(Y, probs).sum(axis=1)
360                )
361
362                w_m /= np.sum(w_m)
363
364                base_learner.set_params(seed=self.seed + (m + 1) * 1000)
365
366                if self.verbose == 1:
367                    pbar.update(m)
368
369            if self.verbose == 1:
370                pbar.update(self.n_estimators)
371
372            self.n_estimators = len(self.base_learners_)
373            self.classes_ = np.unique(y)
374
375            return self
376
377    def predict(self, X, **kwargs):
378        """Predict test data X.
379
380        Parameters:
381
382            X: {array-like}, shape = [n_samples, n_features]
383                Training vectors, where n_samples is the number
384                of samples and n_features is the number of features.
385
386            **kwargs: additional parameters to be passed to
387                  self.cook_test_set
388
389        Returns:
390
391            model predictions: {array-like}
392        """
393        return self.predict_proba(X, **kwargs).argmax(axis=1)
394
395    def predict_proba(self, X, **kwargs):
396        """Predict probabilities for test data X.
397
398        Parameters:
399
400            X: {array-like}, shape = [n_samples, n_features]
401                Training vectors, where n_samples is the number
402                of samples and n_features is the number of features.
403
404            **kwargs: additional parameters to be passed to
405                  self.cook_test_set
406
407        Returns:
408
409            probability estimates for test data: {array-like}
410
411        """
412
413        n_iter = len(self.base_learners_)
414
415        if self.method == "SAMME":
416            ensemble_learner = np.zeros((X.shape[0], self.n_classes))
417
418            # if self.verbose == 1:
419            #    pbar = Progbar(n_iter)
420
421            for idx, base_learner in self.base_learners_.items():
422                preds = base_learner.predict(X, **kwargs)
423
424                ensemble_learner += self.alpha_[idx] * mo.one_hot_encode2(
425                    preds, self.n_classes
426                )
427
428                # if self.verbose == 1:
429                #    pbar.update(idx)
430
431            # if self.verbose == 1:
432            #    pbar.update(n_iter)
433
434            expit_ensemble_learner = expit(ensemble_learner)
435
436            sum_ensemble = expit_ensemble_learner.sum(axis=1)
437
438            return expit_ensemble_learner / sum_ensemble[:, None]
439
440        # if self.method == "SAMME.R":
441        ensemble_learner = 0
442
443        # if self.verbose == 1:
444        #    pbar = Progbar(n_iter)
445
446        for idx, base_learner in self.base_learners_.items():
447            probs = base_learner.predict_proba(X, **kwargs)
448
449            np.clip(a=probs, a_min=2.220446049250313e-16, a_max=1.0, out=probs)
450
451            log_preds_proba = np.log(probs)
452
453            ensemble_learner += (
454                log_preds_proba - log_preds_proba.mean(axis=1)[:, None]
455            )
456
457            # if self.verbose == 1:
458            #    pbar.update(idx)
459
460        ensemble_learner *= self.n_classes - 1
461
462        # if self.verbose == 1:
463        #    pbar.update(n_iter)
464
465        expit_ensemble_learner = expit(ensemble_learner)
466
467        sum_ensemble = expit_ensemble_learner.sum(axis=1)
468
469        return expit_ensemble_learner / sum_ensemble[:, None]
470
471    @property
472    def _estimator_type(self):
473        return "classifier"

  48class Base(BaseEstimator):
  49    """Base model from which all the other classes inherit.
  50
  51    This class contains the most important data preprocessing/feature engineering methods.
  52
  53    Parameters:
  54
  55        n_hidden_features: int
  56            number of nodes in the hidden layer
  57
  58        activation_name: str
  59            activation function: 'relu', 'tanh', 'sigmoid', 'prelu' or 'elu'
  60
  61        a: float
  62            hyperparameter for 'prelu' or 'elu' activation function
  63
  64        nodes_sim: str
  65            type of simulation for hidden layer nodes: 'sobol', 'hammersley', 'halton',
  66            'uniform'
  67
  68        bias: boolean
  69            indicates if the hidden layer contains a bias term (True) or
  70            not (False)
  71
  72        dropout: float
  73            regularization parameter; (random) percentage of nodes dropped out
  74            of the training
  75
  76        direct_link: boolean
  77            indicates if the original features are included (True) in model's
  78            fitting or not (False)
  79
  80        n_clusters: int
  81            number of clusters for type_clust='kmeans' or type_clust='gmm'
  82            clustering (could be 0: no clustering)
  83
  84        cluster_encode: bool
  85            defines how the variable containing clusters is treated (default is one-hot);
  86            if `False`, then labels are used, without one-hot encoding
  87
  88        type_clust: str
  89            type of clustering method: currently k-means ('kmeans') or Gaussian
  90            Mixture Model ('gmm')
  91
  92        type_scaling: a tuple of 3 strings
  93            scaling methods for inputs, hidden layer, and clustering respectively
  94            (and when relevant).
  95            Currently available: standardization ('std') or MinMax scaling ('minmax') or robust scaling ('robust') or  max absolute scaling ('maxabs')
  96
  97        col_sample: float
  98            percentage of features randomly chosen for training
  99
 100        row_sample: float
 101            percentage of rows chosen for training, by stratified bootstrapping
 102
 103        seed: int
 104            reproducibility seed for nodes_sim=='uniform', clustering and dropout
 105
 106        backend: str
 107            "cpu" or "gpu" or "tpu"
 108
 109    """
 110
 111    # construct the object -----
 112
 113    def __init__(
 114        self,
 115        n_hidden_features=5,
 116        activation_name="relu",
 117        a=0.01,
 118        nodes_sim="sobol",
 119        bias=True,
 120        dropout=0,
 121        direct_link=True,
 122        n_clusters=2,
 123        cluster_encode=True,
 124        type_clust="kmeans",
 125        type_scaling=("std", "std", "std"),
 126        col_sample=1,
 127        row_sample=1,
 128        seed=123,
 129        backend="cpu",
 130    ):
 131        # input checks -----
 132
 133        sys_platform = platform.system()
 134
 135        if (sys_platform == "Windows") and (backend in ("gpu", "tpu")):
 136            warnings.warn(
 137                "No GPU/TPU computing on Windows yet, backend set to 'cpu'"
 138            )
 139            backend = "cpu"
 140
 141        assert activation_name in (
 142            "relu",
 143            "tanh",
 144            "sigmoid",
 145            "prelu",
 146            "elu",
 147        ), "'activation_name' must be in ('relu', 'tanh', 'sigmoid','prelu', 'elu')"
 148
 149        assert nodes_sim in (
 150            "sobol",
 151            "hammersley",
 152            "uniform",
 153            "halton",
 154        ), "'nodes_sim' must be in ('sobol', 'hammersley', 'uniform', 'halton')"
 155
 156        assert type_clust in (
 157            "kmeans",
 158            "gmm",
 159        ), "'type_clust' must be in ('kmeans', 'gmm')"
 160
 161        assert (len(type_scaling) == 3) & all(
 162            type_scaling[i] in ("minmax", "std", "robust", "maxabs")
 163            for i in range(len(type_scaling))
 164        ), "'type_scaling' must have length 3, and available scaling methods are 'minmax' scaling, standardization ('std'), robust scaling ('robust') and max absolute ('maxabs')"
 165
 166        assert (col_sample >= 0) & (
 167            col_sample <= 1
 168        ), "'col_sample' must be comprised between 0 and 1 (both included)"
 169
 170        assert backend in (
 171            "cpu",
 172            "gpu",
 173            "tpu",
 174        ), "must have 'backend' in ('cpu', 'gpu', 'tpu')"
 175
 176        self.n_hidden_features = n_hidden_features
 177        self.activation_name = activation_name
 178        self.a = a
 179        self.nodes_sim = nodes_sim
 180        self.bias = bias
 181        self.seed = seed
 182        self.backend = backend
 183        self.dropout = dropout
 184        self.direct_link = direct_link
 185        self.cluster_encode = cluster_encode
 186        self.type_clust = type_clust
 187        self.type_scaling = type_scaling
 188        self.col_sample = col_sample
 189        self.row_sample = row_sample
 190        self.n_clusters = n_clusters
 191        if isinstance(self, RegressorMixin):
 192            self.type_fit = "regression"
 193        elif isinstance(self, ClassifierMixin):
 194            self.type_fit = "classification"
 195        self.subsampler_ = None
 196        self.index_col_ = None
 197        self.index_row_ = True
 198        self.clustering_obj_ = None
 199        self.clustering_scaler_ = None
 200        self.nn_scaler_ = None
 201        self.scaler_ = None
 202        self.encoder_ = None
 203        self.W_ = None
 204        self.X_ = None
 205        self.y_ = None
 206        self.y_mean_ = None
 207        self.beta_ = None
 208
 209        # activation function -----
 210        if sys_platform in ("Linux", "Darwin"):
 211            activation_options = {
 212                "relu": ac.relu if (self.backend == "cpu") else jnn.relu,
 213                "tanh": np.tanh if (self.backend == "cpu") else jnp.tanh,
 214                "sigmoid": (
 215                    ac.sigmoid if (self.backend == "cpu") else jnn.sigmoid
 216                ),
 217                "prelu": partial(ac.prelu, a=a),
 218                "elu": (
 219                    partial(ac.elu, a=a)
 220                    if (self.backend == "cpu")
 221                    else partial(jnn.elu, a=a)
 222                ),
 223            }
 224        else:  # on Windows currently, no JAX
 225            activation_options = {
 226                "relu": (
 227                    ac.relu if (self.backend == "cpu") else NotImplementedError
 228                ),
 229                "tanh": (
 230                    np.tanh if (self.backend == "cpu") else NotImplementedError
 231                ),
 232                "sigmoid": (
 233                    ac.sigmoid
 234                    if (self.backend == "cpu")
 235                    else NotImplementedError
 236                ),
 237                "prelu": partial(ac.prelu, a=a),
 238                "elu": (
 239                    partial(ac.elu, a=a)
 240                    if (self.backend == "cpu")
 241                    else NotImplementedError
 242                ),
 243            }
 244        self.activation_func = activation_options[activation_name]
 245
 246    # "preprocessing" methods to be inherited -----
 247
 248    def encode_clusters(self, X=None, predict=False, scaler=None, **kwargs):  #
 249        """Create new covariates with kmeans or GMM clustering
 250
 251        Parameters:
 252
 253            X: {array-like}, shape = [n_samples, n_features]
 254                Training vectors, where n_samples is the number
 255                of samples and n_features is the number of features.
 256
 257            predict: boolean
 258                is False on training set and True on test set
 259
 260            scaler: {object} of class StandardScaler, MinMaxScaler, RobustScaler or MaxAbsScaler
 261                if scaler has already been fitted on training data (online training), it can be passed here
 262
 263            **kwargs:
 264                additional parameters to be passed to the
 265                clustering method
 266
 267        Returns:
 268
 269            Clusters' matrix, one-hot encoded: {array-like}
 270
 271        """
 272
 273        np.random.seed(self.seed)
 274
 275        if X is None:
 276            X = self.X_
 277
 278        if isinstance(X, pd.DataFrame):
 279            X = copy.deepcopy(X.values.astype(float))
 280
 281        if len(X.shape) == 1:
 282            X = X.reshape(1, -1)
 283
 284        if predict is False:  # encode training set
 285
 286            # scale input data before clustering
 287            self.clustering_scaler_, scaled_X = mo.scale_covariates(
 288                X, choice=self.type_scaling[2], scaler=self.clustering_scaler_
 289            )
 290
 291            self.clustering_obj_, X_clustered = mo.cluster_covariates(
 292                scaled_X,
 293                self.n_clusters,
 294                self.seed,
 295                type_clust=self.type_clust,
 296                **kwargs
 297            )
 298
 299            if self.cluster_encode:
 300                return mo.one_hot_encode(X_clustered, self.n_clusters).astype(
 301                    np.float16
 302                )
 303
 304            return X_clustered.astype(np.float16)
 305
 306        # if predict == True, encode test set
 307        X_clustered = self.clustering_obj_.predict(
 308            self.clustering_scaler_.transform(X)
 309        )
 310
 311        if self.cluster_encode == True:
 312            return mo.one_hot_encode(X_clustered, self.n_clusters).astype(
 313                np.float16
 314            )
 315
 316        return X_clustered.astype(np.float16)
 317
 318    def create_layer(self, scaled_X, W=None):
 319        """Create hidden layer.
 320
 321        Parameters:
 322
 323            scaled_X: {array-like}, shape = [n_samples, n_features]
 324                Training vectors, where n_samples is the number
 325                of samples and n_features is the number of features
 326
 327            W: {array-like}, shape = [n_features, hidden_features]
 328                if provided, constructs the hidden layer with W; otherwise computed internally
 329
 330        Returns:
 331
 332            Hidden layer matrix: {array-like}
 333
 334        """
 335
 336        n_features = scaled_X.shape[1]
 337
 338        # hash_sim = {
 339        #         "sobol": generate_sobol,
 340        #         "hammersley": generate_hammersley,
 341        #         "uniform": generate_uniform,
 342        #         "halton": generate_halton
 343        #     }
 344
 345        if self.bias is False:  # no bias term in the hidden layer
 346            if W is None:
 347                if self.nodes_sim == "sobol":
 348                    self.W_ = generate_sobol(
 349                        n_dims=n_features,
 350                        n_points=self.n_hidden_features,
 351                        seed=self.seed,
 352                    )
 353                elif self.nodes_sim == "hammersley":
 354                    self.W_ = generate_hammersley(
 355                        n_dims=n_features,
 356                        n_points=self.n_hidden_features,
 357                        seed=self.seed,
 358                    )
 359                elif self.nodes_sim == "uniform":
 360                    self.W_ = generate_uniform(
 361                        n_dims=n_features,
 362                        n_points=self.n_hidden_features,
 363                        seed=self.seed,
 364                    )
 365                else:
 366                    self.W_ = generate_halton(
 367                        n_dims=n_features,
 368                        n_points=self.n_hidden_features,
 369                        seed=self.seed,
 370                    )
 371
 372                assert (
 373                    scaled_X.shape[1] == self.W_.shape[0]
 374                ), "check dimensions of covariates X and matrix W"
 375
 376                return mo.dropout(
 377                    x=self.activation_func(
 378                        mo.safe_sparse_dot(
 379                            a=scaled_X, b=self.W_, backend=self.backend
 380                        )
 381                    ),
 382                    drop_prob=self.dropout,
 383                    seed=self.seed,
 384                )
 385
 386            # W is not none
 387            assert (
 388                scaled_X.shape[1] == W.shape[0]
 389            ), "check dimensions of covariates X and matrix W"
 390
 391            # self.W_ = W
 392            return mo.dropout(
 393                x=self.activation_func(
 394                    mo.safe_sparse_dot(a=scaled_X, b=W, backend=self.backend)
 395                ),
 396                drop_prob=self.dropout,
 397                seed=self.seed,
 398            )
 399
 400        # with bias term in the hidden layer
 401        if W is None:
 402            n_features_1 = n_features + 1
 403
 404            if self.nodes_sim == "sobol":
 405                self.W_ = generate_sobol(
 406                    n_dims=n_features_1,
 407                    n_points=self.n_hidden_features,
 408                    seed=self.seed,
 409                )
 410            elif self.nodes_sim == "hammersley":
 411                self.W_ = generate_hammersley(
 412                    n_dims=n_features_1,
 413                    n_points=self.n_hidden_features,
 414                    seed=self.seed,
 415                )
 416            elif self.nodes_sim == "uniform":
 417                self.W_ = generate_uniform(
 418                    n_dims=n_features_1,
 419                    n_points=self.n_hidden_features,
 420                    seed=self.seed,
 421                )
 422            else:
 423                self.W_ = generate_halton(
 424                    n_dims=n_features_1,
 425                    n_points=self.n_hidden_features,
 426                    seed=self.seed,
 427                )
 428
 429            # self.W_ = hash_sim[self.nodes_sim](
 430            #         n_dims=n_features_1,
 431            #         n_points=self.n_hidden_features,
 432            #         seed=self.seed,
 433            #     )
 434
 435            return mo.dropout(
 436                x=self.activation_func(
 437                    mo.safe_sparse_dot(
 438                        a=mo.cbind(
 439                            np.ones(scaled_X.shape[0]),
 440                            scaled_X,
 441                            backend=self.backend,
 442                        ),
 443                        b=self.W_,
 444                        backend=self.backend,
 445                    )
 446                ),
 447                drop_prob=self.dropout,
 448                seed=self.seed,
 449            )
 450
 451        # W is not None
 452        # self.W_ = W
 453        return mo.dropout(
 454            x=self.activation_func(
 455                mo.safe_sparse_dot(
 456                    a=mo.cbind(
 457                        np.ones(scaled_X.shape[0]),
 458                        scaled_X,
 459                        backend=self.backend,
 460                    ),
 461                    b=W,
 462                    backend=self.backend,
 463                )
 464            ),
 465            drop_prob=self.dropout,
 466            seed=self.seed,
 467        )
 468
 469    def _jax_create_layer(
 470        self, scaled_X: jnp.ndarray, W: Optional[jnp.ndarray] = None
 471    ) -> jnp.ndarray:
 472        """JAX-compatible version of create_layer that exactly matches the original functionality."""
 473        key = jax.random.PRNGKey(self.seed)
 474        n_features = scaled_X.shape[1]
 475
 476        # Generate weights if not provided
 477        if W is None:
 478            if self.bias:
 479                n_features_1 = n_features + 1
 480                shape = (n_features_1, self.n_hidden_features)
 481            else:
 482                shape = (n_features, self.n_hidden_features)
 483
 484            # JAX-compatible weight generation matching original behavior
 485            if self.nodes_sim == "sobol":
 486                W_np = generate_sobol(
 487                    n_dims=n_features_1,
 488                    n_points=self.n_hidden_features,
 489                    seed=self.seed,
 490                )
 491                W = jnp.asarray(W_np)
 492            elif self.nodes_sim == "hammersley":
 493                W_np = generate_hammersley(
 494                    n_dims=n_features_1,
 495                    n_points=self.n_hidden_features,
 496                    seed=self.seed,
 497                )
 498                W = jnp.asarray(W_np)
 499            elif self.nodes_sim == "uniform":
 500                key, subkey = jax.random.split(key)
 501                W = jax.random.uniform(
 502                    subkey, shape=shape, minval=-1.0, maxval=1.0
 503                )
 504            else:  # halton
 505                W_np = generate_halton(
 506                    n_dims=n_features_1,
 507                    n_points=self.n_hidden_features,
 508                    seed=self.seed,
 509                )
 510                W = jnp.asarray(W_np)
 511
 512            self.W_ = np.array(W)  # Store as numpy for original methods
 513
 514        # Prepare input with bias if needed
 515        if self.bias:
 516            X_with_bias = jnp.hstack(
 517                [jnp.ones((scaled_X.shape[0], 1)), scaled_X]
 518            )
 519            print("X_with_bias shape:", X_with_bias.shape)
 520            print("W shape:", W.shape)
 521            linear_output = jnp.dot(X_with_bias, W)
 522        else:
 523            linear_output = jnp.dot(scaled_X, W)
 524
 525        # Apply activation function
 526        if self.activation_name == "relu":
 527            activated = jax.nn.relu(linear_output)
 528        elif self.activation_name == "tanh":
 529            activated = jnp.tanh(linear_output)
 530        elif self.activation_name == "sigmoid":
 531            activated = jax.nn.sigmoid(linear_output)
 532        else:  # leaky relu
 533            activated = jax.nn.leaky_relu(linear_output, negative_slope=self.a)
 534
 535        # Apply dropout
 536        if self.dropout > 0:
 537            key, subkey = jax.random.split(key)
 538            mask = jax.random.bernoulli(
 539                subkey, p=1 - self.dropout, shape=activated.shape
 540            )
 541            activated = jnp.where(mask, activated / (1 - self.dropout), 0)
 542
 543        return activated
 544
 545    def cook_training_set(self, y=None, X=None, W=None, **kwargs):
 546        """Create new hidden features for training set, with hidden layer, center the response.
 547
 548        Parameters:
 549
 550            y: array-like, shape = [n_samples]
 551                Target values
 552
 553            X: {array-like}, shape = [n_samples, n_features]
 554                Training vectors, where n_samples is the number
 555                of samples and n_features is the number of features
 556
 557            W: {array-like}, shape = [n_features, hidden_features]
 558                if provided, constructs the hidden layer via W
 559
 560        Returns:
 561
 562            (centered response, direct link + hidden layer matrix): {tuple}
 563
 564        """
 565
 566        # either X and y are stored or not
 567        # assert ((y is None) & (X is None)) | ((y is not None) & (X is not None))
 568        if self.n_hidden_features > 0:  # has a hidden layer
 569            assert (
 570                len(self.type_scaling) >= 2
 571            ), "must have len(self.type_scaling) >= 2 when self.n_hidden_features > 0"
 572
 573        if X is None:
 574
 575            if self.col_sample == 1:
 576                input_X = self.X_
 577            else:
 578                n_features = self.X_.shape[1]
 579                new_n_features = int(np.ceil(n_features * self.col_sample))
 580                assert (
 581                    new_n_features >= 1
 582                ), "check class attribute 'col_sample' and the number of covariates provided for X"
 583                np.random.seed(self.seed)
 584                index_col = np.random.choice(
 585                    range(n_features), size=new_n_features, replace=False
 586                )
 587                self.index_col_ = index_col
 588                input_X = self.X_[:, self.index_col_]
 589
 590        else:  # X is not None # keep X vs self.X_
 591
 592            if isinstance(X, pd.DataFrame):
 593                X = copy.deepcopy(X.values.astype(float))
 594
 595            if self.col_sample == 1:
 596                input_X = X
 597            else:
 598                n_features = X.shape[1]
 599                new_n_features = int(np.ceil(n_features * self.col_sample))
 600                assert (
 601                    new_n_features >= 1
 602                ), "check class attribute 'col_sample' and the number of covariates provided for X"
 603                np.random.seed(self.seed)
 604                index_col = np.random.choice(
 605                    range(n_features), size=new_n_features, replace=False
 606                )
 607                self.index_col_ = index_col
 608                input_X = X[:, self.index_col_]
 609
 610        if self.n_clusters <= 0:
 611            # data without any clustering: self.n_clusters is None -----
 612
 613            if self.n_hidden_features > 0:  # with hidden layer
 614
 615                self.nn_scaler_, scaled_X = mo.scale_covariates(
 616                    input_X, choice=self.type_scaling[1], scaler=self.nn_scaler_
 617                )
 618                Phi_X = (
 619                    self.create_layer(scaled_X)
 620                    if W is None
 621                    else self.create_layer(scaled_X, W=W)
 622                )
 623                Z = (
 624                    mo.cbind(input_X, Phi_X, backend=self.backend)
 625                    if self.direct_link is True
 626                    else Phi_X
 627                )
 628                self.scaler_, scaled_Z = mo.scale_covariates(
 629                    Z, choice=self.type_scaling[0], scaler=self.scaler_
 630                )
 631            else:  # no hidden layer
 632                Z = input_X
 633                self.scaler_, scaled_Z = mo.scale_covariates(
 634                    Z, choice=self.type_scaling[0], scaler=self.scaler_
 635                )
 636
 637        else:
 638
 639            # data with clustering: self.n_clusters is not None ----- # keep
 640
 641            augmented_X = mo.cbind(
 642                input_X,
 643                self.encode_clusters(input_X, **kwargs),
 644                backend=self.backend,
 645            )
 646
 647            if self.n_hidden_features > 0:  # with hidden layer
 648
 649                self.nn_scaler_, scaled_X = mo.scale_covariates(
 650                    augmented_X,
 651                    choice=self.type_scaling[1],
 652                    scaler=self.nn_scaler_,
 653                )
 654                Phi_X = (
 655                    self.create_layer(scaled_X)
 656                    if W is None
 657                    else self.create_layer(scaled_X, W=W)
 658                )
 659                Z = (
 660                    mo.cbind(augmented_X, Phi_X, backend=self.backend)
 661                    if self.direct_link is True
 662                    else Phi_X
 663                )
 664                self.scaler_, scaled_Z = mo.scale_covariates(
 665                    Z, choice=self.type_scaling[0], scaler=self.scaler_
 666                )
 667            else:  # no hidden layer
 668                Z = augmented_X
 669                self.scaler_, scaled_Z = mo.scale_covariates(
 670                    Z, choice=self.type_scaling[0], scaler=self.scaler_
 671                )
 672
 673        # Returning model inputs -----
 674        if mx.is_factor(y) is False:  # regression
 675            # center y
 676            if y is None:
 677                self.y_mean_, centered_y = mo.center_response(self.y_)
 678            else:
 679                self.y_mean_, centered_y = mo.center_response(y)
 680
 681            # y is subsampled
 682            if self.row_sample < 1:
 683                n, p = Z.shape
 684
 685                self.subsampler_ = (
 686                    SubSampler(
 687                        y=self.y_, row_sample=self.row_sample, seed=self.seed
 688                    )
 689                    if y is None
 690                    else SubSampler(
 691                        y=y, row_sample=self.row_sample, seed=self.seed
 692                    )
 693                )
 694
 695                self.index_row_ = self.subsampler_.subsample()
 696
 697                n_row_sample = len(self.index_row_)
 698                # regression
 699                return (
 700                    centered_y[self.index_row_].reshape(n_row_sample),
 701                    self.scaler_.transform(
 702                        Z[self.index_row_, :].reshape(n_row_sample, p)
 703                    ),
 704                )
 705            # y is not subsampled
 706            # regression
 707            return (centered_y, self.scaler_.transform(Z))
 708
 709        # classification
 710        # y is subsampled
 711        if self.row_sample < 1:
 712            n, p = Z.shape
 713
 714            self.subsampler_ = (
 715                SubSampler(
 716                    y=self.y_, row_sample=self.row_sample, seed=self.seed
 717                )
 718                if y is None
 719                else SubSampler(y=y, row_sample=self.row_sample, seed=self.seed)
 720            )
 721
 722            self.index_row_ = self.subsampler_.subsample()
 723
 724            n_row_sample = len(self.index_row_)
 725            # classification
 726            return (
 727                y[self.index_row_].reshape(n_row_sample),
 728                self.scaler_.transform(
 729                    Z[self.index_row_, :].reshape(n_row_sample, p)
 730                ),
 731            )
 732        # y is not subsampled
 733        # classification
 734        return (y, self.scaler_.transform(Z))
 735
 736    def cook_test_set(self, X, **kwargs):
 737        """Transform data from test set, with hidden layer.
 738
 739        Parameters:
 740
 741            X: {array-like}, shape = [n_samples, n_features]
 742                Training vectors, where n_samples is the number
 743                of samples and n_features is the number of features
 744
 745            **kwargs: additional parameters to be passed to self.encode_cluster
 746
 747        Returns:
 748
 749            Transformed test set : {array-like}
 750        """
 751
 752        if isinstance(X, pd.DataFrame):
 753            X = copy.deepcopy(X.values.astype(float))
 754
 755        if len(X.shape) == 1:
 756            X = X.reshape(1, -1)
 757
 758        if (
 759            self.n_clusters == 0
 760        ):  # data without clustering: self.n_clusters is None -----
 761            if self.n_hidden_features > 0:
 762                # if hidden layer
 763                scaled_X = (
 764                    self.nn_scaler_.transform(X)
 765                    if (self.col_sample == 1)
 766                    else self.nn_scaler_.transform(X[:, self.index_col_])
 767                )
 768                Phi_X = self.create_layer(scaled_X, self.W_)
 769                if self.direct_link:
 770                    return self.scaler_.transform(
 771                        mo.cbind(scaled_X, Phi_X, backend=self.backend)
 772                    )
 773                # when self.direct_link == False
 774                return self.scaler_.transform(Phi_X)
 775            # if no hidden layer # self.n_hidden_features == 0
 776            return self.scaler_.transform(X)
 777
 778        # data with clustering: self.n_clusters > 0 -----
 779        if self.col_sample == 1:
 780            predicted_clusters = self.encode_clusters(
 781                X=X, predict=True, **kwargs
 782            )
 783            augmented_X = mo.cbind(X, predicted_clusters, backend=self.backend)
 784        else:
 785            predicted_clusters = self.encode_clusters(
 786                X=X[:, self.index_col_], predict=True, **kwargs
 787            )
 788            augmented_X = mo.cbind(
 789                X[:, self.index_col_], predicted_clusters, backend=self.backend
 790            )
 791
 792        if self.n_hidden_features > 0:  # if hidden layer
 793            scaled_X = self.nn_scaler_.transform(augmented_X)
 794            Phi_X = self.create_layer(scaled_X, self.W_)
 795            if self.direct_link:
 796                return self.scaler_.transform(
 797                    mo.cbind(augmented_X, Phi_X, backend=self.backend)
 798                )
 799            return self.scaler_.transform(Phi_X)
 800
 801        # if no hidden layer
 802        return self.scaler_.transform(augmented_X)
 803
 804    def cook_training_set_jax(self, y=None, X=None, W=None, **kwargs):
 805        """JAX-compatible version of cook_training_set that maintains side effects."""
 806        # Initialize random key
 807        key = jax.random.PRNGKey(self.seed)
 808
 809        # Convert inputs to JAX arrays
 810        X = jnp.asarray(X) if X is not None else jnp.asarray(self.X_)
 811        y = jnp.asarray(y) if y is not None else jnp.asarray(self.y_)
 812
 813        # Handle column sampling
 814        if self.col_sample < 1:
 815            n_features = X.shape[1]
 816            new_n_features = int(jnp.ceil(n_features * self.col_sample))
 817            assert new_n_features >= 1, "Invalid col_sample"
 818
 819            key, subkey = jax.random.split(key)
 820            index_col = jax.random.choice(
 821                subkey, n_features, shape=(new_n_features,), replace=False
 822            )
 823            self.index_col_ = np.array(
 824                index_col
 825            )  # Store as numpy for original methods
 826            input_X = X[:, index_col]
 827            n_features = (
 828                new_n_features  # Update n_features after column sampling
 829            )
 830        else:
 831            input_X = X
 832            n_features = X.shape[1]
 833
 834        augmented_X = input_X
 835
 836        # JAX-compatible scaling
 837        def jax_scale(data, mean=None, std=None):
 838            if mean is None:
 839                mean = jnp.mean(data, axis=0)
 840            if std is None:
 841                std = jnp.std(data, axis=0)
 842            return (data - mean) / (std + 1e-10), mean, std
 843
 844        # Hidden layer processing
 845        if self.n_hidden_features > 0:
 846            # Initialize weights if not provided
 847            if W is None:
 848                shape = (n_features, self.n_hidden_features)
 849
 850                # JAX-compatible weight generation
 851                if self.nodes_sim == "uniform":
 852                    key, subkey = jax.random.split(key)
 853                    W = jax.random.uniform(
 854                        subkey, shape=shape, minval=-1.0, maxval=1.0
 855                    ) * (1 / jnp.sqrt(n_features))
 856                else:
 857                    # For other sequences, use numpy generation then convert to JAX
 858                    if self.nodes_sim == "sobol":
 859                        W_np = generate_sobol(
 860                            n_dims=shape[0],
 861                            n_points=shape[1],
 862                            seed=self.seed,
 863                        )
 864                    elif self.nodes_sim == "hammersley":
 865                        W_np = generate_hammersley(
 866                            n_dims=shape[0],
 867                            n_points=shape[1],
 868                            seed=self.seed,
 869                        )
 870                    elif self.nodes_sim == "halton":
 871                        W_np = generate_halton(
 872                            n_dims=shape[0],
 873                            n_points=shape[1],
 874                            seed=self.seed,
 875                        )
 876                    else:  # default to uniform
 877                        key, subkey = jax.random.split(key)
 878                        W = jax.random.uniform(
 879                            subkey, shape=shape, minval=-1.0, maxval=1.0
 880                        ) * (1 / jnp.sqrt(n_features))
 881
 882                    if self.nodes_sim in ["sobol", "hammersley", "halton"]:
 883                        W = jnp.asarray(W_np) * (1 / jnp.sqrt(n_features))
 884
 885                self.W_ = np.array(W)  # Store as numpy for original methods
 886
 887            # Scale features
 888            scaled_X, self.nn_mean_, self.nn_std_ = jax_scale(
 889                augmented_X,
 890                getattr(self, "nn_mean_", None),
 891                getattr(self, "nn_std_", None),
 892            )
 893
 894            # Create hidden layer with proper bias handling
 895            linear_output = jnp.dot(scaled_X, W)
 896
 897            # Apply activation
 898            if self.activation_name == "relu":
 899                Phi_X = jax.nn.relu(linear_output)
 900            elif self.activation_name == "tanh":
 901                Phi_X = jnp.tanh(linear_output)
 902            elif self.activation_name == "sigmoid":
 903                Phi_X = jax.nn.sigmoid(linear_output)
 904            else:  # leaky relu
 905                Phi_X = jax.nn.leaky_relu(linear_output, negative_slope=self.a)
 906
 907            # Apply dropout
 908            if self.dropout > 0:
 909                key, subkey = jax.random.split(key)
 910                mask = jax.random.bernoulli(
 911                    subkey, p=1 - self.dropout, shape=Phi_X.shape
 912                )
 913                Phi_X = jnp.where(mask, Phi_X / (1 - self.dropout), 0)
 914
 915            Z = jnp.hstack([scaled_X, Phi_X]) if self.direct_link else Phi_X
 916        else:
 917            Z = augmented_X
 918
 919        # Final scaling
 920        scaled_Z, self.scale_mean_, self.scale_std_ = jax_scale(
 921            Z,
 922            getattr(self, "scale_mean_", None),
 923            getattr(self, "scale_std_", None),
 924        )
 925
 926        # Center response for regression
 927        if not hasattr(mx, "is_factor") or not mx.is_factor(
 928            y
 929        ):  # regression case
 930            self.y_mean_ = float(
 931                jnp.mean(y)
 932            )  # Convert to Python float for compatibility
 933            centered_y = y - self.y_mean_
 934        else:
 935            centered_y = y
 936
 937        # Handle row sampling
 938        if self.row_sample < 1:
 939            key, subkey = jax.random.split(key)
 940            n_samples = Z.shape[0]
 941            n_row_sample = int(jnp.ceil(n_samples * self.row_sample))
 942            index_row = jax.random.choice(
 943                subkey, n_samples, shape=(n_row_sample,), replace=False
 944            )
 945            self.index_row_ = np.array(
 946                index_row
 947            )  # Store as numpy for original methods
 948            return (centered_y[index_row], scaled_Z[index_row])
 949
 950        return (centered_y, scaled_Z)
 951
 952    def cook_test_set_jax(self, X, **kwargs):
 953        """JAX-compatible test set processing with matching dimension handling."""
 954        X = jnp.asarray(X)
 955
 956        if len(X.shape) == 1:
 957            X = X.reshape(1, -1)
 958
 959        # Handle column sampling
 960        input_X = (
 961            X if self.col_sample == 1 else X[:, jnp.asarray(self.index_col_)]
 962        )
 963
 964        augmented_X = input_X
 965
 966        # JAX-compatible scaling
 967        scaled_X = (augmented_X - self.nn_mean_) / (self.nn_std_ + 1e-10)
 968
 969        # Process hidden layer if needed
 970        if self.n_hidden_features > 0:
 971            Phi_X = self._jax_create_layer(scaled_X, jnp.asarray(self.W_))
 972            Z = jnp.hstack([scaled_X, Phi_X]) if self.direct_link else Phi_X
 973        else:
 974            Z = augmented_X
 975
 976        # Final scaling
 977        scaled_Z = (Z - self.scale_mean_) / (self.scale_std_ + 1e-10)
 978
 979        return scaled_Z
 980
 981    def _jax_create_layer(self, X, W):
 982        """JAX-compatible hidden layer creation."""
 983        # print("X", X.shape)
 984        # print("W", W.shape)
 985        # print("self.W_", self.W_.shape)
 986        linear_output = jnp.dot(X, W)
 987
 988        if self.activation_name == "relu":
 989            return jax.nn.relu(linear_output)
 990        elif self.activation_name == "tanh":
 991            return jnp.tanh(linear_output)
 992        elif self.activation_name == "sigmoid":
 993            return jax.nn.sigmoid(linear_output)
 994        else:  # leaky relu
 995            return jax.nn.leaky_relu(linear_output, negative_slope=self.a)
 996
 997    def cross_val_score(
 998        self,
 999        X,
1000        y,
1001        cv=5,
1002        scoring="accuracy",
1003        random_state=42,
1004        n_jobs=-1,
1005        epsilon=0.5,
1006        penalized=True,
1007        objective="abs",
1008        **kwargs
1009    ):
1010        """
1011        Penalized Cross-validation score for a model.
1012
1013        Parameters:
1014
1015            X: {array-like}, shape = [n_samples, n_features]
1016                Training vectors, where n_samples is the number
1017                of samples and n_features is the number of features
1018
1019            y: array-like, shape = [n_samples]
1020                Target values
1021
1022            X_test: {array-like}, shape = [n_samples, n_features]
1023                Test vectors, where n_samples is the number
1024                of samples and n_features is the number of features
1025
1026            y_test: array-like, shape = [n_samples]
1027                Target values
1028
1029            cv: int
1030                Number of folds
1031
1032            scoring: str
1033                Scoring metric
1034
1035            random_state: int
1036                Random state
1037
1038            n_jobs: int
1039                Number of jobs to run in parallel
1040
1041            epsilon: float
1042                Penalty parameter
1043
1044            penalized: bool
1045                Whether to obtain penalized cross-validation score or not
1046
1047            objective: str
1048                'abs': Minimize the absolute difference between cross-validation score and validation score
1049                'relative': Minimize the relative difference between cross-validation score and validation score
1050        Returns:
1051
1052            A namedtuple with the following fields:
1053                - cv_score: float
1054                    cross-validation score
1055                - val_score: float
1056                    validation score
1057                - penalized_score: float
1058                    penalized cross-validation score: cv_score / val_score + epsilon*(1/val_score + 1/cv_score)
1059                    If higher scoring metric is better, minimize the function result.
1060                    If lower scoring metric is better, maximize the function result.
1061        """
1062        if scoring == "accuracy":
1063            scoring_func = accuracy_score
1064        elif scoring == "balanced_accuracy":
1065            scoring_func = balanced_accuracy_score
1066        elif scoring == "f1":
1067            scoring_func = f1_score
1068        elif scoring == "roc_auc":
1069            scoring_func = roc_auc_score
1070        elif scoring == "r2":
1071            scoring_func = r2_score
1072        elif scoring == "mse":
1073            scoring_func = mean_squared_error
1074        elif scoring == "mae":
1075            scoring_func = mean_absolute_error
1076        elif scoring == "mape":
1077            scoring_func = mean_absolute_percentage_error
1078        elif scoring == "rmse":
1079
1080            def scoring_func(y_true, y_pred):
1081                return np.sqrt(mean_squared_error(y_true, y_pred))
1082
1083        X_train, X_val, y_train, y_val = train_test_split(
1084            X, y, test_size=0.2, random_state=random_state
1085        )
1086
1087        res = cross_val_score(
1088            self, X_train, y_train, cv=cv, scoring=scoring, n_jobs=n_jobs
1089        )  # cross-validation error
1090
1091        if penalized == False:
1092            return res
1093
1094        DescribeResult = namedtuple(
1095            "DescribeResult", ["cv_score", "val_score", "penalized_score"]
1096        )
1097
1098        numerator = res.mean()
1099
1100        # Evaluate on the (cv+1)-th fold
1101        preds_val = self.fit(X_train, y_train).predict(X_val)
1102        try:
1103            denominator = scoring(y_val, preds_val)  # validation error
1104        except Exception as e:
1105            denominator = scoring_func(y_val, preds_val)
1106
1107        # if higher is better
1108        if objective == "abs":
1109            penalized_score = np.abs(numerator - denominator) + epsilon * (
1110                1 / denominator + 1 / numerator
1111            )
1112        elif objective == "relative":
1113            ratio = numerator / denominator
1114            penalized_score = np.abs(ratio - 1) + epsilon * (
1115                1 / denominator + 1 / numerator
1116            )
1117
1118        return DescribeResult(
1119            cv_score=numerator,
1120            val_score=denominator,
1121            penalized_score=penalized_score,
1122        )