From dee3752ee82f0f07a449d850dd42466468e99a3c Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Juan=20Jos=C3=A9=20Rodr=C3=ADguez=20Aldavero?= Date: Mon, 13 May 2024 12:20:54 +0200 Subject: [PATCH 1/7] Optimize Lagrange interpolation --- src/seemps/analysis/lagrange.py | 16 +++++++++------- src/seemps/analysis/mesh.py | 20 +++++++++++++------- 2 files changed, 22 insertions(+), 14 deletions(-) diff --git a/src/seemps/analysis/lagrange.py b/src/seemps/analysis/lagrange.py index b224da65..a3c865e6 100644 --- a/src/seemps/analysis/lagrange.py +++ b/src/seemps/analysis/lagrange.py @@ -5,9 +5,13 @@ from ..state import MPS, Strategy from ..state.schmidt import _destructive_svd +from ..state._contractions import _contract_last_and_first from ..state.core import destructively_truncate_vector from ..truncate import simplify, SIMPLIFICATION_STRATEGY -from .mesh import affine_transformation +from .mesh import array_affine + + +# TODO: Implement multivariate Lagrange interpolation and multirresolution constructions DEFAULT_LAGRANGE_STRATEGY = SIMPLIFICATION_STRATEGY.replace(normalize=False) @@ -99,7 +103,7 @@ def lagrange_rank_revealing( U_L, R = np.linalg.qr(Al.reshape((2, order + 1))) tensors = [U_L.reshape(1, 2, 2)] for _ in range(sites - 2): - B = np.tensordot(R, Ac, axes=1) + B = _contract_last_and_first(R, Ac) r1, s, r2 = B.shape ## SVD U, S, V = _destructive_svd(B.reshape(r1 * s, r2)) @@ -109,7 +113,7 @@ def lagrange_rank_revealing( R = S.reshape(D, 1) * V[:D, :] ## tensors.append(U.reshape(r1, s, -1)) - U_R = np.tensordot(R, Ar, axes=1) + U_R = _contract_last_and_first(R, Ar) tensors.append(U_R) return MPS(tensors) @@ -246,9 +250,7 @@ def local_chebyshev_cardinal(self, x: float, j: int) -> float: gamma_res = ( -gamma if gamma < 0 - else self.d - (gamma - self.d) - if gamma > self.d - else gamma + else self.d - (gamma - self.d) if gamma > self.d else gamma ) if j == gamma_res: P += self.local_angular_cardinal(theta, gamma) @@ -275,7 +277,7 @@ def A_L(self, func: Callable, start: float, stop: float) -> np.ndarray: """ def affine_func(u): - return func(affine_transformation(u, orig=(0, 1), dest=(start, stop))) + return func(array_affine(u, orig=(0, 1), dest=(start, stop))) A = np.zeros((1, 2, self.D)) for s in range(2): diff --git a/src/seemps/analysis/mesh.py b/src/seemps/analysis/mesh.py index 232d736e..fb124576 100644 --- a/src/seemps/analysis/mesh.py +++ b/src/seemps/analysis/mesh.py @@ -113,7 +113,7 @@ def __getitem__(self, idx: int) -> float: ... def __getitem__(self, idx: Union[int, np.ndarray]) -> Union[float, np.ndarray]: super()._validate_index(idx) zero = np.cos(np.pi * (2 * idx + 1) / (2 * self.size)) - return affine_transformation(zero, orig=(-1, 1), dest=(self.stop, self.start)) + return array_affine(zero, orig=(-1, 1), dest=(self.stop, self.start)) class ChebyshevExtremaInterval(Interval): @@ -132,7 +132,7 @@ def __getitem__(self, idx: int) -> float: ... def __getitem__(self, idx: Union[int, np.ndarray]) -> Union[float, np.ndarray]: super()._validate_index(idx) maxima = np.cos(np.pi * idx / (self.size - 1)) - return affine_transformation(maxima, orig=(-1, 1), dest=(self.stop, self.start)) + return array_affine(maxima, orig=(-1, 1), dest=(self.stop, self.start)) class Mesh: @@ -210,7 +210,11 @@ def to_tensor(self): ) -def affine_transformation(x: np.ndarray, orig: tuple, dest: tuple) -> np.ndarray: +def array_affine( + x: np.ndarray, + orig: tuple, + dest: tuple, +) -> np.ndarray: """ Performs an affine transformation of x as u = a*x + b from orig=(x0, x1) to dest=(u0, u1). """ @@ -220,13 +224,13 @@ def affine_transformation(x: np.ndarray, orig: tuple, dest: tuple) -> np.ndarray a = (u1 - u0) / (x1 - x0) b = 0.5 * ((u1 + u0) - a * (x0 + x1)) x_affine = a * x - if np.abs(b) > np.finfo(np.float64).eps: + if abs(b) > np.finfo(np.float64).eps: x_affine = x_affine + b return x_affine def mps_to_mesh_matrix( - sites_per_dimension: list[int], mps_order: str = "A" + sites_per_dimension: list[int], mps_order: str = "A", base: int = 2 ) -> np.ndarray: """ Returns a matrix that transforms an array of MPS indices @@ -236,13 +240,15 @@ def mps_to_mesh_matrix( T = np.zeros((sum(sites_per_dimension), len(sites_per_dimension)), dtype=int) start = 0 for m, n in enumerate(sites_per_dimension): - T[start : start + n, m] = 2 ** np.arange(n)[::-1] + T[start : start + n, m] = base ** np.arange(n)[::-1] start += n return T elif mps_order == "B": T = np.vstack( [ - np.diag([2 ** (n - i - 1) if n > i else 0 for n in sites_per_dimension]) + np.diag( + [base ** (n - i - 1) if n > i else 0 for n in sites_per_dimension] + ) for i in range(max(sites_per_dimension)) ] ) From 255ac0c2f6a4025f467e12a9371f5bd43e7450cf Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Juan=20Jos=C3=A9=20Rodr=C3=ADguez=20Aldavero?= Date: Mon, 13 May 2024 12:21:43 +0200 Subject: [PATCH 2/7] Add COMPUTER_PRECISION to factories --- src/seemps/analysis/factories.py | 47 ++++++++++++++++++++------------ 1 file changed, 29 insertions(+), 18 deletions(-) diff --git a/src/seemps/analysis/factories.py b/src/seemps/analysis/factories.py index e5e8a205..3d559321 100644 --- a/src/seemps/analysis/factories.py +++ b/src/seemps/analysis/factories.py @@ -6,24 +6,23 @@ MPS, MPSSum, Strategy, - DEFAULT_STRATEGY, Truncation, - DEFAULT_TOLERANCE, Simplification, ) -from ..truncate import simplify +from ..truncate import simplify, SIMPLIFICATION_STRATEGY from .mesh import ( Interval, RegularClosedInterval, RegularHalfOpenInterval, ChebyshevZerosInterval, + ChebyshevExtremaInterval, ) -DEFAULT_FACTORY_STRATEGY = Strategy( +COMPUTER_PRECISION = SIMPLIFICATION_STRATEGY.replace( + tolerance=np.finfo(np.double).eps, + simplification_tolerance=np.finfo(np.double).eps, + simplify=Simplification.DO_NOT_SIMPLIFY, method=Truncation.RELATIVE_SINGULAR_VALUE, - tolerance=DEFAULT_TOLERANCE, - simplify=Simplification.VARIATIONAL, - simplification_tolerance=DEFAULT_TOLERANCE, normalize=False, ) @@ -96,7 +95,7 @@ def mps_exponential(start: float, stop: float, sites: int, c: complex = 1) -> MP def mps_sin( - start: float, stop: float, sites: int, strategy: Strategy = DEFAULT_FACTORY_STRATEGY + start: float, stop: float, sites: int, strategy: Strategy = COMPUTER_PRECISION ) -> MPS: """ Returns an MPS representing a sine function discretized over an interval. @@ -124,7 +123,7 @@ def mps_sin( def mps_cos( - start: float, stop: float, sites: int, strategy: Strategy = DEFAULT_FACTORY_STRATEGY + start: float, stop: float, sites: int, strategy: Strategy = COMPUTER_PRECISION ) -> MPS: """ Returns an MPS representing a cosine function discretized over an interval. @@ -154,7 +153,7 @@ def mps_cos( _State = TypeVar("_State", bound=Union[MPS, MPSSum]) -def mps_affine_transformation(mps: _State, orig: tuple, dest: tuple) -> _State: +def mps_affine(mps: _State, orig: tuple, dest: tuple) -> _State: """ Applies an affine transformation to an MPS, mapping it from one interval [x0, x1] to another [u0, u1]. This is a transformation u = a * x + b, with u0 = a * x0 + b and and u1 = a * x1 + b. @@ -188,7 +187,7 @@ def mps_affine_transformation(mps: _State, orig: tuple, dest: tuple) -> _State: return mps_affine -def mps_interval(interval: Interval, strategy: Strategy = DEFAULT_FACTORY_STRATEGY): +def mps_interval(interval: Interval, strategy: Strategy = COMPUTER_PRECISION): """ Returns an MPS corresponding to a specific type of interval (open, closed, or Chebyshev zeros). @@ -216,7 +215,9 @@ def mps_interval(interval: Interval, strategy: Strategy = DEFAULT_FACTORY_STRATE start_zeros = np.pi / (2 ** (sites + 1)) stop_zeros = np.pi + start_zeros mps_zeros = -1.0 * mps_cos(start_zeros, stop_zeros, sites, strategy=strategy) - return mps_affine_transformation(mps_zeros, (-1, 1), (start, stop)) + return mps_affine(mps_zeros, (-1, 1), (start, stop)) + elif isinstance(interval, ChebyshevExtremaInterval): + raise NotImplementedError() else: raise ValueError(f"Unsupported interval type {type(interval)}") @@ -285,7 +286,7 @@ def extend_mps(mps_id: int, mps_map: list[tuple[int, Tensor3]]) -> MPS: def mps_tensor_product( mps_list: list[MPS], mps_order: str = "A", - strategy: Strategy = DEFAULT_FACTORY_STRATEGY, + strategy: Strategy = COMPUTER_PRECISION, ) -> MPS: """ Returns the tensor product of a list of MPS, with the sites arranged @@ -309,16 +310,20 @@ def mps_tensor_product( nested_sites = [mps._data for mps in mps_list] flattened_sites = [site for sites in nested_sites for site in sites] result = MPS(flattened_sites) - else: + elif mps_order == "B": terms = mps_tensor_terms(mps_list, mps_order) result = terms[0] - for idx, mps in enumerate(terms[1:]): + for _, mps in enumerate(terms[1:]): result = result * mps + else: + raise ValueError(f"Invalid mps order {mps_order}") return simplify(result, strategy=strategy) def mps_tensor_sum( - mps_list: list[MPS], mps_order: str = "A", strategy: Strategy = DEFAULT_STRATEGY + mps_list: list[MPS], + mps_order: str = "A", + strategy: Strategy = COMPUTER_PRECISION, ) -> MPS: """ Returns the tensor sum of a list of MPS, with the sites arranged @@ -340,16 +345,22 @@ def mps_tensor_sum( """ if mps_order == "A": result = _mps_tensor_sum_serial_order(mps_list) - else: + elif mps_order == "B": result = MPSSum( [1.0] * len(mps_list), mps_tensor_terms(mps_list, mps_order) ).join() + else: + raise ValueError(f"Invalid mps order {mps_order}") if strategy.get_simplify_flag(): return simplify(result, strategy=strategy) return result def _mps_tensor_sum_serial_order(mps_list: list[MPS]) -> MPS: + """ + Computes the MPS tensor sum in serial order in an optimized manner. + """ + def extend_tensor(A: Tensor3, first: bool, last: bool) -> Tensor3: a, d, b = A.shape output = np.zeros((a + 2, d, b + 2), dtype=A.dtype) @@ -368,7 +379,7 @@ def extend_tensor(A: Tensor3, first: bool, last: bool) -> Tensor3: output = [ extend_tensor(Ai, i == 0, i == len(A) - 1) - for n, A in enumerate(mps_list) + for _, A in enumerate(mps_list) for i, Ai in enumerate(A) ] output[0] = output[0][[0], :, :] From 9aa8da629ab86a62749b4e5f038aa2c5b4beaa15 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Juan=20Jos=C3=A9=20Rodr=C3=ADguez=20Aldavero?= Date: Mon, 13 May 2024 12:22:37 +0200 Subject: [PATCH 3/7] Add mps_as_mpo to simplify_mpo --- src/seemps/truncate/simplify_mpo.py | 47 +++++++++++++++++++---------- 1 file changed, 31 insertions(+), 16 deletions(-) diff --git a/src/seemps/truncate/simplify_mpo.py b/src/seemps/truncate/simplify_mpo.py index 3bb60ae9..d7a18a06 100644 --- a/src/seemps/truncate/simplify_mpo.py +++ b/src/seemps/truncate/simplify_mpo.py @@ -1,18 +1,36 @@ -from typing import Optional +from typing import Optional, Union +from math import isqrt -from ..operators import MPO +from ..operators import MPO, MPOList, MPOSum from ..state import DEFAULT_STRATEGY, MPS, Strategy from ..truncate import SIMPLIFICATION_STRATEGY, simplify -def mpo_as_mps(mpo): +def mpo_as_mps(mpo: MPO) -> MPS: """Recast MPO as MPS.""" _, i, j, _ = mpo[0].shape return MPS([t.reshape(t.shape[0], i * j, t.shape[-1]) for t in mpo._data]) +def mps_as_mpo( + mps: MPS, + mpo_strategy: Strategy = DEFAULT_STRATEGY, +) -> MPO: + """Recast MPS as MPO.""" + _, S, _ = mps[0].shape + s = isqrt(S) + if s**2 != S: + raise ValueError("The physical dimensions of the MPS must be a perfect square") + return MPO( + [t.reshape(t.shape[0], s, s, t.shape[-1]) for t in mps._data], + strategy=mpo_strategy, + ) + + +# TODO: As opposed to MPS, the MPO class does not have an error attribute to keep track +# of the simplification errors def simplify_mpo( - operator: MPO, + operator: Union[MPO, MPOList, MPOSum], strategy: Strategy = SIMPLIFICATION_STRATEGY, direction: int = +1, guess: Optional[MPS] = None, @@ -23,26 +41,23 @@ def simplify_mpo( Parameters ---------- - operator : MPO - MPO to approximate. + operator : Union[MPO, MPOList, MPOSum] + Operator to approximate. strategy : Strategy - Truncation strategy. Defaults to `SIMPLIFICATION_STRATEGY`. + Truncation strategy. Defaults to `SIMPLIFICATION_STRATEGY`. direction : { +1, -1 } - Initial direction for the sweeping algorithm. Defaults to +1. + Initial direction for the sweeping algorithm. Defaults to +1. guess : MPS - A guess for the new state, to ease the optimization. Defaults to None. + A guess for the new state, to ease the optimization. Defaults to None. mpo_strategy : Strategy - Strategy of the resulting MPO. Defaults to `DEFAULT_STRATEGY`. + Strategy of the resulting MPO. Defaults to `DEFAULT_STRATEGY`. Returns ------- MPO Approximation O to the operator. """ - _, i, j, _ = operator[0].shape + if isinstance(operator, MPOList) or isinstance(operator, MPOSum): + operator = operator.join() mps = simplify(mpo_as_mps(operator), strategy, direction, guess) - [t.reshape(t.shape[0], i, j, t.shape[-1]) for t in mps._data] - return MPO( - [t.reshape(t.shape[0], i, j, t.shape[-1]) for t in mps._data], - strategy=mpo_strategy, - ) + return mps_as_mpo(mps, mpo_strategy=mpo_strategy) From cd82f0bffc1ab6880c063c9057e4df883a71f881 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Juan=20Jos=C3=A9=20Rodr=C3=ADguez=20Aldavero?= Date: Mon, 13 May 2024 12:23:04 +0200 Subject: [PATCH 4/7] Add mpo_affine transformation to operators --- src/seemps/analysis/operators.py | 18 +++++++++++++++++- 1 file changed, 17 insertions(+), 1 deletion(-) diff --git a/src/seemps/analysis/operators.py b/src/seemps/analysis/operators.py index 99bf5884..a87d9285 100644 --- a/src/seemps/analysis/operators.py +++ b/src/seemps/analysis/operators.py @@ -1,6 +1,6 @@ from __future__ import annotations import numpy as np -from ..operators import MPO, MPOList +from ..operators import MPO, MPOList, MPOSum from ..state import Strategy, DEFAULT_STRATEGY from typing import Union @@ -292,3 +292,19 @@ def sin_mpo(n: int, a: float, dx: float, strategy=DEFAULT_STRATEGY): exp2 = exponential_mpo(n, a, dx, c=-1j, strategy=strategy) sin_mpo = (-1j) * 0.5 * (exp1 - exp2) return sin_mpo.join(strategy=strategy) + + +def mpo_affine( + mpo: MPO, + orig: tuple, + dest: tuple, +): + x0, x1 = orig + u0, u1 = dest + a = (u1 - u0) / (x1 - x0) + b = 0.5 * ((u1 + u0) - a * (x0 + x1)) + mpo_affine = a * mpo + if abs(b) > np.finfo(np.float64).eps: + I = MPO([np.ones((1, 2, 2, 1))] * len(mpo_affine)) + mpo_affine = MPOSum(mpos=[mpo_affine, I], weights=[1, b]).join() + return mpo_affine From 36f342b670ecbb824f95015e2d69ed7b1cb49c97 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Juan=20Jos=C3=A9=20Rodr=C3=ADguez=20Aldavero?= Date: Mon, 13 May 2024 12:23:24 +0200 Subject: [PATCH 5/7] Update MPS Chebyshev approximation --- src/seemps/analysis/chebyshev.py | 425 ++++++++++++++++++-------- tests/test_analysis/test_chebyshev.py | 185 +++++++++-- tests/test_analysis/test_mesh.py | 12 + 3 files changed, 454 insertions(+), 168 deletions(-) diff --git a/src/seemps/analysis/chebyshev.py b/src/seemps/analysis/chebyshev.py index 4e22bdd2..f4cb3c1d 100644 --- a/src/seemps/analysis/chebyshev.py +++ b/src/seemps/analysis/chebyshev.py @@ -3,211 +3,366 @@ from math import sqrt import numpy as np from scipy.fft import dct # type: ignore -from ..tools import make_logger -from ..state import CanonicalMPS, MPS, MPSSum, Strategy, Truncation, Simplification -from ..truncate import simplify -from ..operators import MPO -from .mesh import ChebyshevZerosInterval, Interval -from .factories import mps_interval, mps_affine_transformation - -# TODO: All the tests have been done using the `RELATIVE_SINGULAR_VALUE` -# truncation method, which is kind of flunky, because it does not measure -# the actual error. This should be migrated to DEFAULT_STRATEGY, maybe -# strengthening the tolerance -DEFAULT_CHEBYSHEV_STRATEGY = Strategy( - method=Truncation.RELATIVE_SINGULAR_VALUE, - tolerance=1e-8, - simplify=Simplification.VARIATIONAL, - simplification_tolerance=1e-8, - normalize=False, + +from .. import tools +from ..state import CanonicalMPS, MPS, MPSSum, Strategy +from ..truncate import simplify, SIMPLIFICATION_STRATEGY +from ..truncate.simplify_mpo import simplify_mpo +from ..operators import MPO, MPOList, MPOSum +from .mesh import ( + Interval, + ChebyshevZerosInterval, + ChebyshevExtremaInterval, + array_affine, ) +from .operators import mpo_affine +from .factories import mps_interval, mps_affine + + +DEFAULT_CHEBYSHEV_STRATEGY = SIMPLIFICATION_STRATEGY.replace(normalize=False) -# TODO: Implement projection coefficients (coming from integration) -def chebyshev_coefficients( - f: Callable, - order: int, +def interpolation_coefficients( + func: Callable, + order: Optional[int] = None, start: float = -1, stop: float = +1, domain: Optional[Interval] = None, + interpolated_nodes: str = "zeros", ) -> np.polynomial.Chebyshev: """ - Returns the Chebyshev interpolation coefficients for a given function - on a specified interval. - - The error of a Chebyshev interpolation is proportiona to all of the neglected - coefficients of higher order. If the decay is exponential, it can be approximated by - the last coefficient computed by this function. + Returns the coefficients for the Chebyshev interpolation of a function on a given set + of nodes and on a specified interval. Parameters ---------- - f : Callable + func : Callable The target function to approximate with Chebyshev polynomials. - order : int + order : Optional[int], default = None The number of Chebyshev coefficients to compute. + If None, estimates an order that results in an error below machine precision. domain : Optional[Interval], default = None The domain on which the function is defined and in which the approximation is desired. start : float, default = -1 stop : float, default = +1 Alternative way to specify the function's domain. + interpolated_nodes : str, default = "zeros" + The nodes on which the function is interpolated. Use "zeros" for + Chebyshev zeros or "extrema" for Chebyshev extrema. Returns ------- coefficients : `numpy.polynomial.Chebyshev` An array of Chebyshev coefficients scaled to the specified interval. """ + if order is None: + order = estimate_order(func, start, stop, domain) if domain is not None: start, stop = domain.start, domain.stop - chebyshev_zeros = np.flip(ChebyshevZerosInterval(start, stop, order).to_vector()) - coefficients = dct(f(chebyshev_zeros), type=2) / order + if interpolated_nodes == "zeros": + nodes = ChebyshevZerosInterval(start, stop, order).to_vector() + coefficients = (1 / order) * dct(np.flip(func(nodes)), type=2) # type: ignore + elif interpolated_nodes == "extrema": + nodes = ChebyshevExtremaInterval(start, stop, order).to_vector() + coefficients = 2 * dct(np.flip(func(nodes)), type=1, norm="forward") + coefficients[0] /= 2 # type: ignore + return np.polynomial.Chebyshev(coefficients, domain=(start, stop)) + + +def projection_coefficients( + func: Callable, + order: Optional[int] = None, + start: float = -1, + stop: float = +1, + domain: Optional[Interval] = None, +) -> np.polynomial.Chebyshev: + """ + Returns the coefficients for the Chebyshev projection of a function using + Chebyshev-Gauss integration. + + Parameters + ---------- + func : Callable + The target function to approximate with Chebyshev polynomials. + order : Optional[int], default = None + The number of Chebyshev projection coefficients to compute. + If None, estimates an order that results in an error below machine precision. + start : float, default = -1 + stop : float, default = +1 + The domain on which the function is defined and in which the + approximation is desired. + domain : Optional[Interval], default = None + Alternative way to specify the function's domain. + + Returns + ------- + coefficients : `numpy.polynomial.Chebyshev` + An array of Chebyshev coefficients scaled to the specified interval. + """ + if order is None: + order = estimate_order(func, start, stop, domain) + if domain is not None: + start, stop = domain.start, domain.stop + quad_order = order # TODO: Check if this order integrates to machine precision + nodes = np.cos(np.pi * np.arange(1, 2 * quad_order, 2) / (2.0 * quad_order)) + nodes_affine = array_affine(nodes, orig=(-1, 1), dest=(start, stop)) + weights = np.ones(quad_order) * (np.pi / quad_order) + T_matrix = np.cos(np.outer(np.arange(order), np.arccos(nodes))) + coefficients = (2 / np.pi) * (T_matrix * func(nodes_affine)) @ weights coefficients[0] /= 2 return np.polynomial.Chebyshev(coefficients, domain=(start, stop)) -# TODO: Implement adaptivity (starting point) for when using projection coefficients +def estimate_order( + func: Callable, + start: float = -1, + stop: float = +1, + domain: Optional[Interval] = None, + tolerance: float = float(np.finfo(np.float64).eps), + initial_order: int = 2, + max_order: int = 2**12, # 4096 +) -> int: + """ + Returns an estimation of the number of Chebyshev coefficients required to achieve a + given accuracy such that the last pair of coefficients fall below a given tolerance, + as they theoretically bound the maximum error of the expansion. + + Notes + ----- + - The coefficients are evaluated in pairs because even and odd functions respectively + have vanishing even and odd coefficients. + """ + if domain is not None: + start, stop = domain.start, domain.stop + order = initial_order + while order <= max_order: + c = projection_coefficients(func, order, start, stop).coef + max_c_in_pairs = np.maximum(abs(c[::2]), abs(c[1::2])) + c_below_tolerance = np.where(max_c_in_pairs < tolerance)[0] + if c_below_tolerance.size > 0 and c_below_tolerance[0] != 0: + return 2 * c_below_tolerance[0] + 1 + order *= 2 + raise ValueError("Order exceeds max_order without achieving tolerance.") + + def cheb2mps( - c: np.polynomial.Chebyshev, + coefficients: np.polynomial.Chebyshev, + initial_mps: Optional[MPS] = None, domain: Optional[Interval] = None, - x: Optional[MPS] = None, strategy: Strategy = DEFAULT_CHEBYSHEV_STRATEGY, + clenshaw: bool = True, + rescale: bool = True, ) -> MPS: """ - Construct an MPS representation of a function, from a Chebyshev expansion. - - This function takes as input an MPS representation of the first order - polynomial `x` in a given `domain`, with values `[x0, x1]`. It also takes - a Chebyshev expansion `c` of a function `c(x)` defined in a domain that - contains this interval `[x0, x1]`. With this information, it constructs - the MPS that approximates `c(x)`. + Constructs a MPS representation of a function by expanding it on the basis + of Chebyshev polynomials. It takes as input the Chebyshev coefficients of a function + `f(x)` defined in an interval `[a, b]` and, optionally, an initial MPS representing a + function `g(x)` that is taken as the first order polynomial of the expansion. + With this information, it constructs the MPS that approximates `f(g(x))`. Parameters ---------- - c : `numpy.polynomial.Chebyshev` - Chebyshev expansion over a given domain. + coefficients : np.polynomial.Chebyshev + The Chebyshev expansion coefficients representing the target function that + is defined on a given interval `[a, b]`. + initial_mps: Optional[MPS], default = None + The initial MPS on which to apply the expansion. + By default (if `rescale` is True), it must have a support inside the domain of + definition of the function `[a, b]`. + If `rescale` is False, it must have a support inside `[-1, 1]`. domain : Optional[Interval], default = None - Interval of definition for the function, which must be contained in the - Chebyshev's series domain. - x : Optional[MPS], default = None - MPS representation of the `x` function in the series' domain. It will - be computed from `domain` if not provided. - strategy : Strategy, default = DEFAULT_CHEBYSHEV_STRATEGY - Simplification strategy for operations between MPS. + An alternative way to specify the initial MPS by constructing it from the given Interval. + strategy : Strategy + The simplification strategy for operations between MPS. + clenshaw : bool, default = True + Whether to use the Clenshaw algorithm for polynomial evaluation. + rescale : bool, default = True + Whether to perform an affine transformation of the initial MPS from the domain + `[a, b]` of the Chebyshev coefficients to the canonical Chebyshev interval `[-1, 1]`. Returns ------- - f : MPS + f_mps : MPS MPS representation of the polynomial expansion. + + Notes + ----- + - The complexity of the expansion depends on the complexity (e.g. maximum bond dimension) + of the intermediate states which are general hard to anticipate and may depend on + the type of coefficients (interpolation, projection) or evaluation method used (Clenshaw, polynomial). + + - The Clenshaw evaluation method has a better performance overall, but performs worse when + the interpolation order is overestimated, meaning that the expected accuracy is below + machine precision. """ - x_mps: MPS - if domain is not None: - x_mps = mps_interval(domain.map_to(-1, 1)) - elif isinstance(x, MPS): - orig, _ = c.linspace(2) - x_mps = mps_affine_transformation(x, orig, (-1, 1)) + if isinstance(initial_mps, MPS): + pass + elif isinstance(domain, Interval): + initial_mps = mps_interval(domain) else: - raise Exception("In cheb2mps, either domain or an MPS must be provided.") + raise ValueError("Either a domain or an initial MPS must be provided.") + if rescale: + orig = tuple(coefficients.linspace(2)[0]) + initial_mps = mps_affine(initial_mps, orig, (-1, 1)) - I_norm = 2 ** (x_mps.size / 2) + c = coefficients.coef + I_norm = 2 ** (initial_mps.size / 2) normalized_I = CanonicalMPS( - [np.ones((1, 2, 1)) / sqrt(2.0)] * x_mps.size, center=0, is_canonical=True + [np.ones((1, 2, 1)) / sqrt(2.0)] * initial_mps.size, + center=0, + is_canonical=True, + ) + x_norm = initial_mps.norm() + normalized_x = CanonicalMPS( + initial_mps, center=0, normalize=True, strategy=strategy ) - y_i = y_i_plus_1 = normalized_I.zero_state() - with make_logger(2) as logger: - logger(f"Clenshaw evaluation started with {len(c)} steps") - for i, c_i in enumerate(reversed(c.coef)): + if clenshaw: + if tools.DEBUG: + steps = len(c) + tools.log("MPS Clenshaw evaluation started") + y_i = y_i_plus_1 = normalized_I.zero_state() + for i, c_i in enumerate(reversed(c)): y_i_plus_1, y_i_plus_2 = y_i, y_i_plus_1 y_i = simplify( # coef[i] * I - y[i + 2] + (2 * x_mps) * y[i + 1], MPSSum( - [I_norm * c_i, -1, 2], - [normalized_I, y_i_plus_2, x_mps * y_i_plus_1], + weights=[c_i * I_norm, -1, 2 * x_norm], + states=[normalized_I, y_i_plus_2, normalized_x * y_i_plus_1], + check_args=False, + ), + strategy=strategy, + ) + if tools.DEBUG: + tools.log( + f"MPS Clenshaw step {i+1}/{steps} with maximum bond dimension {max(y_i.bond_dimensions())} and error {y_i.error():6e}" + ) + f_mps = simplify( + MPSSum( + weights=[1, -x_norm], + states=[y_i, normalized_x * y_i_plus_1], + check_args=False, + ), + strategy=strategy, + ) + else: + if tools.DEBUG: + steps = len(c) + tools.log("MPS Chebyshev expansion started") + + f_mps = simplify( + MPSSum( + weights=[c[0] * I_norm, c[1] * x_norm], + states=[normalized_I, normalized_x], + check_args=False, + ), + strategy=strategy, + ) + T_i, T_i_plus_1 = I_norm * normalized_I, x_norm * normalized_x + for i, c_i in enumerate(c[2:], start=2): + T_i_plus_2 = simplify( + MPSSum( + weights=[2 * x_norm, -1], + states=[normalized_x * T_i_plus_1, T_i], check_args=False, ), strategy=strategy, ) - if logger: - logger( - f"Clenshaw step {i} with maximum bond dimension {max(y_i.bond_dimensions())} and error {y_i.error():6e}" + f_mps = simplify( + MPSSum(weights=[1, c_i], states=[f_mps, T_i_plus_2], check_args=False), + strategy=strategy, + ) + if tools.DEBUG: + tools.log( + f"MPS expansion step {i+1}/{steps} with maximum bond dimension {max(f_mps.bond_dimensions())} and error {f_mps.error():6e}" ) - return simplify(y_i - x_mps * y_i_plus_1, strategy=strategy) + T_i, T_i_plus_1 = T_i_plus_1, T_i_plus_2 + return f_mps -# TODO: Implement def cheb2mpo( - c: np.polynomial.Chebyshev, - domain: Optional[Interval] = None, - x: Optional[MPO] = None, + coefficients: np.polynomial.Chebyshev, + initial_mpo: MPO, strategy: Strategy = DEFAULT_CHEBYSHEV_STRATEGY, + clenshaw=True, + rescale=True, ) -> MPO: """ - *NOT IMPLEMENTED* - Construct an MPO representation of a function, from a Chebyshev expansion. - - This function takes as input an MPO representation of the first order - polynomial `x` in a given `domain`, with values `[x0,x1]`. It also takes - a Chebyshev expansion `c` of a function `c(x)` defined in a domain that - contains this interval `[x0,x1]`. With this information, it constructs - the MPO that approximates `c(x)`. + Constructs a MPO representation of a function by expanding it on the basis + of Chebyshev polynomials. Parameters ---------- - c : `numpy.polynomial.Chebyshev` - Chebyshev expansion over a given domain. - domain : Optional[Interval], default = None - Interval of definition for the function, whose domain must be included - in that of `c`. - x : Optional[MPO], default = None - MPS representation of the `x` function in the series' domain. It will - be computed from `domain` if not provided. - strategy : Strategy, default = DEFAULT_CHEBYSHEV_STRATEGY - Simplification strategy for operations between MPOs. + coefficients : np.polynomial.Chebyshev + The Chebyshev expansion coefficients representing the target function that + is defined on a given interval `[a, b]`. + initial_mpo: MPO + The initial MPO on which to apply the expansion. + By default (if `rescale` is True), it must have a support inside the domain of + definition of the function `[a, b]`. + If `rescale` is False, it must have a support inside `[-1, 1]`. + strategy : Strategy + The simplification strategy for operations between MPS. + clenshaw : bool, default = True + Whether to use the Clenshaw algorithm for polynomial evaluation. + rescale : bool, default = True + Whether to perform an affine transformation of the initial MPO from the domain + `[a, b]` of the Chebyshev coefficients to the canonical Chebyshev interval `[-1, 1]`. Returns ------- - f : MPO + f_mpo : MPO MPO representation of the polynomial expansion. """ - raise Exception("cheb2mpo not implemented") - - -# TODO: Consider if this helper function is necessary -def chebyshev_approximation( - f: Callable, - order: int, - domain: Interval, - differentiation_order: int = 0, - strategy: Strategy = DEFAULT_CHEBYSHEV_STRATEGY, -) -> MPS: - """ - Load a function as an MPS using Chebyshev expansions. - - This function constructs a Chebyshev series that approximates `f` over - the given `domain`, and uses that expansion to construct an MPS - representation via `cheb2mps`. - - Parameters - ---------- - func : Callable - A univariate scalar function. - order : int - Order of the Chebyshev expansion - domain : Interval - The domain over which the function is to be approximated. - differentiation_order : int, default = 0 - If positive or negative value `N`, integrate or differentiate the - function a total of `abs(N)` times prior to computing the expansion. - strategy : Strategy, default Strategy() - The strategy used for simplifying the MPS or MPO during computation. + if rescale: + orig = tuple(coefficients.linspace(2)[0]) + initial_mpo = mpo_affine(initial_mpo, orig, (-1, 1)) - Returns - ------- - mps : MPS - MPS approximation to the function. - """ - c = chebyshev_coefficients(f, order, domain.start, domain.stop) - if differentiation_order < 0: - c = c.integ(-differentiation_order, lbnd=domain.start) - elif differentiation_order > 0: - c = c.deriv(differentiation_order) - return cheb2mps(c, domain, strategy=strategy) + c = coefficients.coef + I = MPO([np.eye(2).reshape(1, 2, 2, 1)] * len(initial_mpo)) + if clenshaw: + if tools.DEBUG: + steps = len(c) + tools.log("MPO Clenshaw evaluation started") + y_i = y_i_plus_1 = MPO([np.zeros((1, 2, 2, 1))] * len(initial_mpo)) + for i, c_i in enumerate(reversed(coefficients.coef)): + y_i_plus_1, y_i_plus_2 = y_i, y_i_plus_1 + y_i = simplify_mpo( + MPOSum( + mpos=[I, y_i_plus_2, MPOList([initial_mpo, y_i_plus_1])], + weights=[c_i, -1, 2], + ), + strategy=strategy, + ) + if tools.DEBUG: + tools.log( + f"MPO Clenshaw step {i+1}/{steps} with maximum bond dimension {max(y_i.bond_dimensions())}" + ) + f_mpo = simplify_mpo( + MPOSum([y_i, MPOList([initial_mpo, y_i_plus_1])], weights=[1, -1]), + strategy=strategy, + ) + else: + if tools.DEBUG: + steps = len(c) + tools.log("MPO Chebyshev expansion started") + T_i, T_i_plus_1 = I, initial_mpo + f_mpo = simplify_mpo( + MPOSum(mpos=[T_i, T_i_plus_1], weights=[c[0], c[1]]), + strategy=strategy, + ) + for i, c_i in enumerate(c[2:], start=2): + T_i_plus_2 = simplify_mpo( + MPOSum(mpos=[MPOList([initial_mpo, T_i_plus_1]), T_i], weights=[2, -1]), + strategy=strategy, + ) + f_mpo = simplify_mpo( + MPOSum(mpos=[f_mpo, T_i_plus_2], weights=[1, c_i]), + strategy=strategy, + ) + if tools.DEBUG: + tools.log( + f"MPO expansion step {i+1}/{steps} with maximum bond dimension {max(f_mpo.bond_dimensions())}" + ) + T_i, T_i_plus_1 = T_i_plus_1, T_i_plus_2 + return f_mpo diff --git a/tests/test_analysis/test_chebyshev.py b/tests/test_analysis/test_chebyshev.py index f0789044..63d09b07 100644 --- a/tests/test_analysis/test_chebyshev.py +++ b/tests/test_analysis/test_chebyshev.py @@ -1,23 +1,27 @@ import numpy as np from numpy.polynomial import Chebyshev from scipy.special import erf -from math import sqrt + +from seemps.state import MPS, NO_TRUNCATION from seemps.analysis.mesh import RegularHalfOpenInterval from seemps.analysis.factories import mps_tensor_sum, mps_interval from seemps.analysis.chebyshev import ( DEFAULT_CHEBYSHEV_STRATEGY, - chebyshev_coefficients, - chebyshev_approximation, + interpolation_coefficients, + projection_coefficients, cheb2mps, + estimate_order, + cheb2mpo, ) +from seemps.analysis.operators import x_mpo from ..tools import TestCase -class TestChebyshevExpansion(TestCase): - def test_chebyshev_coefficients_exponential(self): +class TestChebyshevCoefficients(TestCase): + def test_interpolation_coefficients_exponential(self): f = lambda x: np.exp(x) - cheb_coeffs = chebyshev_coefficients(f, 15, -1, 1) + cheb_coeffs = interpolation_coefficients(f, 15, -1, 1) correct_coeffs = [ 1.266065877752008, 1.130318207984970, @@ -37,6 +41,21 @@ def test_chebyshev_coefficients_exponential(self): ] assert np.allclose(list(cheb_coeffs), correct_coeffs) + def test_estimate_order(self): + """Assert that the estimated coefficients and accuracy in norm-inf are below a tolerance.""" + tolerance = 1e-12 + f = lambda x: np.exp(-(x**2)) + proj_coeffs = projection_coefficients( + f, order=estimate_order(f, tolerance=tolerance) + ) + n = 6 + domain = RegularHalfOpenInterval(-1, 1, 2**n) + mps = cheb2mps(proj_coeffs, domain=domain, strategy=NO_TRUNCATION) + y_vec = f(domain.to_vector()) + y_mps = mps.to_vector() + assert proj_coeffs.coef[-1] <= tolerance + assert max(abs(y_mps - y_vec)) <= tolerance + def assertSimilarSeries(self, s1, s2, tol=1e-15): """Ensure two Chebyshev series are close up to tolerance.""" if s1.has_sametype(s2) and s1.has_samedomain(s2): @@ -53,45 +72,98 @@ def assertSimilarSeries(self, s1, s2, tol=1e-15): raise self.failureException(f"Not similar series:\nA={s1}\nB={s2}") def test_chebyshev_coefficients_T0(self): + # At low orders, the three types of coefficients (zeros, extrema and projections) are almost similar T0 = Chebyshev([1, 0, 0, 0]) - self.assertSimilarSeries(chebyshev_coefficients(T0, 4, -1, 1), T0) + self.assertSimilarSeries(interpolation_coefficients(T0, 4, -1, 1), T0) + self.assertSimilarSeries( + interpolation_coefficients(T0, 4, -1, 1, interpolated_nodes="extrema"), + T0, + ) + self.assertSimilarSeries(projection_coefficients(T0, 4, -1, 1), T0) def test_chebyshev_coefficients_T1(self): T1 = Chebyshev([0, 1, 0, 0]) - self.assertSimilarSeries(chebyshev_coefficients(T1, 4, -1, 1), T1) + self.assertSimilarSeries(interpolation_coefficients(T1, 4, -1, 1), T1) + self.assertSimilarSeries( + interpolation_coefficients(T1, 4, -1, 1, interpolated_nodes="extrema"), + T1, + ) + self.assertSimilarSeries(projection_coefficients(T1, 4, -1, 1), T1) def test_chebyshev_coefficients_T4(self): T4 = Chebyshev([0, 0, 0, 0, 1]) - self.assertSimilarSeries(chebyshev_coefficients(T4, 5, -1, 1), T4) + self.assertSimilarSeries(interpolation_coefficients(T4, 5, -1, 1), T4) + self.assertSimilarSeries( + interpolation_coefficients(T4, 5 + 1, -1, 1, interpolated_nodes="extrema"), + T4, # The extrema are computed for a polynomial of degree d-1 so we need d+1 + ) + self.assertSimilarSeries(projection_coefficients(T4, 5, -1, 1), T4) def test_chebyshev_coefficients_T0_other_domain(self): T0 = Chebyshev([1, 0, 0, 0], domain=(-2, 4)) - self.assertSimilarSeries(chebyshev_coefficients(T0, 4, -2, 4), T0) + self.assertSimilarSeries(interpolation_coefficients(T0, 4, -2, 4), T0) + self.assertSimilarSeries( + interpolation_coefficients(T0, 4, -2, 4, interpolated_nodes="extrema"), + T0, + ) + self.assertSimilarSeries(projection_coefficients(T0, 4, -2, 4), T0) def test_chebyshev_coefficients_T1_other_domain(self): T1 = Chebyshev([0, 1, 0, 0], domain=(-2, 4)) - self.assertSimilarSeries(chebyshev_coefficients(T1, 4, -2, 4), T1) + self.assertSimilarSeries(interpolation_coefficients(T1, 4, -2, 4), T1) + self.assertSimilarSeries( + interpolation_coefficients(T1, 4, -2, 4, interpolated_nodes="extrema"), + T1, + ) + self.assertSimilarSeries(projection_coefficients(T1, 4, -2, 4), T1) def test_chebyshev_coefficients_T4_other_domain(self): T4 = Chebyshev([0, 0, 0, 0, 1], domain=(-2, 4)) - self.assertSimilarSeries(chebyshev_coefficients(T4, 5, -2, 4), T4) + self.assertSimilarSeries(interpolation_coefficients(T4, 5, -2, 4), T4) + self.assertSimilarSeries( + interpolation_coefficients(T4, 5 + 1, -2, 4, interpolated_nodes="extrema"), + T4, + ) + self.assertSimilarSeries(projection_coefficients(T4, 5, -2, 4), T4) def test_chebyshev_coefficients_gaussian_derivative(self): f = lambda x: np.exp(-x * x) df = lambda x: -2 * x * np.exp(-x * x) self.assertSimilarSeries( - chebyshev_coefficients(f, 22, -1, 2).deriv(), - chebyshev_coefficients(df, 22, -1, 2), + interpolation_coefficients(f, 22, -1, 2).deriv(), + interpolation_coefficients(df, 22, -1, 2), + ) + self.assertSimilarSeries( + interpolation_coefficients( + f, 22, -1, 2, interpolated_nodes="extrema" + ).deriv(), + interpolation_coefficients(df, 22, -1, 2, interpolated_nodes="extrema"), + ) + self.assertSimilarSeries( + projection_coefficients(f, 22, -1, 2).deriv(), + projection_coefficients(df, 22, -1, 2), ) def test_chebyshev_coefficients_gaussian_integral(self): start = -1 stop = 2 f = lambda x: np.exp(-x * x) - f_intg = lambda x: (sqrt(np.pi) / 2) * (erf(x) - erf(start)) + f_intg = lambda x: (np.sqrt(np.pi) / 2) * (erf(x) - erf(start)) + self.assertSimilarSeries( + interpolation_coefficients(f, 22, start, stop).integ(1, lbnd=start), + interpolation_coefficients(f_intg, 22, start, stop), + ) + self.assertSimilarSeries( + interpolation_coefficients( + f, 22, start, stop, interpolated_nodes="extrema" + ).integ(1, lbnd=start), + interpolation_coefficients( + f_intg, 22, start, stop, interpolated_nodes="extrema" + ), + ) self.assertSimilarSeries( - chebyshev_coefficients(f, 22, start, stop).integ(1, lbnd=start), - chebyshev_coefficients(f_intg, 22, start, stop), + projection_coefficients(f, 22, start, stop).integ(1, lbnd=start), + projection_coefficients(f_intg, 22, start, stop), ) @@ -99,47 +171,94 @@ class TestChebyshevMPS(TestCase): def test_gaussian_1d(self): f = lambda x: np.exp(-(x**2)) interval = RegularHalfOpenInterval(-1, 2, 2**5) - mps_cheb = chebyshev_approximation(f, 30, interval) - self.assertSimilar(f(interval.to_vector()), mps_cheb.to_vector()) + mps_cheb_clen = cheb2mps( + interpolation_coefficients(f, 30, domain=interval), + domain=interval, + clenshaw=True, + ) + mps_cheb_poly = cheb2mps( + interpolation_coefficients(f, 30, domain=interval), + domain=interval, + clenshaw=False, + ) + self.assertSimilar(f(interval.to_vector()), mps_cheb_clen.to_vector()) + self.assertSimilar(f(interval.to_vector()), mps_cheb_poly.to_vector()) def test_gaussian_derivative_1d(self): f = lambda x: np.exp(-(x**2)) f_diff = lambda x: -2 * x * np.exp(-(x**2)) interval = RegularHalfOpenInterval(-1, 2, 2**5) - mps_cheb = chebyshev_approximation(f, 30, interval, differentiation_order=1) - self.assertSimilar(f_diff(interval.to_vector()), mps_cheb.to_vector()) + c = interpolation_coefficients(f, 30, domain=interval) + mps_cheb_clen = cheb2mps(c.deriv(1), domain=interval, clenshaw=True) + mps_cheb_poly = cheb2mps(c.deriv(1), domain=interval, clenshaw=False) + self.assertSimilar(f_diff(interval.to_vector()), mps_cheb_clen.to_vector()) + self.assertSimilar(f_diff(interval.to_vector()), mps_cheb_poly.to_vector()) def test_gaussian_integral_1d(self): - f_intg = lambda x: (sqrt(np.pi) / 2) * (erf(x) - erf(-1)) + f_intg = lambda x: (np.sqrt(np.pi) / 2) * (erf(x) - erf(-1)) interval = RegularHalfOpenInterval(-1, 2, 2**5) - mps_cheb = chebyshev_approximation(f_intg, 30, interval) - self.assertSimilar(f_intg(interval.to_vector()), mps_cheb.to_vector()) + c = interpolation_coefficients(f_intg, 30, domain=interval) + mps_cheb_clen = cheb2mps(c, domain=interval, clenshaw=True) + mps_cheb_poly = cheb2mps(c, domain=interval, clenshaw=False) + self.assertSimilar(f_intg(interval.to_vector()), mps_cheb_clen.to_vector()) + self.assertSimilar(f_intg(interval.to_vector()), mps_cheb_poly.to_vector()) def test_gaussian_integral_1d_b(self): f = lambda x: np.exp(-(x**2)) - f_intg = lambda x: (sqrt(np.pi) / 2) * (erf(x) - erf(-1)) + f_intg = lambda x: (np.sqrt(np.pi) / 2) * (erf(x) - erf(-1)) interval = RegularHalfOpenInterval(-1, 2, 2**5) - mps_cheb = chebyshev_approximation(f, 30, interval, differentiation_order=-1) - self.assertSimilar(f_intg(interval.to_vector()), mps_cheb.to_vector()) + c = interpolation_coefficients(f, 30, domain=interval) + mps_cheb_clen = cheb2mps( + c.integ(1, lbnd=interval.start), domain=interval, clenshaw=True + ) + mps_cheb_poly = cheb2mps( + c.integ(1, lbnd=interval.start), domain=interval, clenshaw=False + ) + self.assertSimilar(f_intg(interval.to_vector()), mps_cheb_clen.to_vector()) + self.assertSimilar(f_intg(interval.to_vector()), mps_cheb_poly.to_vector()) def test_gaussian_2d(self): f = lambda z: np.exp(-(z**2)) - c = chebyshev_coefficients(f, 30, -1, 5) + c = interpolation_coefficients(f, 30, -1, 5) sites = 6 interval_x = RegularHalfOpenInterval(-0.5, 2, 2**sites) interval_y = RegularHalfOpenInterval(-0.5, 3, 2**sites) mps_x_plus_y = mps_tensor_sum( [mps_interval(interval_y), mps_interval(interval_x)] ) + tol = 1e-10 strategy = DEFAULT_CHEBYSHEV_STRATEGY.replace( - tolerance=1e-15, simplification_tolerance=1e-15 + tolerance=tol**2, simplification_tolerance=tol**2 + ) + mps_cheb_clen = cheb2mps( + c, initial_mps=mps_x_plus_y, strategy=strategy, clenshaw=True + ) + mps_cheb_poly = cheb2mps( + c, initial_mps=mps_x_plus_y, strategy=strategy, clenshaw=False ) - mps_cheb = cheb2mps(c, x=mps_x_plus_y, strategy=strategy) X, Y = np.meshgrid(interval_x.to_vector(), interval_y.to_vector()) Z_vector = f(X + Y) - Z_mps = mps_cheb.to_vector().reshape([2**sites, 2**sites]) - self.assertSimilar(Z_vector, Z_mps) + Z_mps_clen = mps_cheb_clen.to_vector().reshape([2**sites, 2**sites]) + Z_mps_poly = mps_cheb_poly.to_vector().reshape([2**sites, 2**sites]) + self.assertSimilar(Z_vector, Z_mps_clen) + self.assertSimilar(Z_vector, Z_mps_poly) class TestChebyshevMPO(TestCase): - pass + def test_gaussian_mpo(self): + """ + Tests the interpolation of a diagonal MPO representing a gaussian. + """ + a, b, n = -1, 1, 5 + dx = (b - a) / 2**n + x = np.linspace(a, b, 2**n, endpoint=False) + + f = lambda x: np.sin(-(x**2)) + coefficients = interpolation_coefficients(f) + + I = MPS([np.ones((1, 2, 1))] * n) + mpo_x = x_mpo(n, a, dx) + mpo_gaussian_clen = cheb2mpo(coefficients, mpo_x, clenshaw=True) + mpo_gaussian_poly = cheb2mpo(coefficients, mpo_x, clenshaw=False) + self.assertSimilar(f(x), mpo_gaussian_clen.apply(I).to_vector()) + self.assertSimilar(f(x), mpo_gaussian_poly.apply(I).to_vector()) diff --git a/tests/test_analysis/test_mesh.py b/tests/test_analysis/test_mesh.py index 493217fd..bdb5c484 100644 --- a/tests/test_analysis/test_mesh.py +++ b/tests/test_analysis/test_mesh.py @@ -42,6 +42,18 @@ def test_regular_chebyshev_zeros_interval_constructor(self): self.assertEqual([I[0], I[1]], list(I)) self.assertEqual([I[0], I[1]], [x for x in I]) + def test_regular_chebyshev_extrema_interval_constructor(self): + I = ChebyshevExtremaInterval(-1, 1, 2) + self.assertEqual(I.start, -1) + self.assertEqual(I.stop, 1) + self.assertEqual(I.size, 2) + self.assertEqual(len(I), 2) + self.assertAlmostEqual(I[0], -1) + self.assertAlmostEqual(I[1], 1) + + self.assertEqual([I[0], I[1]], list(I)) + self.assertEqual([I[0], I[1]], [x for x in I]) + def test_rescaled_chebyshev_zeros_interval(self): f = lambda x: 2 * x + 2 # Affine transformation I = ChebyshevZerosInterval(f(-1), f(1), 2) From 563b7925eb5190dc870ea4efaa036238f7006483 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Juan=20Jos=C3=A9=20Rodr=C3=ADguez=20Aldavero?= Date: Mon, 13 May 2024 12:30:03 +0200 Subject: [PATCH 6/7] Fix integrals --- src/seemps/analysis/integrals.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/seemps/analysis/integrals.py b/src/seemps/analysis/integrals.py index 71c88e29..55b8fe18 100644 --- a/src/seemps/analysis/integrals.py +++ b/src/seemps/analysis/integrals.py @@ -6,7 +6,7 @@ from ..truncate import simplify from ..qft import iqft, qft_flip from .mesh import RegularHalfOpenInterval, Mesh -from .factories import mps_tensor_product, mps_affine_transformation +from .factories import mps_tensor_product, mps_affine from .cross import cross_interpolation, CrossStrategy @@ -284,7 +284,7 @@ def func(k): mps_v = mps_k2 * mps_phase mps = (1 / sqrt(2) ** sites) * qft_flip(iqft(mps_v, strategy=strategy)) - return mps_affine_transformation(mps, (-1, 1), (start, stop)).as_mps() + return mps_affine(mps, (-1, 1), (start, stop)).as_mps() # TODO: Consider if this helper function is necessary From 0e6264a660a4948188d8460541fce8a6054be264 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Juan=20Jos=C3=A9=20Rodr=C3=ADguez=20Aldavero?= Date: Mon, 13 May 2024 13:33:25 +0200 Subject: [PATCH 7/7] Fix mypy issues --- src/seemps/analysis/factories.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/seemps/analysis/factories.py b/src/seemps/analysis/factories.py index 3d559321..7698cb64 100644 --- a/src/seemps/analysis/factories.py +++ b/src/seemps/analysis/factories.py @@ -19,8 +19,8 @@ ) COMPUTER_PRECISION = SIMPLIFICATION_STRATEGY.replace( - tolerance=np.finfo(np.double).eps, - simplification_tolerance=np.finfo(np.double).eps, + tolerance=float(np.finfo(np.double).eps), + simplification_tolerance=float(np.finfo(np.double).eps), simplify=Simplification.DO_NOT_SIMPLIFY, method=Truncation.RELATIVE_SINGULAR_VALUE, normalize=False,