improve ar efficiency
3 changed files+326−292
python/ephys_compression_tests/algorithms/ans/__init__.pymodified+48−3View file
@@ -1,8 +1,53 @@
11 import numpy as np
22 import os
3-from .ar import encode_ar, encode_ar_lossy, decode_ar
3+from . import ar_numba
44 from ...types import Algorithm
55
6+
7+# Adapter functions to match the old ar.py API
8+def encode_ar(data: np.ndarray, order: int):
9+ """Encode using AR model - adapter for ar_numba."""
10+ coeffs, initial_points = ar_numba.fit_ar_model(data, k=order)
11+ residuals_full = ar_numba.compute_residuals(data, coeffs, initial_points)
12+ # Extract residuals excluding the initial points (first 'order' rows)
13+ residuals = residuals_full[order:, :]
14+ # Transpose initial_points to match old API: (order, channels)
15+ initial_values = initial_points.T
16+ return coeffs, residuals, initial_values
17+
18+
19+def encode_ar_lossy(data: np.ndarray, order: int, step: int):
20+ """Encode using AR model with lossy quantization - adapter for ar_numba."""
21+ # Fit the AR model
22+ coeffs, initial_points = ar_numba.fit_ar_model(data, k=order)
23+
24+ # Compute residuals with quantization
25+ residuals_full = ar_numba.compute_residuals_lossy(data, coeffs, initial_points, step=step)
26+
27+ # Extract residuals excluding the initial points (first 'order' rows)
28+ residuals = residuals_full[order:, :]
29+
30+ # Transpose initial_points to match old API: (order, channels)
31+ initial_values = initial_points.T
32+ return coeffs, residuals, initial_values
33+
34+
35+def decode_ar(coeffs: np.ndarray, residuals: np.ndarray, initial_values: np.ndarray):
36+ """Decode AR encoded data - adapter for ar_numba."""
37+ # Transpose initial_values from (order, channels) to (channels, order)
38+ initial_points = initial_values.T
39+
40+ # Create full residuals array including initial points
41+ order = coeffs.shape[1]
42+ n_residuals, n_channels = residuals.shape
43+ n_timepoints = n_residuals + order
44+
45+ residuals_full = np.zeros((n_timepoints, n_channels), dtype=np.int16)
46+ residuals_full[:order, :] = initial_points.T
47+ residuals_full[order:, :] = residuals
48+
49+ return ar_numba.reconstruct_from_residuals(residuals_full, coeffs, initial_points)
50+
651 SOURCE_FILE = "ans/__init__.py"
752
853
@@ -273,7 +318,7 @@ for a in algorithm_dicts_base:
273318 return reconstructed
274319 algorithm_dicts.append({
275320 "name": a["name"] + f"-ar{order}",
276- "version": "2",
321+ "version": a["version"] + f".1",
277322 "encode": encode0_ar,
278323 "decode": decode0_ar,
279324 "description": a["description"] + f" with auto-regressive prediction encoding of order {order}",
@@ -316,7 +361,7 @@ for ar_order in [2, 8]:
316361 return decode0_ar_lossy
317362 algorithm_dicts.append({
318363 "name": f"ans-ar{ar_order}-lossy-tol{tolerance}",
319- "version": "3",
364+ "version": "10",
320365 "encode": make_encode_ar_lossy(),
321366 "decode": make_decode_ar_lossy(),
322367 "description": f"ANS with lossy auto-regressive prediction encoding of order {ar_order} and tolerance {tolerance}",
python/ephys_compression_tests/algorithms/ans/ar.pydeleted+0−289View file
@@ -1,289 +0,0 @@
1-"""Auto-regressive model utilities for ANS compression."""
2-
3-import numpy as np
4-from typing import Tuple
5-from numba import njit
6-
7-
8-def _warmup_numba_functions():
9- """Warmup numba JIT compilation with small test data."""
10- print("Warming up numba functions for AR model...")
11- # Create small test data for 2D arrays (timepoints x channels)
12- test_data = np.array([[1, 2], [2, 3], [3, 4], [4, 5], [5, 6], [6, 7], [7, 8], [8, 9]], dtype=np.int16)
13- test_coeffs = np.array([[0.5, 0.3], [0.4, 0.2]], dtype=np.float32) # channels x order
14- test_residuals = np.array([[1, 2], [2, 3], [3, 4], [4, 5]], dtype=np.int16)
15- test_initial = np.array([[1, 2], [2, 3]], dtype=np.int16)
16- test_step = 2
17-
18- # Warmup each numba function
19- _create_design_matrix_channel(test_data[:, 0], 2)
20- _apply_ar_residuals_kernel_channel(test_data[:, 0], test_coeffs[0])
21- _apply_ar_residuals_lossy_kernel_channel(test_data[:, 0], test_coeffs[0], test_step)
22- _decode_ar_kernel_channel(test_coeffs[0], test_residuals[:, 0], test_initial[:, 0])
23-
24-
25-@njit
26-def _create_design_matrix_channel(data: np.ndarray, order: int) -> Tuple[np.ndarray, np.ndarray]:
27- """Numba-optimized design matrix creation for AR model (single channel)."""
28- n = len(data)
29- X_design = np.zeros((n - order, order))
30- y_target = data[order:]
31-
32- for i in range(n - order):
33- for j in range(order):
34- X_design[i, j] = data[i + order - j - 1]
35-
36- return X_design, y_target
37-
38-
39-def fit_ar_model(data: np.ndarray, order: int) -> np.ndarray:
40- """
41- Fit an autoregressive model of given order using least squares.
42-
43- Args:
44- data: Input data array (timepoints x channels)
45- order: AR model order
46-
47- Returns:
48- AR coefficients as numpy array (channels x order)
49- """
50- n_timepoints, n_channels = data.shape
51- if order >= n_timepoints:
52- raise ValueError(f"AR order {order} must be less than data length {n_timepoints}")
53-
54- # Fit AR model for each channel separately
55- coeffs = np.zeros((n_channels, order), dtype=np.float32)
56-
57- for ch in range(n_channels):
58- # Create design matrix using numba-optimized function
59- X_design, y_target = _create_design_matrix_channel(data[:, ch], order)
60-
61- # Use faster solve via normal equations: (X^T X) coeffs = X^T y
62- # This is faster than lstsq for overdetermined systems
63- XtX = X_design.T @ X_design
64- Xty = X_design.T @ y_target
65- coeffs[ch] = np.linalg.solve(XtX, Xty)
66-
67- return coeffs
68-
69-
70-@njit
71-def _apply_ar_residuals_kernel_channel(data: np.ndarray, coeffs: np.ndarray) -> np.ndarray:
72- """Numba-optimized kernel for computing AR residuals (single channel)."""
73- order = len(coeffs)
74- n = len(data)
75- residuals = np.empty(n - order, dtype=data.dtype)
76-
77- for i in range(order, n):
78- # Predict using previous 'order' samples
79- # Use float32 accumulation
80- prediction = np.float32(0.0)
81- for j in range(order):
82- prediction += coeffs[j] * np.float32(data[i - j - 1])
83-
84- # Round to nearest integer using numpy's round (banker's rounding)
85- prediction_int = np.int16(np.round(prediction))
86- residual = data[i] - prediction_int
87- residuals[i - order] = residual
88-
89- return residuals
90-
91-
92-@njit
93-def _apply_ar_residuals_lossy_kernel_channel(data: np.ndarray, coeffs: np.ndarray, step: int) -> np.ndarray:
94- """Numba-optimized kernel for computing AR residuals with lossy quantization (single channel)."""
95- order = len(coeffs)
96- n = len(data)
97- residuals = np.empty(n - order, dtype=data.dtype)
98-
99- # Pre-allocate reconstructed array for efficiency
100- reconstructed = np.empty(n, dtype=np.int16)
101- reconstructed[:order] = data[:order]
102-
103- # Convert step to float32 for consistent float arithmetic
104- step_f32 = np.float32(step)
105-
106- for i in range(order, n):
107- # Predict using previous 'order' samples from reconstructed data
108- # Use float32 accumulation
109- prediction = np.float32(0.0)
110- for j in range(order):
111- prediction += coeffs[j] * np.float32(reconstructed[i - j - 1])
112-
113- # Round to nearest integer
114- prediction_int = np.int16(np.round(prediction))
115-
116- # Compute residual from original data
117- residual = data[i] - prediction_int
118-
119- # Quantize residual to nearest multiple of step
120- quantized_residual = np.int16(np.round(np.float32(residual) / step_f32) * step_f32)
121- residuals[i - order] = quantized_residual
122-
123- # Reconstruct sample using quantized residual for future predictions
124- reconstructed[i] = prediction_int + quantized_residual
125-
126- return residuals
127-
128-
129-def apply_ar_residuals(data: np.ndarray, coeffs: np.ndarray) -> np.ndarray:
130- """
131- Apply AR model with given coefficients and return residuals.
132-
133- Args:
134- data: Input data array (timepoints x channels)
135- coeffs: AR coefficients (channels x order)
136-
137- Returns:
138- Residuals array (timepoints x channels)
139- """
140- # Ensure coeffs is float32
141- coeffs = np.array(coeffs, dtype=np.float32)
142-
143- n_timepoints, n_channels = data.shape
144- order = coeffs.shape[1]
145- residuals = np.zeros((n_timepoints - order, n_channels), dtype=data.dtype)
146-
147- for ch in range(n_channels):
148- residuals[:, ch] = _apply_ar_residuals_kernel_channel(data[:, ch], coeffs[ch])
149-
150- return residuals
151-
152-
153-def apply_ar_residuals_lossy(data: np.ndarray, coeffs: np.ndarray, step: int) -> np.ndarray:
154- """
155- Apply AR model with given coefficients and return lossy residuals.
156-
157- Args:
158- data: Input data array (timepoints x channels)
159- coeffs: AR coefficients (channels x order)
160- step: Quantization step size
161-
162- Returns:
163- Residuals array (timepoints x channels)
164- """
165- coeffs = np.array(coeffs, dtype=np.float32)
166-
167- n_timepoints, n_channels = data.shape
168- order = coeffs.shape[1]
169- residuals = np.zeros((n_timepoints - order, n_channels), dtype=data.dtype)
170-
171- for ch in range(n_channels):
172- residuals[:, ch] = _apply_ar_residuals_lossy_kernel_channel(data[:, ch], coeffs[ch], step)
173-
174- return residuals
175-
176-
177-def encode_ar(data: np.ndarray, order: int) -> Tuple[np.ndarray, np.ndarray, np.ndarray]:
178- """
179- Encode data using AR model - returns coefficients, residuals, and initial values.
180-
181- Args:
182- data: Input data array (timepoints x channels, int16)
183- order: AR model order
184-
185- Returns:
186- Tuple of (coefficients, residuals, initial_values)
187- - coefficients: channels x order (float32)
188- - residuals: (timepoints - order) x channels (int16)
189- - initial_values: order x channels (int16)
190- """
191- # Fit AR model
192- coeffs = fit_ar_model(data, order)
193-
194- # Convert coefficients to float32 to match what will be deserialized
195- coeffs = coeffs.astype(np.float32)
196-
197- # Compute residuals using float32 coefficients
198- residuals = apply_ar_residuals(data, coeffs)
199-
200- # Store initial values
201- initial_values = data[:order, :]
202-
203- return coeffs, residuals, initial_values
204-
205-
206-def encode_ar_lossy(data: np.ndarray, order: int, step: int) -> Tuple[np.ndarray, np.ndarray, np.ndarray]:
207- """
208- Encode data using AR model with lossy quantization.
209-
210- Args:
211- data: Input data array (timepoints x channels, int16)
212- order: AR model order
213- step: Quantization step size
214-
215- Returns:
216- Tuple of (coefficients, residuals, initial_values)
217- - coefficients: channels x order (float32)
218- - residuals: (timepoints - order) x channels (int16)
219- - initial_values: order x channels (int16)
220- """
221- # Fit AR model
222- coeffs = fit_ar_model(data, order)
223-
224- # Convert coefficients to float32 to match what will be deserialized
225- coeffs = coeffs.astype(np.float32)
226-
227- # Compute residuals using float32 coefficients
228- residuals = apply_ar_residuals_lossy(data, coeffs, step=step)
229-
230- # Store initial values
231- initial_values = data[:order, :]
232-
233- return coeffs, residuals, initial_values
234-
235-
236-@njit
237-def _decode_ar_kernel_channel(coeffs: np.ndarray, residuals: np.ndarray, initial_values: np.ndarray) -> np.ndarray:
238- """Numba-optimized kernel for AR decoding (single channel)."""
239- order = len(coeffs)
240- n = len(residuals) + order
241- reconstructed = np.empty(n, dtype=np.int16)
242- reconstructed[:order] = initial_values
243-
244- for i in range(order, n):
245- # Predict using AR model
246- # Use float32 accumulation
247- prediction = np.float32(0.0)
248- for j in range(order):
249- prediction += coeffs[j] * np.float32(reconstructed[i - j - 1])
250-
251- # Round to nearest integer using numpy's round (banker's rounding)
252- prediction_int = np.int16(np.round(prediction))
253-
254- # Add residual
255- reconstructed[i] = prediction_int + residuals[i - order]
256-
257- return reconstructed
258-
259-
260-def decode_ar(coeffs: np.ndarray, residuals: np.ndarray, initial_values: np.ndarray) -> np.ndarray:
261- """
262- Decode AR encoded data.
263-
264- Args:
265- coeffs: AR coefficients (channels x order, float32)
266- residuals: Residuals array ((timepoints - order) x channels)
267- initial_values: Initial values (order x channels)
268-
269- Returns:
270- Reconstructed data array (timepoints x channels)
271- """
272- # Ensure coeffs is float32
273- coeffs = np.array(coeffs, dtype=np.float32)
274-
275- n_channels = coeffs.shape[0]
276- order = coeffs.shape[1]
277- n_residuals = residuals.shape[0]
278- n_timepoints = n_residuals + order
279-
280- reconstructed = np.zeros((n_timepoints, n_channels), dtype=np.int16)
281-
282- for ch in range(n_channels):
283- reconstructed[:, ch] = _decode_ar_kernel_channel(coeffs[ch], residuals[:, ch], initial_values[:, ch])
284-
285- return reconstructed
286-
287-
288-# Warmup numba functions on module import
289-_warmup_numba_functions()
python/ephys_compression_tests/algorithms/ans/ar_numba.pyadded+278−0View file
@@ -0,0 +1,278 @@
1+"""
2+Numba-accelerated implementation of autoregressive model operations.
3+All operations work with int16 data.
4+"""
5+
6+import numpy as np
7+from numba import jit, prange
8+
9+
10+@jit(nopython=True, parallel=False, fastmath=True)
11+def _fit_ar_model_channel(channel_data: np.ndarray, k: int, subsample_factor: int,
12+ min_samples: int) -> np.ndarray:
13+ """
14+ Fit AR model for a single channel using least squares with subsampling.
15+
16+ Args:
17+ channel_data: 1D array for a single channel (int16)
18+ k: AR model order
19+ subsample_factor: Use every Nth sample for fitting
20+ min_samples: Minimum number of samples to use
21+
22+ Returns:
23+ coefficients: 1D array of k coefficients (float32)
24+ """
25+ n = len(channel_data)
26+
27+ # Determine subsampling stride
28+ max_samples = n - k
29+ stride = max(1, max_samples // min_samples, subsample_factor)
30+
31+ # Number of samples we'll actually use
32+ n_samples = (max_samples + stride - 1) // stride
33+
34+ # Build design matrix X and target vector y using subsampled data
35+ # X has shape (n_samples, k), y has shape (n_samples,)
36+ X = np.zeros((n_samples, k), dtype=np.float32)
37+ y = np.zeros(n_samples, dtype=np.float32)
38+
39+ sample_idx = 0
40+ for t in range(k, n, stride):
41+ if sample_idx >= n_samples:
42+ break
43+ for i in range(k):
44+ X[sample_idx, i] = np.float32(channel_data[t - k + i])
45+ y[sample_idx] = np.float32(channel_data[t])
46+ sample_idx += 1
47+
48+ # Truncate if needed
49+ actual_samples = sample_idx
50+ if actual_samples < n_samples:
51+ X = X[:actual_samples, :]
52+ y = y[:actual_samples]
53+
54+ # Solve normal equations: X.T @ X @ coef = X.T @ y
55+ XtX = np.zeros((k, k), dtype=np.float32)
56+ Xty = np.zeros(k, dtype=np.float32)
57+
58+ for i in range(k):
59+ for j in range(k):
60+ for t in range(actual_samples):
61+ XtX[i, j] += X[t, i] * X[t, j]
62+
63+ for i in range(k):
64+ for t in range(actual_samples):
65+ Xty[i] += X[t, i] * y[t]
66+
67+ # Solve the system using numpy's solver
68+ coefficients = np.linalg.solve(XtX, Xty)
69+
70+ return coefficients
71+
72+
73+def fit_ar_model(data: np.ndarray, k: int, subsample_factor: int = 100,
74+ min_samples: int = 1000) -> tuple[np.ndarray, np.ndarray]:
75+ """
76+ Fit an autoregressive model of order k to multi-channel time series data.
77+
78+ Args:
79+ data: 2D array of shape (timepoints, channels) with dtype int16
80+ k: Order of the autoregressive model
81+ subsample_factor: Use every Nth sample for fitting (default: 100)
82+ min_samples: Minimum number of samples to use for fitting (default: 1000)
83+
84+ Returns:
85+ coefficients: Array of shape (channels, k) with dtype float32
86+ initial_points: Array of shape (channels, k) with dtype int16 (first k samples per channel)
87+ """
88+ n_timepoints, n_channels = data.shape
89+
90+ if n_timepoints <= k:
91+ raise ValueError(f"Need at least {k+1} timepoints for AR({k}) model")
92+
93+ # Store initial k points for each channel
94+ initial_points = data[:k, :].T.copy() # (channels, k)
95+
96+ # Fit coefficients for each channel
97+ coefficients = np.zeros((n_channels, k), dtype=np.float32)
98+
99+ for ch in range(n_channels):
100+ coefficients[ch, :] = _fit_ar_model_channel(data[:, ch], k, subsample_factor, min_samples)
101+
102+ return coefficients, initial_points
103+
104+
105+@jit(nopython=True, parallel=True, fastmath=True)
106+def _compute_residuals_jit(data: np.ndarray, coefficients: np.ndarray,
107+ initial_points: np.ndarray, k: int) -> np.ndarray:
108+ """
109+ JIT-compiled residuals computation.
110+ """
111+ n_timepoints, n_channels = data.shape
112+ residuals = np.zeros((n_timepoints, n_channels), dtype=np.int16)
113+
114+ # First k points are copied as-is
115+ residuals[:k, :] = initial_points.T
116+
117+ # Compute residuals for each channel in parallel
118+ for ch in prange(n_channels):
119+ coef = coefficients[ch, :]
120+
121+ for t in range(k, n_timepoints):
122+ # Predict from previous k samples
123+ predicted = np.float32(0.0)
124+ for i in range(k):
125+ predicted += coef[i] * np.float32(data[t - k + i, ch])
126+
127+ # Residual = actual - predicted (rounded)
128+ residuals[t, ch] = data[t, ch] - np.int16(np.round(predicted))
129+
130+ return residuals
131+
132+
133+@jit(nopython=True, parallel=True, fastmath=True)
134+def _compute_residuals_lossy_jit(data: np.ndarray, coefficients: np.ndarray,
135+ initial_points: np.ndarray, k: int, step: int) -> np.ndarray:
136+ """
137+ JIT-compiled lossy residuals computation with quantization feedback.
138+ """
139+ n_timepoints, n_channels = data.shape
140+ residuals = np.zeros((n_timepoints, n_channels), dtype=np.int16)
141+ reconstructed = np.zeros((n_timepoints, n_channels), dtype=np.int16)
142+
143+ # First k points are copied as-is
144+ residuals[:k, :] = initial_points.T
145+ reconstructed[:k, :] = initial_points.T
146+
147+ step_f32 = np.float32(step)
148+
149+ # Compute residuals for each channel in parallel
150+ for ch in prange(n_channels):
151+ coef = coefficients[ch, :]
152+
153+ for t in range(k, n_timepoints):
154+ # Predict from previous k reconstructed samples
155+ predicted = np.float32(0.0)
156+ for i in range(k):
157+ predicted += coef[i] * np.float32(reconstructed[t - k + i, ch])
158+
159+ prediction_int = np.int16(np.round(predicted))
160+
161+ # Compute residual from original data
162+ residual = data[t, ch] - prediction_int
163+
164+ # Quantize residual to nearest multiple of step
165+ quantized_residual = np.int16(np.round(np.float32(residual) / step_f32) * step_f32)
166+ residuals[t, ch] = quantized_residual
167+
168+ # Reconstruct sample using quantized residual for future predictions
169+ reconstructed[t, ch] = prediction_int + quantized_residual
170+
171+ return residuals
172+
173+
174+def compute_residuals(data: np.ndarray, coefficients: np.ndarray,
175+ initial_points: np.ndarray) -> np.ndarray:
176+ """
177+ Compute residuals given data and AR model coefficients.
178+
179+ Args:
180+ data: 2D array of shape (timepoints, channels) with dtype int16
181+ coefficients: Array of shape (channels, k) with dtype float32
182+ initial_points: Array of shape (channels, k) with dtype int16
183+
184+ Returns:
185+ residuals: Array of shape (timepoints, channels) with dtype int16
186+ """
187+ k = coefficients.shape[1]
188+ return _compute_residuals_jit(data, coefficients, initial_points, k)
189+
190+
191+def compute_residuals_lossy(data: np.ndarray, coefficients: np.ndarray,
192+ initial_points: np.ndarray, step: int) -> np.ndarray:
193+ """
194+ Compute lossy residuals with quantization given data and AR model coefficients.
195+
196+ Args:
197+ data: 2D array of shape (timepoints, channels) with dtype int16
198+ coefficients: Array of shape (channels, k) with dtype float32
199+ initial_points: Array of shape (channels, k) with dtype int16
200+ step: Quantization step size
201+
202+ Returns:
203+ residuals: Array of shape (timepoints, channels) with dtype int16
204+ """
205+ k = coefficients.shape[1]
206+ return _compute_residuals_lossy_jit(data, coefficients, initial_points, k, step)
207+
208+
209+@jit(nopython=True, parallel=True, fastmath=True)
210+def _reconstruct_from_residuals_jit(residuals: np.ndarray, coefficients: np.ndarray,
211+ initial_points: np.ndarray, k: int) -> np.ndarray:
212+ """
213+ JIT-compiled reconstruction.
214+ """
215+ n_timepoints, n_channels = residuals.shape
216+ reconstructed = np.zeros((n_timepoints, n_channels), dtype=np.int16)
217+
218+ # First k points are copied from initial_points
219+ reconstructed[:k, :] = initial_points.T
220+
221+ # Reconstruct each channel in parallel
222+ for ch in prange(n_channels):
223+ coef = coefficients[ch, :]
224+
225+ for t in range(k, n_timepoints):
226+ # Predict from previous k reconstructed samples
227+ predicted = np.float32(0.0)
228+ for i in range(k):
229+ predicted += coef[i] * np.float32(reconstructed[t - k + i, ch])
230+
231+ # Reconstruct: actual = predicted (rounded) + residual
232+ reconstructed[t, ch] = np.int16(np.round(predicted)) + residuals[t, ch]
233+
234+ return reconstructed
235+
236+
237+def reconstruct_from_residuals(residuals: np.ndarray, coefficients: np.ndarray,
238+ initial_points: np.ndarray) -> np.ndarray:
239+ """
240+ Reconstruct original data from residuals and AR model coefficients.
241+
242+ Args:
243+ residuals: 2D array of shape (timepoints, channels) with dtype int16
244+ coefficients: Array of shape (channels, k) with dtype float32
245+ initial_points: Array of shape (channels, k) with dtype int16
246+
247+ Returns:
248+ reconstructed: Array of shape (timepoints, channels) with dtype int16
249+ """
250+ k = coefficients.shape[1]
251+ return _reconstruct_from_residuals_jit(residuals, coefficients, initial_points, k)
252+
253+
254+def warmup(n_channels: int = 10, k: int = 10):
255+ """
256+ Warm up Numba JIT compilation for all functions.
257+
258+ Args:
259+ n_channels: Number of channels for warmup data
260+ k: AR model order for warmup
261+ """
262+ print("Warming up JIT...", end="", flush=True)
263+ # Create small warmup data
264+ warmup_data = np.random.randint(-1000, 1000, size=(1000, n_channels), dtype=np.int16)
265+
266+ # Warm up fit_ar_model
267+ coefficients, initial_points = fit_ar_model(warmup_data, k)
268+
269+ # Warm up compute_residuals
270+ residuals = compute_residuals(warmup_data, coefficients, initial_points)
271+
272+ # Warm up compute_residuals_lossy
273+ _ = compute_residuals_lossy(warmup_data, coefficients, initial_points, step=2)
274+
275+ # Warm up reconstruct_from_residuals
276+ _ = reconstruct_from_residuals(residuals, coefficients, initial_points)
277+
278+ print(" done")