tcad-bodeplot/dea/csv_processing.py

import io
import os
import uuid
import re
import tempfile
import time

import numpy as np
import pandas as pd
from fastapi import HTTPException
from scipy import signal
import math

from .config import MAX_CSV_BYTES, MAX_PLOT_POINTS, MAX_ROWS

TEMP_CSV_DIR = os.environ.get(
    "DEA_TEMP_CSV_DIR",
    os.path.join(tempfile.gettempdir(), "diff-eq-analyzer-csv"),
)
TEMP_CSV_TTL_SECONDS = int(os.environ.get("DEA_TEMP_CSV_TTL_SECONDS", str(24 * 60 * 60)))
os.makedirs(TEMP_CSV_DIR, exist_ok=True)


def prune_expired_csv_uploads(now=None):
    if TEMP_CSV_TTL_SECONDS <= 0:
        return
    os.makedirs(TEMP_CSV_DIR, exist_ok=True)
    cutoff = (now or time.time()) - TEMP_CSV_TTL_SECONDS
    for filename in os.listdir(TEMP_CSV_DIR):
        if not filename.endswith(".csv"):
            continue
        path = os.path.join(TEMP_CSV_DIR, filename)
        try:
            if os.path.isfile(path) and os.path.getmtime(path) < cutoff:
                os.remove(path)
        except OSError:
            pass


def get_cached_file_path(file_id):
    if not re.match(r'^[0-9a-f\-]{36}$', file_id):
        raise HTTPException(status_code=400, detail="無效的檔案ID")
    file_path = os.path.join(TEMP_CSV_DIR, f"{file_id}.csv")
    if not os.path.exists(file_path):
        raise HTTPException(status_code=404, detail="找不到快取的 CSV 檔案，請重新上傳")
    return file_path


def downsample_for_plot(index, original, filtered_float, filtered_fixed):
    total = len(index)
    if total <= MAX_PLOT_POINTS:
        return index, original, filtered_float, filtered_fixed, 1
    step = int(np.ceil(total / MAX_PLOT_POINTS))
    return index[::step], original[::step], filtered_float[::step], filtered_fixed[::step], step


async def save_csv_upload(file):
    prune_expired_csv_uploads()
    filename = (file.filename or "").lower()
    if filename and not filename.endswith(".csv"):
        raise HTTPException(status_code=400, detail="請上傳 CSV 檔案")
    contents = await file.read(MAX_CSV_BYTES + 1)
    if len(contents) > MAX_CSV_BYTES:
        raise HTTPException(status_code=413, detail=f"CSV 檔案不可超過 {MAX_CSV_BYTES // (1024 * 1024)}MB")
    if not contents.strip():
        raise HTTPException(status_code=400, detail="CSV 不可為空")

    file_id = str(uuid.uuid4())
    file_path = os.path.join(TEMP_CSV_DIR, f"{file_id}.csv")
    with open(file_path, "wb") as f:
        f.write(contents)
    return file_id


def integer_lfilter(b_int, a_int, x_float, shift_in, shift_out, shift_b, shift_a, use_round=False):
    """
    全整數差分方程式模擬 (Mimic DSP hardware)
    高精度狀態變數架構：前饋不位移，保留小數精度於狀態變數中。
    支援硬體 Rounding (四捨五入) vs Floor (無條件捨去向下取整)。
    """
    b_int = np.asarray(b_int, dtype=np.int64)
    a_int = np.asarray(a_int, dtype=np.int64)
    x_int = np.round(x_float * (2**shift_in)).astype(np.int64)

    # y_hist 將保留在 Q_{in + b} 格式以降低 Truncation Error
    y_hist = np.zeros(len(x_int), dtype=np.int64)
    y_out = np.zeros(len(x_int), dtype=np.int64)

    nb = len(b_int)
    na = len(a_int)

    A0 = int(a_int[0])
    if A0 == 0: A0 = 1

    # 輸出所需的總位移量
    out_shift = shift_in + shift_b - shift_out

    # 預先計算四捨五入的補償值 (+0.5)
    # 韌體開發提示 (C Implementation / RISC-V):
    # 1. 標準 RISC-V (RV32I/IMAC) 的 SRA 指令是純 Floor (無條件捨去)，沒有硬體 Rounding shift。
    #    (除非具備 'P' DSP Extension 才可能有 1-cycle 的硬體 rounding shift)。
    # 2. 演算法秘技：在 C 語言中要實現 1-clock 的四捨五入，不要呼叫 float 的 round()。
    #    請使用 `y = (acc + (1 << (shift - 1))) >> shift`。
    #    編譯器會將 (1 << (shift - 1)) 編譯為常數，整體只消耗 1 個 ADD 指令，極大消除了 DC Bias。
    round_offset_a = (A0 >> 1) if use_round else 0
    round_offset_out = (1 << (out_shift - 1)) if (use_round and out_shift > 0) else 0

    for n in range(len(x_int)):
        sum_b = 0
        # Feedforward: 前饋完全不位移，結果為 Q_{in + b}
        for i in range(nb):
            if n - i >= 0:
                sum_b += b_int[i] * x_int[n - i]

        # Feedback: 歷史紀錄為 Q_{in + b}，係數為 Q_a，乘積為 Q_{in + b + a}
        sum_a = 0
        for j in range(1, na):
            if n - j >= 0:
                sum_a += a_int[j] * y_hist[n - j]

        # Feedback 縮放：前提 A0 = 1 (或者代表放大了 Q_a 倍的常數)，除以 A0 (即 >> shift_a) 歸一化
        # Python 的 // 等同於硬體的 SRA (Arithmetic Right Shift)，會向負無窮大 Floor
        sum_a_scaled = (sum_a + round_offset_a) // A0

        acc = sum_b - sum_a_scaled

        # 將超高精度的 acc 直接存入歷史變數
        y_hist[n] = acc

        # 最終輸出再針對 Q_out 進行位移縮放
        if out_shift > 0:
            y_out[n] = (acc + round_offset_out) >> out_shift
        elif out_shift < 0:
            y_out[n] = acc << (-out_shift)
        else:
            y_out[n] = acc

    return y_out.astype(float) / (2**shift_out)


def filter_preview_response(file_id, b_vals, a_vals, col_idx, b_int=None, a_int=None, shift_in=14, shift_out=14, shift_b=14, shift_a=14, use_round=False):
    path = get_cached_file_path(file_id)

    # 預先讀取欄位名稱，避免用 usecols 讀取後找不到原始索引
    cols = pd.read_csv(path, nrows=0).columns.tolist()
    if col_idx < 0 or col_idx >= len(cols):
        raise HTTPException(status_code=400, detail="欄位索引超出範圍")
    col_to_filter = cols[col_idx]

    # 精準讀取單一欄位，並加上筆數限制 (極大降低記憶體用量與時間)
    df = pd.read_csv(path, usecols=[col_idx], nrows=MAX_ROWS)
    x_signal = pd.to_numeric(df[col_to_filter], errors="coerce")

    if x_signal.isna().any():
        raise HTTPException(status_code=400, detail=f"欄位 {col_to_filter} 含有非數值資料")
    x_values = x_signal.to_numpy(dtype=float)
    if not np.isfinite(x_values).all():
        raise HTTPException(status_code=400, detail=f"欄位 {col_to_filter} 含有非有限數值")

    # 路徑 1: 理想浮點數路徑
    y_float = signal.lfilter(b_vals, a_vals, x_values)

    # 路徑 2: 整數模擬路徑
    if b_int is not None and a_int is not None:
        y_fixed = integer_lfilter(b_int, a_int, x_values, shift_in, shift_out, shift_b, shift_a, use_round)
    else:
        y_fixed = y_float

    index, original, filtered_float, filtered_fixed, step = downsample_for_plot(df.index.to_numpy(), x_signal, y_float, y_fixed)

    return {
        "index": index.tolist(),
        "original": original.tolist(),
        "filtered": filtered_float.tolist(),
        "filtered_fixed": filtered_fixed.tolist(),
        "col_name": col_to_filter,
        "total_points": int(len(df.index)),
        "plot_points": int(len(index)),
        "downsample_step": int(step),
    }


def filtered_csv_text(file_id, b_vals, a_vals, col_idx, b_int=None, a_int=None, shift_in=14, shift_out=14, shift_b=14, shift_a=14, use_round=False):
    path = get_cached_file_path(file_id)
    # 匯出時需要原始所有欄位，但仍受限於 MAX_ROWS
    df = pd.read_csv(path, nrows=MAX_ROWS)

    if col_idx < 0 or col_idx >= len(df.columns):
        raise HTTPException(status_code=400, detail="欄位索引超出範圍")
    col_to_filter = df.columns[col_idx]
    x_signal = pd.to_numeric(df[col_to_filter], errors="coerce")
    x_values = x_signal.to_numpy(dtype=float)

    y_float = signal.lfilter(b_vals, a_vals, x_values)
    if b_int is not None and a_int is not None:
        y_fixed = integer_lfilter(b_int, a_int, x_values, shift_in, shift_out, shift_b, shift_a, use_round)
    else:
        y_fixed = y_float

    df[f"{col_to_filter}_filtered_ideal"] = y_float
    df[f"{col_to_filter}_filtered_fixed"] = y_fixed
    csv_buffer = io.StringIO()
    df.to_csv(csv_buffer, index=False)
    return csv_buffer.getvalue()