Coverage for python/gsfit/tensioned_cubic_splines_regularisations.py: 11%
55 statements
« prev ^ index » next coverage.py v7.16.2, created at 2026-10-14 15:25 +0000
1import typing
3import gsfit_rs
4import numpy as np
5import numpy.typing as npt
8def make_tensioned_cubic_b_spline_regularisations(
9 interior_knots: npt.NDArray[np.float64],
10 interval_tensions: npt.NDArray[np.float64],
11 regularisation_builder: dict[str, typing.Any],
12) -> npt.NDArray[np.float64]:
13 """
14 Build the regularisation matrix for a tensioned cubic B-spline source function.
16 This gathers the boiler-plate that used to be copied into the examples (see the previous
17 `examples/scan_gsfit__st40__real_data_pressure_constrained.py` and
18 `examples/example_07__mastu__freegsnke_data__tensioned_cubic_b_spline_test.ipynb`) into a single
19 reusable function which is driven purely by settings.
21 :param interior_knots: Interior knot locations in normalised poloidal flux [dimensionless].
22 The boundary knots at psi_n = 0 and psi_n = 1 are added automatically and must not be included here.
23 :param interval_tensions: Tension over each interval between knots [dimensionless].
24 :param regularisation_builder: Dictionary describing how to build the regularisation matrix. Keys:
25 * `regularisation_scale` [dimensionless]: Multiplies the "bulk" regularisation rows.
26 * `regularisation_mode` [str]: How to regularise the bulk of the profile. Presently only
27 `"second_derivative_zero"` is supported, which drives the second derivative of the source
28 function to zero at every knot.
29 * `left_boundary_condition` / `right_boundary_condition` [str]: One of `"free"`, `"dirichlet"`,
30 or `"neumann"`. `"free"` adds no boundary row, `"dirichlet"` drives the source function value to
31 zero at the boundary, and `"neumann"` drives the first derivative of the source function to zero
32 at the boundary.
33 * `left_boundary_weight` / `right_boundary_weight` [dimensionless]: Weight applied to the boundary
34 rows (only used when the corresponding boundary condition is not `"free"`).
36 :return: The regularisation matrix, shape = (n_regularisation, n_dof).
38 The analytic formulae for the value, first derivative, and second derivative of a tensioned cubic
39 B-spline at a knot are derived in the following LaTeX document:
40 https://www.overleaf.com/read/vbydpqkjtmds#910342
41 which builds on results from P. E. Koch and T. Lyche, "Interpolation with Exponential B-Splines in
42 Tension", in *Geometric Modelling*, pp. 173-190.
43 """
45 # Ensure numpy arrays (the values may arrive as Python lists from the JSON settings)
46 interior_knots = np.asarray(interior_knots, dtype=np.float64)
47 interval_tensions = np.asarray(interval_tensions, dtype=np.float64)
49 # Number of degrees of freedom for a cubic spline = number of interior knots plus 4
50 n_dof: int = len(interior_knots) + 4
52 # A dummy source function is used purely to access the `sigma1_array` and `sigma2_array` values,
53 # which are needed to evaluate the analytic derivative formulae below.
54 regularisations_dummy: npt.NDArray[np.float64] = np.zeros((1, n_dof), dtype=np.float64)
55 source_function_dummy = gsfit_rs.TensionedCubicBSpline(regularisations_dummy, interior_knots, interval_tensions)
56 sigma1_array: npt.NDArray[np.float64] = source_function_dummy.get_array1(["sigma1_array"])
57 sigma2_array: npt.NDArray[np.float64] = source_function_dummy.get_array1(["sigma2_array"])
59 # All knots, including the boundary knots at psi_n = 0 and psi_n = 1
60 knots: npt.NDArray[np.float64] = np.concatenate(([0.0], interior_knots, [1.0]))
61 n_knots: int = len(knots)
63 regularisation_scale: float = regularisation_builder["regularisation_scale"]
64 regularisation_mode: str = regularisation_builder["regularisation_mode"]
66 # "Bulk" regularisation rows, one per knot (including the boundary knots)
67 bulk_regularisations: npt.NDArray[np.float64] = np.zeros((n_knots, n_dof), dtype=np.float64)
68 if regularisation_mode == "second_derivative_zero":
69 # Impose d^2(source_function)/d(psi_n)^2 = 0 at each knot location
70 for i_knot in range(n_knots):
71 j: int = i_knot + 3
72 bulk_regularisations[i_knot, j - 3] = 1.0 / (sigma2_array[j - 2] * sigma1_array[j - 1])
73 bulk_regularisations[i_knot, j - 2] = -(1.0 / sigma2_array[j - 2] + 1.0 / sigma2_array[j - 1]) / sigma1_array[j - 1]
74 bulk_regularisations[i_knot, j - 1] = 1.0 / (sigma2_array[j - 1] * sigma1_array[j - 1])
75 else:
76 raise ValueError(f"Unknown regularisation_mode: '{regularisation_mode}'. Supported values are: 'second_derivative_zero'.")
78 bulk_regularisations *= regularisation_scale
80 # Boundary rows are stacked before and after the bulk rows
81 left_boundary_regularisation: npt.NDArray[np.float64] | None = _make_boundary_regularisation(
82 side="left",
83 boundary_condition=regularisation_builder["left_boundary_condition"],
84 boundary_weight=regularisation_builder["left_boundary_weight"],
85 n_dof=n_dof,
86 )
87 right_boundary_regularisation: npt.NDArray[np.float64] | None = _make_boundary_regularisation(
88 side="right",
89 boundary_condition=regularisation_builder["right_boundary_condition"],
90 boundary_weight=regularisation_builder["right_boundary_weight"],
91 n_dof=n_dof,
92 )
94 regularisations_list: list[npt.NDArray[np.float64]] = []
95 if left_boundary_regularisation is not None:
96 regularisations_list.append(left_boundary_regularisation)
97 regularisations_list.append(bulk_regularisations)
98 if right_boundary_regularisation is not None:
99 regularisations_list.append(right_boundary_regularisation)
101 regularisations: npt.NDArray[np.float64] = np.concatenate(regularisations_list, axis=0)
103 return regularisations
106def _make_boundary_regularisation(
107 side: str,
108 boundary_condition: str,
109 boundary_weight: float,
110 n_dof: int,
111) -> npt.NDArray[np.float64] | None:
112 """
113 Build a single boundary regularisation row.
115 :param side: `"left"` (psi_n = 0) or `"right"` (psi_n = 1).
116 :param boundary_condition: One of `"free"`, `"dirichlet"`, or `"neumann"`.
117 :param boundary_weight: Weight applied to the boundary row [dimensionless].
118 :param n_dof: Number of degrees of freedom (columns of the regularisation matrix).
120 :return: A row of shape (1, n_dof), or `None` when the boundary condition is `"free"`.
122 At a clamped boundary only a single B-spline is non-zero, so the source function value equals the
123 first (psi_n = 0) or last (psi_n = 1) degree of freedom. The first derivative at the boundary is
124 proportional to the difference of the two degrees of freedom nearest the boundary (see the LaTeX
125 document referenced in `make_tensioned_cubic_b_spline_regularisations`).
126 """
128 if boundary_condition == "free":
129 return None
131 boundary_regularisation: npt.NDArray[np.float64] = np.zeros((1, n_dof), dtype=np.float64)
133 if side == "left":
134 if boundary_condition == "dirichlet":
135 # Impose source_function = 0 at psi_n = 0
136 boundary_regularisation[0, 0] = boundary_weight
137 elif boundary_condition == "neumann":
138 # Impose d(source_function)/d(psi_n) = 0 at psi_n = 0
139 boundary_regularisation[0, 0] = boundary_weight
140 boundary_regularisation[0, 1] = -boundary_weight
141 else:
142 raise ValueError(f"Unknown boundary_condition: '{boundary_condition}'. Supported values are: 'free', 'dirichlet', 'neumann'.")
143 elif side == "right":
144 if boundary_condition == "dirichlet":
145 # Impose source_function = 0 at psi_n = 1
146 boundary_regularisation[0, -1] = boundary_weight
147 elif boundary_condition == "neumann":
148 # Impose d(source_function)/d(psi_n) = 0 at psi_n = 1
149 boundary_regularisation[0, -1] = boundary_weight
150 boundary_regularisation[0, -2] = -boundary_weight
151 else:
152 raise ValueError(f"Unknown boundary_condition: '{boundary_condition}'. Supported values are: 'free', 'dirichlet', 'neumann'.")
153 else:
154 raise ValueError(f"Unknown side: '{side}'. Supported values are: 'left', 'right'.")
156 return boundary_regularisation