From 3e9c499a5946fe2df44f40ed6787816bcef8c06d Mon Sep 17 00:00:00 2001 From: Christina Hedges <14965634+christinahedges@users.noreply.github.com> Date: Thu, 8 Dec 2022 14:47:44 -0500 Subject: [PATCH 1/6] Update spline to be slightly more robust Extracted out of #54, this is just a slightly more robust version of this part of the code. --- src/psfmachine/utils.py | 44 +++++++++++++++++++++++++++++++++-------- 1 file changed, 36 insertions(+), 8 deletions(-) diff --git a/src/psfmachine/utils.py b/src/psfmachine/utils.py index 38d2c5a..660a6dc 100644 --- a/src/psfmachine/utils.py +++ b/src/psfmachine/utils.py @@ -199,15 +199,43 @@ def spline1d(x, knots, degree=3): return X -def _make_A_cartesian(x, y, n_knots=10, radius=3.0, knot_spacing_type="sqrt", degree=3): - if knot_spacing_type == "sqrt": - knots = np.linspace(-np.sqrt(radius), np.sqrt(radius), n_knots) - knots = np.sign(knots) * knots ** 2 +def _make_A_cartesian(x, y, n_knots=10, radius=3.0, spacing="sqrt", degree=3): + # Must be odd + n_odd_knots = n_knots if n_knots % 2 == 1 else n_knots + 1 + if spacing == "sqrt": + x_knots = np.linspace(-np.sqrt(radius), np.sqrt(radius), n_odd_knots) + x_knots = np.sign(x_knots) * x_knots ** 2 + y_knots = np.linspace(-np.sqrt(radius), np.sqrt(radius), n_odd_knots) + y_knots = np.sign(y_knots) * y_knots ** 2 else: - knots = np.linspace(-radius, radius, n_knots) - x_spline = spline1d(x, knots=knots, degree=degree) - y_spline = spline1d(y, knots=knots, degree=degree) - + x_knots = np.linspace(-radius, radius, n_odd_knots) + y_knots = np.linspace(-radius, radius, n_odd_knots) + x_spline = sparse.csr_matrix( + np.asarray( + dmatrix( + "bs(x, knots=knots, degree=degree, include_intercept=True)", + { + "x": list(np.hstack([x_knots.min(), x, x_knots.max()])), + "degree": degree, + "knots": x_knots, + }, + ) + )[1:-1] + ) + y_spline = sparse.csr_matrix( + np.asarray( + dmatrix( + "bs(x, knots=knots, degree=degree, include_intercept=True)", + { + "x": list(np.hstack([y_knots.min(), y, y_knots.max()])), + "degree": degree, + "knots": y_knots, + }, + ) + )[1:-1] + ) + x_spline = x_spline[:, np.asarray(x_spline.sum(axis=0))[0] != 0] + y_spline = y_spline[:, np.asarray(y_spline.sum(axis=0))[0] != 0] X = sparse.hstack( [x_spline.multiply(y_spline[:, idx]) for idx in range(y_spline.shape[1])], format="csr", From c9416749f3c087b8a21b72f16ef5fb8d2ca39439 Mon Sep 17 00:00:00 2001 From: Jorge Date: Mon, 12 Dec 2022 09:57:34 -0800 Subject: [PATCH 2/6] kwargs --- src/psfmachine/utils.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/psfmachine/utils.py b/src/psfmachine/utils.py index 660a6dc..0b50cc9 100644 --- a/src/psfmachine/utils.py +++ b/src/psfmachine/utils.py @@ -199,10 +199,10 @@ def spline1d(x, knots, degree=3): return X -def _make_A_cartesian(x, y, n_knots=10, radius=3.0, spacing="sqrt", degree=3): +def _make_A_cartesian(x, y, n_knots=10, radius=3.0, knot_spacing_type="sqrt", degree=3): # Must be odd n_odd_knots = n_knots if n_knots % 2 == 1 else n_knots + 1 - if spacing == "sqrt": + if knot_spacing_type == "sqrt": x_knots = np.linspace(-np.sqrt(radius), np.sqrt(radius), n_odd_knots) x_knots = np.sign(x_knots) * x_knots ** 2 y_knots = np.linspace(-np.sqrt(radius), np.sqrt(radius), n_odd_knots) From c2bfb7bb236f9eb51f9b07cb48a21e47957484ad Mon Sep 17 00:00:00 2001 From: Jorge Date: Mon, 12 Dec 2022 10:43:58 -0800 Subject: [PATCH 3/6] odd n knots in cartesian A --- tests/test_perturbation.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/tests/test_perturbation.py b/tests/test_perturbation.py index 532fcea..db1d95d 100644 --- a/tests/test_perturbation.py +++ b/tests/test_perturbation.py @@ -93,7 +93,7 @@ def test_perturbation_matrix3d(): p3 = PerturbationMatrix3D( time=time, dx=dx, dy=dy, nknots=4, radius=5, resolution=5, poly_order=1 ) - assert p3.cartesian_matrix.shape == (169, 64) + assert p3.cartesian_matrix.shape == (169, 81) assert p3.vectors.shape == (10, 2) assert p3.shape == ( p3.cartesian_matrix.shape[0] * p3.ntime, From 7b1f69d58ff64fd4b5bac52d533f4c09d892432f Mon Sep 17 00:00:00 2001 From: Jorge Date: Mon, 12 Dec 2022 13:29:16 -0800 Subject: [PATCH 4/6] limiting numpy v --- pyproject.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/pyproject.toml b/pyproject.toml index ed2dfac..27f23e5 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -20,7 +20,7 @@ source = ["src"] [tool.poetry.dependencies] python = "^3.7.1" -numpy = "^1.19.4" +numpy = ">=1.19.4, < 1.23" scipy = "^1.5.4" astropy = "^4.2" tqdm = "^4.54.0" From aa7b0ab372b0bde5922ca0253d32b866dc102c09 Mon Sep 17 00:00:00 2001 From: Jorge Date: Mon, 12 Dec 2022 13:42:05 -0800 Subject: [PATCH 5/6] undo numpy v limit --- pyproject.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/pyproject.toml b/pyproject.toml index 27f23e5..d6d2cc9 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -20,7 +20,7 @@ source = ["src"] [tool.poetry.dependencies] python = "^3.7.1" -numpy = ">=1.19.4, < 1.23" +numpy = ">=1.19.4" scipy = "^1.5.4" astropy = "^4.2" tqdm = "^4.54.0" From 379f9e445fb174dc6d032235536f0bb455457ac9 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Jorge=20Mart=C3=ADnez-Palomera?= Date: Thu, 8 Jun 2023 13:28:42 -0700 Subject: [PATCH 6/6] simplify version of robust cartesian DM and knots --- src/psfmachine/utils.py | 35 ++++++++++------------------------- 1 file changed, 10 insertions(+), 25 deletions(-) diff --git a/src/psfmachine/utils.py b/src/psfmachine/utils.py index 0b50cc9..0901b96 100644 --- a/src/psfmachine/utils.py +++ b/src/psfmachine/utils.py @@ -183,8 +183,10 @@ def _make_A_polar(phi, r, cut_r=6, rmin=1, rmax=18, n_r_knots=12, n_phi_knots=15 return X1 -def spline1d(x, knots, degree=3): +def spline1d(x, knots, degree=3, include_knots=False): """Make a spline in a 1D variable `x`""" + if include_knots: + x = np.hstack([knots.min(), x, knots.max()]) X = sparse.csr_matrix( np.asarray( dmatrix( @@ -193,6 +195,9 @@ def spline1d(x, knots, degree=3): ) ) ) + if include_knots: + X = X[1:-1] + x = x[1:-1] if not X.shape[0] == x.shape[0]: raise ValueError("`patsy` has made the wrong matrix.") X = X[:, np.asarray(X.sum(axis=0) != 0)[0]] @@ -210,30 +215,10 @@ def _make_A_cartesian(x, y, n_knots=10, radius=3.0, knot_spacing_type="sqrt", de else: x_knots = np.linspace(-radius, radius, n_odd_knots) y_knots = np.linspace(-radius, radius, n_odd_knots) - x_spline = sparse.csr_matrix( - np.asarray( - dmatrix( - "bs(x, knots=knots, degree=degree, include_intercept=True)", - { - "x": list(np.hstack([x_knots.min(), x, x_knots.max()])), - "degree": degree, - "knots": x_knots, - }, - ) - )[1:-1] - ) - y_spline = sparse.csr_matrix( - np.asarray( - dmatrix( - "bs(x, knots=knots, degree=degree, include_intercept=True)", - { - "x": list(np.hstack([y_knots.min(), y, y_knots.max()])), - "degree": degree, - "knots": y_knots, - }, - ) - )[1:-1] - ) + + x_spline = spline1d(x, knots=x_knots, degree=degree, include_knots=True) + y_spline = spline1d(y, knots=y_knots, degree=degree, include_knots=True) + x_spline = x_spline[:, np.asarray(x_spline.sum(axis=0))[0] != 0] y_spline = y_spline[:, np.asarray(y_spline.sum(axis=0))[0] != 0] X = sparse.hstack(