LTE model for tidal stations, example for Warnemunde (11)

11a.csv

11a.csv.json

{
  "Aliased": [
    0.422362756,
    0.368617497,
    0.255621397,
    0.209019747,
    0.322015847,
    0.155274488,
    0.262765007,
    0.375761106,
    0.05374526,
    0.225992198,
    0.172246939,
    0.488757205,
    0.112996099,
    0.00714361,
    0.10034691,
    0.04660165,
    0.481613596,
    2,
    1
  ],
  "AliasedAmp": [
    0.162168508,
    0.018437804,
    0.080427203,
    0.100193077,
    -0.13758768,
    -0.237395871,
    -0.177315945,
    -0.112108433,
    -0.090758587,
    -0.124693774,
    0.129941321,
    0.05887287,
    0.047920509,
    -0.02670324,
    0.079661958,
    -0.170565617,
    -0.013872906,
    0.031258073,
    0.073714152
  ],
  "AliasedPhase": [
    12.26665875,
    11.41281403,
    17.90437492,
    9.159760748,
    24.87352569,
    14.65428835,
    10.7419008,
    14.19058454,
    10.14360775,
    8.237883753,
    5.232416953,
    20.48980669,
    5.299072807,
    24.21263636,
    6.574752375,
    6.740362898,
    10.45183646,
    6.226353584,
    3.299100693
  ],
  "Hold": 0.001538141,
  "Imp_Amp": 36.63771512,
  "Imp_Stride": 12,
  "Initial": 0.0212,
  "LTE_Amp": 0.752732079,
  "LTE_Freq": 188.2119915,
  "LTE_Phase": -1.845171154,
  "Periods": [
    27.2122,
    27.3216,
    27.5545,
    13.63339513,
    13.69114014,
    13.6061,
    13.6608,
    13.71877789,
    6795.985773,
    1616.215172,
    2120.513853,
    13.77725,
    3232.430344,
    9.132931547,
    9.108450374,
    9.120674533,
    27.0926041
  ],
  "PeriodsAmp": [
    0.306411813,
    0.105333615,
    0.046157526,
    0.052010303,
    -0.159651577,
    0.116360596,
    0.168727625,
    0.182906578,
    0.041646376,
    0.170302166,
    0.156154628,
    0.066553408,
    0.025133677,
    0.005700373,
    0.048782882,
    0.006749513,
    -0.182128222
  ],
  "PeriodsPhase": [
    13.56978148,
    7.87504352,
    18.32973237,
    7.156452097,
    21.01885068,
    9.702972954,
    6.011334823,
    10.17573262,
    9.25752208,
    7.37083659,
    6.278543069,
    15.92705949,
    5.758181895,
    8.262567708,
    1.540785851,
    6.178219982,
    7.925674634
  ],
  "Year": 365.2520198,
  "final_state": {
    "D_prev": 0.04294
  },
  "metrics": {
    "mse": null,
    "pearson_p": null,
    "pearson_r": null,
    "var_squared_error": null
  },
  "model_params": {
    "Hold": 0.001538141,
    "Imp_Amp": 36.63771512,
    "Imp_Stride": 12,
    "LTE_Amp": 0.752732079,
    "LTE_Freq": 188.219915,
    "LTE_Phase": -1.845171154,
    "Periods": [
      27.2122
    ],
    "PeriodsAmp": [
      0.306
    ],
    "PeriodsPhase": [
      13.57
    ],
    "Year": 365.2520198
  },
  "model_type": "spreadsheet",
  "n_samples": 0,
  "notes": "Default .json file for timeseries_modeler.py. Update 'source_csv' to match your CSV filename. The script will overwrite/augment this file with actual n_samples, final_state and metrics after processing.",
  "source_csv": "data.csv"
}

timeseries_modeler.py

#!/usr/bin/env python3
"""
timeseries_modeler.py

Read a two-column CSV (time, value). Read a JSON file with the same filename
appended by ".json" (e.g. data.csv.json) into a data dictionary. Create a clone of the
time-series and run a loop that steps through timestamps, calling a user-
customizable model-step function to produce a modeled series. At the end the
script computes Pearson's correlation coefficient and the variance of the
squared errors (and also MSE), using library routines.

Usage:
  python timeseries_modeler.py data.csv [--out fitted.csv] [--params '{"key":value}'] [--model spreadsheet|example|clone] [--plot]
  python3 timeseries_modeler.py 11a.csv --plot  --model spreadsheet

Notes:
- The CSV must have at least two columns (time, value). Header is tolerated.
- Default model is a simple example; an optional "spreadsheet" model is
  included that follows the formulas described in your earlier message:
    - impulse every Imp_Stride rows
    - C_t = Imp_Amp * impulse_mask * sum_k PeriodsAmp[k] * sin(2*pi*time/Period_k + PeriodsPhase_k)
    - D_t = (1-Hold)*D_{t-1} + Hold*C_t
    - E_t = LTE_Amp * sin(D_t*LTE_Freq + LTE_Phase)
  Provide model parameters via --params as a JSON string (see defaults below).
- The implementation loops over timestamps (explicit for-loop) so the model can
  be stateful (sample-and-hold needs previous D).
- Pearson correlation computed using scipy.stats.pearsonr if available,
  otherwise a fallback via numpy.corrcoef is used.
- Outputs:
    - CSV with columns: time, observed, model, residual
    - JSON .json file updated with metadata (metrics, parameters used)

Dependencies:
  numpy, pandas, scipy (optional but recommended for Pearson), matplotlib (optional for plotting)
"""

from __future__ import annotations
import argparse
import json
import math
import os
from typing import Any, Dict, Tuple, Callable

import numpy as np
import pandas as pd

# Try to import scipy.stats for pearsonr, else fallback to numpy
try:
    from scipy.stats import pearsonr
    _HAS_SCIPY = True
except Exception:
    pearsonr = None  # type: ignore
    _HAS_SCIPY = False

# Optional plotting
try:
    import matplotlib.pyplot as plt
    _HAS_MPL = True
except Exception:
    _HAS_MPL = False

TWOPI = 2.0 * math.pi


def read_two_column_csv(path: str) -> Tuple[np.ndarray, np.ndarray]:
    """Read CSV-like file with two columns: time, value. Tolerant to header rows."""
    df = pd.read_csv(path, header=0)
    if df.shape[1] < 2:
        raise ValueError("Input CSV must have at least two columns (time, value).")
    t = pd.to_numeric(df.iloc[:, 0], errors='coerce')
    y = pd.to_numeric(df.iloc[:, 1], errors='coerce')
    mask = (~t.isna()) & (~y.isna())
    t = t[mask].to_numpy(dtype=float)
    y = y[mask].to_numpy(dtype=float)
    return t, y


def read_json_p(path: str) -> Dict[str, Any]:
    """Read JSON file at path (if exists), else return empty dict."""
    if not os.path.exists(path):
        return {}
    with open(path, 'r', encoding='utf-8') as fh:
        return json.load(fh)


def write_json_p(path: str, data: Dict[str, Any]) -> None:
    """Write JSON dict to file path."""
    with open(path, 'w', encoding='utf-8') as fh:
        json.dump(data, fh, indent=2, sort_keys=True)


def clone_series(series: np.ndarray) -> np.ndarray:
    """Return a deep copy (clone) of the series."""
    return np.array(series, copy=True)


# ---- Model-step interface ---------------------------------------------------
# A model-step function implements:
#   model_value, state = model_step(i, t_i, clone_i, state, params)
# Where:
#   i: integer time index (0..N-1)
#   t_i: timestamp value (float)
#   clone_i: cloned original series value at index i (float)
#   state: dict carrying model internal state (mutable or replaced)
#   params: dict of user parameters
#
# The runner will set up initial state = {} and call model_step for each i.
# ---------------------------------------------------------------------------


def model_step_clone(i: int, t_i: float, clone_i: float, state: Dict[str, Any], params: Dict[str, Any]):
    """Trivial clone model: model value equals the cloned series value."""
    return float(clone_i), state


def model_step_example(i: int, t_i: float, clone_i: float, state: Dict[str, Any], params: Dict[str, Any]):
    """
    Simple example model: low-pass sample-and-hold of the clone plus a small sinusoidal
    modulation. Demonstrates a stateful update.

    Parameters in params:
      - alpha: exponential smoothing parameter in [0,1] (default 0.9)
      - sin_amp: amplitude of additive sine (default 0.0)
      - sin_freq: frequency in cycles per year (default 1.0)
      - sin_phase: phase in radians (default 0.0)
    """
    alpha = float(params.get('alpha', 0.9))
    sin_amp = float(params.get('sin_amp', 0.0))
    sin_freq = float(params.get('sin_freq', 1.0))
    sin_phase = float(params.get('sin_phase', 0.0))

    # stateful accumulator
    prev = state.get('acc', 0.0)
    # simple exponential smoothing toward the clone value
    acc = alpha * prev + (1.0 - alpha) * float(clone_i)
    state['acc'] = acc

    mod = sin_amp * math.sin(TWOPI * t_i * sin_freq + sin_phase)
    model_value = acc + mod
    return float(model_value), state


def model_step_spreadsheet(i: int, t_i: float, clone_i: float, state: Dict[str, Any], params: Dict[str, Any]):
    """
    Implements the spreadsheet-style model described:
      - Impulse mask every Imp_Stride rows (row numbering starts at 1)
      - C_t = (Imp_Amp if (row_idx % Imp_Stride == 1) else 0) * SUM_k PeriodsAmp[k]*sin(2*pi * t / Period_k + PeriodsPhase[k])
      - D_t = (1-Hold)*D_{t-1} + Hold*C_t    (initialize D_0 = Hold*C_0)
      - E_t = LTE_Amp * sin(D_t * LTE_Freq + LTE_Phase)

    Required params (defaults provided):
      - Imp_Stride: int (default 1)
      - Imp_Amp: float (default 1.0)
      - Periods: list of floats (default [1.0])
      - PeriodsAmp: list of floats (len match Periods, default [1.0])
      - PeriodsPhase: list of floats (len match Periods, default [0.0])
      - Hold: float in [0,1] (default 0.5)
      - LTE_Amp: float (default 1.0)
      - LTE_Freq: float (default 1.0)
      - LTE_Phase: float (default 0.0)
    """
    # Prepare and normalize parameters
    Initial = float(params.get('Initial', 0.04294))
    Year = float(params.get('Year', 365.242))
    Imp_Stride = int(params.get('Imp_Stride', 12))
    Imp_Amp = float(params.get('Imp_Amp', 1.0))
    Periods = list(params.get('Periods', [1.0]))
    PeriodsAmp = list(params.get('PeriodsAmp', [1.0] * len(Periods)))
    PeriodsPhase = list(params.get('PeriodsPhase', [0.0] * len(Periods)))
    Aliased = list(params.get('Aliased', [1.0]))
    AliasedAmp = list(params.get('AliasedAmp', [1.0] * len(Aliased)))
    AliasedPhase = list(params.get('AliasedPhase', [0.0] * len(Aliased)))
    Hold = float(params.get('Hold', 0.5))
    LTE_Amp = float(params.get('LTE_Amp', 1.0))
    LTE_Freq = float(params.get('LTE_Freq', 1.0))
    LTE_Phase = float(params.get('LTE_Phase', 0.0))

    # compute sum of sinusoids
    ssum = 0.0
    for amp, per, ph in zip(PeriodsAmp, Periods, PeriodsPhase):
        per = float(per)
        if per == 0.0:
            continue
        ssum += float(amp) * math.sin(TWOPI * float(t_i) * Year / per + float(ph))

    # row index as 1..N (i is 0..N-1)
    row_idx = i + 2
    impulse = 1.0 if (row_idx % max(1, Imp_Stride) == 1) else 0.0
    C_t = impulse * Imp_Amp * ssum

    # stateful D
    if 'D_prev' not in state:
        # initialize as Hold * C_0 (spreadsheet initialization)
        state['D_prev'] = Initial

    D_prev = float(state['D_prev'])
    D_t = (1.0 - Hold) * D_prev + Hold * C_t
    state['D_prev'] = D_t

    # compute sum of aliases
    ssum = 0.0
    for amp, freq, ph in zip(AliasedAmp, Aliased, AliasedPhase):
        freq = float(freq)
        ssum += float(amp) * math.sin(TWOPI * float(t_i) * freq + float(ph))

    E_t = LTE_Amp * math.sin(D_t * LTE_Freq + LTE_Phase) + ssum

    return float(E_t), state


# Runner that loops through timestamps and produces modeled series ----------------
def run_loop_time_series(time: np.ndarray,
                         observed: np.ndarray,
                         model_step_fn: Callable,
                         params: Dict[str, Any]) -> Tuple[np.ndarray, Dict[str, Any]]:
    """
    Loop over timestamps and compute model values using model_step_fn.

    Returns (model_array, final_state).
    """
    N = time.size
    clone = clone_series(observed)
    state: Dict[str, Any] = {}
    model = np.zeros_like(observed, dtype=float)
    for i in range(N):
        t_i = float(time[i])
        clone_i = float(clone[i])
        v, state = model_step_fn(i, t_i, clone_i, state, params)
        model[i] = float(v)
    return model, state


# Metrics -----------------------------------------------------------------------
def compute_metrics(observed: np.ndarray, model: np.ndarray) -> Dict[str, Any]:
    """Compute Pearson correlation, residuals, MSE, and variance of squared errors."""
    residuals = model - observed
    mse = float(np.mean(residuals ** 2))
    # variance of squared errors (i.e., variance of residuals**2)
    var_sq_err = float(np.var(residuals ** 2))
    # Pearson r (use scipy if available else numpy)
    try:
        if _HAS_SCIPY:
            r_val, p_val = pearsonr(observed, model)  # type: ignore
            pearson_r = float(r_val)
            pearson_p = float(p_val)
        else:
            # fallback
            cov = np.cov(observed, model, ddof=0)
            denom = math.sqrt(cov[0, 0] * cov[1, 1]) if cov[0, 0] > 0 and cov[1, 1] > 0 else 0.0
            if denom == 0:
                pearson_r = float(np.nan)
                pearson_p = float(np.nan)
            else:
                pearson_r = float(cov[0, 1] / denom)
                pearson_p = float(np.nan)
    except Exception:
        pearson_r = float(np.nan)
        pearson_p = float(np.nan)

    return {
        'residuals': residuals,
        'mse': mse,
        'var_squared_error': var_sq_err,
        'pearson_r': pearson_r,
        'pearson_p': pearson_p
    }


# Utilities ---------------------------------------------------------------------
def parse_json_like_arg(s: str | None) -> Dict[str, Any]:
    """Parse a JSON-like string from CLI --params. Returns dict."""
    if s is None:
        return {}
    try:
        return json.loads(s)
    except Exception as e:
        raise argparse.ArgumentTypeError(f"Invalid JSON for --params: {e}")


def main():
    ap = argparse.ArgumentParser(description="Time-series modeler with per-timestamp loop and metrics.")
    ap.add_argument('csv', help='Input CSV file (two columns: time, value).')
    ap.add_argument('--out', default='fitted_out.csv', help='Output CSV file with time, observed, model, residual.')
    ap.add_argument('--model', choices=['clone', 'example', 'spreadsheet'], default='example',
                    help='Which built-in model-step to use. Default "example".')
    ap.add_argument('--plot', action='store_true', help='Show plot of observed vs model (requires matplotlib).')
    args = ap.parse_args()

    # Read CSV
    time, y = read_two_column_csv(args.csv)
    if time.size == 0:
        raise SystemExit("No data read from CSV.")

    # Read JSON file with same name appended by ".json"
    json_path = args.csv + '.json'
    data_dict = read_json_p(json_path)
    params = data_dict

    # Choose model function
    if args.model == 'clone':
        model_fn = model_step_clone
    elif args.model == 'example':
        model_fn = model_step_example
    else:
        model_fn = model_step_spreadsheet

    # Make a clone of the series (not strictly necessary but requested)
    cloned = clone_series(y)
    # y = [0.0] * len(y)

    # Run the time-stepping loop (explicit loop per timestamp)
    model_vals, final_state = run_loop_time_series(time, cloned, model_fn, params)

    # Compute metrics
    metrics = compute_metrics(y, model_vals)

    # Prepare output DataFrame
    out_df = pd.DataFrame({
        'time': time,
        'observed': y,
        'model': model_vals,
        'residual': metrics['residuals']
    })
    out_df.to_csv(args.out, index=False)

    # Update data dictionary with metadata and save to .json file
    data_dict_out = dict(data_dict)  # copy existing
    data_dict_out.setdefault('source_csv', os.path.basename(args.csv))
    data_dict_out.setdefault('n_samples', int(time.size))
    data_dict_out.setdefault('model_type', args.model)
    data_dict_out.setdefault('model_params', params)
    data_dict_out.setdefault('final_state', _serialize_state_for_json(final_state))
    data_dict_out.setdefault('metrics', {
        'mse': metrics['mse'],
        'var_squared_error': metrics['var_squared_error'],
        'pearson_r': metrics['pearson_r'],
        'pearson_p': metrics['pearson_p']
    })

    write_json_p(json_path, data_dict_out)

    # Print summary
    print(f"Processed {time.size} samples from: {args.csv}")
    print(f"Model type: {args.model}")
    print(f"Output CSV: {args.out}")
    print(f"Updated JSON data file: {json_path}")
    print("Metrics:")
    print(f"  MSE = {metrics['mse']:.6e}")
    print(f"  Variance of squared errors = {metrics['var_squared_error']:.6e}")
    print(f"  Pearson r = {metrics['pearson_r']}, p-value = {metrics.get('pearson_p', None)}")

    # Optional plotting
    if args.plot:
        if not _HAS_MPL:
            print("matplotlib not available; cannot plot.")
        else:
            plt.figure(figsize=(12, 6))
            plt.plot(time, y, label='observed', lw=1)
            plt.plot(time, model_vals, label='model', lw=1)
            # plt.plot(time, metrics['residuals'], label='residual', lw=0.8, alpha=0.7)
            plt.legend()
            plt.xlabel('time')
            plt.ylabel('value')
            plt.title(f"Model: {args.model}  Pearson r={metrics['pearson_r']:.4g}  MSE={metrics['mse']:.4g}")
            plt.grid(True)
            plt.show()


def _serialize_state_for_json(state: Dict[str, Any]) -> Dict[str, Any]:
    """Helper to convert final state to JSON-serializable form."""
    out = {}
    for k, v in state.items():
        try:
            json.dumps({k: v})
            out[k] = v
        except Exception:
            # try to coerce numpy scalars
            if isinstance(v, (np.generic,)):
                out[k] = np.asscalar(v)  # type: ignore[attr-defined]
            else:
                out[k] = str(v)
    return out


if __name__ == '__main__':
    main()