Merge pull request #40 from fdtomasi/py310

Support Python 3.10

.github/workflows/python-package.yml

@@ -16,7 +16,7 @@ jobs:
     strategy:
       fail-fast: false
       matrix:
-        python-version: [ "3.7", "3.8" ]
+        python-version: [ "3.8", "3.9", "3.10" ]
 
     steps:
       - uses: actions/checkout@v3

regain/__init__.py

@@ -28,4 +28,4 @@
 # OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
 # OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
 
-__version__ = "0.3.8"
+__version__ = "0.3.9"

regain/bayesian/gwishart_inference.py

@@ -39,8 +39,7 @@
 from scipy import linalg
 from scipy.optimize import minimize
 from scipy.special import comb
-from sklearn.covariance import empirical_covariance
-from sklearn.covariance.empirical_covariance_ import log_likelihood
+from sklearn.covariance import empirical_covariance, log_likelihood
 from sklearn.linear_model import LassoLars
 from sklearn.utils import Bunch, check_array
 from sklearn.utils.extmath import fast_logdet
@@ -51,9 +50,11 @@
 def mk_all_ugs(n_dim):
     """Utility for generating all possible graphs."""
     nedges = int(comb(n_dim, 2))
-    m = 2 ** nedges
+    m = 2**nedges
 
-    ind = np.array([list(binary_repr(x, width=len(binary_repr(m - 1)))) for x in range(m)]).astype(int)
+    ind = np.array(
+        [list(binary_repr(x, width=len(binary_repr(m - 1)))) for x in range(m)]
+    ).astype(int)
     ord = np.argsort(ind.sum(axis=1))
     ind = ind[ord]
 
@@ -98,7 +99,12 @@ def score_blankets(blankets, X, alphas=(0.01, 0.5, 1)):
                 X_mb = np.zeros((X.shape[0], 1))
 
             y_mb = X[:, i]
-            score = np.sum([LassoLars(alpha=alpha).fit(X_mb, y_mb).score(X_mb, y_mb) for alpha in alphas])
+            score = np.sum(
+                [
+                    LassoLars(alpha=alpha).fit(X_mb, y_mb).score(X_mb, y_mb)
+                    for alpha in alphas
+                ]
+            )
 
             scores.append(score)
         scores_all.append(scores)
@@ -109,7 +115,12 @@ def score_blankets(blankets, X, alphas=(0.01, 0.5, 1)):
 
 
 def _get_graphs(blankets, scores, n_dim, n_resampling=200):
-    idx = np.array([np.random.choice(scores.shape[1], p=scores[i], size=n_resampling) for i in range(n_dim)])
+    idx = np.array(
+        [
+            np.random.choice(scores.shape[1], p=scores[i], size=n_resampling)
+            for i in range(n_dim)
+        ]
+    )
 
     graphs_ = np.array([blankets[i][idx[i]] for i in range(n_dim)]).transpose(1, 0, 2)
     # symmetrise with AND operator -> product
@@ -268,7 +279,11 @@ def compute_score(X, G, P, S, GWprior=None, score_method="bic"):
 
         logdetHdiag = sum(np.log(-diagH))
         lognormconst = dof * np.log(2 * np.pi) / 2 + logh - logdetHdiag / 2.0
-        score = lognormconst - GWprior.lognormconst - n_samples * n_dim * np.log(2 * np.pi) / 2
+        score = (
+            lognormconst
+            - GWprior.lognormconst
+            - n_samples * n_dim * np.log(2 * np.pi) / 2
+        )
         GWpost.lognormconst = lognormconst
 
     elif score_method == "laplace":
@@ -294,7 +309,11 @@ def compute_score(X, G, P, S, GWprior=None, score_method="bic"):
         # neg Hessian will be posdef
         logdetH = 2 * sum(np.log(np.diag(linalg.cholesky(-H))))
         lognormconst = dof * np.log(2 * np.pi) / 2 + logh - logdetH / 2.0
-        score = lognormconst - GWprior.lognormconst - n_samples * n_dim * np.log(2 * np.pi) / 2
+        score = (
+            lognormconst
+            - GWprior.lognormconst
+            - n_samples * n_dim * np.log(2 * np.pi) / 2
+        )
         GWpost.lognormconst = lognormconst
 
     GWpost.score = score
@@ -315,7 +334,9 @@ def GWishartScore(X, G, d0=3, S0=None, score_method="bic", mode="covsel"):
         noData = np.zeros((0, n_dim))
 
         P0, S_noData = GWishartFit(noData, G, GWprior)
-        GWtemp = compute_score(noData, G, P0, S_noData, GWprior=GWprior, score_method=score_method)
+        GWtemp = compute_score(
+            noData, G, P0, S_noData, GWprior=GWprior, score_method=score_method
+        )
         GWprior.lognormconst = GWtemp.lognormconst
 
     # Compute the map precision matrix P
@@ -344,13 +365,17 @@ def bayesian_graphical_lasso(
         alphas = np.logspace(-2, 0, 20)
 
     # get a series of Markov blankets for vaiours alphas
-    mdl = GraphicalLasso(assume_centered=assume_centered, tol=tol, max_iter=max_iter, verbose=False)
+    mdl = GraphicalLasso(
+        assume_centered=assume_centered, tol=tol, max_iter=max_iter, verbose=False
+    )
     precisions = [mdl.set_params(alpha=a).fit(X).precision_ for a in alphas]
     mblankets = markov_blankets(precisions, tol=tol, unique=True)
 
     normalized_scores = score_blankets(mblankets, X=X, alphas=[0.01, 0.5, 1])
 
-    graphs = _get_graphs(mblankets, normalized_scores, n_dim=n_dim, n_resampling=n_resampling)
+    graphs = _get_graphs(
+        mblankets, normalized_scores, n_dim=n_dim, n_resampling=n_resampling
+    )
 
     nonzeros_all = [np.triu(g, 1) + np.eye(n_dim, dtype=bool) for g in graphs]
 
@@ -361,7 +386,10 @@ def bayesian_graphical_lasso(
     # Find non-zero elements of upper triangle of G
     # make sure diagonal is non-zero
     # G = nonzeros_all[1] # probably can discard if all zeros?
-    res = [GWishartScore(X, G, d0=d0, S0=S0, mode=mode, score_method=scoring) for G in nonzeros_all]
+    res = [
+        GWishartScore(X, G, d0=d0, S0=S0, mode=mode, score_method=scoring)
+        for G in nonzeros_all
+    ]
 
     top_n = [x.P for x in sorted(res, key=lambda x: x.score)[::-1][:top_n]]
     return np.mean(top_n, axis=0)
@@ -439,7 +467,12 @@ def __init__(
         top_n=1,
     ):
         super(GraphicalLasso, self).__init__(
-            alpha=alpha, tol=tol, max_iter=max_iter, verbose=verbose, assume_centered=assume_centered, mode=mode
+            alpha=alpha,
+            tol=tol,
+            max_iter=max_iter,
+            verbose=verbose,
+            assume_centered=assume_centered,
+            mode=mode,
         )
         self.alphas = alphas
         self.n_resampling = n_resampling

regain/bayesian/wishart_process_.py

@@ -34,9 +34,9 @@
 import numpy as np
 from scipy import linalg, stats
 from sklearn.covariance import empirical_covariance
-from sklearn.datasets.base import Bunch
 from sklearn.gaussian_process import kernels
 from sklearn.metrics.pairwise import rbf_kernel
+from sklearn.utils import Bunch
 from sklearn.utils.validation import check_X_y
 from tqdm import trange
 
@@ -333,17 +333,17 @@ def fit(self, X, y):
             samples_u, loglikes, lps = out
 
         # Burn in
-        self.lps_after_burnin = lps[self.burn_in :]
-        self.samples_u_after_burnin = samples_u[self.burn_in :]
-        self.loglikes_after_burnin = loglikes[self.burn_in :]
+        self.lps_after_burnin = lps[self.burn_in:]
+        self.samples_u_after_burnin = samples_u[self.burn_in:]
+        self.loglikes_after_burnin = loglikes[self.burn_in:]
 
         # % Select the best hyperparameters based on the loglikes_after_burnin
         pos = np.argmax(self.loglikes_after_burnin)
         self.lmap = self.lps_after_burnin[pos]
         self.u_map = self.samples_u_after_burnin[pos]
 
         if self.learn_ell:
-            self.Ls_after_burnin = Ls[self.burn_in :]
+            self.Ls_after_burnin = Ls[self.burn_in:]
             self.Lmap = self.Ls_after_burnin[pos]
         else:
             self.Lmap = L

regain/covariance/graphical_lasso_.py

@@ -44,7 +44,7 @@
 from regain.norm import l1_od_norm
 from regain.prox import prox_logdet, soft_thresholding_off_diagonal
 from regain.update_rules import update_rho
-from regain.utils import convergence
+from regain.utils import Convergence
 
 
 def logl(emp_cov, precision):
@@ -169,7 +169,7 @@ def graphical_lasso(
         obj = objective(emp_cov, K, Z, alpha) if compute_objective else np.nan
         rnorm = np.linalg.norm(K - Z, "fro")
         snorm = rho * np.linalg.norm(Z - Z_old, "fro")
-        check = convergence(
+        check = Convergence(
             obj=obj,
             rnorm=rnorm,
             snorm=snorm,
@@ -180,10 +180,7 @@ def graphical_lasso(
 
         Z_old = Z.copy()
         if verbose:
-            print(
-                "obj: %.4f, rnorm: %.4f, snorm: %.4f,"
-                "eps_pri: %.4f, eps_dual: %.4f" % check[:5]
-            )
+            print(check)
 
         checks.append(check)
         if check.rnorm <= check.e_pri and check.snorm <= check.e_dual:

regain/covariance/infimal_convolution_.py

@@ -40,7 +40,7 @@
 from regain.norm import l1_od_norm
 from regain.prox import prox_laplacian, prox_trace_indicator, soft_thresholding
 from regain.update_rules import update_rho
-from regain.utils import convergence
+from regain.utils import Convergence
 
 
 def objective(S, R, K, L, alpha, tau):
@@ -136,24 +136,27 @@ def infimal_convolution(
         obj = objective(S, R, K, L, alpha, tau) if compute_objective else np.nan
         rnorm = np.linalg.norm(R - K + L)
         snorm = rho * np.linalg.norm(R - R_old)
-        check = convergence(
+        check = Convergence(
             obj=obj,
             rnorm=rnorm,
             snorm=snorm,
-            e_pri=np.sqrt(R.size) * tol + rtol * max(np.linalg.norm(R), np.linalg.norm(K - L)),
+            e_pri=np.sqrt(R.size) * tol
+            + rtol * max(np.linalg.norm(R), np.linalg.norm(K - L)),
             e_dual=np.sqrt(R.size) * tol + rtol * rho * np.linalg.norm(U),
         )
         R_old = R.copy()
 
         if verbose:
-            print("obj: %.4f, rnorm: %.4f, snorm: %.4f," "eps_pri: %.4f, eps_dual: %.4f" % check[:5])
+            print(check)
 
         checks.append(check)
         if check.rnorm <= check.e_pri and check.snorm <= check.e_dual:
             break
         if check.obj == np.inf:
             break
-        rho_new = update_rho(rho, rnorm, snorm, iteration=iteration_, **(update_rho_options or {}))
+        rho_new = update_rho(
+            rho, rnorm, snorm, iteration=iteration_, **(update_rho_options or {})
+        )
         # scaled dual variables should be also rescaled
         U *= rho / rho_new
         rho = rho_new

regain/covariance/kernel_latent_time_graphical_lasso_.py

@@ -49,7 +49,7 @@
 )
 from regain.prox import prox_logdet, prox_trace_indicator, soft_thresholding
 from regain.update_rules import update_rho
-from regain.utils import convergence
+from regain.utils import Convergence
 from regain.validation import check_norm_prox
 
 
@@ -311,7 +311,7 @@ def kernel_latent_time_graphical_lasso(
             else np.nan
         )
 
-        check = convergence(
+        check = Convergence(
             obj=obj,
             rnorm=rnorm,
             snorm=snorm,
@@ -360,10 +360,7 @@ def kernel_latent_time_graphical_lasso(
             W_M_old[m] = (W_M[m][0].copy(), W_M[m][1].copy())
 
         if verbose:
-            print(
-                "obj: %.4f, rnorm: %.4f, snorm: %.4f,"
-                "eps_pri: %.4f, eps_dual: %.4f" % check[:5]
-            )
+            print(check)
 
         checks.append(check)
         if check.rnorm <= check.e_pri and check.snorm <= check.e_dual:

regain/covariance/kernel_time_graphical_lasso_.py

@@ -50,7 +50,7 @@
 from regain.norm import l1_od_norm
 from regain.prox import prox_logdet, soft_thresholding
 from regain.update_rules import update_rho
-from regain.utils import convergence
+from regain.utils import Convergence
 from regain.validation import check_norm_prox
 
 
@@ -163,7 +163,7 @@ def kernel_time_graphical_lasso(
         n_samples = np.ones(n_times)
 
     checks = [
-        convergence(
+        Convergence(
             obj=objective(n_samples, emp_cov, Z_0, Z_0, Z_M, alpha, kernel, psi)
         )
     ]
@@ -245,7 +245,7 @@ def kernel_time_graphical_lasso(
             else np.nan
         )
 
-        check = convergence(
+        check = Convergence(
             obj=obj,
             rnorm=rnorm,
             snorm=snorm,
@@ -283,10 +283,7 @@ def kernel_time_graphical_lasso(
             Z_M_old[m] = (Z_M[m][0].copy(), Z_M[m][1].copy())
 
         if verbose:
-            print(
-                "obj: %.4f, rnorm: %.4f, snorm: %.4f,"
-                "eps_pri: %.4f, eps_dual: %.4f" % check[:5]
-            )
+            print(check)
 
         checks.append(check)
         if stop_at is not None:

regain/covariance/latent_graphical_lasso_.py

@@ -34,13 +34,15 @@
 
 import numpy as np
 from scipy import linalg
-from six.moves import range
 
-from regain.covariance.graphical_lasso_ import GraphicalLasso, init_precision
-from regain.covariance.graphical_lasso_ import objective as obj_gl
+from regain.covariance.graphical_lasso_ import (
+    GraphicalLasso,
+    init_precision,
+    objective as obj_gl,
+)
 from regain.prox import prox_logdet, prox_trace_indicator, soft_thresholding
 from regain.update_rules import update_rho
-from regain.utils import convergence
+from regain.utils import Convergence
 
 
 def objective(emp_cov, R, K, L, alpha, tau):
@@ -144,24 +146,27 @@ def latent_graphical_lasso(
         obj = objective(emp_cov, R, K, L, alpha, tau) if compute_objective else np.nan
         rnorm = np.linalg.norm(R - K + L)
         snorm = rho * np.linalg.norm(R - R_old)
-        check = convergence(
+        check = Convergence(
             obj=obj,
             rnorm=rnorm,
             snorm=snorm,
-            e_pri=np.sqrt(R.size) * tol + rtol * max(np.linalg.norm(R), np.linalg.norm(K - L)),
+            e_pri=np.sqrt(R.size) * tol
+            + rtol * max(np.linalg.norm(R), np.linalg.norm(K - L)),
             e_dual=np.sqrt(R.size) * tol + rtol * rho * np.linalg.norm(U),
         )
         R_old = R.copy()
 
         if verbose:
-            print("obj: %.4f, rnorm: %.4f, snorm: %.4f," "eps_pri: %.4f, eps_dual: %.4f" % check[:5])
+            print(check)
 
         checks.append(check)
         if check.rnorm <= check.e_pri and check.snorm <= check.e_dual:
             break
         if check.obj == np.inf:
             break
-        rho_new = update_rho(rho, rnorm, snorm, iteration=iteration_, **(update_rho_options or {}))
+        rho_new = update_rho(
+            rho, rnorm, snorm, iteration=iteration_, **(update_rho_options or {})
+        )
         # scaled dual variables should be also rescaled
         U *= rho / rho_new
         rho = rho_new
@@ -293,7 +298,12 @@ def _fit(self, emp_cov):
             Empirical covariance of data.
 
         """
-        self.precision_, self.latent_, self.covariance_, self.n_iter_ = latent_graphical_lasso(
+        (
+            self.precision_,
+            self.latent_,
+            self.covariance_,
+            self.n_iter_,
+        ) = latent_graphical_lasso(
             emp_cov,
             alpha=self.alpha,
             tau=self.tau,