diff --git a/.github/workflows/CI.yml b/.github/workflows/CI.yml new file mode 100644 index 00000000..83d85491 --- /dev/null +++ b/.github/workflows/CI.yml @@ -0,0 +1,29 @@ +name: CI + +on: [push, pull_request] + +jobs: + test: + runs-on: ${{ matrix.os }} + strategy: + fail-fast: true + matrix: + os: ["ubuntu-latest", "macos-latest"] + python-version: ["3.5", "3.6", "3.7", "3.8", "3.9"] + + steps: + - name: Checkout source + uses: actions/checkout@v2 + + - name: Setup python + uses: actions/setup-python@v2 + with: + python-version: ${{ matrix.python-version }} + architecture: x64 + + - name: Install + run: | + pip install -e . + pip install pytest + - name: Run tests + run: pytest diff --git a/pyerrors/correlators.py b/pyerrors/correlators.py index 98e20e6b..33809267 100644 --- a/pyerrors/correlators.py +++ b/pyerrors/correlators.py @@ -1,11 +1,11 @@ import numpy as np import autograd.numpy as anp +import matplotlib.pyplot as plt import scipy.linalg from .pyerrors import Obs, dump_object from .fits import standard_fit from .linalg import eigh, mat_mat_op from .roots import find_root -import matplotlib.pyplot as plt class Corr: @@ -24,7 +24,7 @@ class Corr: def __init__(self, data_input, padding_front=0, padding_back=0, prange=None): # All data_input should be a list of things at different timeslices. This needs to be verified - if not (isinstance(data_input, list)): + if not isinstance(data_input, list): raise TypeError('Corr__init__ expects a list of timeslices.') # data_input can have multiple shapes. The simplest one is a list of Obs. # We check, if this is the case @@ -115,9 +115,9 @@ class Corr: def plottable(self): if self.N != 1: raise Exception("Can only make Corr[N=1] plottable") # We could also autoproject to the groundstate or expect vectors, but this is supposed to be a super simple function. - x_list = [x for x in range(self.T) if (not self.content[x] is None)] - y_list = [y[0].value for y in self.content if (y is not None)] - y_err_list = [y[0].dvalue for y in self.content if (y is not None)] + x_list = [x for x in range(self.T) if not self.content[x] is None] + y_list = [y[0].value for y in self.content if y is not None] + y_err_list = [y[0].dvalue for y in self.content if y is not None] return x_list, y_list, y_err_list diff --git a/tests/test_linalg.py b/tests/test_linalg.py new file mode 100644 index 00000000..9d7f6377 --- /dev/null +++ b/tests/test_linalg.py @@ -0,0 +1,56 @@ +import sys +sys.path.append('..') +import autograd.numpy as np +import os +import random +import math +import string +import copy +import scipy.optimize +from scipy.odr import ODR, Model, Data, RealData +import pyerrors as pe +import pytest + +def test_matrix_functions(): + dim = 3 + int(4 * np.random.rand()) + print(dim) + matrix = [] + for i in range(dim): + row = [] + for j in range(dim): + row.append(pe.pseudo_Obs(np.random.rand(), 0.2 + 0.1 * np.random.rand(), 'e1')) + matrix.append(row) + matrix = np.array(matrix) @ np.identity(dim) + + # Check inverse of matrix + inv = pe.linalg.mat_mat_op(np.linalg.inv, matrix) + check_inv = matrix @ inv + + for (i, j), entry in np.ndenumerate(check_inv): + entry.gamma_method() + if(i == j): + assert math.isclose(entry.value, 1.0, abs_tol=1e-9), 'value ' + str(i) + ',' + str(j) + ' ' + str(entry.value) + else: + assert math.isclose(entry.value, 0.0, abs_tol=1e-9), 'value ' + str(i) + ',' + str(j) + ' ' + str(entry.value) + assert math.isclose(entry.dvalue, 0.0, abs_tol=1e-9), 'dvalue ' + str(i) + ',' + str(j) + ' ' + str(entry.dvalue) + + # Check Cholesky decomposition + sym = np.dot(matrix, matrix.T) + cholesky = pe.linalg.mat_mat_op(np.linalg.cholesky, sym) + check = cholesky @ cholesky.T + + for (i, j), entry in np.ndenumerate(check): + diff = entry - sym[i, j] + diff.gamma_method() + assert math.isclose(diff.value, 0.0, abs_tol=1e-9), 'value ' + str(i) + ',' + str(j) + assert math.isclose(diff.dvalue, 0.0, abs_tol=1e-9), 'dvalue ' + str(i) + ',' + str(j) + + # Check eigh + e, v = pe.linalg.eigh(sym) + for i in range(dim): + tmp = sym @ v[:, i] - v[:, i] * e[i] + for j in range(dim): + tmp[j].gamma_method() + assert math.isclose(tmp[j].value, 0.0, abs_tol=1e-9), 'value ' + str(i) + ',' + str(j) + assert math.isclose(tmp[j].dvalue, 0.0, abs_tol=1e-9), 'dvalue ' + str(i) + ',' + str(j) + diff --git a/tests/test_pyerrors.py b/tests/test_pyerrors.py index caf36760..0f61e8d2 100644 --- a/tests/test_pyerrors.py +++ b/tests/test_pyerrors.py @@ -1,18 +1,11 @@ -import sys -sys.path.append('..') import autograd.numpy as np import os import random -import math import string import copy -import scipy.optimize -from scipy.odr import ODR, Model, Data, RealData import pyerrors as pe import pytest -test_iterations = 100 - def test_dump(): value = np.random.normal(5, 10) dvalue = np.abs(np.random.normal(0, 1)) @@ -32,9 +25,9 @@ def test_comparison(): assert (value1 < value2) == (test_obs1 < test_obs2) -def test_man_grad(): - a = pe.pseudo_Obs(17,2.9,'e1') - b = pe.pseudo_Obs(4,0.8,'e1') +def test_function_overloading(): + a = pe.pseudo_Obs(17, 2.9, 'e1') + b = pe.pseudo_Obs(4, 0.8, 'e1') fs = [lambda x: x[0] + x[1], lambda x: x[1] + x[0], lambda x: x[0] - x[1], lambda x: x[1] - x[0], lambda x: x[0] * x[1], lambda x: x[1] * x[0], lambda x: x[0] / x[1], lambda x: x[1] / x[0], @@ -51,8 +44,8 @@ def test_man_grad(): def test_overloading_vectorization(): - a = np.array([5, 4, 8]) - b = pe.pseudo_Obs(4,0.8,'e1') + a = np.random.randint(0, 100, 10) + b = pe.pseudo_Obs(4, 0.8, 'e1') assert [o.value for o in a * b] == [o.value for o in b * a] assert [o.value for o in a + b] == [o.value for o in b + a] @@ -61,8 +54,7 @@ def test_overloading_vectorization(): assert [o.value for o in b / a] == [o.value for o in [b / p for p in a]] -@pytest.mark.parametrize("n", np.arange(test_iterations // 10)) -def test_covariance_is_variance(n): +def test_covariance_is_variance(): value = np.random.normal(5, 10) dvalue = np.abs(np.random.normal(0, 1)) test_obs = pe.pseudo_Obs(value, dvalue, 't') @@ -73,8 +65,7 @@ def test_covariance_is_variance(n): assert np.abs(test_obs.dvalue ** 2 - pe.covariance(test_obs, test_obs)) <= 10 * np.finfo(np.float).eps -@pytest.mark.parametrize("n", np.arange(test_iterations // 10)) -def test_fft(n): +def test_fft(): value = np.random.normal(5, 100) dvalue = np.abs(np.random.normal(0, 5)) test_obs1 = pe.pseudo_Obs(value, dvalue, 't', int(500 + 1000 * np.random.rand())) @@ -85,106 +76,7 @@ def test_fft(n): assert np.abs(test_obs1.dvalue - test_obs2.dvalue) <= 10 * max(test_obs1.dvalue, test_obs2.dvalue) * np.finfo(np.float).eps -@pytest.mark.parametrize('n', np.arange(test_iterations // 10)) -def test_standard_fit(n): - dim = 10 + int(30 * np.random.rand()) - x = np.arange(dim) - y = 2 * np.exp(-0.06 * x) + np.random.normal(0.0, 0.15, dim) - yerr = 0.1 + 0.1 * np.random.rand(dim) - - oy = [] - for i, item in enumerate(x): - oy.append(pe.pseudo_Obs(y[i], yerr[i], str(i))) - - def f(x, a, b): - return a * np.exp(-b * x) - - popt, pcov = scipy.optimize.curve_fit(f, x, y, sigma=[o.dvalue for o in oy], absolute_sigma=True) - - def func(a, x): - y = a[0] * np.exp(-a[1] * x) - return y - - beta = pe.fits.standard_fit(x, oy, func) - - pe.Obs.e_tag_global = 5 - for i in range(2): - beta[i].gamma_method(e_tag=5, S=1.0) - assert math.isclose(beta[i].value, popt[i], abs_tol=1e-5) - assert math.isclose(pcov[i, i], beta[i].dvalue ** 2, abs_tol=1e-3) - assert math.isclose(pe.covariance(beta[0], beta[1]), pcov[0, 1], abs_tol=1e-3) - pe.Obs.e_tag_global = 0 - - chi2_pyerrors = np.sum(((f(x, *[o.value for o in beta]) - y) / yerr) ** 2) / (len(x) - 2) - chi2_scipy = np.sum(((f(x, *popt) - y) / yerr) ** 2) / (len(x) - 2) - assert math.isclose(chi2_pyerrors, chi2_scipy, abs_tol=1e-10) - - -@pytest.mark.parametrize('n', np.arange(test_iterations // 10)) -def test_odr_fit(n): - dim = 10 + int(30 * np.random.rand()) - x = np.arange(dim) + np.random.normal(0.0, 0.15, dim) - xerr = 0.1 + 0.1 * np.random.rand(dim) - y = 2 * np.exp(-0.06 * x) + np.random.normal(0.0, 0.15, dim) - yerr = 0.1 + 0.1 * np.random.rand(dim) - - ox = [] - for i, item in enumerate(x): - ox.append(pe.pseudo_Obs(x[i], xerr[i], str(i))) - - oy = [] - for i, item in enumerate(x): - oy.append(pe.pseudo_Obs(y[i], yerr[i], str(i))) - - def f(x, a, b): - return a * np.exp(-b * x) - - def func(a, x): - y = a[0] * np.exp(-a[1] * x) - return y - - data = RealData([o.value for o in ox], [o.value for o in oy], sx=[o.dvalue for o in ox], sy=[o.dvalue for o in oy]) - model = Model(func) - odr = ODR(data, model, [0,0], partol=np.finfo(np.float).eps) - odr.set_job(fit_type=0, deriv=1) - output = odr.run() - - beta = pe.fits.odr_fit(ox, oy, func) - - pe.Obs.e_tag_global = 5 - for i in range(2): - beta[i].gamma_method(e_tag=5, S=1.0) - assert math.isclose(beta[i].value, output.beta[i], rel_tol=1e-5) - assert math.isclose(output.cov_beta[i,i], beta[i].dvalue**2, rel_tol=2.5e-1), str(output.cov_beta[i,i]) + ' ' + str(beta[i].dvalue**2) - assert math.isclose(pe.covariance(beta[0], beta[1]), output.cov_beta[0,1], rel_tol=2.5e-1) - pe.Obs.e_tag_global = 0 - - -@pytest.mark.parametrize('n', np.arange(test_iterations // 10)) -def test_odr_derivatives(n): - x = [] - y = [] - x_err = 0.01 - y_err = 0.01 - - for n in np.arange(1, 9, 2): - loc_xvalue = n + np.random.normal(0.0, x_err) - x.append(pe.pseudo_Obs(loc_xvalue, x_err, str(n))) - y.append(pe.pseudo_Obs((lambda x: x ** 2 - 1)(loc_xvalue) + - np.random.normal(0.0, y_err), y_err, str(n))) - - def func(a, x): - return a[0] + a[1] * x ** 2 - - fit1 = pe.fits.odr_fit(x, y, func) - - tfit = pe.fits.fit_general(x, y, func, base_step=0.1, step_ratio=1.1, num_steps=20) - assert np.abs(np.max(np.array(list(fit1[1].deltas.values())) - - np.array(list(tfit[1].deltas.values())))) < 10e-8 - - -@pytest.mark.parametrize('n', np.arange(test_iterations)) -def test_covariance_symmetry(n): +def test_covariance_symmetry(): value1 = np.random.normal(5, 10) dvalue1 = np.abs(np.random.normal(0, 1)) test_obs1 = pe.pseudo_Obs(value1, dvalue1, 't') @@ -199,8 +91,7 @@ def test_covariance_symmetry(n): assert np.abs(cov_ab) < test_obs1.dvalue * test_obs2.dvalue * (1 + 10 * np.finfo(np.float).eps) -@pytest.mark.parametrize('n', np.arange(test_iterations)) -def test_gamma_method(n): +def test_gamma_method(): # Construct pseudo Obs with random shape value = np.random.normal(5, 10) dvalue = np.abs(np.random.normal(0, 1)) @@ -213,28 +104,7 @@ def test_gamma_method(n): assert abs(test_obs.dvalue - dvalue) < 1e-10 * dvalue -@pytest.mark.parametrize('n', np.arange(test_iterations)) -def test_overloading(n): - # Construct pseudo Obs with random shape - obs_list = [] - for i in range(5): - value = np.abs(np.random.normal(5, 2)) + 2.0 - dvalue = np.abs(np.random.normal(0, 0.1)) + 1e-5 - obs_list.append(pe.pseudo_Obs(value, dvalue, 't', 2000)) - - # Test if the error is processed correctly - def f(x): - return x[0] * x[1] + np.sin(x[2]) * np.exp(x[3] / x[1] / x[0]) - np.sqrt(2) / np.cosh(x[4] / x[0]) - - o_obs = f(obs_list) - d_obs = pe.derived_observable(f, obs_list) - - assert np.max(np.abs((o_obs.deltas['t'] - d_obs.deltas['t']) / o_obs.deltas['t'])) < 1e-7, str(obs_list) - assert np.abs((o_obs.value - d_obs.value) / o_obs.value) < 1e-10 - - -@pytest.mark.parametrize('n', np.arange(test_iterations)) -def test_derived_observables(n): +def test_derived_observables(): # Construct pseudo Obs with random shape test_obs = pe.pseudo_Obs(2, 0.1 * (1 + np.random.rand()), 't', int(1000 * (1 + np.random.rand()))) @@ -257,8 +127,7 @@ def test_derived_observables(n): assert i_am_one.e_ddvalue['t'] <= 2 * np.finfo(np.float).eps -@pytest.mark.parametrize('n', np.arange(test_iterations // 10)) -def test_multi_ens_system(n): +def test_multi_ens_system(): names = [] for i in range(100 + int(np.random.rand() * 50)): tmp_string = '' @@ -276,8 +145,7 @@ def test_multi_ens_system(n): assert sorted(x for y in sorted(new_obs.e_content.values()) for x in y) == sorted(new_obs.names) -@pytest.mark.parametrize('n', np.arange(test_iterations)) -def test_overloaded_functions(n): +def test_overloaded_functions(): funcs = [np.exp, np.log, np.sin, np.cos, np.tan, np.sinh, np.cosh, np.arcsinh, np.arccosh] deriv = [np.exp, lambda x: 1 / x, np.cos, lambda x: -np.sin(x), lambda x: 1 / np.cos(x) ** 2, np.cosh, np.sinh, lambda x: 1 / np.sqrt(x ** 2 + 1), lambda x: 1 / np.sqrt(x ** 2 - 1)] val = 3 + 0.5 * np.random.rand() @@ -291,49 +159,3 @@ def test_overloaded_functions(n): assert np.max((ad_obs.deltas['t'] - fd_obs.deltas['t']) / ad_obs.deltas['t']) < 1e-8, item.__name__ assert np.abs((ad_obs.value - item(val)) / ad_obs.value) < 1e-10, item.__name__ assert np.abs(ad_obs.dvalue - dval * np.abs(deriv[i](val))) < 1e-6, item.__name__ - - -@pytest.mark.parametrize('n', np.arange(test_iterations // 10)) -def test_matrix_functions(n): - dim = 3 + int(4 * np.random.rand()) - print(dim) - matrix = [] - for i in range(dim): - row = [] - for j in range(dim): - row.append(pe.pseudo_Obs(np.random.rand(), 0.2 + 0.1 * np.random.rand(), 'e1')) - matrix.append(row) - matrix = np.array(matrix) @ np.identity(dim) - - # Check inverse of matrix - inv = pe.linalg.mat_mat_op(np.linalg.inv, matrix) - check_inv = matrix @ inv - - for (i, j), entry in np.ndenumerate(check_inv): - entry.gamma_method() - if(i == j): - assert math.isclose(entry.value, 1.0, abs_tol=1e-9), 'value ' + str(i) + ',' + str(j) + ' ' + str(entry.value) - else: - assert math.isclose(entry.value, 0.0, abs_tol=1e-9), 'value ' + str(i) + ',' + str(j) + ' ' + str(entry.value) - assert math.isclose(entry.dvalue, 0.0, abs_tol=1e-9), 'dvalue ' + str(i) + ',' + str(j) + ' ' + str(entry.dvalue) - - # Check Cholesky decomposition - sym = np.dot(matrix, matrix.T) - cholesky = pe.linalg.mat_mat_op(np.linalg.cholesky, sym) - check = cholesky @ cholesky.T - - for (i, j), entry in np.ndenumerate(check): - diff = entry - sym[i, j] - diff.gamma_method() - assert math.isclose(diff.value, 0.0, abs_tol=1e-9), 'value ' + str(i) + ',' + str(j) - assert math.isclose(diff.dvalue, 0.0, abs_tol=1e-9), 'dvalue ' + str(i) + ',' + str(j) - - # Check eigh - e, v = pe.linalg.eigh(sym) - for i in range(dim): - tmp = sym @ v[:, i] - v[:, i] * e[i] - for j in range(dim): - tmp[j].gamma_method() - assert math.isclose(tmp[j].value, 0.0, abs_tol=1e-9), 'value ' + str(i) + ',' + str(j) - assert math.isclose(tmp[j].dvalue, 0.0, abs_tol=1e-9), 'dvalue ' + str(i) + ',' + str(j) -