Irregularly Sampled Time Series Interpolation for Binary Evolution Simulations Using Dynamic Time Warping

Ugur Demir; Philipp M. Srivastava; Aggelos Katsaggelos; Vicky Kalogera; Santiago L. Tapia; Manuel Ballester; Shamal Lalvani; Patrick Koller

astro-ph.SRApr 15, 2026

Irregularly Sampled Time Series Interpolation for Binary Evolution Simulations Using Dynamic Time Warping

Authors: Ugur Demir, Philipp M. Srivastava, Aggelos Katsaggelos, Vicky Kalogera, Santiago L. Tapia, Manuel Ballester, Shamal Lalvani, Patrick Koller,

Organizations: Electrical and Computer Engineering, Northwestern University, 2145 Sheridan Road, Evanston, IL 60208, USA · Center for Interdisciplinary Exploration and Research in Astrophysics (CIERA), Northwestern University, 1800 Sherman Ave, Evanston, IL 60201, USA · NSF-Simons AI Institute for the Sky (SkAI),172 E. Chestnut St., Chicago, IL 60611, USA · Department of Physics and Astronomy, Northwestern University, 2145 Sheridan Road, Evanston, IL 60208, USA ·

Abstract

Binary stellar evolution simulations are computationally expensive. Stellar population synthesis relies on these detailed evolution models at a fundamental level. Producing thousands of such models requires hundreds of CPU hours, but stellar track interpolation provides one approach to significantly reduce this computational cost. Although single-star track interpolation is straightforward, stellar interactions in binary systems introduce significant complexity to binary evolution, making traditional single-track interpolation methods inapplicable. Binary tracks present fundamentally different challenges compared to single stars, which possess relatively straightforward evolutionary phases identifiable through distinct physical properties. Binary systems are complicated by mutual interactions that can dramatically alter evolutionary trajectories and introduce discontinuities difficult to capture through standard interpolation. In this work, we introduce a novel approach for track alignment and iterative track averaging based on Dynamic Time Warping to address misalignments between neighboring tracks. Our method computes a single shared warping path across all physical parameters simultaneously, placing them on a consistent temporal grid that preserves the causal relationships between parameters. We demonstrate that this joint-alignment strategy maintains key physical relationships such as the Stefan-Boltzmann law in the interpolated tracks. Our comprehensive evaluation across multiple binary configurations demonstrates that proper temporal alignment is crucial for track interpolation methods. The proposed method consistently outperforms existing approaches and enables the efficient generation of more accurate binary population samples for astrophysical studies.

Figures & tables

Figure 1: (a) A slice from the 3D HMS-HMS grid, displaying primary star mass (

M_{1}

) on the x-axis and orbital period (

P

) on the y-axis. Both axes are in logarithmic space, resulting in uniform point separation. The secondary star mass dimension (

M_{2}

) is fixed at

M_{2}/M_{1}=0.7

for visualization clarity. Each colored point corresponds to a simulated initial value configuration, with color coding indicating the evolutionary mass transfer class. Uncolored regions represent parameter combinations that were not simulated and are excluded from the dataset. (b), (c) and (d) illustrate example evolutionary tracks for No Mass Transfer, Unstable Mass Transfer, and Stable Mass Transfer classes, respectively.

Figure 2: Overview of the proposed DTW-based iterative interpolation method. For a given target initial condition

i_{*}

shown in (a), we identify the

K

nearest neighbors (here

K=4

) from the grid

G

i_{n_{1}}

i_{n_{2}}

i_{n_{3}}

, and

i_{n_{4}}

. The grid shown in (a) is plotted without a

\log

transformation. We retrieve their corresponding evolutionary tracks

H_{n_{1}}

H_{n_{2}}

H_{n_{3}}

, and

H_{n_{4}}

shown in (b). The interpolation proceeds through iterative alignment and averaging. First, tracks

H_{n_{1}}

and

H_{n_{2}}

are aligned using DTW as shown in (c), then averaged to produce

H_{\text{avg}}^{(2)}

(d). Next,

H_{n_{3}}

is aligned with

H_{\text{avg}}^{(2)}

(e) and averaged to create

H_{\text{avg}}^{(3)}

(f). Finally,

H_{n_{4}}

is aligned with

H_{\text{avg}}^{(3)}

(g) and averaged to generate the final interpolated track

H_{\text{avg}}^{(4)}

depicted in (h). For visual clarity, panels (b)–(h) display of a single parameter, but the alignment is computed over the alignment parameter subset

\mathcal{A}

and the resulting warping path is applied to all parameters simultaneously.

Figure 3: Dynamic Time Warping alignment process. (Left) The pairwise distance matrix between each point of tracks

H_{A}

and

H_{B}

, with color-coding indicating distance values (lighter colors represent higher values). DTW finds the optimal warping path (red points) that minimizes the cumulative distance function (Equation 7 ). The shaded regions illustrate the global constraints that exclude the constrained region from search areas. The global constraint is used in Constrained DTW (CDTW) to prevent excessive warping. (Middle) The aligned tracks show corresponding points connected by dashed lines. (Right) The close-up panel demonstrates the precise alignment of signal features, particularly at edges where there is a significant drop in parameter values.

Grid	# Samples	# Samples (filtered)	# NMT	# SMT	# UMT	# Test Samples
HMS-HMS	60533	55736	15497	18712	21527	2833
CO-HeMS	39477	32868	23611	8095	1162	2296
CO-HMS	23684	18777	6025	8338	4414	2368

Table 1: Sample Distribution Across Binary Evolution Grids

Figure 4: Comparison of interpolation methods using Relative Absolute Error distributions. Violin plots show the distribution of RAE values for the Change Point, Nearest Neighbor, and our DTW-based iterative averaging methods across all binary evolution parameters. The y-axis is scaled logarithmically (

\log_{10}

) to visualize the full error range, while median values are presented as percentage Relative Absolute Error (shown in parentheses on the x-axis). The purple dashed line indicates the 10% error threshold, above which predictions are considered unacceptable for astrophysical applications. Our DTW-based method consistently achieves lower median errors across most of the parameters for all 3 grids.

Figure 5: Sample interpolation results for the

\dot{M}_{transfer}

parameter from the CO-HMS grid with initial conditions

M_{1}

= 4.31

M_{\odot}

M_{2}

= 1.58

M_{\odot}

and

P

= 129.4 days. The ground truth (green solid line) is compared against predictions from the Change Point (blue dashed) and DTW (orange dashed) algorithms. Despite the DTW prediction being visually closer to the ground truth, it exhibits larger relative absolute error due to slight temporal shifts. At the first spike, the absolute error

e_{DTW}

measures the distance from the spike peak to the interpolated value, while

e_{CP}

is smaller despite the Change Point prediction missing the spike morphology. At the second spike, both methods are severely penalized due to temporal misalignment: the error is calculated from signal peaks (near 0) to minima (around

-100

), despite close visual agreement. The dense clustering of ground truth sampling points during rapid changes results in inflated RAE values, motivating the area-based error metric that better captures track fidelity.

Figure 6: Error calculation methodologies for evaluating interpolation accuracy. (a) Relative Absolute Error (RAE) computation: RAE evaluates deviations at specific ground truth time points. Since interpolated tracks are irregularly sampled and may not contain values at exact ground truth timestamps, we create a continuous representation using linear interpolation between adjacent points. The interpolated prediction is then sampled at ground truth time points to calculate relative deviations. Padding is applied at boundaries when predictions are shorter than ground truth tracks. (b) Area Error calculation: Both predicted time points (blue dots) and ground truth time points (green dots) are linearly interpolated to treat them as continuous signals. Area error quantifies prediction deviations across the entire temporal domain by computing the relative absolute distance between tracks at infinitesimally small time steps

\Delta a

. The calculation is approximated by measuring the area between the predicted and ground truth tracks (blue shaded region) and normalizing by the area under the ground truth curve (green shaded region).

Figure 7: Comparison of interpolation methods using Area Error distributions. Violin plots display the area error distributions for the Nearest Neighbor approach, Change Point algorithm, and our DTW-based iterative averaging method across all binary evolution parameters and three grids (HMS-HMS, CO-HeMS, and CO-HMS). The y-axis is scaled logarithmically, with median values marked within each distribution. Our DTW-based method demonstrates superior performance across all three grids, achieving consistently lower median area errors.

Figure 8: Interpolation results comparing methods across different binary parameters. Initial conditions (

M_{1},M_{2},P

) for each sample are displayed in the panel titles. Each plot displays the neighboring tracks

H_{n_{i}}

(used for interpolation), predictions from the Change Point algorithm ( [ 32 ] ), predictions from our DTW-based method, and the corresponding ground truth evolutionary track. The Change Point algorithm shows significant deviations from the ground truth, producing morphologically inconsistent predictions despite the similar shapes of neighboring tracks. In contrast, our DTW-based alignment approach achieves near-perfect overlap with the ground truth. The visual results demonstrate that proper temporal alignment is crucial for accurate interpolation.

Figure 9: Example interpolation results for

\log_{10}(T_{\text{eff},1})

parameter with initial conditions

M_{1}

= 19.39

M_{\odot}

M_{2}

= 2.90

M_{\odot}

and

P

= 278.12 days. The close-up section highlights DTW-based approach capturing fine-scale variations in the signal, while the Change Point algorithm fails to preserve these critical details.

Method	Class-Aware	$M_{1}$	$M_{2}$	$P$	$\dot{M}_{transfer}$	$T_{eff,1}$	$T_{eff,2}$	$L_{1}$	$L_{2}$	$R_{1}$	$R_{2}$
Change Points (4)	✓	1.8122	0.3381	3.2128	33.675	1.4933	0.3972	1.1133	0.4695	9.716	4.036
Change Points (8)	✓	0.9217	0.2486	1.8003	22.1917	0.9608	0.2417	0.5851	0.2386	8.0558	2.0624
Change Points (16)	✓	0.596	0.2277	1.2938	23.2439	0.7345	0.1712	0.4104	0.1846	5.9181	1.5825
Change Points (32)	✓	0.5846	0.2327	1.2103	36.3708	0.6268	0.1464	0.3267	0.1192	4.8179	1.4481
Change Points (64)	✓	0.6882	0.2332	1.3281	39.5013	0.6803	0.1358	0.3285	0.1103	5.2436	1.3505
Change Points (128)	✓	0.7585	0.247	1.4474	38.4068	0.7273	0.1462	0.3645	0.134	6.0697	1.492

Table 2: Ablation Study Results on HMS-HMS Grid

Method	Neighbor Order	$M_{1}$	$M_{2}$	$P$	$\dot{M}_{transfer}$	$T_{eff,1}$	$T_{eff,2}$	$L_{1}$	$L_{2}$	$R_{1}$	$R_{2}$
DTW (0-1 norm)	Descending	0.4588	0.2016	1.0146	5.9789	0.3935	0.0898	0.1948	0.0616	3.2208	0.9075
DTW (0-1 norm)	Ascending	0.4554	0.2005	1.0168	5.7059	0.3913	0.09	0.1958	0.0613	3.2217	0.9033
DTW (0-1 norm)	Random	0.4615	0.2017	1.0177	5.7902	0.3959	0.0905	0.1958	0.0613	3.2371	0.912
DTW (max norm)	Descending	0.4742	0.2058	1.0458	6.4124	0.4206	0.0935	0.2116	0.0646	3.3949	0.9241
DTW (max norm)	Ascending	0.4906	0.2107	1.0516	6.8531	0.437	0.0971	0.2261	0.0662	3.4889	0.9666
DTW (max norm)	Random	0.4919	0.2091	1.0521	6.4917	0.4273	0.0981	0.2175	0.0651	3.4644	0.9659

Table 3: Effect of Neighbor Ordering on Iterative Weighted Averaging Performance on HMS-HMS Grid

Input	Alignment Method	$M_{1}$	$M_{2}$	$P$	$\dot{M}_{transfer}$	$T_{eff,1}$	$T_{eff,2}$	$L_{1}$	$L_{2}$	$R_{1}$	$R_{2}$
Change Points (4)	-	1.8122	0.3381	3.2128	33.675	1.4933	0.3972	1.1133	0.4695	9.716	4.036
Change Points (4)	DTW (0-1 norm)	1.7594	0.332	3.1172	30.914	1.3727	0.4062	1.1257	0.4635	8.8358	3.966
Change Points (4)	DTW (max norm)	1.6931	0.3318	3.1442	31.3186	1.3926	0.3744	1.0952	0.457	8.7349	3.9666
Change Points (4)	CDTW (radius:3)	1.7594	0.332	3.1172	30.914	1.3727	0.4062	1.1257	0.4635	8.8358	3.966
Change Points (4)	SDTW	1.7836	0.3339	3.2036	30.7565	1.4427	0.3835	1.0913	0.4654	8.9835	3.9371
Change Points (8)	-	0.9217	0.2486	1.8003	22.1917	0.9608	0.2417	0.5851	0.2386	8.0558	2.0624

Table 4: Performance Improvement of Change Point Algorithm When Enhanced with DTW alignment on HMS-HMS Grid

Figure 10: Hertzsprung-Russell diagram analysis for two representative HMS-HMS systems. Each row corresponds to a different binary system with initial conditions shown above the panels. The first two panels display the effective temperature and luminosity as functions of age for both the primary and secondary stars, comparing ground truth (solid) and DTW-based predictions (dashed). The third and fourth panels show the corresponding HR diagrams for the primary and secondary stars, respectively. The top row presents a well-predicted system where the interpolated tracks closely follow the ground truth in both the time domain and the HR diagram. The bottom row shows a more challenging case where deviations occur during rapid evolutionary phases, resulting in visible differences in the HR diagram trajectories.

Figure 11: Stefan-Boltzmann consistency analysis for a representative HMS-HMS system with initial conditions

M_{1}=11.12\,M_{\odot}

M_{2}=8.86\,M_{\odot}

, and

P=14.17

days. The top half shows results for the primary star and the bottom half for the secondary star. For each star, the upper row displays luminosity, radius, and effective temperature as functions of age, comparing ground truth (solid green), DTW-based predictions (dashed red), and Change Point predictions (dotted blue). Below these, the Stefan-Boltzmann ratio

\mathcal{R}_{\mathrm{SB}}=L/(4\pi R^{2}\sigma T_{\mathrm{eff}}^{4})

is shown for each method, with a zoomed-in view highlighting the rapid evolutionary phase. While the individual parameter predictions appear visually similar across all three methods, the Change Point approach produces spurious spikes in

\mathcal{R}_{\mathrm{SB}}

during rapid evolutionary phases due to independent per-parameter sampling. The DTW-based joint alignment maintains

\mathcal{R}_{\mathrm{SB}}\approx 1

throughout the evolution, consistent with the ground truth.

Dynamic Time Warping

Dynamic Time Warping aligns sequences with nonlinear timing differences to compare, classify, or transform time-series data.

See all 14 papers

Generative Refinement Learning for Temporal Interpolation in Scientific Machine Learning
Oct 27, 2026Pu Ren, Michael W. Mahoney, N. Benjamin ErichsonScientific ML Residual Learning
Lawrence Berkeley National Lab · International Computer Science Institute · University of California, Berkeley
ReLaTS: a Reinforcement Learning-based method for dynamically determining the coupling Time Step in multi-scale simulations of self-gravitating systems
Jun 18, 2026Veronica Saz Ulibarrena, Simon Portegies ZwartReinforcement Learning Scientific ML
Leiden Observatory, Leiden University, Einsteinweg 55, 2333 CC, Leiden, The Netherlands
Trajectory inference via Acceleration Matching
Aug 4, 2026Bartolo Dazzini, Giovanni Conforti, Alain Durmus +1Trajectory Generation Trajectory Inference
Department of Mathematics, University of Padova. · CMAP, ´Ecole polytechnique. · Institute for Foundations of Data Science, Yale University.

Oct 27, 2026J/K move · Enter open · S save

Generative Refinement Learning for Temporal Interpolation in Scientific Machine Learning

Pu Ren, Michael W. Mahoney, N. Benjamin Erichson

Lawrence Berkeley National Lab · International Computer Science Institute · University of California, Berkeley

Scientific ML Residual Learning Diffusion Models

Open paper AskPDF arXiv

Computationally expensive simulations and real-world observations often provide data at a lower temporal resolution than downstream tasks require. Averaging observations at fixed grid points can introduce errors, and motion estimates alone do not determine how the physical field values change. In this work, we study whether learning to correct a simple numerical interpolation improves the reconstruction of intermediate scientific fields. Our model, interpolative FLow EXpert (iFLEX), first estimates the missing field by taking a weighted average of the two observations at each grid point, giving more weight to the observation closer to the requested time. During training, subtracting this estimate from an observed intermediate field gives the residual that a conditional diffusion model learns to generate. Both observations condition the model, together with the target time, the interval between observations, and available physical parameters. To reconstruct an unobserved intermediate field, the generated residual is added to the weighted average. We evaluate iFLEX on Navier-Stokes Kraichnan turbulence, Shanghai radar precipitation, and sea surface temperature. Across these datasets, iFLEX achieves lower reconstruction error than the evaluated baselines, including optical-flow and diffusion models developed for video frame interpolation in computer vision. In a controlled turbulence comparison at a Reynolds number of 36,000, learning the residual reduces relative Frobenius norm error from about 13.7 for direct field generation to 10.5, with the diffusion backbone and training setup held fixed. Code is available at https://github.com/paulpuren/iflex.

Abstract

Dynamic Time Warping

Explore similar work

Generative Refinement Learning for Temporal Interpolation in Scientific Machine Learning

ReLaTS: a Reinforcement Learning-based method for dynamically determining the coupling Time Step in multi-scale simulations of self-gravitating systems

Trajectory inference via Acceleration Matching