diff --git a/examples/05_matrix_operations.ipynb b/examples/05_matrix_operations.ipynb index 926d1734..1f08951b 100644 --- a/examples/05_matrix_operations.ipynb +++ b/examples/05_matrix_operations.ipynb @@ -7,8 +7,7 @@ "outputs": [], "source": [ "import pyerrors as pe\n", - "import numpy as np\n", - "import scipy" + "import numpy as np" ] }, { @@ -27,8 +26,8 @@ "name": "stdout", "output_type": "stream", "text": [ - "[[Obs[4.10(20)] Obs[-1.00(10)]]\n", - " [Obs[-1.00(10)] Obs[1.000(10)]]]\n" + "[[4.10(20) -1.00(10)]\n", + " [-1.00(10) 1.000(10)]]\n" ] } ], @@ -63,8 +62,8 @@ "name": "stdout", "output_type": "stream", "text": [ - "[[Obs[17.81] Obs[-5.1]]\n", - " [Obs[-5.1] Obs[2.0]]]\n" + "[[17.81 -5.1]\n", + " [-5.1 2.0]]\n" ] } ], @@ -88,8 +87,8 @@ "name": "stdout", "output_type": "stream", "text": [ - "[[Obs[4.1] Obs[-1.0]]\n", - " [Obs[-1.0] Obs[1.0]]]\n" + "[[4.1 -1.0]\n", + " [-1.0 1.0]]\n" ] } ], @@ -113,8 +112,8 @@ "name": "stdout", "output_type": "stream", "text": [ - "[[Obs[78.12099999999998] Obs[-22.909999999999997]]\n", - " [Obs[-22.909999999999997] Obs[7.1]]]\n" + "[[78.12099999999998 -22.909999999999997]\n", + " [-22.909999999999997 7.1]]\n" ] } ], @@ -138,8 +137,8 @@ "name": "stdout", "output_type": "stream", "text": [ - "[[Obs[30.161857460980094] Obs[-1.1752011936438014]]\n", - " [Obs[-1.1752011936438014] Obs[1.1752011936438014]]]\n" + "[[30.16185746098009 -1.1752011936438014]\n", + " [-1.1752011936438014 1.1752011936438014]]\n" ] } ], @@ -163,7 +162,7 @@ "name": "stdout", "output_type": "stream", "text": [ - "[Obs[2.00(40)] Obs[1.00(10)]]\n" + "[2.00(40) 1.00(10)]\n" ] } ], @@ -171,8 +170,7 @@ "vec1 = pe.pseudo_Obs(2, 0.4, 'e1')\n", "vec2 = pe.pseudo_Obs(1, 0.1, 'e1')\n", "vector = np.asarray([vec1, vec2])\n", - "for (i), entry in np.ndenumerate(vector):\n", - " entry.gamma_method()\n", + "pe.gm(vector)\n", "print(vector)" ] }, @@ -192,14 +190,13 @@ "name": "stdout", "output_type": "stream", "text": [ - "[Obs[7.2(1.7)] Obs[-1.00(46)]]\n" + "[7.2(1.7) -1.00(45)]\n" ] } ], "source": [ "product = matrix @ vector\n", - "for (i), entry in np.ndenumerate(product):\n", - " entry.gamma_method()\n", + "pe.gm(product)\n", "print(product)" ] }, @@ -245,15 +242,14 @@ "name": "stdout", "output_type": "stream", "text": [ - "[[Obs[2.025(49)] Obs[0.0]]\n", - " [Obs[-0.494(50)] Obs[0.870(29)]]]\n" + "[[2.025(49) 0.0]\n", + " [-0.494(50) 0.870(29)]]\n" ] } ], "source": [ "cholesky = pe.linalg.cholesky(matrix)\n", - "for (i, j), entry in np.ndenumerate(cholesky):\n", - " entry.gamma_method()\n", + "pe.gm(cholesky)\n", "print(cholesky)" ] }, @@ -273,8 +269,8 @@ "name": "stdout", "output_type": "stream", "text": [ - "[[Obs[-8.881784197001252e-16] Obs[0.0]]\n", - " [Obs[0.0] Obs[0.0]]]\n" + "[[-8.881784197001252e-16 0.0]\n", + " [0.0 0.0]]\n" ] } ], @@ -299,18 +295,17 @@ "name": "stdout", "output_type": "stream", "text": [ - "[[Obs[0.494(12)] Obs[0.0]]\n", - " [Obs[0.280(40)] Obs[1.150(39)]]]\n", + "[[0.494(12) 0.0]\n", + " [0.280(40) 1.150(38)]]\n", "Check:\n", - "[[Obs[1.0] Obs[0.0]]\n", - " [Obs[0.0] Obs[1.0]]]\n" + "[[1.0 0.0]\n", + " [0.0 1.0]]\n" ] } ], "source": [ "inv = pe.linalg.inv(cholesky)\n", - "for (i, j), entry in np.ndenumerate(inv):\n", - " entry.gamma_method()\n", + "pe.gm(inv)\n", "print(inv)\n", "print('Check:')\n", "check_inv = cholesky @ inv\n", @@ -335,21 +330,19 @@ "output_type": "stream", "text": [ "Eigenvalues:\n", - "[Obs[0.705(57)] Obs[4.39(19)]]\n", + "[0.705(56) 4.39(20)]\n", "Eigenvectors:\n", - "[[Obs[-0.283(26)] Obs[-0.9592(75)]]\n", - " [Obs[-0.9592(75)] Obs[0.283(26)]]]\n" + "[[-0.283(26) -0.9592(75)]\n", + " [-0.9592(75) 0.283(26)]]\n" ] } ], "source": [ "e, v = pe.linalg.eigh(matrix)\n", - "for (i), entry in np.ndenumerate(e):\n", - " entry.gamma_method()\n", + "pe.gm(e)\n", "print('Eigenvalues:')\n", "print(e)\n", - "for (i, j), entry in np.ndenumerate(v):\n", - " entry.gamma_method()\n", + "pe.gm(v)\n", "print('Eigenvectors:')\n", "print(v)" ] @@ -371,9 +364,9 @@ "output_type": "stream", "text": [ "Check eigenvector 1\n", - "[Obs[-5.551115123125783e-17] Obs[0.0]]\n", + "[-5.551115123125783e-17 0.0]\n", "Check eigenvector 2\n", - "[Obs[0.0] Obs[-2.220446049250313e-16]]\n" + "[0.0 -2.220446049250313e-16]\n" ] } ], @@ -407,7 +400,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.8.10" + "version": "3.11.0" } }, "nbformat": 4, diff --git a/pyerrors/obs.py b/pyerrors/obs.py index a25dceb8..31ed84af 100644 --- a/pyerrors/obs.py +++ b/pyerrors/obs.py @@ -1032,6 +1032,22 @@ class CObs: return f"({self.real:{format_type}}{self.imag:+{significance}}j)" +def gamma_method(x, **kwargs): + """Vectorized version of the gamma_method applicable to lists or arrays of Obs. + + See docstring of pe.Obs.gamma_method for details. + """ + return np.vectorize(lambda o: o.gm(**kwargs))(x) + + +def gm(x, **kwargs): + """Short version of the vectorized gamma_method. + + See docstring of pe.Obs.gamma_method for details + """ + return gamma_method(x, **kwargs) + + def _format_uncertainty(value, dvalue, significance=2): """Creates a string of a value and its error in paranthesis notation, e.g., 13.02(45)""" if dvalue == 0.0 or (not np.isfinite(dvalue)): diff --git a/tests/obs_test.py b/tests/obs_test.py index 924d98cd..cc13c845 100644 --- a/tests/obs_test.py +++ b/tests/obs_test.py @@ -1,4 +1,5 @@ -import autograd.numpy as np +import numpy as np +import autograd.numpy as anp import os import copy import matplotlib.pyplot as plt @@ -126,9 +127,9 @@ def test_function_overloading(): 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], - lambda x: np.exp(x[0]), lambda x: np.sin(x[0]), lambda x: np.cos(x[0]), lambda x: np.tan(x[0]), - lambda x: np.log(x[0]), lambda x: np.sqrt(np.abs(x[0])), - lambda x: np.sinh(x[0]), lambda x: np.cosh(x[0]), lambda x: np.tanh(x[0])] + lambda x: anp.exp(x[0]), lambda x: anp.sin(x[0]), lambda x: anp.cos(x[0]), lambda x: anp.tan(x[0]), + lambda x: anp.log(x[0]), lambda x: anp.sqrt(anp.abs(x[0])), + lambda x: anp.sinh(x[0]), lambda x: anp.cosh(x[0]), lambda x: anp.tanh(x[0])] for i, f in enumerate(fs): t1 = f([a, b]) @@ -306,9 +307,9 @@ def test_derived_observables(): test_obs = pe.pseudo_Obs(2, 0.1 * (1 + np.random.rand()), 't', int(1000 * (1 + np.random.rand()))) # Check if autograd and numgrad give the same result - d_Obs_ad = pe.derived_observable(lambda x, **kwargs: x[0] * x[1] * np.sin(x[0] * x[1]), [test_obs, test_obs]) + d_Obs_ad = pe.derived_observable(lambda x, **kwargs: x[0] * x[1] * anp.sin(x[0] * x[1]), [test_obs, test_obs]) d_Obs_ad.gamma_method() - d_Obs_fd = pe.derived_observable(lambda x, **kwargs: x[0] * x[1] * np.sin(x[0] * x[1]), [test_obs, test_obs], num_grad=True) + d_Obs_fd = pe.derived_observable(lambda x, **kwargs: x[0] * x[1] * anp.sin(x[0] * x[1]), [test_obs, test_obs], num_grad=True) d_Obs_fd.gamma_method() assert d_Obs_ad == d_Obs_fd @@ -1283,12 +1284,14 @@ def test_format_uncertainty(): pe.obs._format_uncertainty(1, np.NaN) pe.obs._format_uncertainty(np.NaN, np.inf) + def test_format(): o1 = pe.pseudo_Obs(0.348, 0.0123, "test") assert o1.__format__("+3") == '+0.3480(123)' assert o1.__format__("+2") == '+0.348(12)' assert o1.__format__(" 2") == ' 0.348(12)' + def test_f_string_obs(): o1 = pe.pseudo_Obs(0.348, 0.0123, "test") print(f"{o1}") @@ -1297,6 +1300,7 @@ def test_f_string_obs(): print(f"{o1:-1}") print(f"{o1: 8}") + def test_f_string_cobs(): o_real = pe.pseudo_Obs(0.348, 0.0123, "test") o_imag = pe.pseudo_Obs(0.348, 0.0123, "test") @@ -1307,7 +1311,24 @@ def test_f_string_cobs(): print(f"{o1:-1}") print(f"{o1: 8}") + def test_compute_drho_fails(): obs = pe.input.json.load_json("tests/data/compute_drho_fails.json.gz") obs.gm() assert np.isclose(obs.dvalue, 0.0022150779611891094) + + +def test_vec_gm(): + obs = pe.misc.gen_correlated_data(np.arange(3), np.array([[0.0364 , 0.03627262, 0.03615699], + [0.03627262, 0.03688438, 0.03674798], + [0.03615699, 0.03674798, 0.03732882]]), "qq", 3.8, 1000) + pe.gm(obs[0], S=0) + assert obs[0].S["qq"] == 0 + pe.gm(obs, S=1.3) + assert np.all(np.vectorize(lambda x: x.S["qq"])(obs) == 1.3) + aa = np.array([obs, obs, obs]) + pe.gamma_method(aa, S=2.2) + assert np.all(np.vectorize(lambda x: x.S["qq"])(aa) == 2.20) + cc = pe.Corr(obs) + pe.gm(cc, S=4.12) + assert np.all(np.vectorize(lambda x: x.S["qq"])(cc.content) == 4.12)