additional code and documentation formatting.

Author: k1nshuk

athena/active.py

@@ -31,6 +31,7 @@ class ActiveSubspaces(Subspaces):
         Hristache, et al.
     :param int n_boot: number of bootstrap samples. Default is 100.
     """
+
     def __init__(self, dim, method='exact', n_boot=100):
         super().__init__(dim, method, n_boot)
 
@@ -313,13 +314,13 @@ def _hit_and_run_inactive(self, reduced_input, n_points):
             f, g = b - np.dot(A, z0), np.dot(A, d)
 
             # find an upper bound on the step
-            min_ind = np.logical_and(g <= 0,
-                                     f < -np.sqrt(np.finfo(np.float64).eps))
+            min_ind = np.logical_and(g <= 0, f
+                                     < -np.sqrt(np.finfo(np.float64).eps))
             eps_max = np.amin(f[min_ind] / g[min_ind])
 
             # find a lower bound on the step
-            max_ind = np.logical_and(g > 0,
-                                     f < -np.sqrt(np.finfo(np.float64).eps))
+            max_ind = np.logical_and(g > 0, f
+                                     < -np.sqrt(np.finfo(np.float64).eps))
             eps_min = np.amax(f[max_ind] / g[max_ind])
 
             # randomly sample eps

athena/compatibility.py

@@ -0,0 +1,177 @@
+"""
+Compatibility layer for handling different package versions.
+This module provides uniform interfaces for functionality that might
+depend on specific versions of packages or alternative implementations.
+"""
+import numpy as np
+import warnings
+from packaging import version
+
+# Check if scikit-learn-extra's KMedoids is usable
+# with the current NumPy version
+SKLEARN_EXTRA_AVAILABLE = False
+try:
+    import sklearn_extra
+    from sklearn_extra.cluster import KMedoids as SklearnExtraKMedoids
+    SKLEARN_EXTRA_AVAILABLE = True
+
+    # Check if NumPy version is compatible with sklearn_extra
+    if version.parse(np.__version__) >= version.parse('2.0.0'):
+        warnings.warn(
+            "You are using NumPy >= 2.0.0 with scikit-learn-extra which may "
+            "cause compatibility issues. If you encounter errors, consider "
+            "using the built-in KMedoids implementation in ATHENA.")
+except ImportError:
+    SklearnExtraKMedoids = None
+
+
+# Implementation based on scikit-learn's KMeans but adapted for KMedoids
+class KMedoids:
+    """
+    K-Medoids clustering.
+    
+    A custom implementation that doesn't rely on scikit-learn-extra, thus
+    ensuring compatibility with NumPy 2.0+.
+    
+    Parameters
+    ----------
+    n_clusters : int, default=8
+        The number of clusters to form as well as the number of medoids to generate.
+    
+    init : {'k-medoids++', 'random'} or array of shape (n_clusters, n_features), default='k-medoids++'
+        Method for initialization.
+    
+    max_iter : int, default=300
+        Maximum number of iterations of the k-medoids algorithm for a single run.
+    
+    random_state : int, RandomState instance or None, default=None
+        Determines random number generation for centroid initialization.
+    """
+
+    def __init__(self,
+                 n_clusters=8,
+                 init='k-medoids++',
+                 max_iter=300,
+                 random_state=None):
+        self.n_clusters = n_clusters
+        self.init = init
+        self.max_iter = max_iter
+        self.random_state = random_state
+        self.cluster_centers_ = None
+        self.labels_ = None
+        self.inertia_ = None
+        self.n_iter_ = 0
+
+    def _init_medoids(self, X):
+        """Initialize the medoids."""
+        rng = np.random.RandomState(self.random_state)
+        n_samples = X.shape[0]
+
+        if isinstance(self.init, str) and self.init == 'random':
+            # Random selection
+            indices = rng.permutation(n_samples)[:self.n_clusters]
+            self.cluster_centers_ = X[indices].copy()
+        elif isinstance(self.init, str) and self.init == 'k-medoids++':
+            # Implementation of k-medoids++ initialization
+            # Choose the first medoid randomly
+            indices = np.zeros(self.n_clusters, dtype=int)
+            indices[0] = rng.randint(n_samples)
+
+            # Calculate distances to the first medoid
+            distances = np.sum((X - X[indices[0]])**2, axis=1)
+
+            # Choose remaining medoids
+            for i in range(1, self.n_clusters):
+                # Choose point with probability proportional to distance squared
+                probs = distances / np.sum(distances)
+                indices[i] = rng.choice(n_samples, p=probs)
+
+                # Update distances
+                new_dist = np.sum((X - X[indices[i]])**2, axis=1)
+                distances = np.minimum(distances, new_dist)
+
+            self.cluster_centers_ = X[indices].copy()
+        else:
+            # Use provided initial medoids
+            self.cluster_centers_ = np.asarray(self.init, dtype=X.dtype)
+
+    def fit(self, X):
+        """Compute k-medoids clustering."""
+        X = np.asarray(X)
+        self._init_medoids(X)
+
+        best_labels = None
+        best_inertia = float('inf')
+        best_centers = None
+
+        for i in range(self.max_iter):
+            # Assign each point to closest medoid
+            distances = np.zeros((X.shape[0], self.n_clusters))
+            for j in range(self.n_clusters):
+                distances[:, j] = np.sum((X - self.cluster_centers_[j])**2,
+                                         axis=1)
+
+            labels = np.argmin(distances, axis=1)
+
+            # Update medoids
+            old_centers = self.cluster_centers_.copy()
+
+            # For each cluster, update medoid to be the point minimizing inertia
+            for j in range(self.n_clusters):
+                cluster_points = X[labels == j]
+                if len(cluster_points) > 0:
+                    # Compute pairwise distances within cluster
+                    inertias = np.zeros(len(cluster_points))
+                    for k, point in enumerate(cluster_points):
+                        inertias[k] = np.sum(
+                            np.sum((cluster_points - point)**2, axis=1))
+
+                    # Choose point with minimal inertia as new medoid
+                    min_idx = np.argmin(inertias)
+                    self.cluster_centers_[j] = cluster_points[min_idx].copy()
+
+            # Compute inertia
+            inertia = 0
+            for j in range(self.n_clusters):
+                cluster_points = X[labels == j]
+                if len(cluster_points) > 0:
+                    inertia += np.sum(
+                        np.sum((cluster_points - self.cluster_centers_[j])**2,
+                               axis=1))
+
+            # Store best result
+            if inertia < best_inertia:
+                best_inertia = inertia
+                best_labels = labels
+                best_centers = self.cluster_centers_.copy()
+
+            # Check for convergence
+            center_shift = np.sum(
+                np.sqrt(np.sum((old_centers - self.cluster_centers_)**2,
+                               axis=1)))
+            if center_shift < 1e-4:
+                break
+
+        self.labels_ = best_labels
+        self.cluster_centers_ = best_centers
+        self.inertia_ = best_inertia
+        self.n_iter_ = i + 1
+
+        return self
+
+    def predict(self, X):
+        """Predict the closest cluster for each sample in X."""
+        X = np.asarray(X)
+        distances = np.zeros((X.shape[0], self.n_clusters))
+        for j in range(self.n_clusters):
+            distances[:, j] = np.sum((X - self.cluster_centers_[j])**2, axis=1)
+
+        return np.argmin(distances, axis=1)
+
+
+# Export the appropriate KMedoids implementation
+if SKLEARN_EXTRA_AVAILABLE and version.parse(
+        np.__version__) < version.parse('2.0.0'):
+    # Use sklearn-extra's implementation when available and NumPy < 2.0
+    KMedoids = SklearnExtraKMedoids
+# Otherwise use our implementation which is compatible with NumPy 2.0+

athena/feature_map.py

@@ -39,6 +39,7 @@ class FeatureMap():
 
     :raises TypeError
     """
+
     def __init__(self, distr, bias, input_dim, n_features, params, sigma_f):
         if callable(distr):
             self.distr = distr
@@ -202,9 +203,11 @@ def tune_pr_matrix(self,
             # Reformat bounds for BayesianOptimization package format
             # BayesianOptimization uses a dictionary of parameter names and their range tuples
             # Unlike GPyOpt which used a list of dictionaries with 'name', 'type', and 'domain' keys
-            bounds_dict = {f'var_{i}': (bound.start, bound.stop) 
-                         for i, bound in enumerate(bounds)}
-            
+            bounds_dict = {
+                f'var_{i}': (bound.start, bound.stop)
+                for i, bound in enumerate(bounds)
+            }
+
             # Create wrapper for the objective function to handle the format difference
             # BayesianOptimization passes parameters as keyword arguments, not as an array
             def bayes_wrapper(**kwargs):
@@ -213,25 +216,24 @@ def bayes_wrapper(**kwargs):
                 # BayesianOptimization maximizes functions by default, but we want to minimize
                 # So we negate the score (lower scores are better in our original function)
                 return -func(x, best, **fn_args)
-            
+
             # Initialize optimizer with our wrapper function and parameter bounds
             optimizer = BayesianOptimization(
                 f=bayes_wrapper,
                 pbounds=bounds_dict,
                 random_state=42  # For reproducible results
             )
-            
+
             # Run optimization
             # init_points: how many steps of random exploration to perform
             # n_iter: how many steps of bayesian optimization to perform
-            optimizer.maximize(
-                init_points=2,
-                n_iter=maxiter
-            )
-            
+            optimizer.maximize(init_points=2, n_iter=maxiter)
+
             # Extract the best parameters found and transform back
             # optimizer.max contains the best score and parameters found
-            best_params = [optimizer.max['params'][f'var_{i}'] for i in range(len(bounds))]
+            best_params = [
+                optimizer.max['params'][f'var_{i}'] for i in range(len(bounds))
+            ]
             # Apply 10^ transformation as done in the original implementation
             self.params = 10**np.array(best_params)
         else:

athena/kas.py

@@ -55,6 +55,7 @@ class KernelActiveSubspaces(Subspaces):
     :cvar numpy.ndarray metric: metric matrix for vectorial active
         subspaces.
     """
+
     def __init__(self,
                  dim,
                  feature_map=None,

athena/local.py

@@ -36,6 +36,7 @@ class MaximumASDimensionReached(Exception):
 class ClusterBase():
     """Local Active Subspaces clustering Base class.
     """
+
     def __init__(self):
 
         self.inputs = None
@@ -210,6 +211,7 @@ def plot_clusters(self, save=False, title='2d_clusters', plot=True):
 
 class KMeansAS(ClusterBase):
     """Clustering with k-means"""
+
     def __init__(self):
         super().__init__()
         self.centers = None
@@ -231,6 +233,7 @@ def _fit_clustering(self):
 
 class KMedoidsAS(ClusterBase):
     """Clustering with k-medoids"""
+
     def __init__(self):
         super().__init__()
         self.centers = None
@@ -257,6 +260,7 @@ def as_metric(self, X, Y):
 
 
 class TopDownHierarchicalAS(ClusterBase):
+
     def __init__(self):
         """TODO check states logic.
         1. 2 and 4 are exclusives
@@ -402,7 +406,9 @@ def _fit_clustering(self, print_states=False, plot=False):
 
     def refine_one_step(self):
         """Increase the dimension of the Active Subspace once, when possible."""
+
         class LeafUpdate(object):
+
             def __init__(self):
                 self.score = 0
                 self.leaves_list = []
@@ -437,6 +443,7 @@ def refine_further(self, minimum_score, plot=False):
         print("Start refining: increasing the as dimension when possible.")
 
         class CallRefine(object):
+
             def __init__(self, minimum_score):
                 self.min = minimum_score
 
@@ -474,7 +481,9 @@ def _print_state_debug(self):
 
     def _print_leaves_score(self):
         """Print the information of every leaf."""
+
         class ComputeScore(object):
+
             def __init__(self):
                 self.n_leaves = 0
                 self.leaves_dim = []
@@ -495,7 +504,9 @@ def __call__(self, node):
 
     def assign_leaf_labels(self):
         """Assign integer labels to the leaves."""
+
         class LeafLabels(object):
+
             def __init__(self):
                 self.labels_counter = 0
 
@@ -508,6 +519,7 @@ def __call__(self, node):
 
     def reset_gprs(self):
         """Reset the GPRs of every leaf and root."""
+
         def reset_gpr(node):
             node.gpr = None
             node.ss = None
@@ -519,7 +531,9 @@ def plot_clusters(self,
                       save_data=True,
                       plot=True,
                       save=True):
+
         class SaveLeafInfo(object):
+
             def __init__(self):
                 self.n_leaves = 0
                 self.n_elems = []
@@ -691,6 +705,7 @@ def void_func(*args, **kwargs):
 
 
 class TopDownNode():
+
     def __init__(self, parent, node_indexes, val_indexes, tree_obj):
         """A TopDownNode is defined by the indexes of the triplets (inputs,
         outputs, gradients) of the training data and the parent node. The root
@@ -822,6 +837,7 @@ def refine_further(self, minimum_score):
     class NormalizeDivisive():
         """Inner class for normalization of inputs, gradients w.r.t. local
         clusters"""
+
         def __init__(self, norm_type, ind, inputs):
             self.type = norm_type
 
@@ -978,7 +994,8 @@ def refine_cluster(self):
                 return state, self.children
 
             # check if clustering is possible
-            if self.ind.shape[0] < self.hierarchical.total_clusters + n_clusters:
+            if self.ind.shape[
+                    0] < self.hierarchical.total_clusters + n_clusters:
                 state.add(5)
                 _log.debug("Refine returns 5 : " + str(state) +
                            " and list length " + str(len(self.children)))

athena/local_classification.py

@@ -33,6 +33,7 @@
 class SpectralClassification(metaclass=abc.ABCMeta):
     """Evaluate the connected components from X, n_neighbours, features and custom
     distance that must be defined in concrete class."""
+
     def __init__(self):
         self.X = None
         self.features = None
@@ -158,6 +159,7 @@ class ClassifyAS(SpectralClassification):
     the AS dimension of the n_neighbours neighbouring samples with a resampling
     of neighbour_resampling. The local_as_criterion can be 'min' or 'average'
     over the batches of neighbouring samples."""
+
     def __init__(self):
         super().__init__()