Skip to content

Python Examples

Complete Python examples demonstrating fastloess capabilities with NumPy and matplotlib.

Batch Smoothing

Process complete datasets with confidence intervals, diagnostics, and cross-validation.

#!/usr/bin/env python3
"""
fastloess Batch Smoothing - Comprehensive Examples

 1. Basic smoothing
 2. Robust smoothing with outliers
 3. Uncertainty quantification (confidence/prediction intervals)
 4. Cross-validation (K-Fold)
 5. Complete diagnostic analysis
 6. Different weight functions (kernels)
 7. Robustness methods comparison
 8. Benchmark
 9. Scaling methods (MAR, MAD, Mean)
10. Boundary policies
11. Zero-weight fallback strategies
12. Polynomial degrees + iterations_used
13. Distance metrics
14. Surface modes and standard errors
15. Additional weight functions
16. LOOCV and auto-converge
17. Interpolation tuning (surface_mode effects)
"""

import time

import numpy as np
from fastloess import Loess


def make_linear(n: int):
    x = np.arange(n, dtype=float)
    y = 2.0 * x + 1.0
    return x, y


# ── Example 1: Basic Smoothing ───────────────────────────────────────────────
def example_1_basic_smoothing():
    print("Example 1: Basic Smoothing")

    x = np.array([1, 2, 3, 4, 5], dtype=float)
    y = np.array([2.0, 4.1, 5.9, 8.2, 9.8])

    result = Loess(fraction=0.5, iterations=3).fit(x, y)

    print(f"  fraction_used={result.fraction_used}")
    print(f"  Smoothed: [{', '.join(f'{v:.3f}' for v in result.y)}]")
    print()


# ── Example 2: Robust Smoothing with Outliers ────────────────────────────────
def example_2_robust_with_outliers():
    print("Example 2: Robust Smoothing with Outliers")

    x = np.array([1, 2, 3, 4, 5, 6, 7, 8], dtype=float)
    y = np.array([2.1, 4.0, 5.9, 25.0, 10.1, 12.0, 14.1, 15.9])  # 25.0 outlier

    result = Loess(
        fraction=0.5,
        iterations=5,
        robustness_method="bisquare",
        return_robustness_weights=True,
        return_residuals=True,
    ).fit(x, y)

    if result.robustness_weights is not None:
        for i, w in enumerate(result.robustness_weights):
            if w < 0.5:
                print(f"  Outlier at index {i} (y={y[i]}): weight={w:.3f}")
    print(f"  Smoothed: [{', '.join(f'{v:.2f}' for v in result.y)}]")
    print()


# ── Example 3: Uncertainty Quantification ───────────────────────────────────
def example_3_uncertainty_quantification():
    print("Example 3: Uncertainty Quantification")

    x = np.array([1, 2, 3, 4, 5, 6, 7, 8], dtype=float)
    y = np.array([2.1, 3.8, 6.2, 7.9, 10.3, 11.8, 14.1, 15.7])

    result = Loess(
        fraction=0.5,
        iterations=3,
        confidence_intervals=0.95,
        prediction_intervals=0.95,
    ).fit(x, y)

    print("  x  y_smooth  conf_low  conf_high  pred_low  pred_high")
    cl = result.confidence_lower
    cu = result.confidence_upper
    pl = result.prediction_lower
    pu = result.prediction_upper
    if cl is not None and cu is not None and pl is not None and pu is not None:
        for i in range(len(result.y)):
            print(
                f"  {result.x[i]:.0f}  {result.y[i]:.4f}  "
                f"{cl[i]:.4f}  {cu[i]:.4f}  "
                f"{pl[i]:.4f}  {pu[i]:.4f}"
            )
    print()


# ── Example 4: Cross-Validation ──────────────────────────────────────────────
def example_4_cross_validation():
    print("Example 4: Cross-Validation for Parameter Selection")

    x = np.arange(1, 21, dtype=float)
    y = 2 * x + 1 + np.sin(x * 0.5)

    result = Loess(
        cv_fractions=[0.2, 0.3, 0.5, 0.7],
        cv_method="kfold",
        cv_k=5,
        iterations=2,
        return_diagnostics=True,
    ).fit(x, y)

    print(f"  Selected fraction: {result.fraction_used}")
    if result.cv_scores is not None:
        fracs = [0.2, 0.3, 0.5, 0.7]
        print("  CV Scores (RMSE per fraction):")
        for frac, score in zip(fracs, result.cv_scores):
            print(f"    fraction={frac}: {score:.4f}")
    print()


# ── Example 5: Complete Diagnostic Analysis ──────────────────────────────────
def example_5_complete_diagnostics():
    print("Example 5: Complete Diagnostic Analysis")

    x = np.array([1, 2, 3, 4, 5, 6, 7, 8], dtype=float)
    y = np.array([2.1, 3.8, 6.2, 7.9, 10.3, 11.8, 14.1, 15.7])

    result = Loess(
        fraction=0.5,
        iterations=3,
        confidence_intervals=0.95,
        prediction_intervals=0.95,
        return_diagnostics=True,
        return_residuals=True,
        return_robustness_weights=True,
    ).fit(x, y)

    if result.diagnostics is not None:
        d = result.diagnostics
        print("  Diagnostics:")
        print(f"    RMSE:        {d.rmse:.6f}")
        print(f"    MAE:         {d.mae:.6f}")
        print(f"    R²:          {d.r_squared:.6f}")
        print(f"    Residual SD: {d.residual_sd:.6f}")
        if d.aic is not None:
            print(f"    AIC:         {d.aic:.2f}")
        if d.aicc is not None:
            print(f"    AICc:        {d.aicc:.2f}")
        if d.effective_df is not None:
            print(f"    Eff. DF:     {d.effective_df:.2f}")

    print(f"  smoothed[0]: {result.y[0]:.5f}")
    if result.residuals is not None:
        print(f"  residuals[0]: {result.residuals[0]:.5f}")
    if result.robustness_weights is not None:
        print(f"  rob_weight[0]: {result.robustness_weights[0]:.4f}")
    print()


# ── Example 6: Different Weight Functions (Kernels) ──────────────────────────
def example_6_different_kernels():
    print("Example 6: Different Weight Functions (Kernels)")

    x = np.array([1, 2, 3, 4, 5], dtype=float)
    y = np.array([2.0, 4.1, 5.9, 8.2, 9.8])

    for kernel in ["tricube", "epanechnikov", "gaussian", "biweight"]:
        result = Loess(fraction=0.5, weight_function=kernel).fit(x, y)
        print(f"  {kernel}: [{', '.join(f'{v:.3f}' for v in result.y)}]")
    print()


# ── Example 7: Robustness Methods Comparison ─────────────────────────────────
def example_7_robustness_methods():
    print("Example 7: Robustness Methods Comparison")

    x = np.array([1, 2, 3, 4, 5], dtype=float)
    y = np.array([2.0, 4.1, 20.0, 8.2, 9.8])  # 20.0 is an outlier

    for method in ["bisquare", "huber", "talwar"]:
        result = Loess(
            fraction=0.5,
            iterations=5,
            robustness_method=method,
            return_robustness_weights=True,
        ).fit(x, y)
        print(f"  {method}:")
        print(f"    Smoothed: [{', '.join(f'{v:.2f}' for v in result.y)}]")
        if result.robustness_weights is not None:
            print(
                f"    Weights:  [{', '.join(f'{w:.3f}' for w in result.robustness_weights)}]"
            )
    print()


# ── Example 8: Benchmark ─────────────────────────────────────────────────────
def example_8_benchmark():
    print("Example 8: Benchmark")

    n = 1000
    x = np.arange(n, dtype=float)
    y = np.sin(x * 0.1) + np.cos(x * 0.01)

    t0 = time.perf_counter()
    result = Loess(parallel=True).fit(x, y)
    ms = (time.perf_counter() - t0) * 1000

    print(f"  {n} points in {ms:.2f}ms")
    print(f"  fraction_used={result.fraction_used}, y[0]={result.y[0]:.4f}")
    print()


# ── Example 9: Scaling Methods (MAR, MAD, Mean) ──────────────────────────────
def example_9_scaling_methods():
    print("Example 9: Scaling Methods")

    x, y = make_linear(20)

    for method in ["mar", "mad", "mean"]:
        result = Loess(fraction=0.5, scaling_method=method).fit(x, y)
        print(f"  {method}: y[0]={result.y[0]:.3f}")
    print()


# ── Example 10: Boundary Policies ────────────────────────────────────────────
def example_10_boundary_policies():
    print("Example 10: Boundary Policies")

    x, y = make_linear(30)

    for policy in ["extend", "reflect", "zero", "noboundary"]:
        result = Loess(fraction=0.5, boundary_policy=policy).fit(x, y)
        print(f"  {policy}: first={result.y[0]:.2f}, last={result.y[-1]:.2f}")
    print()


# ── Example 11: Zero-Weight Fallback Strategies ───────────────────────────────
def example_11_zero_weight_fallback():
    print("Example 11: Zero-Weight Fallback Strategies")

    x, y = make_linear(20)

    for fb in ["use_local_mean", "return_original", "return_none"]:
        result = Loess(fraction=0.5, zero_weight_fallback=fb).fit(x, y)
        print(f"  {fb}: y[0]={result.y[0]:.3f}")
    print()


# ── Example 12: Polynomial Degrees + iterations_used ──────────────────────────
def example_12_polynomial_degrees():
    print("Example 12: Polynomial Degrees")

    x, y = make_linear(30)

    for deg in ["constant", "linear", "quadratic", "cubic", "quartic"]:
        result = Loess(fraction=0.5, iterations=2, degree=deg).fit(x, y)
        print(
            f"  {deg}: y[0]={result.y[0]:.3f}, iterations_used={result.iterations_used}"
        )
    print()


# ── Example 13: Distance Metrics ─────────────────────────────────────────────
def example_13_distance_metrics():
    print("Example 13: Distance Metrics")

    x, y = make_linear(20)

    for metric in ["euclidean", "normalized", "manhattan", "chebyshev"]:
        result = Loess(fraction=0.5, distance_metric=metric).fit(x, y)
        print(f"  {metric}: y[0]={result.y[0]:.3f}")

    # Minkowski with custom p via the "minkowski:p" string format
    result_mink = Loess(fraction=0.5, distance_metric="minkowski:3").fit(x, y)
    print(f"  minkowski(p=3): y[0]={result_mink.y[0]:.3f}")
    print()


# ── Example 14: Surface Modes and Standard Errors ────────────────────────────
def example_14_surface_modes_and_se():
    print("Example 14: Surface Modes and Standard Errors")

    x, y = make_linear(30)

    # Direct surface — fits every point exactly; SE fields fully populated
    r_direct = Loess(
        fraction=0.5,
        surface_mode="direct",
        return_se=True,
        confidence_intervals=0.95,
        prediction_intervals=0.95,
    ).fit(x, y)

    print("  surface_mode=direct:")
    print(f"    confidence_lower non-null: {r_direct.confidence_lower is not None}")
    print(f"    prediction_lower non-null: {r_direct.prediction_lower is not None}")
    if r_direct.standard_errors is not None:
        print(f"    standard_errors[0]: {r_direct.standard_errors[0]:.4f}")
    if r_direct.enp is not None:
        print(f"    enp: {r_direct.enp:.3f}")
    if r_direct.trace_hat is not None:
        print(f"    trace_hat: {r_direct.trace_hat:.3f}")
    if r_direct.delta1 is not None:
        print(f"    delta1: {r_direct.delta1:.3f}")
    if r_direct.delta2 is not None:
        print(f"    delta2: {r_direct.delta2:.3f}")
    if r_direct.residual_scale is not None:
        print(f"    residual_scale: {r_direct.residual_scale:.4f}")
    if r_direct.leverage is not None:
        print(f"    leverage[0]: {r_direct.leverage[0]:.4f}")

    # Interpolation surface — faster, approximate
    r_interp = Loess(fraction=0.5, surface_mode="interpolation", return_se=True).fit(
        x, y
    )

    print("  surface_mode=interpolation:")
    print(f"    y[0]: {r_interp.y[0]:.3f}")
    if r_interp.standard_errors is not None:
        print(f"    standard_errors[0]: {r_interp.standard_errors[0]:.4f}")
    print()


# ── Example 15: Additional Weight Functions (Uniform, Triangle, Cosine) ───────
def example_15_additional_kernels():
    print("Example 15: Additional Weight Functions (Uniform, Triangle, Cosine)")

    x = np.array([1, 2, 3, 4, 5], dtype=float)
    y = np.array([2.0, 4.1, 5.9, 8.2, 9.8])

    for kernel in ["uniform", "triangle", "cosine"]:
        result = Loess(fraction=0.5, weight_function=kernel).fit(x, y)
        print(f"  {kernel}: [{', '.join(f'{v:.3f}' for v in result.y)}]")
    print()


# ── Example 16: LOOCV, K-Fold, and Auto-Converge ─────────────────────────────
def example_16_loocv_and_auto_converge():
    print("Example 16: LOOCV, K-Fold, and Auto-Converge")

    x = np.arange(1, 21, dtype=float)
    y = 2 * x + 1 + np.sin(x * 0.5)

    # Leave-one-out cross-validation
    r_loocv = Loess(cv_fractions=[0.3, 0.5, 0.7], cv_method="loocv").fit(x, y)
    print(f"  LOOCV selected fraction: {r_loocv.fraction_used}")
    if r_loocv.cv_scores is not None:
        print(f"  LOOCV scores: [{', '.join(f'{s:.4f}' for s in r_loocv.cv_scores)}]")

    # K-Fold cross-validation
    r_kfold = Loess(cv_fractions=[0.2, 0.4, 0.6], cv_method="kfold", cv_k=5).fit(x, y)
    print(f"  KFold(k=5) selected fraction: {r_kfold.fraction_used}")
    if r_kfold.cv_scores is not None:
        print(f"  KFold scores: [{', '.join(f'{s:.4f}' for s in r_kfold.cv_scores)}]")

    # Auto-converge: stop robustness iterations when change < tolerance
    r_ac = Loess(fraction=0.5, auto_converge=1e-4).fit(x, y)
    print(f"  auto_converge=1e-4: iterations_used={r_ac.iterations_used}")
    print()


# ── Example 17: Interpolation Tuning (surface_mode effects) ──────────────────
def example_17_interpolation_tuning():
    print("Example 17: Interpolation Tuning (surface_mode effects)")

    n = 50
    x, y = make_linear(n)

    # Default (interpolation) — fastest, uses a spatial grid
    r_interp = Loess(fraction=0.5, surface_mode="interpolation").fit(x, y)
    print(f"  interpolation: y[0]={r_interp.y[0]:.3f}, y[-1]={r_interp.y[-1]:.3f}")

    # Direct — fits every point exactly, more accurate but slower
    r_direct = Loess(fraction=0.5, surface_mode="direct").fit(x, y)
    print(f"  direct:        y[0]={r_direct.y[0]:.3f}, y[-1]={r_direct.y[-1]:.3f}")

    # Fraction sweep with direct surface
    for frac in [0.2, 0.5, 0.8]:
        r = Loess(fraction=frac, surface_mode="direct").fit(x, y)
        print(f"  direct fraction={frac}: y[0]={r.y[0]:.3f}")

    # Interpolation + SE for hat-matrix statistics
    r_se = Loess(fraction=0.5, surface_mode="interpolation", return_se=True).fit(x, y)
    if r_se.enp is not None:
        print(f"  interpolation+SE enp: {r_se.enp:.3f}")
    print()


# ── Main ──────────────────────────────────────────────────────────────────────
def main():
    print("=" * 60)
    print("fastloess Batch Smoothing - Comprehensive Examples")
    print("=" * 60)
    print()

    example_1_basic_smoothing()
    example_2_robust_with_outliers()
    example_3_uncertainty_quantification()
    example_4_cross_validation()
    example_5_complete_diagnostics()
    example_6_different_kernels()
    example_7_robustness_methods()
    example_8_benchmark()
    example_9_scaling_methods()
    example_10_boundary_policies()
    example_11_zero_weight_fallback()
    example_12_polynomial_degrees()
    example_13_distance_metrics()
    example_14_surface_modes_and_se()
    example_15_additional_kernels()
    example_16_loocv_and_auto_converge()
    example_17_interpolation_tuning()

    print("=== Batch Smoothing Examples Complete ===")


if __name__ == "__main__":
    main()

Download batch_smoothing.py


Streaming Smoothing

Process large datasets in memory-efficient chunks with overlap merging.

#!/usr/bin/env python3
"""
fastloess Streaming Smoothing - Comprehensive Examples

 1. Basic chunked processing
 2. Chunk size comparison
 3. Overlap strategies
 4. Large dataset processing
 5. Outlier handling in streaming mode
 6. File-based streaming simulation
 7. Benchmark (sequential streaming)
 8. Merge strategies
 9. Advanced streaming options
"""

import time

import numpy as np
from fastloess import StreamingLoess


def make_linear(n: int):
    x = np.arange(n, dtype=float)
    y = 2.0 * x + 1.0
    return x, y


def process_all(
    model: StreamingLoess, x: np.ndarray, y: np.ndarray, chunk_size: int, overlap: int
):
    """Feed only full-size chunks; finalize() handles remaining data."""
    step = chunk_size - overlap
    result_x: list = []
    result_y: list = []
    n = len(x)
    start = 0
    while start + chunk_size <= n:
        res = model.process_chunk(
            x[start : start + chunk_size], y[start : start + chunk_size]
        )
        result_x.extend(res.x)
        result_y.extend(res.y)
        start += step
    fin = model.finalize()
    result_x.extend(fin.x)
    result_y.extend(fin.y)
    return np.array(result_x), np.array(result_y)


# ── Example 1: Basic Chunked Processing ─────────────────────────────────────
def example_1_basic_chunked_processing():
    print("Example 1: Basic Chunked Processing")

    n = 50
    x, y = make_linear(n)
    chunk_size, overlap = 15, 5

    model = StreamingLoess(
        fraction=0.5,
        iterations=2,
        chunk_size=chunk_size,
        overlap=overlap,
        return_residuals=True,
    )

    print(f"  Dataset: {n} pts, chunk={chunk_size}, overlap={overlap}")
    total_x: list = []
    total_y: list = []
    ci = 0
    start = 0
    while start + chunk_size <= n:
        res = model.process_chunk(
            x[start : start + chunk_size], y[start : start + chunk_size]
        )
        if len(res.x) > 0:
            total_x.extend(res.x)
            total_y.extend(res.y)
            print(
                f"  Chunk {ci}: {len(res.x)} pts (x: {res.x[0]:.0f}..{res.x[-1]:.0f})"
            )
        start += chunk_size - overlap
        ci += 1
    fin = model.finalize()
    if len(fin.x) > 0:
        total_x.extend(fin.x)
        total_y.extend(fin.y)
        print(f"  Finalize: {len(fin.x)} remaining pts")
    print(f"  Total: {len(total_y)}/{n}")
    print()


# ── Example 2: Chunk Size Comparison ─────────────────────────────────────────
def example_2_chunk_size_comparison():
    print("Example 2: Chunk Size Comparison")

    n = 100
    x, y = make_linear(n)

    for cs, ov, label in [(20, 5, "Small"), (50, 10, "Medium"), (80, 15, "Large")]:
        model = StreamingLoess(fraction=0.5, iterations=1, chunk_size=cs, overlap=ov)
        chunks = 0
        total = 0
        start = 0
        while start + cs <= n:
            res = model.process_chunk(x[start : start + cs], y[start : start + cs])
            if len(res.x) > 0:
                chunks += 1
                total += len(res.x)
            start += cs - ov
        fin = model.finalize()
        if len(fin.x) > 0:
            chunks += 1
            total += len(fin.x)
        print(f"  {label} (size={cs}, overlap={ov}): chunks={chunks}, total={total}")
    print()


# ── Example 3: Overlap Strategies ────────────────────────────────────────────
def example_3_overlap_strategies():
    print("Example 3: Overlap Strategies")

    n = 100
    x, y = make_linear(n)
    cs = 40

    for overlap, label in [
        (0, "No overlap"),
        (10, "10-pt overlap"),
        (20, "20-pt overlap"),
    ]:
        model = StreamingLoess(fraction=0.5, chunk_size=cs, overlap=overlap)
        total = 0
        step = cs - overlap
        start = 0
        while start + cs <= n:
            total += len(
                model.process_chunk(x[start : start + cs], y[start : start + cs]).x
            )
            start += step
        total += len(model.finalize().x)
        print(f"  {label}: total output={total}")
    print()


# ── Example 4: Large Dataset Processing ──────────────────────────────────────
def example_4_large_dataset_processing():
    print("Example 4: Large Dataset Processing")

    n = 10_000
    x = np.arange(n, dtype=float)
    y = np.sin(x * 0.01) + x * 0.001

    cs, ov = 500, 50
    model = StreamingLoess(fraction=0.05, iterations=2, chunk_size=cs, overlap=ov)

    total = 0
    step = cs - ov
    start = 0
    while start + cs <= n:
        total += len(
            model.process_chunk(x[start : start + cs], y[start : start + cs]).x
        )
        if total > 0 and total % 2000 < step:
            print(f"  Progress: ~{total} pts smoothed")
        start += step
    total += len(model.finalize().x)
    print(f"  Total: {total}/{n}, memory: constant (chunk={cs})")
    print()


# ── Example 5: Outlier Handling in Streaming Mode ─────────────────────────────
def example_5_outlier_handling():
    print("Example 5: Outlier Handling in Streaming Mode")

    n = 100
    x = np.arange(n, dtype=float)
    y = 2 * x + 1 + np.sin(x * 0.2) * 2
    y[[25, 50, 75]] += 50  # Outliers

    for method in ["bisquare", "huber", "talwar"]:
        model = StreamingLoess(
            fraction=0.5,
            iterations=5,
            robustness_method=method,
            chunk_size=30,
            overlap=10,
            return_residuals=True,
        )
        large = 0
        start = 0
        while start + 30 <= n:
            res = model.process_chunk(x[start : start + 30], y[start : start + 30])
            if res.residuals is not None:
                large += int(np.sum(np.abs(res.residuals) > 10))
            start += 20
        fin = model.finalize()
        if fin.residuals is not None:
            large += int(np.sum(np.abs(fin.residuals) > 10))
        print(f"  {method}: pts with |residual|>10: {large}")
    print()


# ── Example 6: File-Based Streaming Simulation ───────────────────────────────
def example_6_file_simulation():
    print("Example 6: File-Based Streaming Simulation")
    print("  Simulating: input.csv -> Smooth -> output.csv")

    total_lines, cs, ov = 200, 50, 10
    model = StreamingLoess(
        fraction=0.5, iterations=2, chunk_size=cs, overlap=ov, return_residuals=True
    )

    out_count = 0
    ci = 0
    start_line = 0
    while start_line < total_lines:
        end_line = min(start_line + cs, total_lines)
        xc = np.arange(start_line, end_line, dtype=float)
        yc = 2 * xc + 1 + np.sin(xc * 0.1) * 3
        print(f"  Reading chunk {ci} (lines {start_line}..{end_line - 1})")
        res = model.process_chunk(xc, yc)
        if len(res.x) > 0:
            out_count += len(res.x)
            print(f"    -> Writing {len(res.x)} smoothed pts (total: {out_count})")
        start_line += cs - ov
        ci += 1
    fin = model.finalize()
    if len(fin.x) > 0:
        out_count += len(fin.x)
        print(f"  Finalizing: {len(fin.x)} remaining pts")
    print(f"  Input: {total_lines}, Output: {out_count}")
    print()


# ── Example 7: Benchmark (Sequential Streaming) ───────────────────────────────
def example_7_benchmark():
    print("Example 7: Benchmark (Sequential Streaming)")

    n, cs, ov = 1000, 100, 10
    model = StreamingLoess(fraction=0.5, iterations=3, chunk_size=cs, overlap=ov)

    t0 = time.perf_counter()
    total = 0
    start = 0
    while start + cs <= n:
        xc = np.arange(start, start + cs, dtype=float)
        yc = np.sin(xc * 0.1) + np.cos(xc * 0.01)
        total += len(model.process_chunk(xc, yc).x)
        start += cs - ov
    total += len(model.finalize().x)
    ms = (time.perf_counter() - t0) * 1000

    print(f"  {total} pts in {ms:.2f}ms (chunk={cs}, overlap={ov})")
    print()


# ── Example 8: Merge Strategies ──────────────────────────────────────────────
def example_8_merge_strategies():
    print("Example 8: Merge Strategies")

    n = 50
    x, y = make_linear(n)

    for strategy in ["average", "weighted_average", "take_first", "take_last"]:
        model = StreamingLoess(
            fraction=0.5,
            iterations=2,
            chunk_size=20,
            overlap=5,
            merge_strategy=strategy,
        )
        total = 0
        start = 0
        while start + 20 <= n:
            total += len(
                model.process_chunk(x[start : start + 20], y[start : start + 20]).x
            )
            start += 15
        total += len(model.finalize().x)
        print(f"  {strategy}: total={total}")
    print()


# ── Example 9: Advanced Streaming Options ─────────────────────────────────────
def example_9_advanced_options():
    print("Example 9: Advanced Streaming Options")

    n = 50
    x, y = make_linear(n)

    model = StreamingLoess(
        fraction=0.5,
        iterations=2,
        degree="quadratic",
        scaling_method="mar",
        boundary_policy="reflect",
        zero_weight_fallback="return_original",
        distance_metric="manhattan",
        surface_mode="direct",
        return_se=True,
        return_diagnostics=True,
        return_robustness_weights=True,
        auto_converge=1e-3,
        chunk_size=20,
        overlap=5,
    )

    total = 0
    start = 0
    while start + 20 <= n:
        total += len(
            model.process_chunk(x[start : start + 20], y[start : start + 20]).x
        )
        start += 15
    fin = model.finalize()
    total += len(fin.x)

    print(f"  total pts: {total}")
    if fin.standard_errors is not None and len(fin.standard_errors) > 0:
        print(f"  standard_errors[0]: {fin.standard_errors[0]:.4f}")
    if fin.diagnostics is not None:
        print(f"  diagnostics.rmse: {fin.diagnostics.rmse:.3f}")
        print(f"  diagnostics.r_squared: {fin.diagnostics.r_squared:.3f}")
        if fin.diagnostics.aic is not None:
            print(f"  diagnostics.aic: {fin.diagnostics.aic:.3f}")
    if fin.robustness_weights is not None and len(fin.robustness_weights) > 0:
        print(f"  robustness_weights[0]: {fin.robustness_weights[0]:.4f}")
    print()


# ── Main ──────────────────────────────────────────────────────────────────────
def main():
    print("=" * 60)
    print("fastloess Streaming Smoothing - Comprehensive Examples")
    print("=" * 60)
    print()

    example_1_basic_chunked_processing()
    example_2_chunk_size_comparison()
    example_3_overlap_strategies()
    example_4_large_dataset_processing()
    example_5_outlier_handling()
    example_6_file_simulation()
    example_7_benchmark()
    example_8_merge_strategies()
    example_9_advanced_options()

    print("=== Streaming Smoothing Examples Complete ===")


if __name__ == "__main__":
    main()

Download streaming_smoothing.py


Online Smoothing

Real-time smoothing with sliding window for streaming data applications.

#!/usr/bin/env python3
"""
fastloess Online Smoothing - Comprehensive Examples

 1. Basic incremental processing
 2. Real-time sensor data simulation
 3. Outlier handling in online mode
 4. Window size comparison
 5. Memory-bounded processing (embedded systems)
 6. Sliding window behavior
 7. Benchmark (sequential online)
 8. Update modes (Full vs Incremental) and min_points
 9. Advanced online options
"""

import time

import numpy as np
from fastloess import OnlineLoess


def add_all_points(model, x, y):
    """Feed all (x, y) pairs through add_point; return array of smoothed values (np.nan when None)."""
    smoothed = []
    for xi, yi in zip(x, y):
        r = model.add_point(float(xi), float(yi))
        smoothed.append(r.smoothed if r is not None else float("nan"))
    return np.array(smoothed)


# ── Example 1: Basic Incremental Processing ──────────────────────────────────
def example_1_basic_streaming():
    print("Example 1: Basic Incremental Processing")

    x = np.array([1, 2, 3, 4, 5, 6, 7, 8, 9, 10], dtype=float)
    y = np.array([3.1, 5.0, 7.2, 8.9, 11.1, 13.0, 15.2, 16.8, 19.1, 21.0])

    model = OnlineLoess(fraction=0.5, iterations=2, window_capacity=5)
    smoothed = add_all_points(model, x, y)

    print(f"  {'X':>8} {'Y_obs':>12} {'Y_smooth':>12}")
    for i in range(len(x)):
        sv = f"{smoothed[i]:12.2f}" if not np.isnan(smoothed[i]) else "  (buffering)"
        print(f"  {x[i]:8.2f} {y[i]:12.2f} {sv}")
    print()


# ── Example 2: Real-Time Sensor Data Simulation ───────────────────────────────
def example_2_sensor_data_simulation():
    print("Example 2: Real-Time Sensor Data Simulation")
    print("  Simulating temperature sensor with noise...")

    hours = np.arange(24, dtype=float)
    temp = 20 + 5 * np.sin(hours * np.pi / 12) + (hours * 7 % 11) * 0.3 - 1.5

    model = OnlineLoess(
        fraction=0.4, iterations=3, robustness_method="bisquare", window_capacity=12
    )
    smoothed = add_all_points(model, hours, temp)

    print(f"  {'Hour':>6} {'Raw':>12} {'Smoothed':>12}")
    for i in range(len(hours)):
        sv = (
            f"{smoothed[i]:10.2f}\u00b0C"
            if not np.isnan(smoothed[i])
            else "  (warming up)"
        )
        print(f"  {hours[i]:6.0f} {temp[i]:10.2f}\u00b0C {sv}")
    print()


# ── Example 3: Outlier Handling in Online Mode ────────────────────────────────
def example_3_outlier_handling():
    print("Example 3: Outlier Handling in Online Mode")

    x = np.array([1, 2, 3, 4, 5, 6, 7, 8, 9, 10], dtype=float)
    y = np.array([2.0, 4.1, 5.9, 25.0, 10.1, 12.0, 14.1, 50.0, 18.0, 20.1])

    for method in ["bisquare", "talwar"]:
        model = OnlineLoess(
            fraction=0.5, iterations=5, robustness_method=method, window_capacity=6
        )
        smoothed = add_all_points(model, x, y)
        valid = smoothed[~np.isnan(smoothed)]
        print(f"  {method}: [{', '.join(f'{v:.1f}' for v in valid)}]")
    print()


# ── Example 4: Window Size Comparison ────────────────────────────────────────
def example_4_window_size_comparison():
    print("Example 4: Window Size Comparison")

    x = np.arange(1, 21, dtype=float)
    y = 2 * x + np.sin(x * 0.5) * 3

    for w in [5, 10, 15]:
        model = OnlineLoess(fraction=0.5, iterations=2, window_capacity=w)
        smoothed = add_all_points(model, x, y)
        valid = smoothed[~np.isnan(smoothed)]
        last5 = valid[-5:]
        print(
            f"  window_capacity={w}: last 5 = [{', '.join(f'{v:.2f}' for v in last5)}]"
        )
    print()


# ── Example 5: Memory-Bounded Processing ──────────────────────────────────────
def example_5_memory_bounded_processing():
    print("Example 5: Memory-Bounded Processing (Embedded Systems)")

    total = 1000
    x = np.arange(total, dtype=float)
    y = 2 * x + np.sin(x * 0.1) * 5 + (np.arange(total) % 7 - 3) * 0.5

    model = OnlineLoess(fraction=0.3, iterations=1, window_capacity=20)

    t0 = time.perf_counter()
    smoothed = add_all_points(model, x, y)
    ms = (time.perf_counter() - t0) * 1000

    valid = smoothed[~np.isnan(smoothed)]
    n_out = len(valid)
    for milestone in [200, 400, 600, 800, 1000]:
        if milestone <= n_out:
            print(
                f"  Processed: {milestone:4d} pts | smoothed={valid[milestone - 1]:.2f}"
            )
    print(f"  Total: {n_out}, final smoothed: {valid[-1]:.2f} ({ms:.1f} ms)")
    print("  Memory: constant (window=20)")
    print()


# ── Example 6: Sliding Window Behavior ───────────────────────────────────────
def example_6_sliding_window_behavior():
    print("Example 6: Sliding Window Behavior")

    x = np.array([1, 2, 3, 4, 5, 6, 7, 8], dtype=float)
    y = np.array([2, 4, 6, 8, 10, 12, 14, 16], dtype=float)

    # Use add_point() directly to show buffering behaviour
    model = OnlineLoess(fraction=0.6, iterations=0, window_capacity=4)
    print(f"  {'Pt':>4} {'X':>4} {'Y':>6} {'Smoothed':>10} Status")
    for i, (xi, yi) in enumerate(zip(x, y), 1):
        r = model.add_point(float(xi), float(yi))
        if r is not None:
            print(
                f"  {i:4d} {xi:4.0f} {yi:6.0f} {r.smoothed:10.2f} Window full (sliding)"
            )
        else:
            print(f"  {i:4d} {xi:4.0f} {yi:6.0f} {'-':>10} Filling ({i}/4)")
    print("  Output starts after window fills (4 pts), then slides.")
    print()


# ── Example 7: Benchmark (Sequential Online) ──────────────────────────────────
def example_7_benchmark():
    print("Example 7: Benchmark (Sequential Online)")

    n = 1000
    x = np.arange(n, dtype=float)
    y = np.sin(x * 0.1) + np.cos(x * 0.01)

    model = OnlineLoess(fraction=0.5, iterations=3, window_capacity=10)

    t0 = time.perf_counter()
    smoothed = add_all_points(model, x, y)
    ms = (time.perf_counter() - t0) * 1000

    valid = smoothed[~np.isnan(smoothed)]
    print(f"  {len(valid)} pts processed in {ms:.2f}ms (window_capacity=10)")
    print()


# ── Example 8: Update Modes (Full vs Incremental) and min_points ───────────────
def example_8_update_modes():
    print("Example 8: Update Modes (Full vs Incremental) and min_points")

    x = np.arange(30, dtype=float)
    y = 2 * x + 1.0

    for mode in ["full", "incremental"]:
        model = OnlineLoess(
            fraction=0.5,
            iterations=2,
            update_mode=mode,
            min_points=5,
            window_capacity=15,
        )
        smoothed = add_all_points(model, x, y)
        valid = smoothed[~np.isnan(smoothed)]
        print(f"  {mode}: {len(valid)} pts emitted (out of {len(x)})")

    # Show last smoothed value from the per-point API
    model = OnlineLoess(fraction=0.5, iterations=2, window_capacity=10, min_points=3)
    smoothed = add_all_points(model, x, y)
    valid = smoothed[~np.isnan(smoothed)]
    if len(valid) > 0:
        print(f"  last smoothed: {valid[-1]:.3f}")
    print()


# ── Example 9: Advanced Online Options ────────────────────────────────────────
def example_9_advanced_online_options():
    print("Example 9: Advanced Online Options")

    x = np.arange(30, dtype=float)
    y = 2 * x + 1.0

    model = OnlineLoess(
        fraction=0.5,
        iterations=2,
        degree="quadratic",
        scaling_method="mar",
        boundary_policy="reflect",
        zero_weight_fallback="return_original",
        distance_metric="chebyshev",
        auto_converge=1e-3,
        return_robustness_weights=True,
        min_points=5,
        window_capacity=15,
    )
    result = add_all_points(model, x, y)
    valid = result[~np.isnan(result)]
    print(f"  emitted: {len(valid)}")
    if len(valid) > 0:
        print(f"  last smoothed: {valid[-1]:.3f}")
    print()


# ── Main ──────────────────────────────────────────────────────────────────────
def main():
    print("=" * 60)
    print("fastloess Online Smoothing - Comprehensive Examples")
    print("=" * 60)
    print()

    example_1_basic_streaming()
    example_2_sensor_data_simulation()
    example_3_outlier_handling()
    example_4_window_size_comparison()
    example_5_memory_bounded_processing()
    example_6_sliding_window_behavior()
    example_7_benchmark()
    example_8_update_modes()
    example_9_advanced_online_options()

    print("=== Online Smoothing Examples Complete ===")


if __name__ == "__main__":
    main()

Download online_smoothing.py


Running the Examples

# Install dependencies
pip install fastloess matplotlib numpy

# Run examples
cd examples/python
python batch_smoothing.py
python streaming_smoothing.py
python online_smoothing.py

Output

The batch smoothing example generates visualization plots in examples/python/plots/:

  • batch_main.png - Main smoothing comparison
  • batch_weights.png - Robustness weights visualization
  • batch_boundary.png - Boundary policy comparison