From 832b6655b47e9ef8d301f059b91c2aa9f6a5815a Mon Sep 17 00:00:00 2001 From: fromor Date: Thu, 19 May 2022 11:48:11 +0200 Subject: [PATCH 01/26] comment tutorial 08, add local tests --- athena/__init__.py | 2 +- athena/local.py | 2 +- athena/local_classification.py | 6 +- tests/test_local.py | 41 +++++++++ tutorials/tutorial08/08_local_AS_2d.ipynb | 100 +++++++++++++++++----- tutorials/tutorial08/08_local_AS_2d.py | 35 ++++---- 6 files changed, 141 insertions(+), 45 deletions(-) create mode 100644 tests/test_local.py diff --git a/athena/__init__.py b/athena/__init__.py index 1f7e0c3..0875fe5 100644 --- a/athena/__init__.py +++ b/athena/__init__.py @@ -15,6 +15,6 @@ from .utils import (Normalizer, initialize_weights, linear_program_ineq, local_linear_gradients, sort_eigpairs, CrossValidation, rrmse, average_rrmse) -from .local import (ClusterBase, KMeansAS, KMedoidsAS, TopDownHierarchical, +from .local import (ClusterBase, KMeansAS, KMedoidsAS, TopDownHierarchicalAS, TopDownNode, plot_scores) from .local_classification import (SpectralClassification, ClassifyAS) diff --git a/athena/local.py b/athena/local.py index b4dd2e1..4d2003d 100644 --- a/athena/local.py +++ b/athena/local.py @@ -255,7 +255,7 @@ def as_metric(self, X, Y): self.full_as.evects.T.dot(X - Y) * self.full_as.evals) -class TopDownHierarchical(ClusterBase): +class TopDownHierarchicalAS(ClusterBase): def __init__(self): """TODO check states logic. 1. 2 and 4 are exclusives diff --git a/athena/local_classification.py b/athena/local_classification.py index b4d0ff8..df30613 100644 --- a/athena/local_classification.py +++ b/athena/local_classification.py @@ -31,7 +31,7 @@ class SpectralClassification(metaclass=abc.ABCMeta): - """Evaluate connected components from X, n_neighbours, features and custom + """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 @@ -202,12 +202,12 @@ def distance(x, y): # if two nodes have different AS dimension the distance is inf if np.abs(x[-1] - y[-1]) != 0: return np.inf - # when two nodes have the same AS dimension the distance is euclidean one + # when two nodes have the same AS dimension the distance is the euclidean one else: return np.linalg.norm(x[:-1] - y[:-1]) def evaluate_minimum_distance(self, X, dX): - """Evaluate local acive subspace dimension depending on the choice of + """Evaluate local active subspace dimension depending on the choice of the number of neighbouring points and of subsampling from the same neighbour in order to evaluate the local active subspace""" diff --git a/tests/test_local.py b/tests/test_local.py new file mode 100644 index 0000000..2d5a61e --- /dev/null +++ b/tests/test_local.py @@ -0,0 +1,41 @@ +from unittest import TestCase +import numpy as np +from athena.utils import Normalizer +from athena.local import TopDownHierarchicalAS, KMeansAS, KMedoidsAS, plot_scores +from athena.local_classification import ClassifyAS +from contextlib import contextmanager +import matplotlib.pyplot as plt + +@contextmanager +def assert_plot_figures_added(): + """ + Assert that the number of figures is higher than when you started the test + """ + num_figures_before = plt.gcf().number + yield + num_figures_after = plt.gcf().number + assert num_figures_before < num_figures_after + +class TestLocalAS(TestCase): + def test_init_local_AS(self): + las = TopDownHierarchicalAS() + self.assertIsNone(las.inputs) + self.assertIsNone(las.outputs) + self.assertIsNone(las.gradients) + self.assertIsNone(las.inputs_test) + self.assertIsNone(las.outputs_test) + self.assertIsNone(las.inputs_val) + self.assertIsNone(las.outputs_val) + + self.assertIsNone(las.outputs_dim) + self.assertIsNone(las.method) + self.assertIsNone(las.as_dim) + self.assertIsNone(las.clustering) + self.assertIsNone(las.labels) + self.assertIsNone(las.unique_labels) + self.assertIsNone(las.full_as) + self.assertIsNone(las.full_gpr) + self.assertIsNone(las.local_ass) + self.assertIsNone(las.local_gprs) + self.assertIsNone(las.max_clusters) + self.assertIsNone(las.random_state) \ No newline at end of file diff --git a/tutorials/tutorial08/08_local_AS_2d.ipynb b/tutorials/tutorial08/08_local_AS_2d.ipynb index 8c40699..cd00e49 100644 --- a/tutorials/tutorial08/08_local_AS_2d.ipynb +++ b/tutorials/tutorial08/08_local_AS_2d.ipynb @@ -11,7 +11,7 @@ "parameter space [Romor, Francesco, Marco Tezzele, and Gianluigi Rozza. \"A local\n", "approach\n", "to parameter space reduction for regression and classification\n", - "tasks.\".](https://arxiv.org/pdf/2107.10867.pdf).\n", + "tasks.\"](https://arxiv.org/pdf/2107.10867.pdf).\n", "\n", "Let $\\mathbf{X}:(\\Omega, \\mathcal{A},\n", "P)\\rightarrow\\mathcal{X}\\subset\\mathbb{R}^{n},\\qquad \\frac{d\\mathbf{X}_{\\#}(P)}{d\\lambda}={\\boldsymbol{\\mu}}$ be the inputs continuous random\n", @@ -24,7 +24,7 @@ " \\boldsymbol{\\mu}):\\mathcal{X}\\subset\\mathbb{R}^{n}\\rightarrow\\mathbb{R}$\n", " that is defined locally for every $S_i\\in\\mathcal{P}$ as\n", " $$\n", - " \\left.g\\right|_{S_i}\\mathbf{X} = \\mathbb{E}_{\\boldsymbol{\\mu}_{i}}[f|P_{r(S_i), i}]\\circ\\mathbf{X}\n", + " \\left.g\\right|_{S_i} = \\mathbb{E}_{\\boldsymbol{\\mu}_{i}}[f|P_{r(S_i), i}]\n", " $$\n", " where\n", " $\\boldsymbol{\\mu}_{i}:=(1/\\boldsymbol{\\mu}({S_{i}}))\\cdot\\left.\\boldsymbol{\\mu}\\right|_{S_{i}}\\mathbb{R}^n$,\n", @@ -50,7 +50,7 @@ "import numpy as np\n", "import matplotlib.pyplot as plt\n", "\n", - "from athena import Normalizer, TopDownHierarchical, KMeansAS, KMedoidsAS, ClassifyAS, plot_scores" + "from athena import Normalizer, TopDownHierarchicalAS, KMeansAS, KMedoidsAS, ClassifyAS, plot_scores" ] }, { @@ -99,6 +99,21 @@ "outputs_test = quartic_2d(inputs_test, scaler)[0]" ] }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "This is the dictionary that the class *TopDownHierarchicalAS* needs to perform a\n", + "hierarchical top-down clustering. The metric we employ is the AS metric:\n", + "$$\n", + "\\| x_i - x_j \\|_{\\Lambda} = \\sqrt{(x_i - x_j)^T \\mathbf{W} \\Lambda^2\n", + " \\mathbf{W}^T (x_i - x_j)}\n", + "$$\n", + "it will cluster the training samples along the AS as direction of preference. In\n", + "the following by *score* it is meant the $R^2$ score, computed on a validation\n", + "set if passed, otherwise on the passed test set." + ] + }, { "cell_type": "code", "execution_count": 3, @@ -108,27 +123,40 @@ "config = {\n", " 'inputs_test': inputs_test,\n", " 'outputs_test': outputs_test,\n", - " 'method': 'kmedoids',\n", - " 'n_clusters': 3,\n", - " 'as_dim': 1,\n", - " 'random_state': 0,\n", - " 'max_clusters': 3,\n", - " 'max_red_dim': 2,\n", + " 'inputs_val': None,\n", + " 'outputs_val': None,\n", + " 'method': 'kmedoids', # methods used for clustering at each refinement level\n", + " 'max_clusters': 3, # number of clusters\n", + " 'as_dim': 1, # global AS dimension\n", + " 'random_state': 0, # random state used for clustering with scipy\n", + " 'max_red_dim': 2, # maximum dimension of local AS\n", " 'max_children': 3,\n", " 'min_children': 3,\n", - " 'min_local': 10,\n", + " 'min_local': 10, # minimum number of elements in each local cluster\n", " 'dim_cut': 0.7,\n", - " 'normalization': 'uniform',\n", - " 'metric': 'as',\n", - " 'refinement_criterion': 'force',\n", - " 'as_dim_criterion': 'residual',\n", - " 'minimum_score': 0.98,\n", - " 'score_tolerance': 0.99,\n", - " 'max_dim_refine_further':2,\n", - " 'verbose': True,\n", + " 'normalization': 'uniform', # normalization of the samples of each cluster\n", + " 'metric': 'as', # metric AS-informed defined above\n", + " 'refinement_criterion': 'force', # when possible, split the clusters\n", + " 'as_dim_criterion': 'residual', # choose the AS dimension minimizing the residual energy\n", + " 'score_tolerance': 0.99, # if the score is above or equal, stop refining the node/local cluster\n", + " 'minimum_score': 0.98,# if the local scores are not below this tolerance, the AS dimension is incerased further until it reaches 'max_dim_refine_further'. Pass 'None' or a value > 1, to avoid further refinements of the result.\n", + " 'max_dim_refine_further':2, # maximum AS dimension reached with the final increasing of the local AS dimensions\n", + " 'verbose': True, # If true, it prints the final state (from 1 to 7) of refinement of the leaves, the number of leaves (final local clusters), the scores on the test set.\n", "}" ] }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "First hierarchical top-down clustering is performed. Afterwards, the dimension\n", + "of the local clusters is raised in case the local score is below the\n", + "'score_tolerance' provided: in this case only the dimension of the first cluster\n", + "is raised from 1 to 2. The 'Test scores' are the $R^2$ scores of the response\n", + "surfaces built with the global AS and local AS and the mean absolute regression\n", + "error corresponding to the global AS and the local AS." + ] + }, { "cell_type": "code", "execution_count": 4, @@ -168,12 +196,32 @@ } ], "source": [ - "# TopDownHierarchical clusters\n", - "cluster_as = TopDownHierarchical()\n", + "# TopDownHierarchicalAS clusters\n", + "cluster_as = TopDownHierarchicalAS()\n", "cluster_as.fit(inputs, outputs, gradients, config)\n", "cluster_as.plot_clusters(with_test=True)" ] }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Returns the number of classes and the labels of the inputs related to a\n", + "classification procedure that groups samples based on a local AS dimension:\n", + "1. for each input a local AS dimension is assigned: we take the min AS dimensions of a number equal to 'neighbour_resampling' points among\n", + " 'n_neighbours'. So in this case we take the minimum of $\\binom 6 5$ AS\n", + " dimensions. Another option is assigning to each sample the local AS itself.\n", + "2. a distance matrix is computed using the previous information. If in the\n", + " previous step the AS itself was assigned, then a Grassmanian distance could\n", + " be employed.\n", + "3. a graph is built, imposing that the inputs with the same local AS dimension belong\n", + " to the same connected components\n", + "4. the connected components of the graph are evlauated and returned as the\n", + " classes of the procedure\n", + "\n", + "Finally we employ a multilayer perceptron to get the decision boundaries." + ] + }, { "cell_type": "code", "execution_count": 5, @@ -229,6 +277,16 @@ "print(\"Classification error on train set {}\".format(score))" ] }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Finally, perform a comparison of hierarchical top-down with a single global\n", + "clustering with k-means or k-medoids (with AS-informed metric). For consistency\n", + "of the implementations, hierarchical top-down clustering with only one refinement\n", + "layer should have the same results as k-medoids (if the clustering method chosen is k-medoids)." + ] + }, { "cell_type": "code", "execution_count": 11, @@ -281,7 +339,7 @@ " config['max_clusters'] = n_clusters\n", " config['min_children'] = n_clusters\n", " config['max_children'] = n_clusters\n", - " cluster_as = TopDownHierarchical()\n", + " cluster_as = TopDownHierarchicalAS()\n", " \n", " cluster_as.fit(inputs, outputs, gradients, config)\n", " score = cluster_as.compute_scores(inputs_test, outputs_test)\n", diff --git a/tutorials/tutorial08/08_local_AS_2d.py b/tutorials/tutorial08/08_local_AS_2d.py index 34deaf6..b887263 100644 --- a/tutorials/tutorial08/08_local_AS_2d.py +++ b/tutorials/tutorial08/08_local_AS_2d.py @@ -2,29 +2,29 @@ import numpy as np import matplotlib.pyplot as plt -from athena import Normalizer, TopDownHierarchical, KMeansAS, KMedoidsAS, ClassifyAS, plot_scores +from athena import Normalizer, TopDownHierarchicalAS, KMeansAS, KMedoidsAS, ClassifyAS, plot_scores # import logging # logging.basicConfig(filename='divisive.log', level=logging.DEBUG) -def plot_scores_(possible_clusters, config): +def plot_scores_(possible_clusters, inputs, outputs, gradients, config): """ method='kmeans' or 'kmedoids' or 'divisive' """ scores = [] for n_clusters in possible_clusters: if config['method'] == 'kmeans': - config['n_clusters'] = n_clusters - cluster_as = KMeansAS(config=config) + config['max_clusters'] = n_clusters + cluster_as = KMeansAS() elif config['method'] == 'kmedoids': - config['n_clusters'] = n_clusters - cluster_as = KMedoidsAS(config=config) + config['max_clusters'] = n_clusters + cluster_as = KMedoidsAS() elif config['method'] == 'top-down': config['max_clusters'] = n_clusters config['min_children'] = n_clusters config['max_children'] = n_clusters - cluster_as = TopDownHierarchical(config=config) - cluster_as.fit() + cluster_as = TopDownHierarchicalAS() + cluster_as.fit(inputs, outputs, gradients, config) score = cluster_as.compute_scores(inputs_test, outputs_test) scores.append(score) scores = np.asarray(scores) @@ -85,9 +85,6 @@ def quartic_2d(X, n): 'max_dim_refine_further':2, 'inputs_test': inputs_test, 'outputs_test': outputs_test, - 'inputs': inputs, - 'outputs': outputs, - 'gradients': gradients, } plot_config = { @@ -102,18 +99,18 @@ def quartic_2d(X, n): 'filename': 'quartic_2d_r2_top-down_2_10.pdf' } -# TopDownHierarchical clusters -cluster_as = TopDownHierarchical(config=config) -cluster_as.fit() +# TopDownHierarchicalAS clusters +cluster_as = TopDownHierarchicalAS() +cluster_as.fit(inputs, outputs, gradients, config) cluster_as.plot_clusters(with_test=True) # classification with local AS -agg = ClassifyAS(inputs, +agg = ClassifyAS() +n_c, labels = agg.fit(inputs, gradients, n_neighbours=6, threshold=0.9999, neighbour_resampling=5) -n_c, labels = agg.fit() print("Number of components: ", n_c) agg.plot() score = agg.plot_decision_boundaries() @@ -127,14 +124,14 @@ def quartic_2d(X, n): # k-means config['method'] = 'kmeans' plot_config['filename'] = 'quartic_2d_r2_kmeans_{}_{}.pdf'.format(cl_min, cl_max) -plot_scores_(possible_clusters, config) +plot_scores_(possible_clusters, inputs, outputs, gradients, config) # k-medoids config['method'] = 'kmedoids' plot_config['filename'] = 'quartic_2d_r2_kmedoids_{}_{}.pdf'.format(cl_min, cl_max) -plot_scores_(possible_clusters, config) +plot_scores_(possible_clusters, inputs, outputs, gradients, config) # top-down clustering config['method'] = 'top-down' plot_config['filename'] = 'quartic_2d_r2_top-down_{}_{}.pdf'.format(cl_min, cl_max) -plot_scores_(possible_clusters, config) \ No newline at end of file +plot_scores_(possible_clusters, inputs, outputs, gradients, config) \ No newline at end of file From 2e965fed36eaf1d41dd3b3339e1bed48dea21c8b Mon Sep 17 00:00:00 2001 From: fromor Date: Tue, 21 Jun 2022 16:11:41 +0200 Subject: [PATCH 02/26] add evaluation of AS dimension based on spectral gap or residual energy. --- athena/meta.py | 13 ++------ athena/nll.py | 10 +++--- athena/subspaces.py | 50 ++++++++++++++++++++++++++---- athena/utils.py | 7 +++-- tests/test_active.py | 65 +++++++++++++++++++++------------------ tests/test_feature_map.py | 6 ++-- tests/test_kas.py | 9 ++---- tests/test_subspaces.py | 18 +++++++++++ 8 files changed, 115 insertions(+), 63 deletions(-) diff --git a/athena/meta.py b/athena/meta.py index 470d2fd..882b268 100644 --- a/athena/meta.py +++ b/athena/meta.py @@ -2,16 +2,9 @@ Package metadata """ __all__ = [ - '__project__', - '__title__', - '__author__', - '__copyright__', - '__license__', - '__version__', - '__mail__', - '__maintainer__', - '__status__', - '__packagename__'] + '__project__', '__title__', '__author__', '__copyright__', '__license__', + '__version__', '__mail__', '__maintainer__', '__status__', '__packagename__' +] __project__ = 'ATHENA' __title__ = "athena" diff --git a/athena/nll.py b/athena/nll.py index 5cf1297..a77627c 100644 --- a/athena/nll.py +++ b/athena/nll.py @@ -93,6 +93,9 @@ def train(self, :raises: ValueError: in interactive mode outputs must be provided for the sufficient summary plot. """ + if inputs.shape[1] % 2 != 0: + raise ValueError('The parameter space\'s dimension must be even.') + if interactive: if outputs is None: raise ValueError( @@ -452,12 +455,11 @@ def customized_loss(self, inputs, mapped_inputs, gradients): for k in range(2 * self.n_params): Jacob[:, j, k] = torch.add(dx[:, k], -1 * inputs[:, k]) - JJ2 = torch.unsqueeze( - torch.sqrt(torch.sum(torch.mul(Jacob, Jacob), 2)), 2) + JJ2 = torch.unsqueeze(torch.sqrt(torch.sum(torch.mul(Jacob, Jacob), 2)), + 2) JJJ = torch.div(Jacob, JJ2.expand(-1, -1, 2 * self.n_params)) ex_data = torch.unsqueeze(gradients, 2) - loss_weights = torch.clone(torch.squeeze(torch.matmul(JJJ, ex_data), - 2)) + loss_weights = torch.clone(torch.squeeze(torch.matmul(JJJ, ex_data), 2)) # anisotropy weigths loss_weights[:, self.omega] = 0.0 loss_anisotropy = torch.sqrt( diff --git a/athena/subspaces.py b/athena/subspaces.py index 37b11cf..f8c3216 100644 --- a/athena/subspaces.py +++ b/athena/subspaces.py @@ -16,6 +16,8 @@ import matplotlib.pyplot as plt from sklearn.utils.extmath import randomized_svd from .utils import sort_eigpairs +import itertools +import operator plt.rcParams.update({'font.size': 16}) @@ -139,8 +141,7 @@ def _compute_bootstrap_ranges(self, gradients, weights, metric=None): sub_dist = np.zeros((range_dim - 1, self.n_boot)) for i in range(self.n_boot): - gradients0, weights0 = self._bootstrap_replicate( - gradients, weights) + gradients0, weights0 = self._bootstrap_replicate(gradients, weights) e0, W0 = self._build_decompose_cov_matrix(gradients=gradients0, weights=weights0, metric=metric) @@ -250,6 +251,43 @@ def _partition(self): 'the eigenvectors cannot have dimension less than dim = {}.'. format(self.dim)) + def partition_spectral_gap(self): + """ + Partition the eigenvectors to define the active and inactive subspaces, + based on the highest spectral gap of the ordered eigenvalues. + + :raises: TypeError + """ + if self.evals is None: + raise TypeError( + 'The method fit has to be called first in order to compute the eigenvalues.' + ) + + spectral_gap_index = np.argmax(self.evals[:-1] - self.evals[1:]) + self.dim = int(spectral_gap_index + 1) + self._partition() + + def partition_residual_energy(self, tol=0.99): + """ + Partition the eigenvectors to define the active and inactive subspaces, + based on the residual energy normalized to 1 i.e. the normalized sum of the eigenvalues + corresponding to the active components must be lower than tol. + + :param float tol: threshold for the residual energy. + + :raises: TypeError, ValueError + """ + if self.evals is None: + raise ValueError( + 'The method fit has to be called first in order to compute the eigenvalues.' + ) + + residual_energies = np.array( + list(itertools.accumulate(self.evals, operator.add))) / np.sum( + self.evals) + self.dim = residual_energies[residual_energies <= tol].shape[0] + self._partition() + def plot_eigenvalues(self, n_evals=None, filename=None, @@ -270,7 +308,7 @@ def plot_eigenvalues(self, """ if self.evals is None: raise TypeError('The eigenvalues have not been computed.' - 'You have to perform the fit method.') + 'You have to perform the fit method.') if n_evals is None: n_evals = self.evals.shape[0] if n_evals > self.evals.shape[0]: @@ -320,7 +358,7 @@ def plot_eigenvalues(self, plt.grid(linestyle='dotted') plt.tight_layout() - + if filename: plt.savefig(filename) else: @@ -347,7 +385,7 @@ def plot_eigenvectors(self, """ if self.evects is None: raise TypeError('The eigenvectors have not been computed.' - 'You have to perform the fit method.') + 'You have to perform the fit method.') if n_evects is None: n_evects = self.dim if n_evects > self.evects.shape[0]: @@ -413,7 +451,7 @@ def plot_sufficient_summary(self, """ if self.evects is None: raise TypeError('The eigenvectors have not been computed.' - 'You have to perform the fit method.') + 'You have to perform the fit method.') plt.figure(figsize=figsize) plt.title(title) diff --git a/athena/utils.py b/athena/utils.py index 78e8547..1ac98ae 100644 --- a/athena/utils.py +++ b/athena/utils.py @@ -5,6 +5,7 @@ from scipy.optimize import linprog import GPy + class Normalizer(): """A class for normalizing and unnormalizing bounded inputs. @@ -209,8 +210,7 @@ def __init__(self, inputs, outputs, gradients, subspace, folds=5, **kwargs): if any([v is None for v in [inputs, outputs, gradients, subspace]]): raise ValueError( - 'Any among inputs, outputs, gradients, subspace is None.' - ) + 'Any among inputs, outputs, gradients, subspace is None.') self.inputs = inputs self.outputs = outputs @@ -344,7 +344,8 @@ def average_rrmse(hyperparams, best, csv, verbose=False, resample=5): :rtype: numpy.float64 """ if isinstance(csv, CrossValidation) is False: - raise ValueError("The argument csv must be of type athena.utils.CrossValidation") + raise ValueError( + "The argument csv must be of type athena.utils.CrossValidation") if len(hyperparams.shape) > 1: hyperparams = np.squeeze(hyperparams) diff --git a/tests/test_active.py b/tests/test_active.py index 4f66bf9..d807de2 100644 --- a/tests/test_active.py +++ b/tests/test_active.py @@ -4,6 +4,7 @@ from contextlib import contextmanager import matplotlib.pyplot as plt + @contextmanager def assert_plot_figures_added(): """ @@ -156,10 +157,10 @@ def test_fit_10(self): def test_fit_11(self): np.random.seed(42) grad = np.array( - [[-0.50183952, 0, 0, 0, 0, 0 ,0, 0, 0, 0, 0, 0, 0, 0, 0], - [-1.26638196, 0, 0, 0, 0, 0 ,0, 0, 0, 0, 0, 0, 0, 0, 0], - [ 0.43017941, 0, 0, 0, 0, 0 ,0, 0, 0, 0, 0, 0, 0, 0, 0], - [ 0.65008914, 0, 0, 0, 0, 0 ,0, 0, 0, 0, 0, 0, 0, 0, 0]]) + [[-0.50183952, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], + [-1.26638196, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], + [0.43017941, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], + [0.65008914, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]]) gradients = (grad[i, :] for i in range(4)) ss = ActiveSubspaces(dim=2, method='exact', n_boot=150) ss.fit(gradients=gradients) @@ -169,25 +170,27 @@ def test_fit_11(self): def test_fit_12(self): np.random.seed(42) grad = np.array( - [[-0.50183952, 0, 0, 0, 0, 0 ,0, 0, 0, 0, 0, 0, 0, 0, 0], - [-1.26638196, 0, 0, 0, 0, 0 ,0, 0, 0, 0, 0, 0, 0, 0, 0], - [ 0.43017941, 0, 0, 0, 0, 0 ,0, 0, 0, 0, 0, 0, 0, 0, 0], - [ 0.65008914, 0, 0, 0, 0, 0 ,0, 0, 0, 0, 0, 0, 0, 0, 0]]) + [[-0.50183952, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], + [-1.26638196, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], + [0.43017941, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], + [0.65008914, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]]) gradients = (grad[i, :] for i in range(4)) ss = ActiveSubspaces(dim=2, method='exact', n_boot=150) ss.fit(gradients=gradients) - true_evects = np.array( - [[1., 0.], [0., 1.], [0., 0.], [0., 0.], [0., 0.], [0., 0.], - [0., 0.], [0., 0.], [0., 0.], [0., 0.], [0., 0.], [0., 0.], - [0., 0.], [0., 0.], [0., 0.]]) + true_evects = np.array([[1., 0.], [0., 1.], [0., 0.], [0., + 0.], [0., 0.], + [0., 0.], [0., 0.], [0., 0.], [0., + 0.], [0., 0.], + [0., 0.], [0., 0.], [0., 0.], [0., 0.], + [0., 0.]]) np.testing.assert_array_almost_equal(true_evects, ss.evects) - + def test_activity_scores_01(self): np.random.seed(42) gradients = np.random.uniform(-1, 1, 180).reshape(15, 3, 4) ss = ActiveSubspaces(dim=1, method='exact', n_boot=150) ss.fit(gradients=gradients) - true_scores = np.array([0.599489, 0.055179, 0.37102 , 0.300374]) + true_scores = np.array([0.599489, 0.055179, 0.37102, 0.300374]) np.testing.assert_array_almost_equal(true_scores, ss.activity_scores) def test_activity_scores_02(self): @@ -211,9 +214,11 @@ def test_activity_scores_04(self): gradients = (inputs[i, :] for i in range(4)) ss = ActiveSubspaces(dim=2, method='exact', n_boot=150) ss.fit(gradients=gradients) - true_scores = np.array([0.061844, 0.914259, 0.868018, 0.830417, 1.081719, 0.652845, - 0.486826, 0.974479, 0.165313, 0.111259, 1.135542, 0.547697, - 1.099458, 0.860475, 0.377612]) + true_scores = np.array([ + 0.061844, 0.914259, 0.868018, 0.830417, 1.081719, 0.652845, + 0.486826, 0.974479, 0.165313, 0.111259, 1.135542, 0.547697, + 1.099458, 0.860475, 0.377612 + ]) np.testing.assert_array_almost_equal(true_scores, ss.activity_scores) def test_transform_01(self): @@ -522,10 +527,10 @@ def test_plot_eigenvalues_04(self): def test_plot_eigenvalues_05(self): np.random.seed(42) grad = np.array( - [[-0.50183952, 0, 0, 0, 0, 0 ,0, 0, 0, 0, 0, 0, 0, 0, 0], - [-1.26638196, 0, 0, 0, 0, 0 ,0, 0, 0, 0, 0, 0, 0, 0, 0], - [ 0.43017941, 0, 0, 0, 0, 0 ,0, 0, 0, 0, 0, 0, 0, 0, 0], - [ 0.65008914, 0, 0, 0, 0, 0 ,0, 0, 0, 0, 0, 0, 0, 0, 0]]) + [[-0.50183952, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], + [-1.26638196, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], + [0.43017941, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], + [0.65008914, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]]) gradients = (grad[i, :] for i in range(4)) ss = ActiveSubspaces(dim=2, method='exact', n_boot=200) ss.fit(gradients=gradients) @@ -568,10 +573,10 @@ def test_plot_eigenvectors_04(self): def test_plot_eigenvectors_05(self): np.random.seed(42) grad = np.array( - [[-0.50183952, 0, 0, 0, 0, 0 ,0, 0, 0, 0, 0, 0, 0, 0, 0], - [-1.26638196, 0, 0, 0, 0, 0 ,0, 0, 0, 0, 0, 0, 0, 0, 0], - [ 0.43017941, 0, 0, 0, 0, 0 ,0, 0, 0, 0, 0, 0, 0, 0, 0], - [ 0.65008914, 0, 0, 0, 0, 0 ,0, 0, 0, 0, 0, 0, 0, 0, 0]]) + [[-0.50183952, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], + [-1.26638196, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], + [0.43017941, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], + [0.65008914, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]]) gradients = (grad[i, :] for i in range(4)) ss = ActiveSubspaces(dim=2, method='exact', n_boot=200) ss.fit(gradients=gradients) @@ -616,12 +621,12 @@ def test_plot_sufficient_summary_04(self): def test_frequent_directions_01(self): np.random.seed(42) grad = np.array( - [[-0.50183952, 0, 0, 0, 0, 0 ,0, 0, 0, 0, 0, 0, 0, 0, 0], - [-1.26638196, 0, 0, 0, 0, 0 ,0, 0, 0, 0, 0, 0, 0, 0, 0], - [ 0.43017941, 0, 0, 0, 0, 0 ,0, 0, 0, 0, 0, 0, 0, 0, 0], - [ 0.65008914, 0, 0, 0, 0, 0 ,0, 0, 0, 0, 0, 0, 0, 0, 0]]) + [[-0.50183952, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], + [-1.26638196, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], + [0.43017941, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], + [0.65008914, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]]) gradients = (grad[i, :] for i in range(4)) ss = ActiveSubspaces(dim=2, method='exact', n_boot=150) evals, v = ss._frequent_directions(gradients=gradients) self.assertEqual(v.shape, (15, 2)) - self.assertEqual(evals.shape, (2,)) + self.assertEqual(evals.shape, (2, )) diff --git a/tests/test_feature_map.py b/tests/test_feature_map.py index 7aacac5..1a13fdd 100644 --- a/tests/test_feature_map.py +++ b/tests/test_feature_map.py @@ -183,7 +183,7 @@ def test_brute(self): subspace=ss) best = fm.tune_pr_matrix(func=average_rrmse, bounds=[slice(-2, 1, 0.2) for i in range(1)], - fn_args={'csv':csv}, + fn_args={'csv': csv}, maxiter=10, save_file=False)[1] true = np.array([[0., 0.], [0., 0.], [0., 0.]]) @@ -208,7 +208,7 @@ def test_dual_annealing(self): subspace=ss) best = fm.tune_pr_matrix(func=average_rrmse, bounds=[slice(-2, 1, 0.2) for i in range(1)], - fn_args={'csv':csv}, + fn_args={'csv': csv}, method='dual_annealing', maxiter=5, save_file=False)[1] @@ -235,7 +235,7 @@ def test_bso(self): subspace=ss) best = fm.tune_pr_matrix(func=average_rrmse, bounds=[slice(-2, 1, 0.2) for i in range(1)], - fn_args={'csv':csv}, + fn_args={'csv': csv}, method='bso', maxiter=10, save_file=False)[1] diff --git a/tests/test_kas.py b/tests/test_kas.py index b58fe32..302c628 100644 --- a/tests/test_kas.py +++ b/tests/test_kas.py @@ -124,9 +124,7 @@ def test_compute_03(self): n_features=4, method='exact', n_boot=49) - ss.fit(inputs=inputs, - gradients=gradients, - weights=weights) + ss.fit(inputs=inputs, gradients=gradients, weights=weights) true_evects = np.array( [[0.74714817, 0.6155644, 0.23414206, 0.08959675], [0.35380297, -0.10917583, -0.91115623, 0.18082704], @@ -236,10 +234,7 @@ def test_compute_11(self): n_features=4, method='local', n_boot=49) - ss.fit(inputs=inputs, - outputs=outputs, - weights=weights, - metric=metric) + ss.fit(inputs=inputs, outputs=outputs, weights=weights, metric=metric) true_evals = np.array( [1.58774145e+05, 2.37399662e+02, 8.73268317e+01, 2.99624379e+01]) np.testing.assert_allclose(true_evals, ss.evals) diff --git a/tests/test_subspaces.py b/tests/test_subspaces.py index ec2f222..81853a4 100644 --- a/tests/test_subspaces.py +++ b/tests/test_subspaces.py @@ -125,3 +125,21 @@ def test_plot_sufficient_summary(self): outputs, figsize=(7, 7), title='Sufficient_summary_plots') + + def test_partition_spectral_gap(self): + np.random.seed(42) + matrix = np.array([[1, 1, 1], [2, -4.5, 2], [1, 1.1, 1]]) + weights = np.ones((3, 1)) + ss = Subspaces(dim=1) + ss.evals, ss.evects = ss._build_decompose_cov_matrix(matrix, weights) + ss.partition_spectral_gap() + self.assertEqual(ss.dim, 1) + + def test_partition_residual_energy(self): + np.random.seed(42) + matrix = np.array([[1, 1, 1], [2, -4.5, 2], [1, 1.1, 1]]) + weights = np.ones((3, 1)) + ss = Subspaces(dim=1) + ss.evals, ss.evects = ss._build_decompose_cov_matrix(matrix, weights) + ss.partition_residual_energy() + self.assertEqual(ss.dim, 1) From 9c9d88d17e92b2ed38166133ae4dc7d407f518b1 Mon Sep 17 00:00:00 2001 From: fromor Date: Tue, 28 Jun 2022 11:14:03 +0200 Subject: [PATCH 03/26] keep API consistent for spectralgap and residual energy options --- athena/active.py | 4 +- athena/subspaces.py | 88 +++++++++++++++++++++++++++-------------- tests/test_active.py | 4 +- tests/test_kas.py | 2 +- tests/test_subspaces.py | 12 +++--- 5 files changed, 69 insertions(+), 41 deletions(-) diff --git a/athena/active.py b/athena/active.py index 88ba47e..642976a 100644 --- a/athena/active.py +++ b/athena/active.py @@ -40,7 +40,7 @@ def activity_scores(self): Return the activity scores as defined in Constantine and Diaz, Global sensitivity metrics from active subspaces, arxiv.org/abs/1510.04361 Equation (21). - + :return: array with the activity score of each parameter. :rtype: numpy.ndarray :raises: TypeError @@ -366,7 +366,7 @@ def _rotate_x(self, reduced_inputs, inactive_inputs): def _frequent_directions(self, gradients): """ - Function that performs the frequent directions algorithm for + Function that performs the frequent directions algorithm for matrix sketching. For more details about the method, see "Frequent directions: Simple and deterministic matrix sketching." Ghashami, Mina, et al. diff --git a/athena/subspaces.py b/athena/subspaces.py index f8c3216..cce84f1 100644 --- a/athena/subspaces.py +++ b/athena/subspaces.py @@ -24,6 +24,11 @@ class Subspaces(): """Active Subspaces base class + :param dim: The dimension of the active subspace. If 0, the method computes + the spectral gap and uses it for truncation; if positive interger, + the method uses the `dim` argument for the truncation; if float between 0 + and 1, a truncation with retained energy.The `dim` + parameter is available using all the available methods. :param str method: method to compute the AS. Possible choices are 'exact' when the gradients are provided, or 'local' to use local linear models. This approach is related to the sufficient dimension reduction @@ -77,13 +82,7 @@ def _build_decompose_cov_matrix(self, X = np.squeeze(gradients * np.sqrt(weights).reshape(-1, 1)) n_samples, n_pars = X.shape - # computational complexity of svd and random svd - svd_complexity = n_samples * n_pars * self.dim - rsvd_complexity = n_samples * n_pars * np.log( - self.dim) + (n_samples + n_pars) * self.dim**2 - - if svd_complexity > rsvd_complexity and (n_samples > 10000 - or n_pars > 10000): + if self._check_rsvd(n_samples, n_pars, self.dim): singular, evects = randomized_svd( M=X, n_components=self.dim, @@ -123,14 +122,9 @@ def _compute_bootstrap_ranges(self, gradients, weights, metric=None): n_pars = gradients.shape[-1] n_samples = gradients.shape[0] - svd_complexity = n_samples * n_pars * self.dim - rsvd_complexity = n_samples * n_pars * np.log( - self.dim) + (n_samples + n_pars) * self.dim**2 - # randomized_svd is not implemented for vectorial as yet if len(gradients.shape) == 2: - if svd_complexity > rsvd_complexity and (n_samples > 10000 - or n_pars > 10000): + if self._check_rsvd(n_samples, n_pars, self.dim): range_dim = self.dim else: range_dim = min(n_pars, n_samples) @@ -231,10 +225,13 @@ def _partition(self): :raises: TypeError, ValueError """ + + self.dim = self._set_dim() + if not isinstance(self.dim, int): raise TypeError('dim should be an integer.') - if self.dim < 1 or self.dim > self.evects.shape[1]: + if self.dim < 0 or self.dim > self.evects.shape[1]: raise ValueError( 'dim must be positive and less than the dimension of the ' ' eigenvectors: dim = {}.'.format(self.dim)) @@ -251,27 +248,49 @@ def _partition(self): 'the eigenvectors cannot have dimension less than dim = {}.'. format(self.dim)) - def partition_spectral_gap(self): + def _set_dim(self): """ - Partition the eigenvectors to define the active and inactive subspaces, - based on the highest spectral gap of the ordered eigenvalues. + Set the active subspace dimension based on the spectral gap if self.dim + is 0, residual energy if 0=1. - :raises: TypeError + :raises: ValueError + """ + + # spectral gap + if isinstance(self.dim, int) and self.dim == 0: + dim = self._set_dim_spectral_gap() + # residual energy + elif self.dim > 0 and self.dim < 1: + dim = self._set_dim_residual_energy() + # manual set of AS dimension + elif isinstance(self.dim, int) and self.dim >=1: + dim = self.dim + else: + raise ValueError( + "The parameter `dim`={} has not a valid value.".format(self.dim) + ) + + return dim + + def _set_dim_spectral_gap(self): + """ + Set the active subspace dimension based on the highest spectral gap of the ordered eigenvalues. + + :raises: ValueError """ if self.evals is None: - raise TypeError( + raise ValueError( 'The method fit has to be called first in order to compute the eigenvalues.' ) spectral_gap_index = np.argmax(self.evals[:-1] - self.evals[1:]) - self.dim = int(spectral_gap_index + 1) - self._partition() + return int(spectral_gap_index + 1) - def partition_residual_energy(self, tol=0.99): + def _set_dim_residual_energy(self, tol=0.99): """ - Partition the eigenvectors to define the active and inactive subspaces, - based on the residual energy normalized to 1 i.e. the normalized sum of the eigenvalues - corresponding to the active components must be lower than tol. + Set the active subspace dimension based on the residual energy + normalized to 1: the normalized sum of the eigenvalues + corresponding to the active components will be lower than tol. :param float tol: threshold for the residual energy. @@ -282,11 +301,20 @@ def partition_residual_energy(self, tol=0.99): 'The method fit has to be called first in order to compute the eigenvalues.' ) - residual_energies = np.array( - list(itertools.accumulate(self.evals, operator.add))) / np.sum( - self.evals) - self.dim = residual_energies[residual_energies <= tol].shape[0] - self._partition() + cumulative_energy = np.cumsum(self.evals**2 / (self.evals**2).sum()) + return np.searchsorted(cumulative_energy, tol) + + def _check_rsvd(self, n_samples, n_pars, dim): + """ Check if random svd is to be applied.""" + if isinstance(dim, int) and dim > 0: + # computational complexity of svd and random svd + svd_complexity = n_samples * n_pars * dim + rsvd_complexity = n_samples * n_pars * np.log( + dim) + (n_samples + n_pars) * dim**2 + + return svd_complexity > rsvd_complexity and (n_samples > 10000 or n_pars > 10000) + else: + return False def plot_eigenvalues(self, n_evals=None, diff --git a/tests/test_active.py b/tests/test_active.py index d807de2..cdce55e 100644 --- a/tests/test_active.py +++ b/tests/test_active.py @@ -391,13 +391,13 @@ def test_partition_03(self): matrix = np.random.uniform(-1, 1, 9).reshape(3, 3) ss = ActiveSubspaces(dim=2.0) ss.evects = matrix - with self.assertRaises(TypeError): + with self.assertRaises(ValueError): ss._partition() def test_partition_04(self): np.random.seed(42) matrix = np.random.uniform(-1, 1, 9).reshape(3, 3) - ss = ActiveSubspaces(dim=0) + ss = ActiveSubspaces(dim=0.) ss.evects = matrix with self.assertRaises(ValueError): ss._partition() diff --git a/tests/test_kas.py b/tests/test_kas.py index 302c628..f4715ef 100644 --- a/tests/test_kas.py +++ b/tests/test_kas.py @@ -309,7 +309,7 @@ def test_partition_03(self): matrix = np.random.uniform(-1, 1, 9).reshape(3, 3) ss = KernelActiveSubspaces(dim=2.0) ss.evects = matrix - with self.assertRaises(TypeError): + with self.assertRaises(ValueError): ss._partition() def test_partition_04(self): diff --git a/tests/test_subspaces.py b/tests/test_subspaces.py index 81853a4..757af49 100644 --- a/tests/test_subspaces.py +++ b/tests/test_subspaces.py @@ -68,7 +68,7 @@ def test_partition_03(self): matrix = np.random.uniform(-1, 1, 9).reshape(3, 3) ss = Subspaces(dim=2.0) ss.evects = matrix - with self.assertRaises(TypeError): + with self.assertRaises(ValueError): ss._partition() def test_partition_04(self): @@ -130,10 +130,10 @@ def test_partition_spectral_gap(self): np.random.seed(42) matrix = np.array([[1, 1, 1], [2, -4.5, 2], [1, 1.1, 1]]) weights = np.ones((3, 1)) - ss = Subspaces(dim=1) + ss = Subspaces(dim=0) ss.evals, ss.evects = ss._build_decompose_cov_matrix(matrix, weights) - ss.partition_spectral_gap() - self.assertEqual(ss.dim, 1) + dim = ss._set_dim_spectral_gap() + self.assertEqual(dim, 1) def test_partition_residual_energy(self): np.random.seed(42) @@ -141,5 +141,5 @@ def test_partition_residual_energy(self): weights = np.ones((3, 1)) ss = Subspaces(dim=1) ss.evals, ss.evects = ss._build_decompose_cov_matrix(matrix, weights) - ss.partition_residual_energy() - self.assertEqual(ss.dim, 1) + dim = ss._set_dim_residual_energy() + self.assertEqual(dim, 1) From c78af47eb3e9d956e0512a6c8f1ce2c851d96dc6 Mon Sep 17 00:00:00 2001 From: mtezzele Date: Thu, 10 Nov 2022 18:25:55 -0600 Subject: [PATCH 04/26] fix typo in axis name --- athena/nll.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/athena/nll.py b/athena/nll.py index a77627c..7cda4a7 100644 --- a/athena/nll.py +++ b/athena/nll.py @@ -215,7 +215,7 @@ def plot_sufficient_summary(self, reduced_inputs = self.forward(inputs)[:, 0] plt.plot(reduced_inputs.detach().numpy(), outputs, 'bo') plt.xlabel('Reduced input') - plt.xlabel('Output') + plt.ylabel('Output') else: raise ValueError( 'Sufficient summary plots cannot be made in more than 1 ' \ From 892fdec1cb4177521ea4ee19d1d67c384b566f2a Mon Sep 17 00:00:00 2001 From: mtezzele Date: Tue, 13 Dec 2022 18:02:37 -0600 Subject: [PATCH 05/26] clean up packages needed in the Coverage ci --- .github/workflows/ci.yml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml index e1530d9..9ccfdbc 100644 --- a/.github/workflows/ci.yml +++ b/.github/workflows/ci.yml @@ -21,7 +21,7 @@ jobs: - name: Install Python dependencies run: | python3 -m pip install --upgrade pip - python3 -m pip install numpy scipy matplotlib pytest pytest-cov future torch nose gpy gpyopt sphinx sklearn + python3 -m pip install numpy scipy matplotlib pytest pytest-cov future torch gpy gpyopt scikit-learn scikit-learn-extra - name: Test with pytest env: From 08903a0a37c47ccbb2d871eadc2fe8f9356a0d6e Mon Sep 17 00:00:00 2001 From: mtezzele Date: Tue, 13 Dec 2022 18:30:25 -0600 Subject: [PATCH 06/26] add test for W2 = None --- tests/test_active.py | 10 ++++++++++ 1 file changed, 10 insertions(+) diff --git a/tests/test_active.py b/tests/test_active.py index cdce55e..c8a66db 100644 --- a/tests/test_active.py +++ b/tests/test_active.py @@ -275,6 +275,16 @@ def test_transform_05(self): [-0.49850367, -0.37146678]]) np.testing.assert_array_almost_equal(true_inactive, inactive) + def test_transform_06(self): + np.random.seed(42) + inputs = np.random.uniform(-1, 1, 60).reshape(15, 4) + outputs = np.random.uniform(0, 5, 15) + ss = ActiveSubspaces(dim=2, method='local', n_boot=250) + ss.fit(inputs=inputs, outputs=outputs) + ss.W2 = None + inactive = ss.transform(np.random.uniform(-1, 1, 8).reshape(2, 4))[1] + self.assertIsNone(inactive) + def test_inverse_transform_01(self): np.random.seed(42) inputs = np.random.uniform(-1, 1, 60).reshape(15, 4) From d1e69c995f6a62a45f3880ba9410cae733b9f945 Mon Sep 17 00:00:00 2001 From: Nicola Demo Date: Wed, 14 Dec 2022 11:20:28 +0100 Subject: [PATCH 07/26] Create monthly-tag.yml --- .github/workflows/monthly-tag.yml | 46 +++++++++++++++++++++++++++++++ 1 file changed, 46 insertions(+) create mode 100644 .github/workflows/monthly-tag.yml diff --git a/.github/workflows/monthly-tag.yml b/.github/workflows/monthly-tag.yml new file mode 100644 index 0000000..134388b --- /dev/null +++ b/.github/workflows/monthly-tag.yml @@ -0,0 +1,46 @@ +name: Monthly Automated Tag + +on: + schedule: + - cron: '20 2 1 * *' + +jobs: + + test: + runs-on: ${{ matrix.os }} + strategy: + matrix: + os: [windows-latest, macos-latest, ubuntu-latest] + python-version: [3.7, 3.8] + steps: + - uses: actions/checkout@v2 + - name: Set up Python + uses: actions/setup-python@v2 + with: + python-version: ${{ matrix.python-version }} + - name: Install Python dependencies + run: | + python3 -m pip install --upgrade pip + python3 -m pip install .[test] + - name: Test with pytest + run: | + python3 -m pytest + + monthly_tag: + runs-on: ubuntu-latest + needs: test + steps: + - uses: actions/checkout@v2 + with: + token: ${{ secrets.NDEMO_PAT_TOKEN }} + + - name: Create and push the tag + run: | + python utils/mathlab_versioning.py set --only-date "post$(date +%y%m)" + VERS=$(python utils/mathlab_versioning.py get) + git config --global user.name 'Monthly Tag bot' + git config --global user.email 'mtbot@noreply.github.com' + git add athena/meta.py + git commit -m "monthly version $VERS" + git tag -a "v$VERS" -m "Monthly version $VERS" + git push origin "v$VERS" From 9a3f906dc62c4d4ba6370405427fbd77400a42a5 Mon Sep 17 00:00:00 2001 From: Nicola Demo Date: Wed, 14 Dec 2022 11:22:27 +0100 Subject: [PATCH 08/26] Create mathlab_versioning.py --- utils/mathlab_versioning.py | 100 ++++++++++++++++++++++++++++++++++++ 1 file changed, 100 insertions(+) create mode 100644 utils/mathlab_versioning.py diff --git a/utils/mathlab_versioning.py b/utils/mathlab_versioning.py new file mode 100644 index 0000000..577a8ec --- /dev/null +++ b/utils/mathlab_versioning.py @@ -0,0 +1,100 @@ +import re +import os +import argparse + + +module = 'pina' +meta_file = os.path.join(module, 'meta.py') +version_line = r'__version__.*=.*"(.+?)"' + + +class Version: + def __init__(self, major, minor, patch, date_patch=None): + self.major = major + self.minor = minor + self.patch = patch + self.date_patch = date_patch + + def __str__(self): + + if self.date_patch: + version_string = '{}.{}.{}.{}'.format( + self.major, + self.minor, + self.patch, + self.date_patch + ) + else: + version_string = '{}.{}.{}'.format( + self.major, + self.minor, + self.patch, + ) + return version_string + + +def get_version(): + with open(meta_file, 'r') as fp: + content = fp.read() + + try: + found = re.search(r'__version__.*=.*"(.+?)"', content).group(1) + except AttributeError: + pass + + version = re.split(r'[-\.]', found) + v = Version(*version) + return v + + +def set_version(version): + with open(meta_file, 'r') as fp: + content = fp.read() + + line_string = '__version__ = "{}"'.format(version) + text_after = re.sub('__version__.*=.*"(.+?)"', line_string, content) + + with open(meta_file, 'w') as fp: + fp.write(text_after) + + +if __name__ == '__main__': + parser = argparse.ArgumentParser(description='Manipulate Version') + + subparsers = parser.add_subparsers(dest='command') + + get_ = subparsers.add_parser('get', + help='Get information about current version') + set_ = subparsers.add_parser('set', help='Set version') + flags = set_.add_mutually_exclusive_group(required=False) + flags.add_argument('--only-major', action='store_true') + flags.add_argument('--only-minor', action='store_true') + flags.add_argument('--only-patch', action='store_true') + flags.add_argument('--only-date', action='store_true') + set_.add_argument('version', nargs='+', action="store") + + args = parser.parse_args() + + if args.command == 'get': + print(get_version()) + elif args.command == 'set': + if args.only_major: + current_version = get_version() + current_version.major = args.version[0] + set_version(current_version) + elif args.only_minor: + current_version = get_version() + current_version.minor = args.version[0] + set_version(current_version) + elif args.only_patch: + current_version = get_version() + current_version.patch = args.version[0] + set_version(current_version) + elif args.only_date: + current_version = get_version() + current_version.date_patch = args.version[0] + set_version(current_version) + elif len(args.version) in [3, 4]: + set_version(Version(*args.version)) + else: + raise RuntimeError From f215c0e6abb0dfa36eec3da3fac73f5fe0fee3df Mon Sep 17 00:00:00 2001 From: Nicola Demo Date: Wed, 14 Dec 2022 11:22:43 +0100 Subject: [PATCH 09/26] Update mathlab_versioning.py --- utils/mathlab_versioning.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/utils/mathlab_versioning.py b/utils/mathlab_versioning.py index 577a8ec..b9da118 100644 --- a/utils/mathlab_versioning.py +++ b/utils/mathlab_versioning.py @@ -3,7 +3,7 @@ import argparse -module = 'pina' +module = 'athena' meta_file = os.path.join(module, 'meta.py') version_line = r'__version__.*=.*"(.+?)"' From db05c6e6a43ffaa5aa28c5a77a1a80466f0c8c38 Mon Sep 17 00:00:00 2001 From: mtezzele Date: Wed, 21 Dec 2022 14:43:49 -0600 Subject: [PATCH 10/26] refactor code, nested ifs --- athena/active.py | 13 ++++--------- 1 file changed, 4 insertions(+), 9 deletions(-) diff --git a/athena/active.py b/athena/active.py index 3a8267c..63ce995 100644 --- a/athena/active.py +++ b/athena/active.py @@ -96,9 +96,8 @@ def fit(self, >>> ss3 = ActiveSubspaces(dim=1, method='local', n_boot=150) >>> ss3.fit(inputs=inputs, outputs=outputs) """ - if self.method == 'exact': - if gradients is None: - raise TypeError('gradients argument is None.') + if self.method == 'exact' and gradients is None: + raise TypeError('gradients argument is None.') # estimate active subspace with local linear models. if self.method == 'local': @@ -146,11 +145,7 @@ def transform(self, inputs): active = np.dot(inputs, self.W1) # allow evaluation of active variables only - if self.W2 is None: - inactive = None - else: - inactive = np.dot(inputs, self.W2) - + inactive = None if self.W2 is None else np.dot(inputs, self.W2) return active, inactive def inverse_transform(self, reduced_inputs, n_points=1): @@ -378,7 +373,7 @@ def _frequent_directions(self, gradients): the reduced matrix. :rtype: numpy.ndarray, numpy.ndarray """ - s = np.array([next(gradients) for i in range(self.dim)]).T + s = np.array([next(gradients) for _ in range(self.dim)]).T for grad in gradients: evects, sigma = np.linalg.svd(s, full_matrices=False)[:2] s = np.dot( From 818bd066b7a220351ca6560dc68dab27d195645e Mon Sep 17 00:00:00 2001 From: mtezzele Date: Wed, 21 Dec 2022 14:50:11 -0600 Subject: [PATCH 11/26] refactor code, f-strings --- athena/subspaces.py | 91 ++++++++++++++++++++++----------------------- 1 file changed, 44 insertions(+), 47 deletions(-) diff --git a/athena/subspaces.py b/athena/subspaces.py index cce84f1..9b49452 100644 --- a/athena/subspaces.py +++ b/athena/subspaces.py @@ -66,35 +66,36 @@ def _build_decompose_cov_matrix(self, :return: the sorted eigenvalues, and the corresponding eigenvectors. :rtype: numpy.ndarray, numpy.ndarray """ - if self.method == 'exact' or self.method == 'local': - if metric is not None: - cov_matrix = np.array( - np.sum([ - weights[i, 0] * - np.dot(gradients[i, :, :].T, - np.dot(metric, gradients[i, :, :])) - for i in range(gradients.shape[0]) - ], - axis=0)) - evals, evects = sort_eigpairs(cov_matrix) - return np.squeeze(evals), evects - - X = np.squeeze(gradients * np.sqrt(weights).reshape(-1, 1)) - n_samples, n_pars = X.shape - - if self._check_rsvd(n_samples, n_pars, self.dim): - singular, evects = randomized_svd( - M=X, - n_components=self.dim, - n_oversamples=10, - n_iter='auto', - power_iteration_normalizer='auto', - transpose='auto')[1:] - else: - singular, evects = np.linalg.svd(X, full_matrices=False)[1:] + if self.method not in ['exact', 'local']: + return + if metric is not None: + cov_matrix = np.array( + np.sum([ + weights[i, 0] * + np.dot(gradients[i, :, :].T, + np.dot(metric, gradients[i, :, :])) + for i in range(gradients.shape[0]) + ], + axis=0)) + evals, evects = sort_eigpairs(cov_matrix) + return np.squeeze(evals), evects + + X = np.squeeze(gradients * np.sqrt(weights).reshape(-1, 1)) + n_samples, n_pars = X.shape + + if self._check_rsvd(n_samples, n_pars, self.dim): + singular, evects = randomized_svd( + M=X, + n_components=self.dim, + n_oversamples=10, + n_iter='auto', + power_iteration_normalizer='auto', + transpose='auto')[1:] + else: + singular, evects = np.linalg.svd(X, full_matrices=False)[1:] - evals = singular**2 - return evals, evects.T + evals = singular**2 + return evals, evects.T def _compute_bootstrap_ranges(self, gradients, weights, metric=None): """Compute bootstrap ranges for eigenvalues and subspaces. @@ -183,8 +184,8 @@ def fit(self, *args, **kwargs): to be implemented in subclasses. """ raise NotImplementedError( - 'Subclass must implement abstract method {}.fit'.format( - self.__class__.__name__)) + f'Subclass must implement abstract method {self.__class__.__name__}.fit' + ) def transform(self, inputs): """ @@ -200,9 +201,9 @@ def transform(self, inputs): variables. :rtype: numpy.ndarray, numpy.ndarray """ - raise NotImplementedError('Subclass must implement abstract method ' - '{}.transform'.format( - self.__class__.__name__)) + raise NotImplementedError( + f'Subclass must implement abstract method {self.__class__.__name__}.transform' + ) def inverse_transform(self, reduced_inputs, n_points): """ @@ -216,8 +217,8 @@ def inverse_transform(self, reduced_inputs, n_points): space that are returned that map to the given active variables. """ raise NotImplementedError( - 'Subclass must implement abstract method {}.inverse_transform'. - format(self.__class__.__name__)) + f'Subclass must implement abstract method {self.__class__.__name__}.inverse_transform' + ) def _partition(self): """ @@ -233,8 +234,8 @@ def _partition(self): if self.dim < 0 or self.dim > self.evects.shape[1]: raise ValueError( - 'dim must be positive and less than the dimension of the ' - ' eigenvectors: dim = {}.'.format(self.dim)) + f'dim must be positive and less than the dimension of the eigenvectors: dim = {self.dim}.' + ) # allow evaluation of active eigenvectors only if self.evects.shape[1] < self.evects.shape[0]: @@ -245,8 +246,8 @@ def _partition(self): self.W2 = self.evects[:, self.dim:] else: raise ValueError( - 'the eigenvectors cannot have dimension less than dim = {}.'. - format(self.dim)) + f'the eigenvectors cannot have dimension less than dim = {self.dim}.' + ) def _set_dim(self): """ @@ -259,16 +260,12 @@ def _set_dim(self): # spectral gap if isinstance(self.dim, int) and self.dim == 0: dim = self._set_dim_spectral_gap() - # residual energy elif self.dim > 0 and self.dim < 1: dim = self._set_dim_residual_energy() - # manual set of AS dimension elif isinstance(self.dim, int) and self.dim >=1: dim = self.dim else: - raise ValueError( - "The parameter `dim`={} has not a valid value.".format(self.dim) - ) + raise ValueError(f"The parameter `dim`={self.dim} has not a valid value.") return dim @@ -397,7 +394,7 @@ def plot_eigenvectors(self, filename=None, figsize=None, labels=None, - title=''): + title=''): # sourcery skip: extract-method """ Plot the eigenvectors. @@ -440,7 +437,7 @@ def plot_eigenvectors(self, if labels: ax.set_xticklabels(labels) - ax.set_ylabel('Active eigenvector {}'.format(i + 1)) + ax.set_ylabel(f'Active eigenvector {i + 1}') ax.grid(linestyle='dotted') ax.axis([0, n_pars + 1, -1 - 0.1, 1 + 0.1]) @@ -459,7 +456,7 @@ def plot_sufficient_summary(self, outputs, filename=None, figsize=(10, 8), - title=''): + title=''): # sourcery skip: extract-method """ Plot the sufficient summary. From a39928180a10df63cf9bcbb23881a3ca8bb0ac15 Mon Sep 17 00:00:00 2001 From: mtezzele Date: Wed, 21 Dec 2022 14:53:03 -0600 Subject: [PATCH 12/26] refactor code, f-strings, returns --- athena/utils.py | 27 +++++++++++---------------- 1 file changed, 11 insertions(+), 16 deletions(-) diff --git a/athena/utils.py b/athena/utils.py index 1ac98ae..eb8113c 100644 --- a/athena/utils.py +++ b/athena/utils.py @@ -28,8 +28,7 @@ def fit_transform(self, inputs): between -1 and 1. :rtype: numpy.ndarray """ - inputs_norm = 2.0 * (inputs - self.lb) / (self.ub - self.lb) - 1.0 - return inputs_norm + return 2.0 * (inputs - self.lb) / (self.ub - self.lb) - 1.0 def inverse_transform(self, inputs): """Return corresponding points shifted and scaled to @@ -42,8 +41,7 @@ def inverse_transform(self, inputs): between `self.lb` and `self.ub`. :rtype: numpy.ndarray """ - inputs_unnorm = (self.ub - self.lb) * (inputs + 1.0) / 2.0 + self.lb - return inputs_unnorm + return (self.ub - self.lb) * (inputs + 1.0) / 2.0 + self.lb def initialize_weights(matrix): @@ -81,11 +79,11 @@ def linear_program_ineq(c, A, b): b = b.reshape(-1, ) # make unbounded bounds - bounds = [(None, None) for i in range(c.shape[0])] + bounds = [(None, None) for _ in range(c.shape[0])] res = linprog(c=c, A_ub=-A, b_ub=-b, bounds=bounds) if res.success: return res.x.reshape(-1, 1) - raise RuntimeError('Scipy did not solve the LP. {}'.format(res.message)) + raise RuntimeError(f'Scipy did not solve the LP. {res.message}') def local_linear_gradients(inputs, outputs, weights=None, n_neighbors=None): @@ -119,14 +117,12 @@ def local_linear_gradients(inputs, outputs, weights=None, n_neighbors=None): if n_neighbors is None: n_neighbors = int(min(np.floor(1.7 * n_pars), n_samples)) elif not isinstance(n_neighbors, int): - raise TypeError( - 'n_neighbors ({}) must be an integer.'.format(n_neighbors)) + raise TypeError(f'n_neighbors ({n_neighbors}) must be an integer.') if n_neighbors <= n_pars or n_neighbors > n_samples: raise ValueError( - 'n_neighbors must be between the number of parameters ' - 'and the number of samples. Unsatisfied: {} < {} < {}.'.format( - n_pars, n_neighbors, n_samples)) + f'n_neighbors must be between the number of parameters and the number of samples. Unsatisfied: {n_pars} < {n_neighbors} < {n_samples}.' + ) if weights is None: weights = initialize_weights(inputs) @@ -208,7 +204,7 @@ class CrossValidation(): """ def __init__(self, inputs, outputs, gradients, subspace, folds=5, **kwargs): - if any([v is None for v in [inputs, outputs, gradients, subspace]]): + if any(v is None for v in [inputs, outputs, gradients, subspace]): raise ValueError( 'Any among inputs, outputs, gradients, subspace is None.') @@ -278,8 +274,7 @@ def predict(self, inputs): :rtype: numpy.ndarray """ x_test = self.ss.transform(inputs)[0] - y = self.gp.predict(np.atleast_2d(x_test))[0] - return y + return self.gp.predict(np.atleast_2d(x_test))[0] def scorer(self, inputs, outputs): """ @@ -343,7 +338,7 @@ def average_rrmse(hyperparams, best, csv, verbose=False, resample=5): a specified number of resamples of the projection matrix. :rtype: numpy.float64 """ - if isinstance(csv, CrossValidation) is False: + if not isinstance(csv, CrossValidation): raise ValueError( "The argument csv must be of type athena.utils.CrossValidation") @@ -371,7 +366,7 @@ def average_rrmse(hyperparams, best, csv, verbose=False, resample=5): # save the best parameters if verbose is True: - print("params {} mean {}, std {}".format(hyperparams, mean, std)) + print(f"params {hyperparams} mean {mean}, std {std}") score_records.append(mean) # skip resampling from the same hyperparam if the error is not below From 87bfe0f763917e0a31e5960cb420496bec14a899 Mon Sep 17 00:00:00 2001 From: mtezzele Date: Wed, 21 Dec 2022 14:54:10 -0600 Subject: [PATCH 13/26] refactor code, cls --- athena/projection_factory.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/athena/projection_factory.py b/athena/projection_factory.py index d21555f..c98bd16 100644 --- a/athena/projection_factory.py +++ b/athena/projection_factory.py @@ -147,11 +147,11 @@ def uniform(input_dim, n_features, params): 'uniform': uniform.__func__ } - def __new__(self, fname): + def __new__(cls, fname): # to make the str callable we have to use a dictionary with all the # implemented projection matrices - if fname in self.projections: - return self.__projections[fname] + if fname in cls.projections: + return cls.__projections[fname] raise NameError( """The name of the projection matrix is not correct or not implemented. Check the documentation for all the available From c7400a8deeef4eff8998f1cd2637258da7add208 Mon Sep 17 00:00:00 2001 From: mtezzele Date: Wed, 21 Dec 2022 14:56:36 -0600 Subject: [PATCH 14/26] reference --- athena/feature_map.py | 26 +++++++++++++------------- athena/kas.py | 7 ++++--- 2 files changed, 17 insertions(+), 16 deletions(-) diff --git a/athena/feature_map.py b/athena/feature_map.py index e31dee4..d36d30d 100644 --- a/athena/feature_map.py +++ b/athena/feature_map.py @@ -4,9 +4,10 @@ :References: - Francesco Romor, Marco Tezzele, Andrea Lario, Gianluigi Rozza. - Kernel-based Active Subspaces with application to CFD problems using - Discontinuous Galerkin Method. 2020. - arxiv: https://arxiv.org/abs/2008.12083 + Kernel-based active subspaces with application to computational fluid + dynamics parametric problems using discontinuous Galerkin method. + International Journal for Numerical Methods in Engineering, + 123(23):6000–6027, 2022. doi:10.1002/nme.7099 """ from functools import partial @@ -176,11 +177,7 @@ def tune_pr_matrix(self, best = [0.8, np.zeros((self.n_features, self.input_dim))] if method is None: - if len(self.params) < 4: - method = 'brute' - else: - method = 'dual_annealing' - + method = 'brute' if len(self.params) < 4 else 'dual_annealing' if method == 'brute': self.params = brute(func=func, ranges=bounds, @@ -200,11 +197,14 @@ def tune_pr_matrix(self, maxiter=maxiter, no_local_search=False).x elif method == 'bso': - bounds = [{ - 'name': 'var_' + str(i), - 'type': 'continuous', - 'domain': [bound.start, bound.stop] - } for i, bound in enumerate(bounds)] + bounds = [ + { + 'name': f'var_{str(i)}', + 'type': 'continuous', + 'domain': [bound.start, bound.stop], + } + for i, bound in enumerate(bounds) + ] func_obj = partial(func, best=best, **fn_args) bopt = GPyOpt.methods.BayesianOptimization(func_obj, domain=bounds, diff --git a/athena/kas.py b/athena/kas.py index ae10b39..80f4cd0 100644 --- a/athena/kas.py +++ b/athena/kas.py @@ -4,9 +4,10 @@ :References: - Francesco Romor, Marco Tezzele, Andrea Lario, Gianluigi Rozza. - Kernel-based Active Subspaces with application to CFD problems using - Discontinuous Galerkin Method. 2020. - arxiv: https://arxiv.org/abs/2008.12083 + Kernel-based active subspaces with application to computational fluid + dynamics parametric problems using discontinuous Galerkin method. + International Journal for Numerical Methods in Engineering, + 123(23):6000–6027, 2022. doi:10.1002/nme.7099 """ import numpy as np From ca84b6ab7efd92a3ad9cbf87b8ecc2ae3bb55a5f Mon Sep 17 00:00:00 2001 From: mtezzele Date: Wed, 21 Dec 2022 14:59:42 -0600 Subject: [PATCH 15/26] refactor code, f-strings --- athena/nll.py | 66 ++++++++++++++++++++++++--------------------------- 1 file changed, 31 insertions(+), 35 deletions(-) diff --git a/athena/nll.py b/athena/nll.py index 7cda4a7..9ba3903 100644 --- a/athena/nll.py +++ b/athena/nll.py @@ -130,7 +130,7 @@ def train(self, gradients) if i % 10 == 0 or i == self.epochs - 1: - print('epoch = {}, loss = {}'.format(i, loss)) + print(f'epoch = {i}, loss = {loss}') if interactive: ax1.cla() ax1.set_title('Sufficient summary plot') @@ -149,8 +149,8 @@ def train(self, self.backward.zero_grad() for j in range(self.backward.n_layers): - name_y = 'fc{}_y'.format(j + 1) - name_z = 'fc{}_z'.format(j + 1) + name_y = f'fc{j + 1}_y' + name_z = f'fc{j + 1}_z' getattr(self.backward, name_y).weight = torch.nn.Parameter( getattr(self.forward, name_y).weight) getattr(self.backward, name_z).weight = torch.nn.Parameter( @@ -329,10 +329,8 @@ def __init__(self, n_params, n_layers, dh, active_dim): self.omega = slice(active_dim) for i in range(self.n_layers): - setattr(self, 'fc{}_y'.format(i + 1), - nn.Linear(self.n_params, 2 * self.n_params)) - setattr(self, 'fc{}_z'.format(i + 1), - nn.Linear(self.n_params, 2 * self.n_params)) + setattr(self, f'fc{i + 1}_y', nn.Linear(self.n_params, 2 * self.n_params)) + setattr(self, f'fc{i + 1}_z', nn.Linear(self.n_params, 2 * self.n_params)) def forward(self, inputs): """ @@ -349,12 +347,12 @@ def forward(self, inputs): vars()['var0_z'] = torch.clone(bb[1]) for i in range(self.n_layers): - name_y = 'fc{}_y'.format(i + 1) - name_z = 'fc{}_z'.format(i + 1) - var_y0 = 'var{}_y'.format(i) - var_z0 = 'var{}_z'.format(i) - var_y1 = 'var{}_y'.format(i + 1) - var_z1 = 'var{}_z'.format(i + 1) + name_y = f'fc{i + 1}_y' + name_z = f'fc{i + 1}_z' + var_y0 = f'var{i}_y' + var_z0 = f'var{i}_z' + var_y1 = f'var{i + 1}_y' + var_z1 = f'var{i + 1}_z' sig_y = torch.unsqueeze( torch.tanh(getattr(self, name_y)(vars()[var_z0])), 2) @@ -386,16 +384,16 @@ def customized_loss(self, inputs, mapped_inputs, gradients): """ # Define the weights and bias of the inverse network for i in range(self.n_layers): - name_y = 'fc{}_y'.format(i + 1) - name_z = 'fc{}_z'.format(i + 1) + name_y = f'fc{i + 1}_y' + name_z = f'fc{i + 1}_z' - inv_name_y_weight = 'inv_fc{}_y_weight'.format(i + 1) - inv_name_z_weight = 'inv_fc{}_z_weight'.format(i + 1) + inv_name_y_weight = f'inv_fc{i + 1}_y_weight' + inv_name_z_weight = f'inv_fc{i + 1}_z_weight' vars()[inv_name_y_weight] = getattr(self, name_y).weight vars()[inv_name_z_weight] = getattr(self, name_z).weight - inv_name_y_bias = 'inv_fc{}_y_bias'.format(i + 1) - inv_name_z_bias = 'inv_fc{}_z_bias'.format(i + 1) + inv_name_y_bias = f'inv_fc{i + 1}_y_bias' + inv_name_z_bias = f'inv_fc{i + 1}_z_bias' vars()[inv_name_y_bias] = getattr(self, name_y).bias vars()[inv_name_z_bias] = getattr(self, name_z).bias @@ -407,20 +405,20 @@ def customized_loss(self, inputs, mapped_inputs, gradients): output_dy[:, j] += 0.001 bb = torch.split(output_dy, self.n_params, dim=1) - var_y0 = 'var{}_y'.format(self.n_layers - 1) - var_z0 = 'var{}_z'.format(self.n_layers - 1) + var_y0 = f'var{self.n_layers - 1}_y' + var_z0 = f'var{self.n_layers - 1}_z' vars()[var_y0] = bb[0] vars()[var_z0] = bb[1] for i in range(self.n_layers - 1, -1, -1): - inv_name_y_weight = 'inv_fc{}_y_weight'.format(i + 1) - inv_name_z_weight = 'inv_fc{}_z_weight'.format(i + 1) - inv_name_y_bias = 'inv_fc{}_y_bias'.format(i + 1) - inv_name_z_bias = 'inv_fc{}_z_bias'.format(i + 1) - var_y0 = 'var{}_y'.format(i) - var_z0 = 'var{}_z'.format(i) - var_y1 = 'var{}_y'.format(i - 1) - var_z1 = 'var{}_z'.format(i - 1) + inv_name_y_weight = f'inv_fc{i + 1}_y_weight' + inv_name_z_weight = f'inv_fc{i + 1}_z_weight' + inv_name_y_bias = f'inv_fc{i + 1}_y_bias' + inv_name_z_bias = f'inv_fc{i + 1}_z_bias' + var_y0 = f'var{i}_y' + var_z0 = f'var{i}_z' + var_y1 = f'var{i - 1}_y' + var_z1 = f'var{i - 1}_z' sig_z = torch.tanh( torch.add( @@ -487,10 +485,8 @@ def __init__(self, n_params, n_layers, dh): self.dh = dh for i in range(self.n_layers): - setattr(self, 'fc{}_y'.format(i + 1), - nn.Linear(self.n_params, 2 * self.n_params)) - setattr(self, 'fc{}_z'.format(i + 1), - nn.Linear(self.n_params, 2 * self.n_params)) + setattr(self, f'fc{i + 1}_y', nn.Linear(self.n_params, 2 * self.n_params)) + setattr(self, f'fc{i + 1}_z', nn.Linear(self.n_params, 2 * self.n_params)) def forward(self, mapped_inputs): """ @@ -505,8 +501,8 @@ def forward(self, mapped_inputs): y, z = torch.split(mapped_inputs, self.n_params, dim=1) for i in range(self.n_layers - 1, -1, -1): - name_y = 'fc{}_y'.format(i + 1) - name_z = 'fc{}_z'.format(i + 1) + name_y = f'fc{i + 1}_y' + name_z = f'fc{i + 1}_z' sig_z = torch.unsqueeze(torch.tanh(getattr(self, name_z)(y)), 2) K_z = torch.transpose(getattr(self, name_z).weight, 0, 1) From c622782ea164d789d641ca0bbd26a9bd08e613bf Mon Sep 17 00:00:00 2001 From: mtezzele Date: Wed, 21 Dec 2022 15:18:55 -0600 Subject: [PATCH 16/26] refactor code, swap if-else, improve readability, formatting --- athena/feature_map.py | 17 ++- athena/local.py | 247 +++++++++++++++------------------ athena/local_classification.py | 72 +++++----- athena/nll.py | 20 ++- athena/subspaces.py | 31 +++-- 5 files changed, 188 insertions(+), 199 deletions(-) diff --git a/athena/feature_map.py b/athena/feature_map.py index d36d30d..e380c45 100644 --- a/athena/feature_map.py +++ b/athena/feature_map.py @@ -110,7 +110,7 @@ def tune_pr_matrix(self, method=None, maxiter=50, save_file=False, - fn_args={}): + fn_args=None): """ Tune the parameters of the spectral distribution. Three methods are available: log-grid-search (brute), annealing (dual_annealing) and @@ -174,6 +174,8 @@ def tune_pr_matrix(self, numpy.ndarray of zeros. :rtype: list """ + if fn_args is None: + fn_args = {} best = [0.8, np.zeros((self.n_features, self.input_dim))] if method is None: @@ -197,14 +199,11 @@ def tune_pr_matrix(self, maxiter=maxiter, no_local_search=False).x elif method == 'bso': - bounds = [ - { - 'name': f'var_{str(i)}', - 'type': 'continuous', - 'domain': [bound.start, bound.stop], - } - for i, bound in enumerate(bounds) - ] + bounds = [{ + 'name': f'var_{str(i)}', + 'type': 'continuous', + 'domain': [bound.start, bound.stop], + } for i, bound in enumerate(bounds)] func_obj = partial(func, best=best, **fn_args) bopt = GPyOpt.methods.BayesianOptimization(func_obj, domain=bounds, diff --git a/athena/local.py b/athena/local.py index 4d2003d..1c804a3 100644 --- a/athena/local.py +++ b/athena/local.py @@ -45,14 +45,14 @@ def __init__(self): self.outputs_val = None self.outputs_dim = None - self.method = None # gradients evaluation method + self.method = None # gradients evaluation method self.as_dim = None - self.clustering = None # implemented KMeans or KMedoids + self.clustering = None # implemented KMeans or KMedoids # cluster labels start from 0 self.labels = None # labels list of each input - self.unique_labels = None # set of labels + self.unique_labels = None # set of labels self.full_as = None self.full_gpr = None @@ -62,15 +62,19 @@ def __init__(self): self.max_clusters = None self.random_state = None - def _init_local_as(self, inputs=None, outputs=None, gradients=None, config=None): + def _init_local_as(self, + inputs=None, + outputs=None, + gradients=None, + config=None): assert inputs is not None and outputs is not None, 'inputs or outputs argument is None' - if isinstance(inputs, np.ndarray) and isinstance(outputs, np.ndarray): - self.inputs = inputs - self.outputs = outputs.squeeze() - else: + if not isinstance(inputs, np.ndarray) or not isinstance( + outputs, np.ndarray): raise TypeError('Inputs and outputs type is not np.ndarray.') + self.inputs = inputs + self.outputs = outputs.squeeze() if gradients is None or isinstance(gradients, np.ndarray): self.gradients = gradients else: @@ -95,10 +99,8 @@ def _init_local_as(self, inputs=None, outputs=None, gradients=None, config=None) self.outputs_val = config['outputs_val'] if config.get( 'outputs_val', None) is not None else self.outputs - if len(self.outputs.shape) == 2: - self.outputs_dim = self.outputs.shape[1] - else: - self.outputs_dim = 1 + self.outputs_dim = self.outputs.shape[1] if len( + self.outputs.shape) == 2 else 1 def fit(self, inputs=None, outputs=None, gradients=None, config=None): """ @@ -112,8 +114,8 @@ def fit(self, inputs=None, outputs=None, gradients=None, config=None): def _fit_clustering(self): """Init unique_labels, labels, method predict""" raise NotImplementedError( - 'Subclass must implement abstract method {}._fit_clustering'. - format(self.__class__.__name__)) + f'Subclass must implement abstract method {self.__class__.__name__}._fit_clustering' + ) def _fit_global(self, plot=False): """Compute global Active Subspace""" @@ -185,15 +187,11 @@ def compute_scores(self, inputs_test, outputs_test): mae_full = mean_absolute_error(full_pred, outputs_test) mae_local = mean_absolute_error(local_pred, outputs_test) - scores = np.array([r2_full, r2_local, mae_full, mae_local]) - return scores + return np.array([r2_full, r2_local, mae_full, mae_local]) @staticmethod def cluster_metric(x, y): - if x[-1] == y[-1]: - return np.linalg.norm(x - y) - else: - return 1000 + return np.linalg.norm(x - y) if x[-1] == y[-1] else 1000 # TODO plot clusters with dimensions higher than 2 with TSNE ? def plot_clusters(self, save=False, title='2d_clusters', plot=True): @@ -203,7 +201,8 @@ def plot_clusters(self, save=False, title='2d_clusters', plot=True): plt.scatter(self.inputs[:, 0], self.inputs[:, 1], c=self.labels) plt.grid(linestyle='dotted') plt.tight_layout() - if save: plt.savefig('{}.pdf'.format(title)) + if save: + plt.savefig(f'{title}.pdf') if plot: plt.show() @@ -221,7 +220,7 @@ def _fit_clustering(self): self.labels = self.clustering.labels_ self.unique_labels = list(set(self.labels)) self.centers = self.clustering.cluster_centers_ - _log.debug("Number of clusters k-means: {}".format(self.max_clusters)) + _log.debug(f"Number of clusters k-means: {self.max_clusters}") silhouette_ = silhouette_score(self.inputs, self.labels) _log.debug( @@ -243,7 +242,7 @@ def _fit_clustering(self): self.labels = self.clustering.labels_ self.unique_labels = list(set(self.labels)) self.centers = self.clustering.cluster_centers_ - _log.debug("Number of clusters k-means: {}".format(self.max_clusters)) + _log.debug(f"Number of clusters k-means: {self.max_clusters}") silhouette_ = silhouette_score(self.inputs, self.labels) _log.debug( @@ -277,12 +276,12 @@ def __init__(self): self.score_tolerance = None self.dim_cut = None - self.normalization = None # uniform, gaussian, root + self.normalization = None # uniform, gaussian, root self.metric = None # as, CI, euclidean - self.refinement_criterion = None #global, mean, force - self.as_dim_criterion = None # spectral_gap, residual + self.refinement_criterion = None #global, mean, force + self.as_dim_criterion = None # spectral_gap, residual - self.minimum_score = None # it increases the cluster dimension until the score threshold is reached + self.minimum_score = None # it increases the cluster dimension until the score threshold is reached self.max_dim_refine_further = None self.verbose = False @@ -311,23 +310,19 @@ def _init_hierarchical(self, config): self.min_local = config['min_local'] self.score_tolerance = config['score_tolerance'] self.dim_cut = config['dim_cut'] - self.normalization = config[ - 'normalization'] + self.normalization = config['normalization'] self.metric = config['metric'] - self.refinement_criterion = config[ - 'refinement_criterion'] - self.as_dim_criterion = config[ - 'as_dim_criterion'] + self.refinement_criterion = config['refinement_criterion'] + self.as_dim_criterion = config['as_dim_criterion'] except KeyError as k: raise KeyError( - "Missing mandatory keyword {} for class TopDownHierarchical". - format(k.args[0])) from k + f"Missing mandatory keyword {k.args[0]} for class TopDownHierarchical" + ) from k - self.minimum_score = config.get( - 'minimum_score', None - ) - self.max_dim_refine_further = config.get('max_dim_refine_further', self.inputs.shape[1]) + self.minimum_score = config.get('minimum_score', None) + self.max_dim_refine_further = config.get('max_dim_refine_further', + self.inputs.shape[1]) self.verbose = config.get('verbose', True) def fit(self, inputs=None, outputs=None, gradients=None, config=None): @@ -358,16 +353,15 @@ def _fit_clustering(self, print_states=False, plot=False): # continue until children_fifo is empty while (children_fifo != []): _log.debug( - "Before pop total queue length is: {},\n total clusters: {}" - .format(str(len(children_fifo)), str(self.total_clusters))) + f"Before pop total queue length is: {len(children_fifo)},\n total clusters: {str(self.total_clusters)}" + ) node = children_fifo.pop(0) state, children = node.refine_cluster() _log.debug( - "After pop {} and len children list added is {} total clusters: {}" - .format(str(state), str(len(children)), - str(self.total_clusters))) + f"After pop {str(state)} and len children list added is {len(children)} total clusters: {str(self.total_clusters)}" + ) if 3 in state: break @@ -375,13 +369,13 @@ def _fit_clustering(self, print_states=False, plot=False): continue elif 4 in state: # refine further the present node self.total_clusters += len(children) - 1 - for child in children: - if child.state == 4: - children_fifo.append(child) + children_fifo.extend(child for child in children + if child.state == 4) elif 2 in state: # tolerance reached, continue self.total_clusters += len(children) - 1 - elif any(x in state for x in - [1, 5]): # abort further clustering conditions + elif any( + x in state + for x in [1, 5]): # abort further clustering conditions continue res = self.compute_scores(self.inputs_test, self.outputs_test) @@ -391,10 +385,8 @@ def _fit_clustering(self, print_states=False, plot=False): print("Hierarchical top-down clustering completed with states") for sta in state: print(sta, " : ", self.state_d[str(sta)]) - print("n_leaves: {:d}".format(n_leaves)) - - print("Test score: {}".format(res)) + print(f"Test score: {res}") if print_states: self._print_state() @@ -475,7 +467,7 @@ def _print_state_debug(self): children_fifo = [self.root] while (children_fifo != []): node = children_fifo.pop(0) - _log.debug(str(node.state) + " : " + self.state_d[str(node.state)]) + _log.debug(f"{str(node.state)} : {self.state_d[str(node.state)]}") children_fifo.extend(node.children) def _print_leaves_score(self): @@ -557,7 +549,7 @@ def __call__(self, node): leaves_info.r_dims[i] * col_ones, leaves_info.leaf_labels[i] * col_ones))) arr = np.vstack(to_be_saved) - _log.debug("Saved data with shape: {}".format(arr.shape)) + _log.debug(f"Saved data with shape: {arr.shape}") np.save("clusters_r2_asdim_labels.npy", to_be_saved) fig = plt.figure(figsize=(12, 10)) @@ -583,8 +575,8 @@ def __call__(self, node): plt.legend() plt.title( - "Hierarchical top-down clustering, total clusters= {} ".format( - self.total_clusters)) + f"Hierarchical top-down clustering, total clusters= {self.total_clusters} " + ) plt.tight_layout() if save: plt.savefig("clusters_top_down.png") @@ -667,9 +659,7 @@ def predict(self, test_inputs): return predictions, labels def _cluster_min_dist(self, test, clusters): - group_list = [] - for node in clusters: - group_list.append(self.inputs[node.ind]) + group_list = [self.inputs[node.ind] for node in clusters] return self.set_min_dist(test, group_list) @staticmethod @@ -734,8 +724,9 @@ def __init__(self, parent, node_indexes, val_indexes, tree_obj): self.state = None self.type = "root" if self.parent is None else self.hierarchical.normalization - _log.debug("Validation set length: {}".format( - self.hierarchical.inputs_val[self.t_ind].shape)) + _log.debug( + f"Validation set length: {self.hierarchical.inputs_val[self.t_ind].shape}" + ) ####### Estabilish invariants r_dim and score self.normalizer = self.NormalizeDivisive(self.type, self.ind, @@ -767,7 +758,7 @@ def __init__(self, parent, node_indexes, val_indexes, tree_obj): # change as dimension ass.dim = self.r_dim ass._partition() - _log.debug("reduced dim is: " + str(self.r_dim)) + _log.debug(f"reduced dim is: {str(self.r_dim)}") # 2. the score is above the score_tolerance gpr = GPy.models.GPRegression(ass.transform(inp_normalized)[0], @@ -787,14 +778,14 @@ def __init__(self, parent, node_indexes, val_indexes, tree_obj): local_pred, multioutput='raw_values'), 0, 1)) - _log.debug("r2 score is: " + str(self.score)) + _log.debug(f"r2 score is: {str(self.score)}") def refine_further(self, minimum_score): - _log.debug("r2 score inside refine_further is: " + str(self.score)) + _log.debug(f"r2 score inside refine_further is: {str(self.score)}") # create self.ss if it does not exist, through property decorator - while (self.score < minimum_score and - self.ss.dim + 1 <= self.hierarchical.max_dim_refine_further): + while (self.score < minimum_score + and self.ss.dim + 1 <= self.hierarchical.max_dim_refine_further): # change as dimension self.ss.dim += 1 self.r_dim += 1 @@ -811,9 +802,9 @@ def refine_further(self, minimum_score): _log.debug("reduced dim inside refine_further loop is: " + str(self.r_dim)) - t_inp_normalized = self.normalizer( - self.t_ind, self.hierarchical.inputs_val, - self.hierarchical.outputs_val)[0] + t_inp_normalized = self.normalizer(self.t_ind, + self.hierarchical.inputs_val, + self.hierarchical.outputs_val)[0] # re-optimize gpr self._gpr = None @@ -824,7 +815,7 @@ def refine_further(self, minimum_score): local_pred, multioutput='uniform_average') - _log.debug("r2 score inside while cycle is: " + str(self.score)) + _log.debug(f"r2 score inside while cycle is: {str(self.score)}") class NormalizeDivisive(): """Inner class for normalization of inputs, gradients w.r.t. local @@ -856,11 +847,7 @@ def __init__(self, norm_type, ind, inputs): "A proper flag was not passed to normalize class: " + self.type) - def __call__(self, - ind=None, - inputs=None, - outputs=None, - gradients=None): + def __call__(self, ind=None, inputs=None, outputs=None, gradients=None): if ind is not None: inp = inputs[ind] out = outputs[ind] @@ -898,21 +885,18 @@ def __call__(self, @property # lazy property def gpr(self): - # TODO check type and behaviour e.g. return only if state is 2 or 4 - if self._gpr is None: - inp_normalized, out, _ = self.normalizer(self.ind, - self.hierarchical.inputs, - self.hierarchical.outputs) - - gpr = GPy.models.GPRegression(self.ss.transform(inp_normalized)[0], - out.reshape(inp_normalized.shape[0], - -1), - normalizer=False) - gpr.optimize_restarts(20, verbose=False) - self._gpr = gpr - return gpr - else: + if self._gpr is not None: return self._gpr + inp_normalized, out, _ = self.normalizer(self.ind, + self.hierarchical.inputs, + self.hierarchical.outputs) + + gpr = GPy.models.GPRegression(self.ss.transform(inp_normalized)[0], + out.reshape(inp_normalized.shape[0], -1), + normalizer=False) + gpr.optimize_restarts(20, verbose=False) + self._gpr = gpr + return gpr @gpr.setter def gpr(self, value): @@ -925,22 +909,21 @@ def gpr(self): @property # lazy property def ss(self): # TODO check type and behaviour e.g. return only if state is 2 or 4 - if self._ss is None: - inp_normalized, out, grad_normalized = self.normalizer( - self.ind, self.hierarchical.inputs, self.hierarchical.outputs, - self.hierarchical.gradients) - - if grad_normalized is not None: - ss = ActiveSubspaces(dim=self.r_dim) - ss.fit(gradients=grad_normalized) - else: - ss = ActiveSubspaces(dim=self.r_dim, method='local') - ss.fit(inputs=inp_normalized, outputs=out) + if self._ss is not None: + return self._ss + inp_normalized, out, grad_normalized = self.normalizer( + self.ind, self.hierarchical.inputs, self.hierarchical.outputs, + self.hierarchical.gradients) - self._ss = ss - return ss + if grad_normalized is not None: + ss = ActiveSubspaces(dim=self.r_dim) + ss.fit(gradients=grad_normalized) else: - return self._ss + ss = ActiveSubspaces(dim=self.r_dim, method='local') + ss.fit(inputs=inp_normalized, outputs=out) + + self._ss = ss + return ss @ss.setter def ss(self, value): @@ -993,8 +976,7 @@ 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))) @@ -1004,10 +986,9 @@ def refine_cluster(self): # cluster with metric of choice if self.hierarchical.metric == 'as': _log.debug("AS metric") - cluster_ = KMedoids( - metric=self.as_metric, - n_clusters=n_clusters, - random_state=self.hierarchical.random_state) + cluster_ = KMedoids(metric=self.as_metric, + n_clusters=n_clusters, + random_state=self.hierarchical.random_state) local_inp = self.normalizer(self.ind, self.hierarchical.inputs, self.hierarchical.outputs)[0] @@ -1021,9 +1002,8 @@ def refine_cluster(self): # C1 norm with outputs included elif self.hierarchical.metric == 'C1': _log.debug("C1 metric") - cluster_ = KMedoids( - n_clusters=n_clusters, - random_state=self.hierarchical.random_state) + cluster_ = KMedoids(n_clusters=n_clusters, + random_state=self.hierarchical.random_state) inp_normalized, out, grad_normalized = self.normalizer( self.ind, self.hierarchical.inputs, @@ -1033,9 +1013,8 @@ def refine_cluster(self): elif self.hierarchical.metric == 'euclidean': _log.debug("Euclidean metric") - cluster_ = KMedoids( - n_clusters=n_clusters, - random_state=self.hierarchical.random_state) + cluster_ = KMedoids(n_clusters=n_clusters, + random_state=self.hierarchical.random_state) local_inp = self.normalizer(self.ind, self.hierarchical.inputs, self.hierarchical.outputs)[0] @@ -1155,14 +1134,12 @@ def global_children_priority(self, children, res_max, res, children_states, multioutput='uniform_average') if res_max < full_children_r2: - _log.debug("Best score bested: {} < {}".format( - res_max, full_children_r2)) + _log.debug(f"Best score bested: {res_max} < {full_children_r2}") res_max = full_children_r2 self.children = children self.cluster = cluster else: - _log.debug("Best score not bested: {} < {}".format( - res_max, full_children_r2)) + _log.debug(f"Best score not bested: {res_max} < {full_children_r2}") if res_max >= self.hierarchical.score_tolerance: return_value = "break" @@ -1182,32 +1159,32 @@ def mean_children_priority(self, children, res_max, res, children_states, # stop refinement: tolerance achieved for the first time w.r.t. # n_clusters loop - if 2 == max_state: + if max_state == 2: _log.debug("Leaf reached") self.children = children self.cluster = cluster return_value = "break" - elif 4 == max_state: + elif max_state == 4: _log.debug("Leaf not reached, children states" + str(children_states)) # if the overall score is bested change children DEBUG changed to # 1.01 instead of * 1 if res_max < res_: - _log.debug("Bested: " + str(res_max) + " < " + str(res_)) + _log.debug(f"Bested: {str(res_max)} < {str(res_)}") res_max = res_ self.children = children self.cluster = cluster else: - _log.debug("Not bested: " + str(res_max) + " " + str(res_)) + _log.debug(f"Not bested: {str(res_max)} {str(res_)}") return_value = "continue" return return_value, res_max, state - def force_refinement_priority(self, children, res_max, res, - children_states, cluster, state): + def force_refinement_priority(self, children, res_max, res, children_states, + cluster, state): """The refinement is chosen if the local r2 scores are above the threshold or if res_max is below the average of the local r2 scores. It can NOT happen that refinement is not performed because the parent score @@ -1218,13 +1195,13 @@ def force_refinement_priority(self, children, res_max, res, # stop refinement: tolerance achieved for the first time w.r.t. # n_clusters loop - if 2 == max_state: + if max_state == 2: _log.debug("Leaf reached") self.children = children self.cluster = cluster return_value = "break" - elif 4 == max_state: + elif max_state == 4: _log.debug("Leaf not reached, children states" + str(children_states)) @@ -1239,12 +1216,12 @@ def force_refinement_priority(self, children, res_max, res, # if the overall score is bested change children if res_max < res_: - _log.debug("Bested: " + str(res_max) + " < " + str(res_)) + _log.debug(f"Bested: {str(res_max)} < {str(res_)}") res_max = res_ self.children = children self.cluster = cluster else: - _log.debug("Not Bested: " + str(res_max) + " > " + str(res_)) + _log.debug(f"Not Bested: {str(res_max)} > {str(res_)}") return_value = "continue" @@ -1267,17 +1244,17 @@ def as_dim_gap(self, max_dim): assert max_dim < self.evals.shape[ 0], "max_red_dim cannot be equal or grater than the input dimension" max_gap = self.evals[0] - self.evals[1] - _log.debug("Spectral gaps: " + str(max_gap)) + _log.debug(f"Spectral gaps: {str(max_gap)}") r_dim = 1 i = 1 while (i < max_dim): gap = self.evals[i] - self.evals[i + 1] - _log.debug("Spectral gaps: " + str(gap)) + _log.debug(f"Spectral gaps: {str(gap)}") if gap > max_gap: max_gap = gap r_dim = i + 1 i += 1 - _log.debug("Returned: " + str(r_dim)) + _log.debug(f"Returned: {str(r_dim)}") return r_dim def as_dim(self, max_dim, threshold=0.95): @@ -1288,13 +1265,13 @@ def as_dim(self, max_dim, threshold=0.95): format(self.cum_evals[r_dim - 1], threshold, r_dim)) if self.cum_evals[r_dim - 1] >= threshold: - _log.debug("break {}".format(r_dim)) + _log.debug(f"break {r_dim}") break r_dim += 1 else: - _log.debug("else {}".format(r_dim)) + _log.debug(f"else {r_dim}") r_dim -= 1 - _log.debug("Returned: " + str(r_dim)) + _log.debug(f"Returned: {r_dim}") return r_dim def as_metric(self, X, Y): diff --git a/athena/local_classification.py b/athena/local_classification.py index df30613..e505aac 100644 --- a/athena/local_classification.py +++ b/athena/local_classification.py @@ -73,10 +73,12 @@ def plot(self, id1=0, id2=1, save=False): for i in range(self.n_components): # dim_mask = self.features[:, -1] == i + 1 dim_mask = self.labels == i - ax1.scatter(self.X[dim_mask, id1], - self.X[dim_mask, id2], - c=colors[i].reshape(1, -1), - label="id " + str(i)) + ax1.scatter( + self.X[dim_mask, id1], + self.X[dim_mask, id2], + c=colors[i].reshape(1, -1), + label=f"id {str(i)}", + ) plt.title("clusters") plt.grid() plt.legend() @@ -85,10 +87,12 @@ def plot(self, id1=0, id2=1, save=False): colors1 = cm.rainbow(np.linspace(0, 1, self.X.shape[1])) for i in range(self.X.shape[1]): dim_mask = self.features[:, -1] == i + 1 - ax1.scatter(self.X[dim_mask, id1], - self.X[dim_mask, id2], - c=colors1[i].reshape(1, -1), - label="dim " + str(i + 1)) + ax1.scatter( + self.X[dim_mask, id1], + self.X[dim_mask, id2], + c=colors1[i].reshape(1, -1), + label=f"dim {str(i + 1)}", + ) plt.title("as dimensions") plt.grid() plt.legend() @@ -102,11 +106,11 @@ def plot(self, id1=0, id2=1, save=False): def make_graph(self): """Use scipy sparse format COO to create adjacency list""" distance_matrix = self.custom_distance(self.features) - _log.debug("distance matrix {}".format(distance_matrix)) + _log.debug(f"distance matrix {distance_matrix}") adjacency_list = np.argsort(distance_matrix, axis=1)[:, :self.n_neighbours] - _log.debug("adjacency list {}".format(adjacency_list)) + _log.debug(f"adjacency list {adjacency_list}") distance_restricted = np.take_along_axis(distance_matrix, adjacency_list, @@ -117,19 +121,19 @@ def make_graph(self): neighbours_count = np.array([self.n_neighbours]) - np.count_nonzero( inf_mask, axis=1) - _log.debug("neighbours count {}".format(neighbours_count)) + _log.debug(f"neighbours count {neighbours_count}") assert (all(neighbours_count >= 1)) row = np.hstack((i * np.ones(neighbours_count[i], dtype=np.int8) for i in range(neighbours_count.shape[0]))) - _log.debug("row {}".format(row)) + _log.debug(f"row {row}") col = np.hstack((adjacency_list[i, :neighbours_count[i]] for i in range(neighbours_count.shape[0]))) - _log.debug("col {}".format(col)) + _log.debug(f"col {col}") data = np.ones(np.sum(neighbours_count), dtype=np.int8) - _log.debug("data {}".format(data)) + _log.debug(f"data {data}") assert (data.shape[0] == row.shape[0] == col.shape[0]) return coo_matrix((data, (row, col)), dtype=np.int8) @@ -218,7 +222,7 @@ def evaluate_minimum_distance(self, X, dX): for i in range(X.shape[0]): D = np.linalg.norm((X - X[i].reshape(1, -1)), axis=1) ind = np.argsort(D) - _log.debug("ind {}, neigh {}, dist {}".format(i, ind, D)) + _log.debug(f"ind {i}, neigh {ind}, dist {D}") res = [] @@ -235,13 +239,13 @@ def evaluate_minimum_distance(self, X, dX): if cumulative > self.threshold: break res.append(as_dim) - _log.debug(" {} evals {}, res {}, {}, \n{}".format( - mask, as_dim, res, evals[:4], - dX[ind[:self.n_neighbours]][mask, :])) + _log.debug( + f" {mask} evals {as_dim}, res {res}, {evals[:4]}, \n{dX[ind[:self.n_neighbours]][mask, :]}" + ) - _log.debug("{} as dimensiones {} approx {}".format( - i, int(as_dim_features[i] / self.n_neighbours), - as_dim_features[i] / self.n_neighbours)) + _log.debug( + f"{i} as dimensiones {int(as_dim_features[i] / self.n_neighbours)} approx {as_dim_features[i] / self.n_neighbours}" + ) if self.local_as_criterion == 'min': as_dim_features[i] = round(min(res)) @@ -249,10 +253,10 @@ def evaluate_minimum_distance(self, X, dX): as_dim_features[i] = round(sum(res) / self.n_neighbours) else: raise ValueError( - "The local_as_criterion must be 'min' or 'mean'. Passed value is {}" - .format(self.local_as_criterion)) + f"The local_as_criterion must be 'min' or 'mean'. Passed value is {self.local_as_criterion}" + ) - _log.debug("as dimensiones {}".format(as_dim_features)) + _log.debug(f"as dimensiones {as_dim_features}") return as_dim_features.reshape(-1, 1) @@ -265,7 +269,7 @@ def decision_boundaries(X, test_size=0.2, classifier=MLPClassifier(alpha=1, max_iter=1000)): - _log.debug("Shapes: {} {} ".format(X.shape, y.shape)) + _log.debug(f"Shapes: {X.shape} {y.shape} ") clf = classifier h = .02 # step size in the mesh @@ -277,15 +281,15 @@ def decision_boundaries(X, test_size=test_size, random_state=42) - _log.debug("Shapes train, test: {}, {}, {}, {}".format( - X_train.shape, X_test.shape, y_train.shape, y_test.shape)) + _log.debug( + f"Shapes train, test: {X_train.shape}, {X_test.shape}, {y_train.shape}, {y_test.shape}" + ) x_min, x_max = X[:, components[0]].min() - .1, X[:, components[0]].max() + .1 y_min, y_max = X[:, components[1]].min() - .1, X[:, components[1]].max() + .1 - xx, yy = np.meshgrid(np.arange(x_min, x_max, h), - np.arange(y_min, y_max, h)) + xx, yy = np.meshgrid(np.arange(x_min, x_max, h), np.arange(y_min, y_max, h)) # train lassifier clf.fit(X_train, y_train) @@ -293,7 +297,7 @@ def decision_boundaries(X, if true_inputs is not None and true_labels is not None: score_true = clf.score(true_inputs, true_labels) - _log.debug("Score {} {}".format(score, score_true)) + _log.debug(f"Score {score} {score_true}") if plot: # Plot the decision boundary. For that, we will assign a color to @@ -303,12 +307,12 @@ def decision_boundaries(X, Z = clf.decision_function( np.hstack((np.c_[xx.ravel(), yy.ravel()], np.zeros((xx.ravel().shape[0], X.shape[1] - 2))))) - _log.debug("Decision function {}".format(Z.shape)) + _log.debug(f"Decision function {Z.shape}") else: Z = clf.predict_proba( np.hstack((np.c_[xx.ravel(), yy.ravel()], np.zeros((xx.ravel().shape[0], X.shape[1] - 2))))) - _log.debug("Predictions probability {}".format(Z.shape[1])) + _log.debug(f"Predictions probability {Z.shape[1]}") n_labels = Z.shape[1] @@ -332,7 +336,7 @@ def decision_boundaries(X, ax[l].set_ylim(yy.min(), yy.max()) ax[l].set_xticks(()) ax[l].set_yticks(()) - ax[l].set_title("Local cluster {}".format(l)) + ax[l].set_title(f"Local cluster {l}") plt.tight_layout() plt.show() @@ -341,4 +345,4 @@ def decision_boundaries(X, if true_inputs is not None and true_labels is not None: return score, score_true else: - return score \ No newline at end of file + return score diff --git a/athena/nll.py b/athena/nll.py index 9ba3903..fb1ca7c 100644 --- a/athena/nll.py +++ b/athena/nll.py @@ -72,8 +72,8 @@ def train(self, outputs=None, interactive=False, target_loss=0.0001, - optim_args={}, - scheduler_args={}): + optim_args=None, + scheduler_args=None): """ Train the whole RevNet. @@ -93,6 +93,10 @@ def train(self, :raises: ValueError: in interactive mode outputs must be provided for the sufficient summary plot. """ + if optim_args is None: + optim_args = {} + if scheduler_args is None: + scheduler_args = {} if inputs.shape[1] % 2 != 0: raise ValueError('The parameter space\'s dimension must be even.') @@ -329,8 +333,10 @@ def __init__(self, n_params, n_layers, dh, active_dim): self.omega = slice(active_dim) for i in range(self.n_layers): - setattr(self, f'fc{i + 1}_y', nn.Linear(self.n_params, 2 * self.n_params)) - setattr(self, f'fc{i + 1}_z', nn.Linear(self.n_params, 2 * self.n_params)) + setattr(self, f'fc{i + 1}_y', + nn.Linear(self.n_params, 2 * self.n_params)) + setattr(self, f'fc{i + 1}_z', + nn.Linear(self.n_params, 2 * self.n_params)) def forward(self, inputs): """ @@ -485,8 +491,10 @@ def __init__(self, n_params, n_layers, dh): self.dh = dh for i in range(self.n_layers): - setattr(self, f'fc{i + 1}_y', nn.Linear(self.n_params, 2 * self.n_params)) - setattr(self, f'fc{i + 1}_z', nn.Linear(self.n_params, 2 * self.n_params)) + setattr(self, f'fc{i + 1}_y', + nn.Linear(self.n_params, 2 * self.n_params)) + setattr(self, f'fc{i + 1}_z', + nn.Linear(self.n_params, 2 * self.n_params)) def forward(self, mapped_inputs): """ diff --git a/athena/subspaces.py b/athena/subspaces.py index 9b49452..47e5b34 100644 --- a/athena/subspaces.py +++ b/athena/subspaces.py @@ -18,6 +18,7 @@ from .utils import sort_eigpairs import itertools import operator + plt.rcParams.update({'font.size': 16}) @@ -71,9 +72,8 @@ def _build_decompose_cov_matrix(self, if metric is not None: cov_matrix = np.array( np.sum([ - weights[i, 0] * - np.dot(gradients[i, :, :].T, - np.dot(metric, gradients[i, :, :])) + weights[i, 0] * np.dot(gradients[i, :, :].T, + np.dot(metric, gradients[i, :, :])) for i in range(gradients.shape[0]) ], axis=0)) @@ -84,13 +84,12 @@ def _build_decompose_cov_matrix(self, n_samples, n_pars = X.shape if self._check_rsvd(n_samples, n_pars, self.dim): - singular, evects = randomized_svd( - M=X, - n_components=self.dim, - n_oversamples=10, - n_iter='auto', - power_iteration_normalizer='auto', - transpose='auto')[1:] + singular, evects = randomized_svd(M=X, + n_components=self.dim, + n_oversamples=10, + n_iter='auto', + power_iteration_normalizer='auto', + transpose='auto')[1:] else: singular, evects = np.linalg.svd(X, full_matrices=False)[1:] @@ -262,10 +261,11 @@ def _set_dim(self): dim = self._set_dim_spectral_gap() elif self.dim > 0 and self.dim < 1: dim = self._set_dim_residual_energy() - elif isinstance(self.dim, int) and self.dim >=1: + elif isinstance(self.dim, int) and self.dim >= 1: dim = self.dim else: - raise ValueError(f"The parameter `dim`={self.dim} has not a valid value.") + raise ValueError( + f"The parameter `dim`={self.dim} has not a valid value.") return dim @@ -306,10 +306,11 @@ def _check_rsvd(self, n_samples, n_pars, dim): if isinstance(dim, int) and dim > 0: # computational complexity of svd and random svd svd_complexity = n_samples * n_pars * dim - rsvd_complexity = n_samples * n_pars * np.log( - dim) + (n_samples + n_pars) * dim**2 + rsvd_complexity = n_samples * n_pars * np.log(dim) + ( + n_samples + n_pars) * dim**2 - return svd_complexity > rsvd_complexity and (n_samples > 10000 or n_pars > 10000) + return svd_complexity > rsvd_complexity and (n_samples > 10000 + or n_pars > 10000) else: return False From fd27594e74f55d62248e43de334ed3e5a1930c51 Mon Sep 17 00:00:00 2001 From: mtezzele Date: Wed, 21 Dec 2022 18:21:40 -0600 Subject: [PATCH 17/26] fix eigenpair sort, replace gpy with scikit-learn --- athena/active.py | 4 +-- athena/subspaces.py | 9 +++--- athena/utils.py | 23 ++++++++------ tests/test_active.py | 28 ++++++++--------- tests/test_feature_map.py | 65 ++++++++++++++++++++------------------- tests/test_nll.py | 8 ++--- tests/test_utils.py | 14 ++++----- 7 files changed, 80 insertions(+), 71 deletions(-) diff --git a/athena/active.py b/athena/active.py index 63ce995..edbabe8 100644 --- a/athena/active.py +++ b/athena/active.py @@ -314,12 +314,12 @@ def _hit_and_run_inactive(self, reduced_input, n_points): # find an upper bound on the step min_ind = np.logical_and(g <= 0, - f < -np.sqrt(np.finfo(np.float).eps)) + 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.float).eps)) + f < -np.sqrt(np.finfo(np.float64).eps)) eps_min = np.amax(f[max_ind] / g[max_ind]) # randomly sample eps diff --git a/athena/subspaces.py b/athena/subspaces.py index 47e5b34..48a0467 100644 --- a/athena/subspaces.py +++ b/athena/subspaces.py @@ -75,9 +75,8 @@ def _build_decompose_cov_matrix(self, weights[i, 0] * np.dot(gradients[i, :, :].T, np.dot(metric, gradients[i, :, :])) for i in range(gradients.shape[0]) - ], - axis=0)) - evals, evects = sort_eigpairs(cov_matrix) + ], axis=0)) + evals, evects = sort_eigpairs(*np.linalg.eigh(cov_matrix)) return np.squeeze(evals), evects X = np.squeeze(gradients * np.sqrt(weights).reshape(-1, 1)) @@ -94,6 +93,8 @@ def _build_decompose_cov_matrix(self, singular, evects = np.linalg.svd(X, full_matrices=False)[1:] evals = singular**2 + + evals, evects = sort_eigpairs(evals, evects) return evals, evects.T def _compute_bootstrap_ranges(self, gradients, weights, metric=None): @@ -318,7 +319,7 @@ def plot_eigenvalues(self, n_evals=None, filename=None, figsize=(8, 8), - title=''): + title=''): # sourcery skip: class-extract-method """ Plot the eigenvalues. diff --git a/athena/utils.py b/athena/utils.py index eb8113c..0ea44e2 100644 --- a/athena/utils.py +++ b/athena/utils.py @@ -4,6 +4,7 @@ import matplotlib.pyplot as plt from scipy.optimize import linprog import GPy +from sklearn.gaussian_process import GaussianProcessRegressor class Normalizer(): @@ -155,10 +156,11 @@ def local_linear_gradients(inputs, outputs, weights=None, n_neighbors=None): return gradients, new_inputs -def sort_eigpairs(matrix): - """Compute eigenpairs and sort. +def sort_eigpairs(evals, evects): + """Sort eigenpairs. - :param numpy.ndarray matrix: matrix whose eigenpairs you want. + :param numpy.ndarray evals: eigenvalues. + :param numpy.ndarray evects: eigenvectors. :return: vector of sorted eigenvalues; orthogonal matrix of corresponding eigenvectors. :rtype: numpy.ndarray, numpy.ndarray @@ -169,7 +171,6 @@ def sort_eigpairs(matrix): the eigenvectors so that the first component of each eigenvector is positive. This normalization is very helpful for the bootstrapping. """ - evals, evects = np.linalg.eigh(matrix) evals = abs(evals) ind = np.argsort(evals) evals = evals[ind[::-1]] @@ -177,7 +178,7 @@ def sort_eigpairs(matrix): s = np.sign(evects[0, :]) s[s == 0] = 1 evects *= s - return evals.reshape(-1, 1), evects + return evals, evects class CrossValidation(): @@ -260,9 +261,12 @@ def fit(self, inputs, gradients, outputs): y = self.ss.transform(inputs)[0] # build response surface - kern = GPy.kern.RBF(input_dim=y.shape[1], ARD=True) - self.gp = GPy.models.GPRegression(y, np.atleast_2d(outputs), kern) - self.gp.optimize_restarts(15, verbose=False) + # kern = GPy.kern.RBF(input_dim=y.shape[1], ARD=True) + # self.gp = GPy.models.GPRegression(y, np.atleast_2d(outputs), kern) + # self.gp.optimize_restarts(15, verbose=False) + + self.gp = GaussianProcessRegressor(n_restarts_optimizer=15) + self.gp.fit(y, np.atleast_2d(outputs)) def predict(self, inputs): """ @@ -274,7 +278,8 @@ def predict(self, inputs): :rtype: numpy.ndarray """ x_test = self.ss.transform(inputs)[0] - return self.gp.predict(np.atleast_2d(x_test))[0] + # return self.gp.predict(np.atleast_2d(x_test))[0] + return self.gp.predict(np.atleast_2d(x_test), return_std=False) def scorer(self, inputs, outputs): """ diff --git a/tests/test_active.py b/tests/test_active.py index c8a66db..aabbfe6 100644 --- a/tests/test_active.py +++ b/tests/test_active.py @@ -65,10 +65,10 @@ def test_fit_03(self): weights = np.ones((15, 1)) / 15 ss = ActiveSubspaces(dim=1, n_boot=200) ss.fit(gradients=gradients, weights=weights) - true_evects = np.array([[-0.019091, -0.408566, 0.861223, -0.301669], - [-0.767799, 0.199069, 0.268823, 0.546434], - [-0.463451, -0.758442, -0.427696, -0.164486], - [-0.441965, 0.467131, -0.055723, -0.763774]]) + true_evects = np.array([[0.019091, -0.408566, -0.861223, -0.301669], + [0.767799, 0.199069, -0.268823, 0.546434], + [0.463451, -0.758442, 0.427696, -0.164486], + [0.441965, 0.467131, 0.055723, -0.763774]]) np.testing.assert_array_almost_equal(true_evects, ss.evects) def test_fit_04(self): @@ -86,10 +86,10 @@ def test_fit_05(self): outputs = np.random.uniform(0, 5, 15) ss = ActiveSubspaces(dim=1, method='local', n_boot=200) ss.fit(inputs=inputs, outputs=outputs) - true_evects = np.array([[-0.164383, 0.717021, -0.237246, -0.634486], - [-0.885808, 0.177628, 0.004112, 0.428691], - [0.255722, 0.558199, 0.734083, 0.290071], - [0.350612, 0.377813, -0.636254, 0.574029]]) + true_evects = np.array([[0.164383, 0.717021, 0.237246, -0.634486], + [0.885808, 0.177628, -0.004112, 0.428691], + [0.255722, -0.558199, 0.734083, -0.290071], + [0.350612, -0.377813, -0.636254, -0.574029]]) np.testing.assert_array_almost_equal(true_evects, ss.evects) def test_fit_06(self): @@ -228,7 +228,7 @@ def test_transform_01(self): ss = ActiveSubspaces(dim=2, method='local', n_boot=250) ss.fit(inputs=inputs, outputs=outputs) active = ss.transform(np.random.uniform(-1, 1, 8).reshape(2, 4))[0] - true_active = np.array([[-0.004748, 0.331107], [-0.949099, 0.347534]]) + true_active = np.array([[0.232762, 0.419052], [0.613532, 1.004439]]) np.testing.assert_array_almost_equal(true_active, active) def test_transform_02(self): @@ -238,8 +238,8 @@ def test_transform_02(self): ss = ActiveSubspaces(dim=2, method='local', n_boot=250) ss.fit(inputs=inputs, outputs=outputs) inactive = ss.transform(np.random.uniform(-1, 1, 8).reshape(2, 4))[1] - true_inactive = np.array([[-1.03574242, -0.04662904], - [-0.49850367, -0.37146678]]) + true_inactive = np.array([[-0.792408, -0.57175], + [-0.144381, 0.043564]]) np.testing.assert_array_almost_equal(true_inactive, inactive) def test_transform_03(self): @@ -271,8 +271,8 @@ def test_transform_05(self): ss = ActiveSubspaces(dim=2, method='local', n_boot=250) ss.fit(inputs=inputs, outputs=outputs) inactive = ss.transform(np.random.uniform(-1, 1, 8).reshape(2, 4))[1] - true_inactive = np.array([[-1.03574242, -0.04662904], - [-0.49850367, -0.37146678]]) + true_inactive = np.array([[-0.792408, -0.57175 ], + [-0.144381, 0.043564]]) np.testing.assert_array_almost_equal(true_inactive, inactive) def test_transform_06(self): @@ -480,7 +480,7 @@ def test_fit_bootstrap_ranges_02(self): weights = np.ones((30, 1)) / 30 ss = ActiveSubspaces(dim=1, method='exact', n_boot=100) ss.fit(gradients=gradients, weights=weights) - true_bounds_subspace = np.array([[0.00261813, 0.58863862, 0.99998352]]) + true_bounds_subspace = np.array([[0.002618, 0.290506, 0.648893]]) np.testing.assert_array_almost_equal(true_bounds_subspace, ss.subs_br) def test_fit_bootstrap_ranges_03(self): diff --git a/tests/test_feature_map.py b/tests/test_feature_map.py index 1a13fdd..88dfedf 100644 --- a/tests/test_feature_map.py +++ b/tests/test_feature_map.py @@ -182,11 +182,13 @@ def test_brute(self): folds=2, subspace=ss) best = fm.tune_pr_matrix(func=average_rrmse, - bounds=[slice(-2, 1, 0.2) for i in range(1)], + bounds=[slice(-2, 1, 0.2) for _ in range(1)], fn_args={'csv': csv}, maxiter=10, save_file=False)[1] - true = np.array([[0., 0.], [0., 0.], [0., 0.]]) + true = np.array([[-0.781768, -1.871064], + [-0.545585, -1.13183], + [1.961803, 0.95774]]) np.testing.assert_array_almost_equal(true, best) def test_dual_annealing(self): @@ -207,38 +209,39 @@ def test_dual_annealing(self): folds=2, subspace=ss) best = fm.tune_pr_matrix(func=average_rrmse, - bounds=[slice(-2, 1, 0.2) for i in range(1)], + bounds=[slice(-2, 1, 0.2) for _ in range(1)], fn_args={'csv': csv}, method='dual_annealing', maxiter=5, save_file=False)[1] - true = np.array([[0.21895786, 0.52423581], [-5.3186504, 14.9646475], - [4.2713126, 11.28870881]]) + true = np.array([[-18.205881, 0.130872], [-4.232873, 1.833328], + [-5.631037, 2.571455]]) np.testing.assert_array_almost_equal(true, best) - def test_bso(self): - np.random.seed(42) - inputs = np.random.uniform(-1, 1, 10).reshape(5, 2) - outputs = np.random.uniform(0, 5, 10).reshape(5, 2) - gradients = np.random.uniform(-1, 1, 20).reshape(5, 2, 2) - fm = FeatureMap(distr='laplace', - bias=np.random.uniform(-1, 1, 3), - input_dim=2, - n_features=3, - params=np.zeros(1), - sigma_f=outputs.var()) - ss = KernelActiveSubspaces(dim=1, feature_map=fm) - csv = CrossValidation(inputs=inputs, - outputs=outputs, - gradients=gradients, - folds=2, - subspace=ss) - best = fm.tune_pr_matrix(func=average_rrmse, - bounds=[slice(-2, 1, 0.2) for i in range(1)], - fn_args={'csv': csv}, - method='bso', - maxiter=10, - save_file=False)[1] - true = np.array([[14.9646475, 4.2713126], [11.28870881, 8.33313971], - [1.16475035, 9.92216877]]) - np.testing.assert_array_almost_equal(true, best) + # TODO: remove GPyOpt dependency with Emukit + # def test_bso(self): + # np.random.seed(42) + # inputs = np.random.uniform(-1, 1, 10).reshape(5, 2) + # outputs = np.random.uniform(0, 5, 10).reshape(5, 2) + # gradients = np.random.uniform(-1, 1, 20).reshape(5, 2, 2) + # fm = FeatureMap(distr='laplace', + # bias=np.random.uniform(-1, 1, 3), + # input_dim=2, + # n_features=3, + # params=np.zeros(1), + # sigma_f=outputs.var()) + # ss = KernelActiveSubspaces(dim=1, feature_map=fm) + # csv = CrossValidation(inputs=inputs, + # outputs=outputs, + # gradients=gradients, + # folds=2, + # subspace=ss) + # best = fm.tune_pr_matrix(func=average_rrmse, + # bounds=[slice(-2, 1, 0.2) for _ in range(1)], + # fn_args={'csv': csv}, + # method='bso', + # maxiter=10, + # save_file=False)[1] + # true = np.array([[14.9646475, 4.2713126], [11.28870881, 8.33313971], + # [1.16475035, 9.92216877]]) + # np.testing.assert_array_almost_equal(true, best) diff --git a/tests/test_nll.py b/tests/test_nll.py index 54777e5..90a11c2 100644 --- a/tests/test_nll.py +++ b/tests/test_nll.py @@ -20,17 +20,17 @@ def assert_plot_figures_added(): def read_data(): - data = np.loadtxt('tests/data/naca0012.txt', skiprows=1, delimiter=',') - real_inputs = data[:, 1:19] + data_in = np.loadtxt('tests/data/naca0012.txt', skiprows=1, delimiter=',') + real_inputs = data_in[:, 1:19] n_params = real_inputs.shape[1] lb = -0.01 * np.ones(n_params) ub = 0.01 * np.ones(n_params) normalizer = Normalizer(lb=lb, ub=ub) # inputs in [-1, 1] inputs = normalizer.fit_transform(real_inputs) - lift = data[:, 19] + lift = data_in[:, 19] # gradients with respect to normalized inputs - grad_lift = data[:, 21:39] + grad_lift = data_in[:, 21:39] return inputs, lift, grad_lift diff --git a/tests/test_utils.py b/tests/test_utils.py index ddf9e6b..1fbab74 100644 --- a/tests/test_utils.py +++ b/tests/test_utils.py @@ -96,14 +96,14 @@ def test_local_linear_gradients_05(self): def test_sort_eigpairs_evals(self): np.random.seed(42) matrix = np.random.uniform(-1, 1, 9).reshape(3, 3) - evals = sort_eigpairs(matrix)[0] - true_evals = np.array([[1.466942], [1.025235], [0.294945]]) + evals = sort_eigpairs(*np.linalg.eigh(matrix))[0] + true_evals = np.array([1.466942, 1.025235, 0.294945]) np.testing.assert_array_almost_equal(true_evals, evals) def test_sort_eigpairs_evects(self): np.random.seed(42) matrix = np.random.uniform(-1, 1, 9).reshape(3, 3) - evects = sort_eigpairs(matrix)[1] + evects = sort_eigpairs(*np.linalg.eigh(matrix))[1] true_evects = np.array([[0.511409, 0.515297, 0.687699], [-0.653819, -0.286001, 0.700517], [0.557657, -0.807881, 0.190647]]) @@ -183,7 +183,7 @@ def test_cross_validation_fit_01(self): folds=3, subspace=ss) csv.fit(inputs, gradients, outputs) - self.assertEqual(csv.gp.X.shape[1], 2) + self.assertEqual(csv.gp.X_train_.shape[1], 2) def test_cross_validation_fit_02(self): np.random.seed(42) @@ -210,7 +210,7 @@ def test_cross_validation_run_01(self): gradients=gradients, folds=2, subspace=ss) - true_value = (8.91180601306311, 6.806790947903309) + true_value = (9.696572, 6.45413) np.testing.assert_array_almost_equal(csv.run(), true_value) def test_cross_validation_run_02(self): @@ -230,7 +230,7 @@ def test_cross_validation_run_02(self): gradients=gradients, folds=2, subspace=ss) - true_value = (2.26333743325053, 0.43902733603381605) + true_value = (2.362181, 0.41159) np.testing.assert_array_almost_equal(csv.run(), true_value) def test_rrmse_01(self): @@ -296,5 +296,5 @@ def test_average_rrmse_02(self): best = [0.1, np.zeros((3, 2))] hyperparams = np.array([-1.]) score = average_rrmse(hyperparams, best, csv, verbose=True, resample=1) - true = 2.1357048222078703 + true = 7.409494 np.testing.assert_array_almost_equal(score, true) From f3eef4214a48d4e6707d97f80e4c615419678e71 Mon Sep 17 00:00:00 2001 From: mtezzele Date: Wed, 21 Dec 2022 18:22:09 -0600 Subject: [PATCH 18/26] rm gpy --- setup.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/setup.py b/setup.py index 0c6d461..bea898f 100644 --- a/setup.py +++ b/setup.py @@ -21,7 +21,7 @@ ) REQUIRED = [ - 'numpy', 'scipy', 'matplotlib', 'torch', 'GPy', 'GPyOpt', 'scikit-learn', 'scikit-learn-extra' + 'numpy', 'scipy', 'matplotlib', 'torch', 'GPyOpt', 'scikit-learn', 'scikit-learn-extra' ] EXTRAS = { From 4fe9d24afd9f12f7890dc56481c464390b476161 Mon Sep 17 00:00:00 2001 From: mtezzele Date: Wed, 21 Dec 2022 18:25:18 -0600 Subject: [PATCH 19/26] rm commented code --- athena/subspaces.py | 6 +++--- athena/utils.py | 6 ------ 2 files changed, 3 insertions(+), 9 deletions(-) diff --git a/athena/subspaces.py b/athena/subspaces.py index 48a0467..f184f97 100644 --- a/athena/subspaces.py +++ b/athena/subspaces.py @@ -319,7 +319,7 @@ def plot_eigenvalues(self, n_evals=None, filename=None, figsize=(8, 8), - title=''): # sourcery skip: class-extract-method + title=''): """ Plot the eigenvalues. @@ -396,7 +396,7 @@ def plot_eigenvectors(self, filename=None, figsize=None, labels=None, - title=''): # sourcery skip: extract-method + title=''): """ Plot the eigenvectors. @@ -458,7 +458,7 @@ def plot_sufficient_summary(self, outputs, filename=None, figsize=(10, 8), - title=''): # sourcery skip: extract-method + title=''): """ Plot the sufficient summary. diff --git a/athena/utils.py b/athena/utils.py index 0ea44e2..569c321 100644 --- a/athena/utils.py +++ b/athena/utils.py @@ -260,11 +260,6 @@ def fit(self, inputs, gradients, outputs): **self.kwargs) y = self.ss.transform(inputs)[0] - # build response surface - # kern = GPy.kern.RBF(input_dim=y.shape[1], ARD=True) - # self.gp = GPy.models.GPRegression(y, np.atleast_2d(outputs), kern) - # self.gp.optimize_restarts(15, verbose=False) - self.gp = GaussianProcessRegressor(n_restarts_optimizer=15) self.gp.fit(y, np.atleast_2d(outputs)) @@ -278,7 +273,6 @@ def predict(self, inputs): :rtype: numpy.ndarray """ x_test = self.ss.transform(inputs)[0] - # return self.gp.predict(np.atleast_2d(x_test))[0] return self.gp.predict(np.atleast_2d(x_test), return_std=False) def scorer(self, inputs, outputs): From c049af610209565a69c5656958435a7f42220f73 Mon Sep 17 00:00:00 2001 From: mtezzele Date: Fri, 23 Dec 2022 11:53:34 -0600 Subject: [PATCH 20/26] rm tests for py 3.7, add py 3.9 --- .github/workflows/testing_pr.yml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/.github/workflows/testing_pr.yml b/.github/workflows/testing_pr.yml index 8b7e185..f3b6960 100644 --- a/.github/workflows/testing_pr.yml +++ b/.github/workflows/testing_pr.yml @@ -13,7 +13,7 @@ jobs: fail-fast: false matrix: os: [windows-latest, macos-latest, ubuntu-latest] - python-version: [3.7, 3.8] + python-version: [3.8, 3.9] steps: - uses: actions/checkout@v2 From 4752bdf693a523cc2ad4dfc030836cbe81fb4aef Mon Sep 17 00:00:00 2001 From: mtezzele Date: Fri, 23 Dec 2022 12:26:17 -0600 Subject: [PATCH 21/26] rm gpy import, f-string in tut1 --- athena/utils.py | 5 ++--- tutorials/tutorial01/01_intro.py | 2 +- 2 files changed, 3 insertions(+), 4 deletions(-) diff --git a/athena/utils.py b/athena/utils.py index 569c321..c3609ed 100644 --- a/athena/utils.py +++ b/athena/utils.py @@ -3,7 +3,6 @@ import numpy as np import matplotlib.pyplot as plt from scipy.optimize import linprog -import GPy from sklearn.gaussian_process import GaussianProcessRegressor @@ -200,8 +199,8 @@ class CrossValidation(): :param dict kwargs: additional paramters organized in a dictionary to pass to subspace.fit method. For example 'weights' or 'metric'. - :cvar `GPy.models.GPRegression` gp: Gaussian process of the response surface - built with GPy library. + :cvar `sklearn.gaussian_process.GaussianProcessRegressor` gp: Gaussian + process of the response surface built with scikit-learn. """ def __init__(self, inputs, outputs, gradients, subspace, folds=5, **kwargs): diff --git a/tutorials/tutorial01/01_intro.py b/tutorials/tutorial01/01_intro.py index 1574715..bbeedc7 100644 --- a/tutorials/tutorial01/01_intro.py +++ b/tutorials/tutorial01/01_intro.py @@ -71,7 +71,7 @@ def output(x, normalizer, r): title = '2D sine' asub.plot_eigenvalues(figsize=(6, 4), title=title) -print("Eigenvalues: {}".format(np.squeeze(asub.evals))) +print(f"Eigenvalues: {np.squeeze(asub.evals)}") asub.plot_eigenvectors(figsize=(6, 4), title=title) asub.plot_sufficient_summary(x, f, figsize=(6, 4), title=title) From b4e8f4db63e99abae1e9f76c1b83d4d0dfaff7cd Mon Sep 17 00:00:00 2001 From: mtezzele Date: Thu, 12 Jan 2023 17:14:24 -0600 Subject: [PATCH 22/26] update to py38 and recent works --- README.md | 18 ++++++++++++++---- athena/meta.py | 2 +- setup.py | 2 +- 3 files changed, 16 insertions(+), 6 deletions(-) diff --git a/README.md b/README.md index 39af1e0..30caa4d 100644 --- a/README.md +++ b/README.md @@ -46,9 +46,9 @@ See the [**Examples and Tutorials**](#examples-and-tutorials) section below and ## Dependencies and installation -**ATHENA** requires `numpy`, `matplotlib`, `scipy`, `torch`, `GPy`, `GPyOpt`, +**ATHENA** requires `numpy`, `matplotlib`, `scipy`, `torch`, `GPyOpt`, `scikit-learn`, `scikit-learn-extra`, `sphinx` (for the documentation) and `pytest` (for local test). -The code is compatible with Python 3.6 and above. It can be installed directly +The code is compatible with Python 3.8 and above. It can be installed directly from the source code or via pip. ### Installing from pip @@ -143,21 +143,31 @@ To implement the numerical methods present in this package we followed these wor * Constantine et al. Python Active-subspaces Utility Library, Journal of Open Source Software, 1(5), 79, 2016. [[DOI](https://doi.org/10.21105/joss.00079)]. -* Romor, Tezzele, Lario, Rozza. *Kernel-based Active Subspaces with application to CFD parametric problems using Discontinuous Galerkin method*. 2020. [[arXiv](https://arxiv.org/abs/2008.12083)] +* Romor, Tezzele, Lario, Rozza. *Kernel-based Active Subspaces with application to CFD parametric problems using Discontinuous Galerkin method*. 2022. [[DOI](https://doi.org/10.1002/nme.7099)] * Zhang, Zhang, Hinkle. *Learning nonlinear level sets for dimensionality reduction in function approximation*. Advances in Neural Information Processing Systems 32: Annual Conference on Neural Information Processing Systems 2019, NeurIPS 2019, 8-14 December 2019, Vancouver, BC, Canada. [[arXiv](https://arxiv.org/abs/1902.10652)]. +* Romor, Tezzele, Rozza. *A local approach to parameter space reduction for regression and classification tasks*. 2021. [[arXiv](https://arxiv.org/abs/2107.10867)]. + ### Recent works with ATHENA Here there is a list of the scientific works involving **ATHENA** you can consult and/or cite. If you want to add one, please open a PR. +* Tezzele, Romor, Rozza. *Reduction in Parameter Space*, in Advanced Reduced Order Methods and Applications in Computational Fluid Dynamics. 2022. [[DOI](https://doi.org/10.1137/1.9781611977257.ch16)]. + +* Tezzele, Fabris, Sidari, Sicchiero, Rozza. *A multi-fidelity approach coupling parameter space reduction and non-intrusive POD with application to structural optimization of passenger ship hulls*. 2022. [[DOI](https://doi.org/10.1002/nme.7159)]. + +* Teng, Wang, Ju, Gruber, Zhang. *Level set learning with pseudo-reversible neural networks for nonlinear dimension reduction in function approximation*. 2021. [[arXiv](https://arxiv.org/abs/2112.01438)]. + * Romor, Tezzele, Rozza. *A local approach to parameter space reduction for regression and classification tasks*. 2021. [[arXiv](https://arxiv.org/abs/2107.10867)]. +* Meneghetti, Demo, Rozza. *A Dimensionality Reduction Approach for Convolutional Neural Networks*. 2021. [[arXiv](https://arxiv.org/abs/2110.09163)]. + * Demo, Tezzele, Mola, Rozza. *Hull shape design optimization with parameter space and model reductions and self-learning mesh morphing*. 2021. [[DOI](https://doi.org/10.3390/jmse9020185)], [[arXiv](https://arxiv.org/abs/2101.03781)]. * Romor, Tezzele, Rozza. *Multi-fidelity data fusion for the approximation of scalar functions with low intrinsic dimensionality through active subspaces*. 2020. [[DOI](https://doi.org/10.1002/pamm.202000349)], [[arXiv](https://arxiv.org/abs/2010.08349)]. -* Romor, Tezzele, Lario, Rozza. *Kernel-based Active Subspaces with application to CFD parametric problems using Discontinuous Galerkin method*. 2020. [[arXiv](https://arxiv.org/abs/2008.12083)]. +* Romor, Tezzele, Lario, Rozza. *Kernel-based Active Subspaces with application to CFD parametric problems using Discontinuous Galerkin method*. 2022. [[DOI](https://doi.org/10.1002/nme.7099)] * Demo, Tezzele, Rozza. *A supervised learning approach involving active subspaces for an efficient genetic algorithm in high-dimensional optimization problems*. 2020. [[DOI](https://doi.org/10.1137/20M1345219)], [[arXiv](https://arxiv.org/abs/2006.07282)]. diff --git a/athena/meta.py b/athena/meta.py index 882b268..f5d4dc8 100644 --- a/athena/meta.py +++ b/athena/meta.py @@ -9,7 +9,7 @@ __project__ = 'ATHENA' __title__ = "athena" __author__ = "Marco Tezzele, Francesco Romor" -__copyright__ = "Copyright 2019-2021, Athena contributors" +__copyright__ = "Copyright 2019-2023, Athena contributors" __license__ = "MIT" __version__ = "0.1.2" __mail__ = 'marcotez@gmail.com, francesco.romor@gmail.com' diff --git a/setup.py b/setup.py index bea898f..03cbfd0 100644 --- a/setup.py +++ b/setup.py @@ -53,7 +53,7 @@ classifiers=[ 'Development Status :: 5 - Production/Stable', 'License :: OSI Approved :: MIT License', - 'Programming Language :: Python :: 3.6', + 'Programming Language :: Python :: 3.8', 'Intended Audience :: Science/Research', 'Topic :: Scientific/Engineering :: Mathematics' ], From 971e00a37c1dbfe78d28df313490bb48b24292d9 Mon Sep 17 00:00:00 2001 From: mtezzele Date: Fri, 13 Jan 2023 17:06:32 -0600 Subject: [PATCH 23/26] update badges and fix CI name --- .github/workflows/ci.yml | 4 ++-- README.md | 13 ++++++------- 2 files changed, 8 insertions(+), 9 deletions(-) diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml index 9ccfdbc..145fa3f 100644 --- a/.github/workflows/ci.yml +++ b/.github/workflows/ci.yml @@ -1,4 +1,4 @@ -name: "Coverage Deploy to Codacy" +name: "CI" on: push: @@ -21,7 +21,7 @@ jobs: - name: Install Python dependencies run: | python3 -m pip install --upgrade pip - python3 -m pip install numpy scipy matplotlib pytest pytest-cov future torch gpy gpyopt scikit-learn scikit-learn-extra + python3 -m pip install numpy scipy matplotlib pytest pytest-cov future torch gpyopt scikit-learn scikit-learn-extra - name: Test with pytest env: diff --git a/README.md b/README.md index 30caa4d..994eb66 100644 --- a/README.md +++ b/README.md @@ -10,14 +10,13 @@ PyPI version - - Build Status + + Build Status - - Coverage Status - - - Codacy Badge + + Coverage Status + + Codacy Badge

From ce8fdcbd45fe3ac75144591099d28822e2adecf4 Mon Sep 17 00:00:00 2001 From: mtezzele Date: Wed, 1 Feb 2023 11:13:26 -0600 Subject: [PATCH 24/26] fix python version matrix --- .github/workflows/monthly-tag.yml | 2 +- .github/workflows/pypi-publish.yml | 4 ++-- 2 files changed, 3 insertions(+), 3 deletions(-) diff --git a/.github/workflows/monthly-tag.yml b/.github/workflows/monthly-tag.yml index 134388b..a2e9ea8 100644 --- a/.github/workflows/monthly-tag.yml +++ b/.github/workflows/monthly-tag.yml @@ -11,7 +11,7 @@ jobs: strategy: matrix: os: [windows-latest, macos-latest, ubuntu-latest] - python-version: [3.7, 3.8] + python-version: [3.8, 3.9] steps: - uses: actions/checkout@v2 - name: Set up Python diff --git a/.github/workflows/pypi-publish.yml b/.github/workflows/pypi-publish.yml index b2f8cdd..b1af411 100644 --- a/.github/workflows/pypi-publish.yml +++ b/.github/workflows/pypi-publish.yml @@ -12,10 +12,10 @@ jobs: steps: - uses: actions/checkout@v2 - - name: Set up Python 3.7 + - name: Set up Python 3.8 uses: actions/setup-python@v1 with: - python-version: 3.7 + python-version: 3.8 - name: Install build run: >- From 1430681a02a23a00491fa06638241d9db5bb0d7e Mon Sep 17 00:00:00 2001 From: mtezzele Date: Thu, 27 Apr 2023 16:56:00 -0500 Subject: [PATCH 25/26] restore old covariance method and sort eigs, fix tests --- athena/subspaces.py | 3 +-- athena/utils.py | 8 +++++--- tests/test_active.py | 28 ++++++++++++++-------------- tests/test_utils.py | 6 +++--- 4 files changed, 23 insertions(+), 22 deletions(-) diff --git a/athena/subspaces.py b/athena/subspaces.py index f184f97..d030c94 100644 --- a/athena/subspaces.py +++ b/athena/subspaces.py @@ -76,7 +76,7 @@ def _build_decompose_cov_matrix(self, np.dot(metric, gradients[i, :, :])) for i in range(gradients.shape[0]) ], axis=0)) - evals, evects = sort_eigpairs(*np.linalg.eigh(cov_matrix)) + evals, evects = sort_eigpairs(cov_matrix) return np.squeeze(evals), evects X = np.squeeze(gradients * np.sqrt(weights).reshape(-1, 1)) @@ -94,7 +94,6 @@ def _build_decompose_cov_matrix(self, evals = singular**2 - evals, evects = sort_eigpairs(evals, evects) return evals, evects.T def _compute_bootstrap_ranges(self, gradients, weights, metric=None): diff --git a/athena/utils.py b/athena/utils.py index c3609ed..40f729c 100644 --- a/athena/utils.py +++ b/athena/utils.py @@ -155,8 +155,8 @@ def local_linear_gradients(inputs, outputs, weights=None, n_neighbors=None): return gradients, new_inputs -def sort_eigpairs(evals, evects): - """Sort eigenpairs. +def sort_eigpairs(matrix): + """Compute eigenpairs and sort. :param numpy.ndarray evals: eigenvalues. :param numpy.ndarray evects: eigenvectors. @@ -170,6 +170,7 @@ def sort_eigpairs(evals, evects): the eigenvectors so that the first component of each eigenvector is positive. This normalization is very helpful for the bootstrapping. """ + evals, evects = np.linalg.eigh(matrix) evals = abs(evals) ind = np.argsort(evals) evals = evals[ind[::-1]] @@ -177,7 +178,8 @@ def sort_eigpairs(evals, evects): s = np.sign(evects[0, :]) s[s == 0] = 1 evects *= s - return evals, evects + + return evals.reshape(-1, 1), evects class CrossValidation(): diff --git a/tests/test_active.py b/tests/test_active.py index aabbfe6..58946b6 100644 --- a/tests/test_active.py +++ b/tests/test_active.py @@ -65,10 +65,10 @@ def test_fit_03(self): weights = np.ones((15, 1)) / 15 ss = ActiveSubspaces(dim=1, n_boot=200) ss.fit(gradients=gradients, weights=weights) - true_evects = np.array([[0.019091, -0.408566, -0.861223, -0.301669], - [0.767799, 0.199069, -0.268823, 0.546434], - [0.463451, -0.758442, 0.427696, -0.164486], - [0.441965, 0.467131, 0.055723, -0.763774]]) + true_evects = np.array([[-0.019091, 0.408566, 0.861223, 0.301669], + [-0.767799, -0.199069, 0.268823, -0.546434], + [-0.463451, 0.758442, -0.427696, 0.164486], + [-0.441965, -0.467131, -0.055723, 0.763774]]) np.testing.assert_array_almost_equal(true_evects, ss.evects) def test_fit_04(self): @@ -86,10 +86,10 @@ def test_fit_05(self): outputs = np.random.uniform(0, 5, 15) ss = ActiveSubspaces(dim=1, method='local', n_boot=200) ss.fit(inputs=inputs, outputs=outputs) - true_evects = np.array([[0.164383, 0.717021, 0.237246, -0.634486], - [0.885808, 0.177628, -0.004112, 0.428691], - [0.255722, -0.558199, 0.734083, -0.290071], - [0.350612, -0.377813, -0.636254, -0.574029]]) + true_evects = np.array([[-0.164383, -0.717021, -0.237246, 0.634486], + [-0.885808, -0.177628, 0.004112, -0.428691], + [ 0.255722, -0.558199, 0.734083, -0.290071], + [ 0.350612, -0.377813, -0.636254, -0.574029]]) np.testing.assert_array_almost_equal(true_evects, ss.evects) def test_fit_06(self): @@ -228,7 +228,7 @@ def test_transform_01(self): ss = ActiveSubspaces(dim=2, method='local', n_boot=250) ss.fit(inputs=inputs, outputs=outputs) active = ss.transform(np.random.uniform(-1, 1, 8).reshape(2, 4))[0] - true_active = np.array([[0.232762, 0.419052], [0.613532, 1.004439]]) + true_active = np.array([[-0.004748, -0.331107], [-0.949099, -0.347534]]) np.testing.assert_array_almost_equal(true_active, active) def test_transform_02(self): @@ -238,8 +238,8 @@ def test_transform_02(self): ss = ActiveSubspaces(dim=2, method='local', n_boot=250) ss.fit(inputs=inputs, outputs=outputs) inactive = ss.transform(np.random.uniform(-1, 1, 8).reshape(2, 4))[1] - true_inactive = np.array([[-0.792408, -0.57175], - [-0.144381, 0.043564]]) + true_inactive = np.array([[-1.035742, 0.046629], + [-0.498504, 0.371467]]) np.testing.assert_array_almost_equal(true_inactive, inactive) def test_transform_03(self): @@ -271,8 +271,8 @@ def test_transform_05(self): ss = ActiveSubspaces(dim=2, method='local', n_boot=250) ss.fit(inputs=inputs, outputs=outputs) inactive = ss.transform(np.random.uniform(-1, 1, 8).reshape(2, 4))[1] - true_inactive = np.array([[-0.792408, -0.57175 ], - [-0.144381, 0.043564]]) + true_inactive = np.array([[-1.035742, 0.046629], + [-0.498504, 0.371467]]) np.testing.assert_array_almost_equal(true_inactive, inactive) def test_transform_06(self): @@ -480,7 +480,7 @@ def test_fit_bootstrap_ranges_02(self): weights = np.ones((30, 1)) / 30 ss = ActiveSubspaces(dim=1, method='exact', n_boot=100) ss.fit(gradients=gradients, weights=weights) - true_bounds_subspace = np.array([[0.002618, 0.290506, 0.648893]]) + true_bounds_subspace = np.array([[0.002618, 0.588639, 0.999984]]) np.testing.assert_array_almost_equal(true_bounds_subspace, ss.subs_br) def test_fit_bootstrap_ranges_03(self): diff --git a/tests/test_utils.py b/tests/test_utils.py index 1fbab74..463cc7a 100644 --- a/tests/test_utils.py +++ b/tests/test_utils.py @@ -96,14 +96,14 @@ def test_local_linear_gradients_05(self): def test_sort_eigpairs_evals(self): np.random.seed(42) matrix = np.random.uniform(-1, 1, 9).reshape(3, 3) - evals = sort_eigpairs(*np.linalg.eigh(matrix))[0] - true_evals = np.array([1.466942, 1.025235, 0.294945]) + evals = sort_eigpairs(matrix)[0] + true_evals = np.array([[1.466942], [1.025235], [0.294945]]) np.testing.assert_array_almost_equal(true_evals, evals) def test_sort_eigpairs_evects(self): np.random.seed(42) matrix = np.random.uniform(-1, 1, 9).reshape(3, 3) - evects = sort_eigpairs(*np.linalg.eigh(matrix))[1] + evects = sort_eigpairs(matrix)[1] true_evects = np.array([[0.511409, 0.515297, 0.687699], [-0.653819, -0.286001, 0.700517], [0.557657, -0.807881, 0.190647]]) From f88152a7f7467cd2ac93d56c9e2177a950acb8c1 Mon Sep 17 00:00:00 2001 From: Giovanni Stabile Date: Mon, 29 Jun 2026 10:01:02 +0200 Subject: [PATCH 26/26] Fix label formatting in subspaces.py There was an extract } that was preventing the plot in tutorial to work --- athena/subspaces.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/athena/subspaces.py b/athena/subspaces.py index d030c94..8e9bff1 100644 --- a/athena/subspaces.py +++ b/athena/subspaces.py @@ -489,7 +489,7 @@ def plot_sufficient_summary(self, s=40, alpha=0.9, edgecolors='k') - plt.xlabel('Active variable ' + r'$W_1^T \mathbf{\mu}}$', + plt.xlabel('Active variable ' + r'$W_1^T \mathbf{\mu}$', fontsize=18) plt.ylabel(r'$f \, (\mathbf{\mu})$', fontsize=18) elif self.dim == 2: