Reconstruction

Reconstruction and forecasting from DMD results.

reconstruction

Reconstruction and forecasting from DMD results.

Two functions replace the seven in the old dmd_reconstruction module:

• reconstruct — rebuild the signal from selected modes/pairs, with optional stabilisation.
• forecast — like reconstruct but designed for extrapolation, with eigenvalue modification support.

Functions:

Name	Description
reconstruct	Rebuild the signal from a DMDResult using selected modes or pairs.
forecast	Extrapolate beyond the training window, optionally modifying eigenvalues for stabilisation or frequency changes.

reconstruct

reconstruct(result: DMDResult, times: ndarray | None = None, pairs: list[int] | None = None, mode_indices: list[int] | None = None, stable: bool = False) -> np.ndarray

Rebuild the signal from a DMD decomposition.

This single function replaces reconstruct_dmd, reconstruct_dmd_complex, reconstruct_specific_modes, and reconstruct_conjugate_pairs from the old API.

Parameters:

Name	Type	Description	Default
result	DMDResult	Output of run_dmd.	required
times	ndarray	Time vector. Defaults to result.times[1:] (the training window minus the first frame, matching BOPDMD convention).	None
pairs	list of int	Conjugate-pair numbers (0-indexed) to include. Converted to mode indices internally. Mutually exclusive with mode_indices.	None
mode_indices	list of int	Raw mode indices to include. Mutually exclusive with pairs.	None
stable	bool	If True, use phasor reconstruction (frequencies only, no growth/decay). If False (default), use full complex eigenvalues.	False

Returns:

Type	Description
ndarray	Reconstructed keypoints, shape (n_times, n_markers, 3), with the mean shape added back.

Notes

With stable=False, uses the full eigenvalue reconstruction:

x̃(t) = Re[ Σⱼ bⱼ φⱼ exp(λⱼ t) ]

With stable=True, strips growth/decay for bounded output:

x̃(t) = Σ_{pairs} 2|b| |φ| cos(|ω|t - θ)

Raises:

Type	Description
ValueError	If both pairs and mode_indices are supplied.

Source code in src/birddmd/reconstruction.py

def reconstruct(
    result: DMDResult,
    times: np.ndarray | None = None,
    pairs: list[int] | None = None,
    mode_indices: list[int] | None = None,
    stable: bool = False,
) -> np.ndarray:
    """Rebuild the signal from a DMD decomposition.

    This single function replaces ``reconstruct_dmd``,
    ``reconstruct_dmd_complex``, ``reconstruct_specific_modes``, and
    ``reconstruct_conjugate_pairs`` from the old API.

    Parameters
    ----------
    result : DMDResult
        Output of `run_dmd`.
    times : np.ndarray, optional
        Time vector.  Defaults to ``result.times[1:]`` (the training
        window minus the first frame, matching BOPDMD convention).
    pairs : list of int, optional
        Conjugate-pair numbers (0-indexed) to include.  Converted to
        mode indices internally.  Mutually exclusive with
        *mode_indices*.
    mode_indices : list of int, optional
        Raw mode indices to include.  Mutually exclusive with *pairs*.
    stable : bool
        If ``True``, use phasor reconstruction (frequencies only,
        no growth/decay).  If ``False`` (default), use full complex
        eigenvalues.

    Returns
    -------
    np.ndarray
        Reconstructed keypoints, shape ``(n_times, n_markers, 3)``,
        with the mean shape added back.

    Notes
    -----
    With ``stable=False``, uses the full eigenvalue reconstruction:

        x̃(t) = Re[ Σⱼ bⱼ φⱼ exp(λⱼ t) ]

    With ``stable=True``, strips growth/decay for bounded output:

        x̃(t) = Σ_{pairs} 2|b| |φ| cos(|ω|t - θ)

    Raises
    ------
    ValueError
        If both *pairs* and *mode_indices* are supplied.
    """
    if pairs is not None and mode_indices is not None:
        msg = "Specify *pairs* or *mode_indices*, not both"
        raise ValueError(msg)

    if times is None:
        times = result.times[1:]

    dmd = result._dmd_object
    # Ensure we only use the coordinate rows, not extra Hankel rows
    n_markers = result.mode_magnitudes.shape[0]
    n_actual_coords = n_markers * 3
    modes_trimmed = dmd.modes[:n_actual_coords, :]

    # Resolve pairs → mode_indices
    if pairs is not None:
        mode_indices = result.pair_indices(pairs)

    if stable:
        recon = _reconstruct_phasor(
            times,
            dmd.eigs,
            modes_trimmed,
            dmd.amplitudes,
            result.conjugate_pairs,
            mode_indices,
        )
    else:
        recon = _reconstruct_complex(
            times,
            dmd.eigs,
            modes_trimmed,
            dmd.amplitudes,
            mode_indices,
        )

    # Reshape and add mean shape
    avg = result.average_shape.reshape(-1, n_markers, 3)
    return recon.reshape(-1, n_markers, 3) + avg

forecast

forecast(result: DMDResult, times: ndarray, pairs: list[int] | None = None, mode_indices: list[int] | None = None, stable: bool = True, modify_eigenvalues: dict | None = None) -> np.ndarray

Extrapolate beyond the training window.

Like reconstruct but with additional controls for eigenvalue modification (e.g. changing frequencies, zeroing specific modes).

Parameters:

Name	Type	Description	Default
result	DMDResult	Output of run_dmd.	required
times	ndarray	Forecast time vector (may extend beyond training window).	required
pairs	list of int	Conjugate-pair numbers to include.	None
mode_indices	list of int	Raw mode indices to include.	None
stable	bool	If True (default), strip growth/decay.	True
modify_eigenvalues	dict	Eigenvalue modifications. Keys: "zero_modes" — list of mode indices whose eigenvalues are set to zero (and their conjugates). "changed_freq" — new frequency (rad/s) to assign to the modes in "zero_modes" instead of zero.	None

Returns:

Type	Description
ndarray	Forecast keypoints, shape (n_times, n_markers, 3).

Notes

This replaces run_forecast, run_forecast_with_modified_modes, and modify_mode_frequencies from the old API.

Source code in src/birddmd/reconstruction.py

def forecast(
    result: DMDResult,
    times: np.ndarray,
    pairs: list[int] | None = None,
    mode_indices: list[int] | None = None,
    stable: bool = True,
    modify_eigenvalues: dict | None = None,
) -> np.ndarray:
    """Extrapolate beyond the training window.

    Like `reconstruct` but with additional controls for
    eigenvalue modification (e.g. changing frequencies, zeroing
    specific modes).

    Parameters
    ----------
    result : DMDResult
        Output of `run_dmd`.
    times : np.ndarray
        Forecast time vector (may extend beyond training window).
    pairs : list of int, optional
        Conjugate-pair numbers to include.
    mode_indices : list of int, optional
        Raw mode indices to include.
    stable : bool
        If ``True`` (default), strip growth/decay.
    modify_eigenvalues : dict, optional
        Eigenvalue modifications.  Keys:

        * ``"zero_modes"`` — list of mode indices whose eigenvalues
          are set to zero (and their conjugates).
        * ``"changed_freq"`` — new frequency (rad/s) to assign to
          the modes in ``"zero_modes"`` instead of zero.

    Returns
    -------
    np.ndarray
        Forecast keypoints, shape ``(n_times, n_markers, 3)``.

    Notes
    -----
    This replaces ``run_forecast``, ``run_forecast_with_modified_modes``,
    and ``modify_mode_frequencies`` from the old API.
    """
    if pairs is not None and mode_indices is not None:
        msg = "Specify *pairs* or *mode_indices*, not both"
        raise ValueError(msg)

    dmd = result._dmd_object
    n_markers = result.mode_magnitudes.shape[0]
    n_coords = n_markers * 3
    modes = dmd.modes[:n_coords, :]
    amplitudes = dmd.amplitudes.copy()
    eigenvalues = dmd.eigs.copy()

    # Apply eigenvalue modifications
    modified_indices: set = set()
    if modify_eigenvalues is not None:
        modified_indices = _apply_eigenvalue_modifications(
            eigenvalues,
            modify_eigenvalues,
            result.conjugate_pairs,
            stable,
        )

    # Resolve pairs → mode_indices
    if pairs is not None:
        mode_indices = result.pair_indices(pairs)

    # Strip growth/decay if stable
    if stable and not modified_indices:
        eigenvalues = 1j * np.imag(eigenvalues)
    elif stable and modified_indices:
        # Strip real parts only from unmodified eigenvalues
        for k in range(len(eigenvalues)):
            if k not in modified_indices:
                eigenvalues[k] = 1j * np.imag(eigenvalues[k])

    # Subset modes if requested
    if mode_indices is not None:
        eigenvalues = eigenvalues[mode_indices]
        modes = modes[:, mode_indices]
        amplitudes = amplitudes[mode_indices]

    # Reconstruct
    dynamics = np.exp(np.outer(eigenvalues, times))
    weighted = modes * amplitudes[np.newaxis, :]
    recon = np.real((weighted @ dynamics).T)

    avg = result.average_shape.reshape(-1, n_markers, 3)
    return recon.reshape(-1, n_markers, 3) + avg