diff --git a/src/seemps/analysis/chebyshev.py b/src/seemps/analysis/chebyshev.py index da5dbb6b..094e63ef 100644 --- a/src/seemps/analysis/chebyshev.py +++ b/src/seemps/analysis/chebyshev.py @@ -1,5 +1,5 @@ from __future__ import annotations -from typing import Callable, Optional, Union +from typing import Callable, Optional import numpy as np from scipy.fft import dct # type: ignore @@ -7,9 +7,8 @@ from .mesh import ChebyshevZerosInterval, Interval from .factories import mps_interval -from .sampling import infinity_norm from ..operators import MPO -from ..state import MPS, Strategy, DEFAULT_STRATEGY, Truncation, Simplification +from ..state import MPS, Strategy, Truncation, Simplification from ..truncate import simplify @@ -24,10 +23,11 @@ def chebyshev_coefficients( Returns the Chebyshev coefficients for a given function on a specified interval using the Discrete Cosine Transform (DCT II). - The accuracy of the Chebyshev approximation is correlated to the magnitude - of the last few coefficients in the series (depending on their periodicity), - with smaller absolute values typically indicating a better approximation of - the function. + The error of the Chebyshev approximation is related to the magnitude of the + last (or second-to-last) coefficient of the series. As a rule of thumb, the + error is given by the sum of all the neglected coefficients, which can be + close to the magnitude of the last coefficient if the series decays sufficiently + fast (for example, exponentially fast for analytical functions). Parameters ---------- @@ -134,7 +134,7 @@ def cheb2mps( L = len(x_mps) I = MPS([np.ones((1, 2, 1))] * L) if np.abs(b) > np.finfo(np.float64).eps: - x_mpo = (x_mps + b * I).join() + x_mps = (x_mps + b * I).join() else: raise Exception("In cheb2mps, either domain or an MPS must be provided.") diff --git a/src/seemps/analysis/factories.py b/src/seemps/analysis/factories.py index 4f802253..3c249d96 100644 --- a/src/seemps/analysis/factories.py +++ b/src/seemps/analysis/factories.py @@ -1,7 +1,7 @@ import numpy as np from typing import List -from ..state import MPS, Strategy +from ..state import MPS, Strategy, DEFAULT_STRATEGY from ..truncate import simplify from .mesh import ( Interval, @@ -10,10 +10,6 @@ ChebyshevZerosInterval, ) -# TODO: These descriptions are wrong. `sites` is not the number of -# discretization points, but rather the number of qubits used for this -# representation. The size is 2**sites - def mps_equispaced(start: float, stop: float, sites: int): """ @@ -26,12 +22,12 @@ def mps_equispaced(start: float, stop: float, sites: int): stop : float The end of the interval. sites : int - The number of discretization points. + The number of sites or qubits for the MPS. Returns ------- MPS - An MPS representing the equispaced discretization of the interval. + An MPS representing an equispaced discretization within [start, stop]. """ step = (stop - start) / 2**sites tensor_1 = np.zeros((1, 2, 2)) @@ -47,7 +43,7 @@ def mps_equispaced(start: float, stop: float, sites: int): return MPS(tensors) -def mps_exponential(start: float, stop: float, sites: int, c: complex) -> MPS: +def mps_exponential(start: float, stop: float, sites: int, c: complex = 1) -> MPS: """ Returns an MPS representing an exponential function discretized over an interval. @@ -58,7 +54,7 @@ def mps_exponential(start: float, stop: float, sites: int, c: complex) -> MPS: stop : float The end of the interval. sites : int - The number of discretization points. + The number of sites or qubits for the MPS. c : complex The coefficient in the exponent of the exponential function. @@ -84,7 +80,37 @@ def mps_exponential(start: float, stop: float, sites: int, c: complex) -> MPS: return MPS(tensors) -def mps_cosine(start: float, stop: float, sites: int) -> MPS: +def mps_sine( + start: float, stop: float, sites: int, strategy: Strategy = DEFAULT_STRATEGY +) -> MPS: + """ + Returns an MPS representing a sine function discretized over an interval. + + Parameters + ---------- + start : float + The start of the interval. + stop : float + The end of the interval. + sites : int + The number of sites or qubits for the MPS. + strategy : Strategy, default = DEFAULT_STRATEGY + The MPS simplification strategy to apply. + + Returns + ------- + MPS + An MPS representing the discretized sine function over the interval. + """ + mps_1 = mps_exponential(start, stop, sites, c=1j) + mps_2 = mps_exponential(start, stop, sites, c=-1j) + + return simplify(-0.5j * (mps_1 - mps_2), strategy=strategy) + + +def mps_cosine( + start: float, stop: float, sites: int, strategy: Strategy = DEFAULT_STRATEGY +) -> MPS: """ Returns an MPS representing a cosine function discretized over an interval. @@ -95,7 +121,9 @@ def mps_cosine(start: float, stop: float, sites: int) -> MPS: stop : float The end of the interval. sites : int - The number of discretization points. + The number of sites or qubits for the MPS. + strategy : Strategy, default = DEFAULT_STRATEGY + The MPS simplification strategy to apply. Returns ------- @@ -105,11 +133,10 @@ def mps_cosine(start: float, stop: float, sites: int) -> MPS: mps_1 = mps_exponential(start, stop, sites, c=1j) mps_2 = mps_exponential(start, stop, sites, c=-1j) - return simplify(0.5 * (mps_1 + mps_2)) + return simplify(0.5 * (mps_1 + mps_2), strategy=strategy) -# TODO: Eliminate the `rescale` argument now that we have `Interval.map_to` -def mps_interval(interval: Interval, rescale: bool = False): +def mps_interval(interval: Interval, strategy: Strategy = DEFAULT_STRATEGY): """ Returns an MPS corresponding to a specific type of interval (open, closed, or Chebyshev zeros). @@ -117,16 +144,16 @@ def mps_interval(interval: Interval, rescale: bool = False): ---------- interval : Interval The interval object containing start and stop points and the interval type. - rescale : bool, optional - Flag to rescale the interval to [-1, 1]. + strategy : Strategy, default = DEFAULT_STRATEGY + The MPS simplification strategy to apply. Returns ------- MPS An MPS representing the interval according to its type. """ - start = interval.start if not rescale else -1 - stop = interval.stop if not rescale else 1 + start = interval.start + stop = interval.stop sites = int(np.log2(interval.size)) if isinstance(interval, RegularHalfOpenInterval): return mps_equispaced(start, stop, sites) @@ -136,7 +163,7 @@ def mps_interval(interval: Interval, rescale: bool = False): elif isinstance(interval, ChebyshevZerosInterval): start_mapped = np.pi / (2 ** (sites + 1)) stop_mapped = np.pi + start_mapped - return -1.0 * mps_cosine(start_mapped, stop_mapped, sites) + return -1.0 * mps_cosine(start_mapped, stop_mapped, sites, strategy=strategy) else: raise ValueError(f"Unsupported interval type {type(interval)}") @@ -160,7 +187,7 @@ def mps_tensor_product(mps_list: List[MPS]) -> MPS: return MPS(flattened_sites) -def mps_tensor_sum(mps_list: List[MPS], strategy: Strategy = Strategy()) -> MPS: +def mps_tensor_sum(mps_list: List[MPS], strategy: Strategy = DEFAULT_STRATEGY) -> MPS: """ Returns the tensor sum of a list of MPS. @@ -168,7 +195,7 @@ def mps_tensor_sum(mps_list: List[MPS], strategy: Strategy = Strategy()) -> MPS: ---------- mps_list : List[MPS] The list of MPS objects to sum. - strategy : Strategy, optional + strategy : Strategy, default = DEFAULT_STRATEGY The MPS simplification strategy to apply. Returns diff --git a/src/seemps/analysis/mesh.py b/src/seemps/analysis/mesh.py index 179d1781..6b09adad 100644 --- a/src/seemps/analysis/mesh.py +++ b/src/seemps/analysis/mesh.py @@ -39,6 +39,9 @@ def to_vector(self) -> np.ndarray: def map_to(self, start: float, stop: float) -> Interval: return type(self)(start, stop, self.size) + def update_size(self, size: int) -> Interval: + return type(self)(self.start, self.stop, size) + def __iter__(self) -> Iterator: return (self[i] for i in range(self.size)) diff --git a/tests/test_analysis_factories.py b/tests/test_analysis_factories.py new file mode 100644 index 00000000..6728050a --- /dev/null +++ b/tests/test_analysis_factories.py @@ -0,0 +1,83 @@ +import numpy as np +from seemps.analysis import ( + mps_equispaced, + mps_exponential, + mps_sine, + mps_cosine, + RegularHalfOpenInterval, + RegularClosedInterval, + ChebyshevZerosInterval, + mps_interval, + mps_tensor_sum, + mps_tensor_product, +) + +from .tools import TestCase + + +class TestMPSFactories(TestCase): + def test_mps_equispaced(self): + self.assertSimilar( + mps_equispaced(-1, 1, 5).to_vector(), + np.linspace(-1, 1, 2**5, endpoint=False), + ) + + def test_mps_exponential(self): + self.assertSimilar( + mps_exponential(-1, 1, 5, c=1).to_vector(), + np.exp(np.linspace(-1, 1, 2**5, endpoint=False)), + ) + self.assertSimilar( + mps_exponential(-1, 1, 5, c=-1).to_vector(), + np.exp(-np.linspace(-1, 1, 2**5, endpoint=False)), + ) + + def test_mps_sine(self): + self.assertSimilar( + mps_sine(-1, 1, 5).to_vector(), + np.sin(np.linspace(-1, 1, 2**5, endpoint=False)), + ) + + def test_mps_cosine(self): + self.assertSimilar( + mps_cosine(-1, 1, 5).to_vector(), + np.cos(np.linspace(-1, 1, 2**5, endpoint=False)), + ) + + def test_mps_interval(self): + start = -1 + stop = 1 + sites = 5 + mps_half_open = mps_interval(RegularHalfOpenInterval(start, stop, 2**sites)) + mps_closed = mps_interval(RegularClosedInterval(start, stop, 2**sites)) + mps_zeros = mps_interval(ChebyshevZerosInterval(start, stop, 2**sites)) + zeros = lambda d: np.array( + [np.cos(np.pi * (2 * k - 1) / (2 * d)) for k in range(d, 0, -1)] + ) + self.assertSimilar( + mps_half_open, np.linspace(start, stop, 2**sites, endpoint=False) + ) + self.assertSimilar( + mps_closed, np.linspace(start, stop, 2**sites, endpoint=True) + ) + self.assertSimilar(mps_zeros, zeros(2**sites)) + + +class TestMPSOperations(TestCase): + def test_tensor_product(self): + sites = 5 + interval = RegularHalfOpenInterval(-1, 2, 2**sites) + mps_x = mps_interval(interval) + mps_x_times_y = mps_tensor_product([mps_x, mps_x]) + Z_mps = mps_x_times_y.to_vector().reshape((2**sites, 2**sites)) + X, Y = np.meshgrid(interval.to_vector(), interval.to_vector()) + self.assertSimilar(Z_mps, X * Y) + + def test_tensor_sum(self): + sites = 5 + interval = RegularHalfOpenInterval(-1, 2, 2**sites) + mps_x = mps_interval(interval) + mps_x_plus_y = mps_tensor_sum([mps_x, mps_x]) + Z_mps = mps_x_plus_y.to_vector().reshape((2**sites, 2**sites)) + X, Y = np.meshgrid(interval.to_vector(), interval.to_vector()) + self.assertSimilar(Z_mps, X + Y) diff --git a/tests/test_chebyshev.py b/tests/test_chebyshev.py index c37bcd7b..c14addad 100644 --- a/tests/test_chebyshev.py +++ b/tests/test_chebyshev.py @@ -7,13 +7,10 @@ chebyshev_approximation, cheb2mps, mps_tensor_sum, - mps_tensor_product, mps_interval, ) from seemps.analysis.chebyshev import DEFAULT_CHEBYSHEV_STRATEGY -from seemps.state import Strategy -from seemps.operators import MPO -from numpy.polynomial import Polynomial, Chebyshev +from numpy.polynomial import Chebyshev from .tools import TestCase @@ -101,105 +98,49 @@ def test_chebyshev_coefficients_gaussian_integral(self): class TestChebyshevMPS(TestCase): def test_gaussian_1d(self): - start = 0 - stop = 2 - sites = 5 - order = 20 f = lambda x: np.exp(-(x**2)) - interval = RegularHalfOpenInterval(start, stop, 2**sites) - mps_cheb = chebyshev_approximation(f, order, interval) + interval = RegularHalfOpenInterval(-1, 2, 2**5) + mps_cheb = chebyshev_approximation(f, 20, interval) self.assertSimilar(f(interval.to_vector()), mps_cheb.to_vector()) def test_gaussian_derivative_1d(self): - start = -1 - stop = 2 - sites = 5 - order = 22 f = lambda x: np.exp(-(x**2)) f_diff = lambda x: -2 * x * np.exp(-(x**2)) - interval = RegularHalfOpenInterval(start, stop, 2**sites) - mps_cheb = chebyshev_approximation(f, order, interval, differentiation_order=1) + interval = RegularHalfOpenInterval(-1, 2, 2**5) + mps_cheb = chebyshev_approximation(f, 22, interval, differentiation_order=1) self.assertSimilar(f_diff(interval.to_vector()), mps_cheb.to_vector()) def test_gaussian_integral_1d(self): - start = -1 - stop = 2 - sites = 5 - order = 20 - f_intg = lambda x: (np.sqrt(np.pi) / 2) * (erf(x) - erf(start)) - interval = RegularHalfOpenInterval(start, stop, 2**sites) - mps_cheb = chebyshev_approximation(f_intg, order, interval) + 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, 20, interval) self.assertSimilar(f_intg(interval.to_vector()), mps_cheb.to_vector()) def test_gaussian_integral_1d_b(self): - start = -1 - stop = 2 - sites = 5 - order = 20 f = lambda x: np.exp(-(x**2)) - f_intg = lambda x: (np.sqrt(np.pi) / 2) * (erf(x) - erf(start)) - interval = RegularHalfOpenInterval(start, stop, 2**sites) - mps_cheb = chebyshev_approximation(f, order, interval, differentiation_order=-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, 20, interval, differentiation_order=-1) self.assertSimilar(f_intg(interval.to_vector()), mps_cheb.to_vector()) - # TODO: This does not have a place here. Move tests to separate file - def test_tensor_product(self): - start = -1 - stop = 2 - sites = 5 - interval = RegularHalfOpenInterval(start, stop, 2**sites) - vector = interval.to_vector() - mps_intv = mps_interval(interval) - mps_domain = mps_tensor_product([mps_intv, mps_intv]) - Z_mps = mps_domain.to_vector().reshape((2**sites, 2**sites)) - X, Y = np.meshgrid(vector, vector) - self.assertSimilar(Z_mps, X * Y) - - # TODO: This does not have a place here. Move tests to separate file - def test_tensor_sum(self): - start = -2 - stop = 2 - sites = 5 - interval = RegularHalfOpenInterval(start, stop, 2**sites) - vector = interval.to_vector() - mps_intv = mps_interval(interval) - mps_domain = mps_tensor_sum([mps_intv, mps_intv]) - Z_mps = mps_domain.to_vector().reshape((2**sites, 2**sites)) - X, Y = np.meshgrid(vector, vector) - self.assertSimilar(Z_mps, X + Y) - def test_gaussian_2d(self): - start = -0.5 - stop = 0.5 - sites = 5 - order = 40 - interval = RegularHalfOpenInterval(start, stop, 2**sites) - vector = interval.to_vector() - func = lambda x: np.exp(-(x**2)) + f = lambda z: np.exp(-(z**2)) + c = chebyshev_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)] + ) strategy = DEFAULT_CHEBYSHEV_STRATEGY.replace( tolerance=1e-15, simplification_tolerance=1e-15 ) - mps_x_plus_y = mps_tensor_sum([mps_interval(interval)] * 2) - mps_cheb = cheb2mps( - # Note that we give `chebyshev_coefficients` an interval - # that covers the smallest and largest values of `X + Y` - chebyshev_coefficients(func, order, 2 * start, 2 * stop), - x=mps_x_plus_y, - strategy=strategy, - ) - Z_mps = mps_cheb.to_vector().reshape((2**sites, 2**sites)) - X, Y = np.meshgrid(vector, vector) - Z_vector = func(X + Y) - self.assertSimilar(Z_mps, Z_vector) + 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) class TestChebyshevMPO(TestCase): pass - # def test_gaussian_on_diagonal_mpo(self): - # start = -1 - # stop = 2 - # sites = 5 - # order = 10 - # f = lambda x: np.exp(-(x**2)) - # mpo_domain = position_mpo(start, stop, sites) - # mpo_cheb = chebyshev_approximation(f, order, mpo_domain)