Merge branch 'main' of https://github.com/PaulaGarciaMolina/seemps2 i…

…nto main

src/seemps/analysis/chebyshev.py

@@ -233,7 +233,7 @@ def cheb2mps(
                 strategy=strategy,
             )
             logger(
-                f"MPS Clenshaw step {i+1}/{steps} with maximum bond dimension {max(y_i.bond_dimensions())} and error {y_i.error():6e}"
+                f"MPS Clenshaw step {i+1}/{steps} with maximum bond dimension {y_i.max_bond_dimension()} and error {y_i.error():6e}"
             )
         f_mps = simplify(
             MPSSum(
@@ -270,9 +270,10 @@ def cheb2mps(
                 strategy=strategy,
             )
             logger(
-                f"MPS expansion step {i+1}/{steps} with maximum bond dimension {max(f_mps.bond_dimensions())} and error {f_mps.error():6e}"
+                f"MPS expansion step {i+1}/{steps} with maximum bond dimension {f_mps.max_bond_dimension()} and error {f_mps.error():6e}"
             )
             T_i, T_i_plus_1 = T_i_plus_1, T_i_plus_2
+    logger.close()
     return f_mps
 
 
@@ -330,7 +331,7 @@ def cheb2mpo(
                 strategy=strategy,
             )
             logger(
-                f"MPO Clenshaw step {i+1}/{steps} with maximum bond dimension {max(y_i.bond_dimensions())}"
+                f"MPO Clenshaw step {i+1}/{steps} with maximum bond dimension {y_i.max_bond_dimension()}"
             )
         f_mpo = simplify_mpo(
             MPOSum([y_i, MPOList([initial_mpo, y_i_plus_1])], weights=[1, -1]),
@@ -354,7 +355,7 @@ def cheb2mpo(
                 strategy=strategy,
             )
             logger(
-                f"MPO expansion step {i+1}/{steps} with maximum bond dimension {max(f_mpo.bond_dimensions())}"
+                f"MPO expansion step {i+1}/{steps} with maximum bond dimension {f_mpo.max_bond_dimension()}"
             )
             T_i, T_i_plus_1 = T_i_plus_1, T_i_plus_2
     return f_mpo

src/seemps/analysis/cross/cross.py

@@ -9,7 +9,7 @@
 from .black_box import BlackBox
 from ..sampling import evaluate_mps, random_mps_indices
 from ...state import MPS, random_mps
-from ...tools import make_logger
+from ...tools import Logger
 from ...typing import VectorLike
 
 
@@ -188,7 +188,8 @@ def _check_convergence(
     cross: CrossInterpolation,
     sweep: int,
     cross_strategy: CrossStrategy,
-) -> tuple[bool, str]:
+    logger: Logger,
+) -> bool:
     allowed_sampling_indices = getattr(
         cross.black_box, "allowed_sampling_indices", None
     )
@@ -198,24 +199,23 @@ def _check_convergence(
         allowed_indices=allowed_sampling_indices,
         rng=cross_strategy.rng,
     )
-    maxbond = max(cross.mps.bond_dimensions())
-    logger = make_logger()
-    logger(
-        f"Cross sweep {1+sweep:3d} with error({cross_strategy.num_samples} samples "
-        f"in norm-{cross_strategy.norm_sampling})={error}, maxbond={maxbond}, evals(cumulative)={cross.black_box.evals}"
-    )
-    converged = False
-    message = f"Maximum number of sweeps {cross_strategy.maxiter} reached"
+    # TODO: Use max_bond_dimension() from MPS
+    maxbond = cross.mps.max_bond_dimension()
+    if logger:
+        logger(
+            f"Cross sweep {1+sweep:3d} with error({cross_strategy.num_samples} samples "
+            f"in norm-{cross_strategy.norm_sampling})={error}, maxbond={maxbond}, evals(cumulative)={cross.black_box.evals}"
+        )
     if cross_strategy.check_norm_2:
         change_norm = cross.norm_2_change()
         logger(f"Norm-2 change {change_norm}")
         if change_norm <= cross_strategy.tol_norm_2:
-            converged = True
-            message = f"Stationary state reached with norm-2 change {change_norm}"
+            logger(f"Stationary state reached with norm-2 change {change_norm}")
+            return True
     if error < cross_strategy.tol_sampling:
-        converged = True
-        message = f"State converged within tolerance {cross_strategy.tol_sampling}"
+        logger(f"State converged within tolerance {cross_strategy.tol_sampling}")
+        return True
     elif maxbond > cross_strategy.maxbond:
-        converged = True
-        message = f"Maxbond reached above the threshold {cross_strategy.maxbond}"
-    return converged, message
+        logger(f"Maxbond reached above the threshold {cross_strategy.maxbond}")
+        return True
+    return False

src/seemps/analysis/cross/cross_dmrg.py

@@ -14,7 +14,7 @@
 from ...state.schmidt import _destructive_svd
 from ...state.core import destructively_truncate_vector
 from ...truncate import SIMPLIFICATION_STRATEGY
-from ...tools import log
+from ...tools import make_logger
 
 DEFAULT_CROSS_STRATEGY = SIMPLIFICATION_STRATEGY.replace(
     normalize=False,
@@ -81,21 +81,21 @@ def cross_dmrg(
         low=0, high=black_box.base, size=black_box.sites
     )
     cross = CrossInterpolationDMRG(black_box, initial_point)
-    for i in range(cross_strategy.maxiter):
-        # Forward sweep
-        for k in range(cross.sites - 1):
-            _update_dmrg(cross, k, True, cross_strategy)
-        converged, message = _check_convergence(cross, i, cross_strategy)
-        if converged:
-            break
-        # Backward sweep
-        for k in reversed(range(cross.sites - 1)):
-            _update_dmrg(cross, k, False, cross_strategy)
-        converged, message = _check_convergence(cross, i, cross_strategy)
-        if converged:
-            break
-
-    log(message)
+    converged = False
+    with make_logger(2) as logger:
+        for i in range(cross_strategy.maxiter):
+            # Forward sweep
+            for k in range(cross.sites - 1):
+                _update_dmrg(cross, k, True, cross_strategy)
+            if converged := _check_convergence(cross, i, cross_strategy, logger):
+                break
+            # Backward sweep
+            for k in reversed(range(cross.sites - 1)):
+                _update_dmrg(cross, k, False, cross_strategy)
+            if converged := _check_convergence(cross, i, cross_strategy, logger):
+                break
+        if not converged:
+            logger("Maximum number of TT-Cross iterations reached")
     return CrossResults(mps=cross.mps, evals=black_box.evals)
 
 

src/seemps/analysis/cross/cross_maxvol.py

@@ -11,7 +11,7 @@
     _check_convergence,
 )
 from ...state._contractions import _contract_last_and_first
-from ...tools import log
+from ...tools import make_logger
 
 
 @dataclasses.dataclass
@@ -105,17 +105,19 @@ def cross_maxvol(
         low=0, high=black_box.base, size=black_box.sites
     )
     cross = CrossInterpolationMaxvol(black_box, initial_point)
-    for i in range(cross_strategy.maxiter):
-        # Forward sweep
-        for k in range(cross.sites):
-            _update_maxvol(cross, k, True, cross_strategy)
-        # Backward sweep
-        for k in reversed(range(cross.sites)):
-            _update_maxvol(cross, k, False, cross_strategy)
-        converged, message = _check_convergence(cross, i, cross_strategy)
-        if converged:
-            break
-    log(message)
+    converged = False
+    with make_logger(2) as logger:
+        for i in range(cross_strategy.maxiter):
+            # Forward sweep
+            for k in range(cross.sites):
+                _update_maxvol(cross, k, True, cross_strategy)
+            # Backward sweep
+            for k in reversed(range(cross.sites)):
+                _update_maxvol(cross, k, False, cross_strategy)
+            if converged := _check_convergence(cross, i, cross_strategy, logger):
+                break
+        if not converged:
+            logger("Maximum number of iterations reached")
     return CrossResults(mps=cross.mps, evals=black_box.evals)
 
 

src/seemps/operators/mpo.py

@@ -134,6 +134,10 @@ def bond_dimensions(self) -> list[int]:
         """Return the bond dimensions of the MPO."""
         return [A.shape[-1] for A in self._data][:-1]
 
+    def max_bond_dimension(self) -> int:
+        """Return the largest bond dimension."""
+        return max(A.shape[0] for A in self._data)
+
     @property
     def T(self) -> MPO:
         output = self.copy()

src/seemps/optimization/arnoldi.py

@@ -3,7 +3,7 @@
 import numpy as np
 import scipy.linalg  # type: ignore
 from numpy.typing import NDArray
-from ..tools import log
+from ..tools import make_logger
 from ..state import (
     MPS,
     CanonicalMPS,
@@ -205,7 +205,8 @@ def arnoldi_eigh(
     arnoldi = MPSArnoldiRepresentation(
         operator, strategy, tol_ill_conditioning=tol_ill, gamma=gamma
     )
-    log(f"Arnoldi expansion with maxiter={maxiter}, relative tolerance={tol}")
+    logger = make_logger()
+    logger(f"Arnoldi expansion with maxiter={maxiter}, relative tolerance={tol}")
     v: MPS = arnoldi.restart_with_vector(guess)
     energy, variance = arnoldi.energy_and_variance()
     results = OptimizeResults(
@@ -216,7 +217,7 @@ def arnoldi_eigh(
         converged=False,
         message=f"Exceeded maximum number of steps {maxiter}",
     )
-    log(f"start, energy={energy}, variance={variance}")
+    logger(f"start, energy={energy}, variance={variance}")
     if callback is not None:
         callback(arnoldi.eigenvector(), results)
     last_energy = energy
@@ -233,7 +234,7 @@ def arnoldi_eigh(
         results.variances.append(variance)
         if energy < results.energy:
             results.energy, results.state = energy, arnoldi.eigenvector()
-        log(f"step={step}, energy={energy}, variance={variance}")
+        logger(f"step={step}, energy={energy}, variance={variance}")
         if callback is not None:
             callback(arnoldi.eigenvector(), results)
         if len(arnoldi.V) == 1:
@@ -252,8 +253,9 @@ def arnoldi_eigh(
                 results.converged = True
                 break
             last_energy = energy
-    log(
+    logger(
         f"Arnoldi finished with {step} iterations:\n"
         f"message = {results.message}\nconverged = {results.converged}"
     )
+    logger.close()
     return results

src/seemps/optimization/power.py

@@ -141,4 +141,5 @@ def cgs_callback(state, residual):
             state = simplify(H_v, strategy=strategy)
             total_steps += 1
     logger(f"power_method() finished with results\n{results}")
+    logger.close()
     return results

src/seemps/state/canonical_mps.py

@@ -5,8 +5,6 @@
 from typing import Optional, Sequence, Iterable
 from ..typing import Vector, Tensor3, Tensor4, VectorLike, Environment
 from .schmidt import (
-    _ortho_left,
-    _ortho_right,
     _vector2mps,
     _schmidt_weights,
     _left_orth_2site,
@@ -17,57 +15,16 @@
     _update_left_environment,
     _update_right_environment,
 )
-from ._contractions import _contract_last_and_first
-from .core import DEFAULT_STRATEGY, Strategy
+from .core import (
+    DEFAULT_STRATEGY,
+    Strategy,
+    _update_in_canonical_form_right,
+    _update_in_canonical_form_left,
+    _canonicalize,
+)
 from .mps import MPS
 
 
-# TODO: Replace einsum by a more efficient form
-def _update_in_canonical_form_right(
-    Ψ: list[Tensor3], A: Tensor3, site: int, truncation: Strategy
-) -> tuple[int, float]:
-    """Insert a tensor in canonical form into the MPS Ψ at the given site.
-    Update the neighboring sites in the process."""
-    if site + 1 == len(Ψ):
-        Ψ[site] = A
-        return site, 0.0
-    Ψ[site], sV, err = _ortho_right(A, truncation)
-    site += 1
-    # np.einsum("ab,bic->aic", sV, Ψ[site])
-    Ψ[site] = _contract_last_and_first(sV, Ψ[site])
-    return site, err
-
-
-# TODO: Replace einsum by a more efficient form
-def _update_in_canonical_form_left(
-    Ψ: list[Tensor3], A: Tensor3, site: int, truncation: Strategy
-) -> tuple[int, float]:
-    """Insert a tensor in canonical form into the MPS Ψ at the given site.
-    Update the neighboring sites in the process."""
-    if site == 0:
-        Ψ[site] = A
-        return site, 0.0
-    Ψ[site], Us, err = _ortho_left(A, truncation)
-    site -= 1
-    # np.einsum("aib,bc->aic", Ψ[site], Us)
-    Ψ[site] = np.matmul(Ψ[site], Us)
-    return site, err
-
-
-def _canonicalize(Ψ: list[Tensor3], center: int, truncation: Strategy) -> float:
-    """Update a list of `Tensor3` objects to be in canonical form
-    with respect to `center`."""
-    # TODO: Revise the cumulative error update. Does it follow update_error()?
-    err = 0.0
-    for i in range(0, center):
-        _, errk = _update_in_canonical_form_right(Ψ, Ψ[i], i, truncation)
-        err += sqrt(errk)
-    for i in range(len(Ψ) - 1, center, -1):
-        _, errk = _update_in_canonical_form_left(Ψ, Ψ[i], i, truncation)
-        err += sqrt(errk)
-    return err
-
-
 class CanonicalMPS(MPS):
     """Canonical MPS class.