diff --git a/src/algo1.py b/src/algo1.py index 3965dd3..6c67732 100644 --- a/src/algo1.py +++ b/src/algo1.py @@ -3,94 +3,316 @@ # Copyright (C) 2026 Association Exergie # SPDX-License-Identifier: GPL-3.0-or-later -import time -import pandas as pd -from pathlib import Path +from __future__ import annotations + import argparse -from tqdm import tqdm +import sys +from bisect import bisect_left +from dataclasses import dataclass +from pathlib import Path -parser = argparse.ArgumentParser() -parser.add_argument("file", type=Path, help="Source data") +import pandas as pd -args = parser.parse_args() -run_data_path: Path = args.file +DEFAULT_DATA_ROOTS = ( + Path("~/projects/Exergie/Engine-Recorder/recordings-rpm").expanduser(), + Path("~/projects/Exergie/Engine-Recorder/Python/recordings_rpm").expanduser(), +) +GP_COLUMNS = ("gp0_falling", "gp1_falling") +SPARK_COLUMNS = ("spark_1", "spark_2", "spark_3") +INTERPOLATED_COLUMNS = ( + "algo1_continuous_angle_deg", + "algo1_revolution_index", + "algo1_crank_angle_deg", +) -df = pd.read_csv(run_data_path).set_index("time_us", drop=False) -sync_ok = False -cycle_ctr = 0 -cam_ctr = 0 -exhaust_timestamp = 0 -intake_timestamp = 0 -spark_release_time_1 = 0 -spark_release_time_2 = 0 -spark_release_time_3 = 0 -spark_adv = 25 +@dataclass(frozen=True) +class EncoderPoint: + time_us: int + continuous_angle_deg: float -rows = [] -for _, row in tqdm(df.iterrows(), total=len(df)): - time_us: int = row["time_us"] - crank: int = row["crank"] - cam: int = row["cam"] +@dataclass(frozen=True) +class SparkPulse: + time_us: int + column: str - if cam == 1: - if cycle_ctr < 4: - sync_ok = False - if not sync_ok and cycle_ctr == 4: - cam_ctr += 1 - if cam_ctr >= 2: - sync_ok = True - cycle_ctr = 0 - if crank == 1: +def resolve_input_path(path: Path) -> Path: + expanded = path.expanduser() + if expanded.exists(): + return expanded + + for root in DEFAULT_DATA_ROOTS: + candidate = root / path + if candidate.exists(): + return candidate + + searched = ", ".join(str(root) for root in DEFAULT_DATA_ROOTS) + raise FileNotFoundError(f"{path} does not exist; also searched {searched}") + + +def read_rpm(path: Path) -> pd.DataFrame: + df = pd.read_csv(path) + missing = {"time_us", *GP_COLUMNS} - set(df.columns) + if missing: + raise ValueError(f"missing columns: {', '.join(sorted(missing))}") + + df["time_us"] = df["time_us"].astype("int64") + return df + + +def gp_times(df: pd.DataFrame, column: str) -> list[int]: + pulse_rows = df.loc[df[column].fillna(0).astype(int) == 1] + return sorted(int(row.time_us) for row in pulse_rows.itertuples()) + + +def generate_sparks( + crank_times: list[int], + cam_times: list[int], + *, + spark_advance_deg: float, +) -> list[SparkPulse]: + events: list[tuple[int, int, str]] = [] + events.extend((time_us, 0, "cam") for time_us in cam_times) + events.extend((time_us, 1, "crank") for time_us in crank_times) + events.sort() + + sync_ok = False + cycle_ctr = 0 + cam_ctr = 0 + exhaust_timestamp = 0 + intake_timestamp = 0 + sparks: list[SparkPulse] = [] + + for time_us, _, event_type in events: + if event_type == "cam": + if cycle_ctr < 4: + sync_ok = False + if not sync_ok and cycle_ctr == 4: + cam_ctr += 1 + if cam_ctr >= 2: + sync_ok = True + cycle_ctr = 0 + continue + cycle_ctr += 1 if cycle_ctr >= 5: sync_ok = False if cycle_ctr == 2: exhaust_timestamp = time_us - - if cycle_ctr == 3: + elif cycle_ctr == 3: intake_timestamp = time_us - - if cycle_ctr == 4 and sync_ok: + elif cycle_ctr == 4 and sync_ok: d1a = time_us - intake_timestamp d1b = intake_timestamp - exhaust_timestamp - spark_release_time_1 = int( - time_us + (2 * d1a - d1b) * (45 - spark_adv) / 180 - ) + predicted_interval = 2 * d1a - d1b + if d1a <= 0 or predicted_interval <= 0: + continue + + spark_1_time = int(time_us + d1a * (45.0 - spark_advance_deg) / 180.0) + spark_2_time = int(time_us + d1b * (45.0 - spark_advance_deg) / 180.0) t_s = ( - d1a * (135 + spark_adv) / 180 + (2 * d1a - d1b) * (45 - spark_adv) / 180 + d1a * (135.0 + spark_advance_deg) / 180.0 + + predicted_interval * (45.0 - spark_advance_deg) / 180.0 + ) + t_a = d1a / 2.0 + t_s / 2.0 + spark_3_time = int(t_a * (45.0 - spark_advance_deg) / 180.0 + time_us) + + sparks.extend( + ( + SparkPulse(spark_1_time, "spark_1"), + SparkPulse(spark_2_time, "spark_2"), + SparkPulse(spark_3_time, "spark_3"), + ) ) - t_a = d1a / 2 + t_s / 2 + return sorted(set(sparks), key=lambda pulse: (pulse.time_us, pulse.column)) - spark_release_time_2 = int(t_a * (45 - spark_adv) / 180 + time_us) - f_a = (315 + spark_adv) / (d1a) + (45 - spark_adv) / (2 * d1a - d1b) +def encoder_points(df: pd.DataFrame) -> list[EncoderPoint]: + required = {"time_us", "revolution_index", "crank_angle_deg"} + if not required <= set(df.columns): + return [] - t_a = 1 / f_a + points: list[EncoderPoint] = [] + angle_df = df.dropna(subset=["revolution_index", "crank_angle_deg"]) + for row in angle_df.itertuples(): + revolution_index = int(row.revolution_index) + if revolution_index < 0: + continue + angle = revolution_index * 360.0 + float(row.crank_angle_deg) + points.append( + EncoderPoint(time_us=int(row.time_us), continuous_angle_deg=angle) + ) - spark_release_time_3 = int(2 * (t_a * (45 - spark_adv)) + time_us) + deduped: dict[int, EncoderPoint] = {} + for point in sorted(points, key=lambda item: item.time_us): + deduped[point.time_us] = point + return list(deduped.values()) - spark_1 = 1 if spark_release_time_1 == time_us and sync_ok else 0 - spark_2 = 1 if spark_release_time_2 == time_us and sync_ok else 0 - spark_3 = 1 if spark_release_time_3 == time_us and sync_ok else 0 - rows.append( - { - "time_us": time_us, - "crank": crank, - "cam": cam, - "spark_1": spark_1, - "spark_2": spark_2, - "spark_3": spark_3, - } + +def interpolate_angle(points: list[EncoderPoint], time_us: int) -> float | None: + if not points: + return None + + times = [point.time_us for point in points] + index = bisect_left(times, time_us) + + if index < len(points) and points[index].time_us == time_us: + return points[index].continuous_angle_deg + if index == 0 or index == len(points): + return None + + before = points[index - 1] + after = points[index] + span = after.time_us - before.time_us + if span <= 0: + return before.continuous_angle_deg + + fraction = (time_us - before.time_us) / span + return before.continuous_angle_deg + fraction * ( + after.continuous_angle_deg - before.continuous_angle_deg ) -out_df = pd.DataFrame(rows) -base_name, _ = run_data_path.stem.rsplit("_", 1) -output = run_data_path.parent / f"{base_name}_output.csv" -out_df.to_csv(output) + +def add_spark_rows(df: pd.DataFrame, sparks: list[SparkPulse]) -> pd.DataFrame: + for column in SPARK_COLUMNS: + df[column] = 0 + for column in INTERPOLATED_COLUMNS: + if column not in df.columns: + df[column] = pd.NA + + rows_by_time = { + int(row.time_us): index for index, row in enumerate(df.itertuples()) + } + new_rows_by_time: dict[int, dict[str, object]] = {} + + for pulse in sparks: + if pulse.time_us in rows_by_time: + df.at[rows_by_time[pulse.time_us], pulse.column] = 1 + continue + + if pulse.time_us not in new_rows_by_time: + row = {column: pd.NA for column in df.columns} + row["time_us"] = pulse.time_us + for column in SPARK_COLUMNS: + row[column] = 0 + new_rows_by_time[pulse.time_us] = row + new_rows_by_time[pulse.time_us][pulse.column] = 1 + + if new_rows_by_time: + new_rows = list(new_rows_by_time.values()) + df = pd.concat( + [df, pd.DataFrame(new_rows, columns=df.columns)], ignore_index=True + ) + + return df.sort_values("time_us", kind="stable").reset_index(drop=True) + + +def add_interpolated_angles( + df: pd.DataFrame, points: list[EncoderPoint] +) -> pd.DataFrame: + for index, row in df.iterrows(): + needs_angle = bool(sum(int(row[column] or 0) for column in SPARK_COLUMNS)) + if not needs_angle: + continue + + continuous_angle = interpolate_angle(points, int(row["time_us"])) + if continuous_angle is None: + continue + + revolution_index = int(continuous_angle // 360.0) + crank_angle = continuous_angle - revolution_index * 360.0 + df.at[index, "algo1_continuous_angle_deg"] = round(continuous_angle, 6) + df.at[index, "algo1_revolution_index"] = revolution_index + df.at[index, "algo1_crank_angle_deg"] = round(crank_angle, 6) + + return df + + +def output_path_for(input_path: Path, output_path: Path | None) -> Path: + if output_path is not None: + return output_path.expanduser() + + stem = input_path.stem + if stem.endswith("_rpm"): + stem = stem[: -len("_rpm")] + return input_path.with_name(f"{stem}_algo1.csv") + + +def parse_args() -> argparse.Namespace: + parser = argparse.ArgumentParser( + description="Run algo1 crank/cam angle prediction against one RPM dataset." + ) + parser.add_argument( + "file", + type=Path, + help="RPM CSV file. Bare names are searched in recorder RPM folders.", + ) + parser.add_argument("-o", "--output", type=Path, help="Output CSV path.") + parser.add_argument( + "--crank-gp", + choices=GP_COLUMNS, + default="gp0_falling", + help="GP falling column to treat as crank pulses.", + ) + parser.add_argument( + "--cam-gp", + choices=GP_COLUMNS, + default="gp1_falling", + help="GP falling column to treat as cam pulses.", + ) + parser.add_argument( + "--spark-advance", + type=float, + default=25.0, + help="Spark advance in crank degrees, matching the original algo1 constant.", + ) + return parser.parse_args() + + +def main() -> int: + args = parse_args() + if args.crank_gp == args.cam_gp: + print( + "--crank-gp and --cam-gp must select different GP columns.", file=sys.stderr + ) + return 1 + + try: + input_path = resolve_input_path(args.file) + df = read_rpm(input_path) + except (OSError, ValueError) as exc: + print(f"Could not read input: {exc}", file=sys.stderr) + return 1 + + crank_times = gp_times(df, args.crank_gp) + cam_times = gp_times(df, args.cam_gp) + sparks = generate_sparks( + crank_times, + cam_times, + spark_advance_deg=args.spark_advance, + ) + + output_df = add_spark_rows(df.copy(), sparks) + output_df = add_interpolated_angles(output_df, encoder_points(df)) + + output_path = output_path_for(input_path, args.output) + output_path.parent.mkdir(parents=True, exist_ok=True) + output_df.to_csv(output_path, index=False) + + print( + f"Wrote {output_path} " + f"crank={args.crank_gp}({len(crank_times)}) " + f"cam={args.cam_gp}({len(cam_times)}) " + f"sparks={len(sparks)}" + ) + return 0 + + +if __name__ == "__main__": + sys.exit(main()) diff --git a/src/plot.py b/src/plot.py index a2b647e..3597f36 100644 --- a/src/plot.py +++ b/src/plot.py @@ -2,73 +2,380 @@ # Copyright (C) 2026 Association Exergie # SPDX-License-Identifier: GPL-3.0-or-later -import pandas as pd -from pathlib import Path +from __future__ import annotations + import argparse -from tqdm import tqdm +import sys +from bisect import bisect_left +from dataclasses import dataclass +from pathlib import Path + import matplotlib.pyplot as plt +import pandas as pd +from matplotlib.ticker import MultipleLocator -parser = argparse.ArgumentParser() -parser.add_argument("directory", type=Path, help="Source data directory") -args = parser.parse_args() +GP_COLUMNS = ("gp0_falling", "gp1_falling") +SPARK_COLUMNS = ("spark_1", "spark_2", "spark_3") +PULSE_STYLES = { + "gp0_falling": ("gp0 falling", "tab:green", 5.0, 0.18), + "gp1_falling": ("gp1 falling", "tab:red", 5.0, 0.18), + "spark_1": ("spark 1", "tab:blue", 2.4, 0.65), + "spark_2": ("spark 2", "tab:orange", 2.4, 0.65), + "spark_3": ("spark 3", "tab:purple", 2.4, 0.65), +} -directory: Path = args.directory -if not directory.is_dir(): - parser.error(f"{directory} is not a valid directory") +@dataclass(frozen=True) +class EncoderPoint: + time_us: int + continuous_angle_deg: float -print(f"Processing data in: {directory}") -files: list[Path] = [] +def read_csv(path: Path) -> pd.DataFrame: + df = pd.read_csv(path) + if "time_us" not in df.columns: + raise ValueError("missing column: time_us") + df["time_us"] = df["time_us"].astype("int64") + return df -for path in directory.glob("*.csv"): - stem = path.stem +def encoder_points(df: pd.DataFrame) -> list[EncoderPoint]: + required = {"time_us", "revolution_index", "crank_angle_deg"} + if not required <= set(df.columns): + return [] + + points: list[EncoderPoint] = [] + angle_df = df.dropna(subset=["revolution_index", "crank_angle_deg"]) + for row in angle_df.itertuples(): + revolution_index = int(row.revolution_index) + if revolution_index < 0: + continue + continuous_angle = revolution_index * 360.0 + float(row.crank_angle_deg) + points.append(EncoderPoint(int(row.time_us), continuous_angle)) + + deduped: dict[int, EncoderPoint] = {} + for point in sorted(points, key=lambda item: item.time_us): + deduped[point.time_us] = point + return list(deduped.values()) + + +def interpolate_angle(points: list[EncoderPoint], time_us: int) -> float | None: + if not points: + return None + + times = [point.time_us for point in points] + index = bisect_left(times, time_us) + if index < len(points) and points[index].time_us == time_us: + return points[index].continuous_angle_deg + if index == 0 or index == len(points): + return None + + before = points[index - 1] + after = points[index] + span = after.time_us - before.time_us + if span <= 0: + return before.continuous_angle_deg + + fraction = (time_us - before.time_us) / span + return before.continuous_angle_deg + fraction * ( + after.continuous_angle_deg - before.continuous_angle_deg + ) + + +def iter_csv_paths(path: Path) -> list[Path]: + if path.is_dir(): + return sorted(path.rglob("*.csv")) + return [path] + + +def has_angle_data(points: list[EncoderPoint]) -> bool: + return len(points) >= 2 + + +def rpm_rows(df: pd.DataFrame) -> pd.DataFrame: + if "rpm_raw" not in df.columns or "rpm_lpf" not in df.columns: + return pd.DataFrame() + return df.dropna(subset=["rpm_raw", "rpm_lpf"]) + + +def plot_time(df: pd.DataFrame, ax: plt.Axes) -> None: + rpm_df = rpm_rows(df) + if not rpm_df.empty: + time_s = rpm_df["time_us"] / 1_000_000.0 + ax.plot(time_s, rpm_df["rpm_raw"], label="raw", linewidth=0.75, alpha=0.35) + ax.plot(time_s, rpm_df["rpm_lpf"], label="lpf", linewidth=1.4) + ax.set_xlabel("time (s)") + + +def plot_angle_overlay(df: pd.DataFrame, ax: plt.Axes) -> None: + rpm_df = rpm_rows(df) + if rpm_df.empty: + return + + raw_label_used = False + lpf_label_used = False + for revolution in sorted(rpm_df["revolution_index"].dropna().astype(int).unique()): + if revolution < 0: + continue + rev_df = rpm_df[rpm_df["revolution_index"].astype(int) == revolution] + ax.plot( + rev_df["crank_angle_deg"], + rev_df["rpm_raw"], + color="0.55", + linewidth=0.5, + alpha=0.12, + label="raw" if not raw_label_used else None, + ) + ax.plot( + rev_df["crank_angle_deg"], + rev_df["rpm_lpf"], + linewidth=0.8, + alpha=0.28, + label="lpf per rev" if not lpf_label_used else None, + ) + raw_label_used = True + lpf_label_used = True + + ax.set_xlim(0, 360) + ax.set_xlabel("crank angle (deg)") + + +def plot_angle_continuous(df: pd.DataFrame, ax: plt.Axes) -> None: + rpm_df = rpm_rows(df) + if rpm_df.empty: + return + + angle_df = rpm_df.dropna(subset=["revolution_index", "crank_angle_deg"]) + angle_df = angle_df[angle_df["revolution_index"].astype(int) >= 0] + if angle_df.empty: + return + + continuous_angle = ( + angle_df["revolution_index"].astype(float) * 360.0 + + angle_df["crank_angle_deg"].astype(float) + ) + ax.plot(continuous_angle, angle_df["rpm_raw"], label="raw", linewidth=0.75, alpha=0.35) + ax.plot(continuous_angle, angle_df["rpm_lpf"], label="lpf", linewidth=1.4) + + max_angle = float(continuous_angle.max()) + turn_angle = 0.0 + while turn_angle <= max_angle: + ax.axvline(turn_angle, color="0.7", linewidth=0.5, alpha=0.35) + turn_angle += 360.0 + + ax.set_xlim(0, max_angle) + ax.xaxis.set_major_locator(MultipleLocator(360.0)) + ax.xaxis.set_minor_locator(MultipleLocator(180.0)) + ax.set_xlabel("continuous crank angle (deg)") + + +def selected_pulse_columns( + df: pd.DataFrame, + gp_mode: str, + spark_columns: list[str], +) -> list[str]: + columns: list[str] = [] + if gp_mode in {"gp0", "both"} and "gp0_falling" in df.columns: + columns.append("gp0_falling") + if gp_mode in {"gp1", "both"} and "gp1_falling" in df.columns: + columns.append("gp1_falling") + columns.extend(column for column in spark_columns if column in df.columns) + return columns + + +def pulse_x_value( + row: pd.Series, + *, + x_axis: str, + points: list[EncoderPoint], +) -> float | None: + time_us = int(row["time_us"]) + if x_axis == "time": + return time_us / 1_000_000.0 + + continuous_angle = interpolate_angle(points, time_us) + if continuous_angle is None: + return None + if x_axis in {"angle", "angle-overlay"}: + return continuous_angle % 360.0 + return continuous_angle + + +def plot_pulse_bars( + df: pd.DataFrame, + ax: plt.Axes, + *, + x_axis: str, + gp_mode: str, + spark_columns: list[str], + points: list[EncoderPoint], +) -> None: + for column in selected_pulse_columns(df, gp_mode, spark_columns): + label, color, linewidth, alpha = PULSE_STYLES[column] + pulse_df = df.loc[df[column].fillna(0).astype(int) == 1] + plotted_label = False + for _, row in pulse_df.iterrows(): + x_value = pulse_x_value(row, x_axis=x_axis, points=points) + if x_value is None: + continue + ax.axvline( + x_value, + color=color, + alpha=alpha, + linewidth=linewidth, + label=label if not plotted_label else None, + ) + plotted_label = True + + +def plot_file( + csv_path: Path, + output: Path | None, + *, + x_axis: str, + gp_mode: str, + spark_columns: list[str], +) -> bool: try: - base_name, channel = stem.rsplit("_", 1) - except ValueError: - print(f"Skipping badly named file: {path}") - continue + df = read_csv(csv_path) + except (OSError, ValueError) as exc: + print(f"Could not read {csv_path}: {exc}", file=sys.stderr) + return False - if channel != "output": - print(f"Skipping unknown file: {path}") - continue + points = encoder_points(df) + if x_axis.startswith("angle") and not has_angle_data(points): + print(f"Could not plot crank angle for {csv_path}: missing angle columns", file=sys.stderr) + return False - files.append(path) + fig, ax = plt.subplots(figsize=(12, 5)) + if x_axis in {"angle", "angle-overlay"}: + plot_angle_overlay(df, ax) + elif x_axis == "angle-continuous": + plot_angle_continuous(df, ax) + else: + plot_time(df, ax) - -for file in files: - output_df = pd.read_csv(file).set_index("time_us", drop=False) - fig, ax = plt.subplots(figsize=(12, 6)) - - # RPM - ax.plot(output_df["time_us"] / 1_000_000, output_df["crank"], label="crank") - ax.plot( - output_df["time_us"] / 1_000_000, output_df["cam"], label="cam", color="red" + plot_pulse_bars( + df, + ax, + x_axis=x_axis, + gp_mode=gp_mode, + spark_columns=spark_columns, + points=points, ) - ax.plot( - output_df["time_us"] / 1_000_000, - output_df["spark_1"], - label="spark_1", - color="green", - ) - ax.plot( - output_df["time_us"] / 1_000_000, - output_df["spark_2"], - label="spark_2", - color="yellow", - ) - ax.plot( - output_df["time_us"] / 1_000_000, - output_df["spark_3"], - label="spark_3", - color="orange", - ) - ax.set_xlabel("Time (s)") - ax.set_ylim(0, 1) - ax.grid(True) + ax.set_ylabel("RPM") + ax.set_title(str(csv_path)) + ax.grid(alpha=0.25) + if x_axis == "angle-continuous": + ax.grid(which="minor", axis="x", alpha=0.15, linewidth=0.6) + ax.legend() fig.tight_layout() - plt.show() - plt.close(fig) + + if output is not None: + output.parent.mkdir(parents=True, exist_ok=True) + fig.savefig(output, dpi=160) + print(f"Wrote {output}") + plt.close(fig) + else: + plt.show() + plt.close(fig) + + return True + + +def output_path_for(csv_path: Path, output: Path | None, multiple: bool) -> Path | None: + if output is None: + return None + if not multiple: + return output + return output / f"{csv_path.stem}.png" + + +def parse_spark_columns(value: str) -> list[str]: + if value == "none": + return [] + if value == "all": + return list(SPARK_COLUMNS) + + columns = [column.strip() for column in value.split(",") if column.strip()] + if not columns: + raise argparse.ArgumentTypeError("must specify at least one spark column") + + invalid = sorted(set(columns) - set(SPARK_COLUMNS)) + if invalid: + valid = ", ".join(("none", "all", *SPARK_COLUMNS)) + raise argparse.ArgumentTypeError( + f"invalid spark column(s): {', '.join(invalid)}; choose from {valid}" + ) + + deduped: list[str] = [] + for column in columns: + if column not in deduped: + deduped.append(column) + return deduped + + +def parse_args() -> argparse.Namespace: + parser = argparse.ArgumentParser(description="Plot algo1 RPM output CSV files.") + parser.add_argument("csv_path", type=Path, help="Algo1 CSV file or directory of CSV files.") + parser.add_argument("-o", "--output", type=Path, help="Optional image output path or directory.") + parser.add_argument( + "--x", + choices=("time", "angle", "angle-overlay", "angle-continuous"), + default="time", + help="Plot RPM against time, overlaid 0-360 crank angle, or continuous crank angle.", + ) + parser.add_argument( + "--gp", + choices=("none", "gp0", "gp1", "both"), + default="both", + help="Overlay GP falling events as vertical bars.", + ) + parser.add_argument( + "--sparks", + type=parse_spark_columns, + default="all", + help="Overlay spark events: none, all, one spark, or a comma-separated subset.", + ) + return parser.parse_args() + + +def main() -> int: + args = parse_args() + + if not args.csv_path.exists(): + print(f"Path does not exist: {args.csv_path}", file=sys.stderr) + return 1 + + csv_paths = iter_csv_paths(args.csv_path) + if not csv_paths: + print(f"No CSV files found under {args.csv_path}", file=sys.stderr) + return 1 + + multiple = len(csv_paths) > 1 + if multiple and args.output is not None and args.output.suffix: + print("When plotting a directory, --output must be a directory.", file=sys.stderr) + return 1 + + spark_columns = ( + args.sparks if isinstance(args.sparks, list) else parse_spark_columns(args.sparks) + ) + failures = 0 + for csv_path in csv_paths: + output = output_path_for(csv_path, args.output, multiple) + if not plot_file( + csv_path, + output, + x_axis=args.x, + gp_mode=args.gp, + spark_columns=spark_columns, + ): + failures += 1 + + return 1 if failures else 0 + + +if __name__ == "__main__": + sys.exit(main()) diff --git a/src/spark_error.py b/src/spark_error.py new file mode 100644 index 0000000..b78bf74 --- /dev/null +++ b/src/spark_error.py @@ -0,0 +1,837 @@ +# Copyright (C) 2026 Hector van der Aa +# Copyright (C) 2026 Association Exergie +# SPDX-License-Identifier: GPL-3.0-or-later + +from __future__ import annotations + +import argparse +import math +import statistics +import sys +from bisect import bisect_left, bisect_right +from dataclasses import dataclass +from pathlib import Path + +import matplotlib.pyplot as plt +import pandas as pd + + +GP_COLUMNS = ("gp0_falling", "gp1_falling") +SPARK_COLUMNS = ("spark_1", "spark_2", "spark_3") +DEFAULT_PPR = 256 +CURVE_STYLES = { + "spark_1": {"linestyle": (0, (1.0, 1.2)), "linewidth": 1.7}, + "spark_2": {"linestyle": (0, (1.0, 2.2)), "linewidth": 2.3}, + "spark_3": {"linestyle": (0, (1.0, 3.2)), "linewidth": 2.9}, +} + + +@dataclass(frozen=True) +class EncoderPoint: + time_us: int + continuous_angle_deg: float + + +@dataclass(frozen=True) +class SparkError: + spark_column: str + spark_time_us: int + reference_crank_time_us: int + expected_delta_deg: float + actual_delta_deg: float + error_deg: float + + +@dataclass(frozen=True) +class AngleEvent: + column: str + time_us: int + continuous_angle_deg: float + + +@dataclass(frozen=True) +class ErrorStats: + count: int + minimum: float + maximum: float + mean: float + median: float + stdev: float + rms: float + abs_mean: float + abs_max: float + + +def read_algo1(path: Path) -> pd.DataFrame: + df = pd.read_csv(path) + missing = {"time_us", "revolution_index", "crank_angle_deg"} - set(df.columns) + if missing: + raise ValueError(f"missing columns: {', '.join(sorted(missing))}") + + spark_missing = set(SPARK_COLUMNS) - set(df.columns) + if spark_missing: + raise ValueError(f"missing spark columns: {', '.join(sorted(spark_missing))}") + + df["time_us"] = df["time_us"].astype("int64") + return df.sort_values("time_us", kind="stable").reset_index(drop=True) + + +def encoder_points(df: pd.DataFrame) -> list[EncoderPoint]: + angle_df = df.dropna(subset=["revolution_index", "crank_angle_deg"]) + points: list[EncoderPoint] = [] + + for row in angle_df.itertuples(): + revolution_index = int(row.revolution_index) + if revolution_index < 0: + continue + continuous_angle = revolution_index * 360.0 + float(row.crank_angle_deg) + points.append(EncoderPoint(int(row.time_us), continuous_angle)) + + deduped: dict[int, EncoderPoint] = {} + for point in sorted(points, key=lambda item: item.time_us): + deduped[point.time_us] = point + return list(deduped.values()) + + +def interpolate_angle(points: list[EncoderPoint], time_us: int) -> float | None: + if not points: + return None + + times = [point.time_us for point in points] + index = bisect_left(times, time_us) + if index < len(points) and points[index].time_us == time_us: + return points[index].continuous_angle_deg + if index == 0 or index == len(points): + return None + + before = points[index - 1] + after = points[index] + span_us = after.time_us - before.time_us + if span_us <= 0: + return before.continuous_angle_deg + + fraction = (time_us - before.time_us) / span_us + return before.continuous_angle_deg + fraction * ( + after.continuous_angle_deg - before.continuous_angle_deg + ) + + +def pulse_times(df: pd.DataFrame, column: str) -> list[int]: + pulse_rows = df.loc[df[column].fillna(0).astype(int) == 1] + return sorted(int(row.time_us) for row in pulse_rows.itertuples()) + + +def spark_events(df: pd.DataFrame, spark_columns: list[str]) -> list[tuple[int, str]]: + events: list[tuple[int, str]] = [] + for column in spark_columns: + for time_us in pulse_times(df, column): + events.append((time_us, column)) + return sorted(events) + + +def previous_time(times: list[int], time_us: int) -> int | None: + index = bisect_right(times, time_us) - 1 + if index < 0: + return None + return times[index] + + +def wrap_signed_180(angle_deg: float) -> float: + return (angle_deg + 180.0) % 360.0 - 180.0 + + +def calculate_errors( + df: pd.DataFrame, + *, + crank_gp: str, + spark_columns: list[str], + spark_advance_deg: float, +) -> list[SparkError]: + points = encoder_points(df) + crank_times = pulse_times(df, crank_gp) + expected_delta_deg = 45.0 - spark_advance_deg + errors: list[SparkError] = [] + + for spark_time_us, spark_column in spark_events(df, spark_columns): + reference_time_us = previous_time(crank_times, spark_time_us) + if reference_time_us is None: + continue + + reference_angle = interpolate_angle(points, reference_time_us) + spark_angle = interpolate_angle(points, spark_time_us) + if reference_angle is None or spark_angle is None: + continue + + actual_delta_deg = spark_angle - reference_angle + error_deg = wrap_signed_180(actual_delta_deg - expected_delta_deg) + errors.append( + SparkError( + spark_column=spark_column, + spark_time_us=spark_time_us, + reference_crank_time_us=reference_time_us, + expected_delta_deg=expected_delta_deg, + actual_delta_deg=actual_delta_deg, + error_deg=error_deg, + ) + ) + + return errors + + +def stats_for(values: list[float]) -> ErrorStats | None: + if not values: + return None + + return ErrorStats( + count=len(values), + minimum=min(values), + maximum=max(values), + mean=statistics.fmean(values), + median=statistics.median(values), + stdev=statistics.stdev(values) if len(values) >= 2 else 0.0, + rms=math.sqrt(statistics.fmean(value * value for value in values)), + abs_mean=statistics.fmean(abs(value) for value in values), + abs_max=max(abs(value) for value in values), + ) + + +def format_stats(label: str, stats: ErrorStats | None) -> list[str]: + if stats is None: + return [f"{label}: no valid spark samples"] + + return [ + f"{label}:", + f" count: {stats.count}", + f" min error deg: {stats.minimum:.6f}", + f" max error deg: {stats.maximum:.6f}", + f" mean error deg: {stats.mean:.6f}", + f" median error deg: {stats.median:.6f}", + f" stdev error deg: {stats.stdev:.6f}", + f" rms error deg: {stats.rms:.6f}", + f" mean abs error deg: {stats.abs_mean:.6f}", + f" max abs error deg: {stats.abs_max:.6f}", + ] + + +def build_report( + input_path: Path, + errors: list[SparkError], + *, + crank_gp: str, + spark_columns: list[str], + spark_advance_deg: float, + ppr: int, + band_size_deg: float, +) -> str: + expected_delta_deg = 45.0 - spark_advance_deg + encoder_step_deg = 360.0 / ppr + lines = [ + f"input: {input_path}", + f"crank reference GP: {crank_gp}", + f"spark columns: {', '.join(spark_columns)}", + f"spark advance deg: {spark_advance_deg:.6f}", + f"expected angle after 45 BTDC crank pulse deg: {expected_delta_deg:.6f}", + f"encoder ppr: {ppr}", + f"encoder step deg: {encoder_step_deg:.6f}", + "", + ] + + all_errors = [error.error_deg for error in errors] + lines.extend(format_stats("all sparks", stats_for(all_errors))) + + for column in spark_columns: + column_errors = [ + error.error_deg for error in errors if error.spark_column == column + ] + lines.append("") + lines.extend(format_stats(column, stats_for(column_errors))) + + lines.append("") + lines.append("samples:") + lines.append( + "spark,time_us,reference_crank_time_us,error_band,expected_delta_deg," + "actual_delta_deg,error_deg" + ) + for error in errors: + lines.append( + f"{error.spark_column},{error.spark_time_us}," + f"{error.reference_crank_time_us}," + f"{band_label(error.error_deg, band_size_deg)}," + f"{error.expected_delta_deg:.6f}," + f"{error.actual_delta_deg:.6f}," + f"{error.error_deg:.6f}" + ) + + return "\n".join(lines) + "\n" + + +def plot_distribution( + errors: list[SparkError], + output_path: Path, + *, + spark_columns: list[str], + band_size_deg: float, +) -> None: + fig, axes = plt.subplots( + len(spark_columns), + 1, + figsize=(10, max(3.0, 2.5 * len(spark_columns))), + sharex=True, + ) + if len(spark_columns) == 1: + axes = [axes] + + all_values = [error.error_deg for error in errors] + if all_values: + minimum = math.floor(min(all_values) / band_size_deg) * band_size_deg + maximum = math.ceil(max(all_values) / band_size_deg) * band_size_deg + bins = [] + edge = minimum + while edge <= maximum + band_size_deg * 0.5: + bins.append(edge) + edge += band_size_deg + if len(bins) < 2: + bins.append(minimum + band_size_deg) + else: + bins = 10 + + for ax, column in zip(axes, spark_columns): + values = [error.error_deg for error in errors if error.spark_column == column] + if values: + ax.hist( + values, + bins=bins, + density=True, + alpha=0.75, + color="tab:blue", + edgecolor="black", + ) + else: + ax.text( + 0.5, + 0.5, + "No valid samples", + ha="center", + va="center", + transform=ax.transAxes, + ) + + ax.axvline(0.0, color="black", linewidth=1.0, alpha=0.65) + ax.text( + 0.0, + 0.92, + "target", + ha="center", + va="top", + transform=ax.get_xaxis_transform(), + fontsize=8, + ) + ax.set_ylabel("density") + ax.set_title(column) + ax.grid(alpha=0.25) + + axes[-1].set_xlabel("angular error (deg)") + fig.suptitle(f"Spark angular error distribution ({band_size_deg:g} deg bins)") + fig.tight_layout() + output_path.parent.mkdir(parents=True, exist_ok=True) + fig.savefig(output_path, dpi=160) + plt.close(fig) + + +def density_curve(values: list[float], points: int = 240) -> tuple[list[float], list[float]]: + minimum = min(values) + maximum = max(values) + count = len(values) + + if count >= 2: + stdev = statistics.stdev(values) + else: + stdev = 0.0 + + span = maximum - minimum + bandwidth = 1.06 * stdev * count ** (-1.0 / 5.0) if stdev > 0.0 else 0.0 + if bandwidth <= 0.0: + bandwidth = max(span / 8.0, 1.0) + + margin = max(span * 0.15, bandwidth * 3.0, 1.0) + start = minimum - margin + stop = maximum + margin + step = (stop - start) / (points - 1) + + x_values = [start + index * step for index in range(points)] + normalizer = count * bandwidth * math.sqrt(2.0 * math.pi) + y_values = [ + sum( + math.exp(-0.5 * ((x_value - value) / bandwidth) ** 2) + for value in values + ) + / normalizer + for x_value in x_values + ] + return x_values, y_values + + +def band_edges(error_deg: float, band_size_deg: float = 1.0) -> tuple[float, float]: + lower = math.floor(error_deg / band_size_deg) * band_size_deg + return lower, lower + band_size_deg + + +def format_band_value(value: float) -> str: + if value == 0: + return "0" + if value.is_integer(): + return f"{int(value):+d}" + return f"{value:+.3f}".rstrip("0").rstrip(".") + + +def band_label(error_deg: float, band_size_deg: float = 1.0) -> str: + lower, upper = band_edges(error_deg, band_size_deg) + return f"{format_band_value(lower)}_to_{format_band_value(upper)}_deg" + + +def band_directory_name(error_deg: float, band_size_deg: float = 1.0) -> str: + return band_label(error_deg, band_size_deg).replace("+", "pos").replace("-", "neg") + + +def angle_events( + df: pd.DataFrame, + points: list[EncoderPoint], + columns: list[str], +) -> list[AngleEvent]: + events: list[AngleEvent] = [] + for column in columns: + if column not in df.columns: + continue + for time_us in pulse_times(df, column): + angle = interpolate_angle(points, time_us) + if angle is None: + continue + events.append(AngleEvent(column, time_us, angle)) + return sorted(events, key=lambda event: (event.continuous_angle_deg, event.column)) + + +def cycle_rows( + df: pd.DataFrame, + center_angle_deg: float, + before_deg: float, + after_deg: float, +) -> pd.DataFrame: + if "rpm_raw" not in df.columns or "rpm_lpf" not in df.columns: + return pd.DataFrame() + + rows = df.dropna(subset=["revolution_index", "crank_angle_deg", "rpm_raw", "rpm_lpf"]).copy() + if rows.empty: + return rows + + rows = rows[rows["revolution_index"].astype(int) >= 0] + rows["continuous_angle_deg"] = ( + rows["revolution_index"].astype(float) * 360.0 + + rows["crank_angle_deg"].astype(float) + ) + rows["relative_angle_deg"] = rows["continuous_angle_deg"] - center_angle_deg + return rows[ + (rows["relative_angle_deg"] >= -before_deg) + & (rows["relative_angle_deg"] <= after_deg) + ] + + +def plot_cycle( + df: pd.DataFrame, + points: list[EncoderPoint], + error: SparkError, + output_path: Path, + *, + before_deg: float, + after_deg: float, + band_size_deg: float, +) -> bool: + spark_angle = interpolate_angle(points, error.spark_time_us) + if spark_angle is None: + return False + + rows = cycle_rows(df, spark_angle, before_deg, after_deg) + events = angle_events(df, points, [*GP_COLUMNS, *SPARK_COLUMNS]) + window_events = [ + event + for event in events + if -before_deg <= event.continuous_angle_deg - spark_angle <= after_deg + ] + + fig, (rpm_ax, angle_ax, pulse_ax) = plt.subplots( + 3, + 1, + figsize=(12, 8), + sharex=True, + gridspec_kw={"height_ratios": [2.0, 1.2, 1.0]}, + ) + + if not rows.empty: + rpm_ax.plot( + rows["relative_angle_deg"], + rows["rpm_raw"], + label="rpm raw", + linewidth=0.75, + alpha=0.35, + ) + rpm_ax.plot( + rows["relative_angle_deg"], + rows["rpm_lpf"], + label="rpm lpf", + linewidth=1.4, + ) + angle_ax.plot( + rows["relative_angle_deg"], + rows["crank_angle_deg"], + label="crank angle", + color="tab:orange", + linewidth=1.2, + ) + + for event in window_events: + relative_angle = event.continuous_angle_deg - spark_angle + if event.column == "gp0_falling": + pulse_ax.vlines(relative_angle, 0.0, 0.8, color="tab:green", linewidth=1.6) + elif event.column == "gp1_falling": + pulse_ax.vlines(relative_angle, 1.0, 1.8, color="tab:red", linewidth=1.6) + elif event.time_us == error.spark_time_us and event.column == error.spark_column: + pulse_ax.vlines(relative_angle, 2.0, 3.0, color="black", linewidth=2.8) + else: + pulse_ax.vlines(relative_angle, 2.0, 2.8, color="tab:blue", linewidth=1.4, alpha=0.6) + + for ax in (rpm_ax, angle_ax, pulse_ax): + ax.axvline(0.0, color="black", linewidth=1.0, alpha=0.65) + ax.grid(alpha=0.25) + + reference_angle = interpolate_angle(points, error.reference_crank_time_us) + if reference_angle is not None: + reference_relative = reference_angle - spark_angle + target_relative = reference_relative + error.expected_delta_deg + for ax in (rpm_ax, angle_ax, pulse_ax): + ax.axvline(reference_relative, color="0.35", linestyle="--", linewidth=1.0, alpha=0.65) + ax.axvline(target_relative, color="tab:purple", linestyle=":", linewidth=2.0, alpha=0.85) + pulse_ax.text( + target_relative, + 3.05, + "target", + color="tab:purple", + ha="center", + va="bottom", + fontsize=8, + ) + + rpm_ax.set_ylabel("RPM") + rpm_ax.legend(loc="upper right") + angle_ax.set_ylabel("crank angle (deg)") + angle_ax.set_ylim(0, 360) + angle_ax.legend(loc="upper right") + pulse_ax.set_ylabel("pulses") + pulse_ax.set_yticks([0.4, 1.4, 2.5]) + pulse_ax.set_yticklabels(["gp0", "gp1", "spark"]) + pulse_ax.set_ylim(-0.2, 3.2) + pulse_ax.set_xlabel("crank angle relative to analyzed spark (deg)") + pulse_ax.set_xlim(-before_deg, after_deg) + fig.suptitle( + f"{error.spark_column} at {error.spark_time_us} us " + f"error={error.error_deg:.3f} deg " + f"band={band_label(error.error_deg, band_size_deg)}" + ) + fig.tight_layout() + output_path.parent.mkdir(parents=True, exist_ok=True) + fig.savefig(output_path, dpi=160) + plt.close(fig) + return True + + +def write_cycle_plots( + df: pd.DataFrame, + errors: list[SparkError], + output_root: Path, + *, + before_deg: float, + after_deg: float, + band_size_deg: float, +) -> int: + points = encoder_points(df) + written = 0 + errors_by_band: dict[str, list[SparkError]] = {} + + for index, error in enumerate(errors, start=1): + band_name = band_directory_name(error.error_deg, band_size_deg) + errors_by_band.setdefault(band_name, []).append(error) + directory = output_root / band_name + output_path = directory / ( + f"{index:04d}_{error.spark_column}_{error.spark_time_us}_" + f"err_{error.error_deg:+.3f}.png" + ) + if not plot_cycle( + df, + points, + error, + output_path, + before_deg=before_deg, + after_deg=after_deg, + band_size_deg=band_size_deg, + ): + continue + written += 1 + + for band_name, band_errors in errors_by_band.items(): + output_path = output_root / band_name / "_stacked_cycles.png" + if plot_stacked_cycles( + df, + points, + band_errors, + output_path, + before_deg=before_deg, + after_deg=after_deg, + band_name=band_name, + ): + written += 1 + + return written + + +def plot_stacked_cycles( + df: pd.DataFrame, + points: list[EncoderPoint], + errors: list[SparkError], + output_path: Path, + *, + before_deg: float, + after_deg: float, + band_name: str, +) -> bool: + if not errors: + return False + + fig, (rpm_ax, angle_ax) = plt.subplots( + 2, + 1, + figsize=(12, 7), + sharex=True, + gridspec_kw={"height_ratios": [2.0, 1.2]}, + ) + plotted = False + + for error in errors: + spark_angle = interpolate_angle(points, error.spark_time_us) + reference_angle = interpolate_angle(points, error.reference_crank_time_us) + if spark_angle is None or reference_angle is None: + continue + + rows = cycle_rows(df, spark_angle, before_deg, after_deg) + if rows.empty: + continue + + rpm_ax.plot( + rows["relative_angle_deg"], + rows["rpm_raw"], + color="0.45", + linewidth=0.55, + alpha=0.18, + ) + rpm_ax.plot( + rows["relative_angle_deg"], + rows["rpm_lpf"], + color="tab:blue", + linewidth=0.8, + alpha=0.28, + ) + angle_ax.plot( + rows["relative_angle_deg"], + rows["crank_angle_deg"], + color="tab:orange", + linewidth=0.7, + alpha=0.22, + ) + target_relative = reference_angle + error.expected_delta_deg - spark_angle + for ax in (rpm_ax, angle_ax): + ax.axvline(target_relative, color="tab:purple", linestyle=":", linewidth=0.8, alpha=0.18) + plotted = True + + if not plotted: + plt.close(fig) + return False + + for ax in (rpm_ax, angle_ax): + ax.axvline(0.0, color="black", linewidth=1.0, alpha=0.65) + ax.grid(alpha=0.25) + + rpm_ax.plot([], [], color="0.45", linewidth=0.8, alpha=0.45, label="rpm raw") + rpm_ax.plot([], [], color="tab:blue", linewidth=1.0, alpha=0.55, label="rpm lpf") + rpm_ax.axvline(float("nan"), color="tab:purple", linestyle=":", linewidth=1.2, label="target") + rpm_ax.set_ylabel("RPM") + rpm_ax.legend(loc="upper right") + angle_ax.set_ylabel("crank angle (deg)") + angle_ax.set_ylim(0, 360) + angle_ax.set_xlabel("crank angle relative to analyzed spark (deg)") + angle_ax.set_xlim(-before_deg, after_deg) + fig.suptitle(f"Stacked cycles for {band_name} ({len(errors)} samples)") + fig.tight_layout() + output_path.parent.mkdir(parents=True, exist_ok=True) + fig.savefig(output_path, dpi=160) + plt.close(fig) + return True + + +def default_output_path(input_path: Path, suffix: str) -> Path: + stem = input_path.stem + return input_path.with_name(f"{stem}_{suffix}") + + +def parse_spark_columns(value: str) -> list[str]: + if value == "all": + return list(SPARK_COLUMNS) + + columns = [column.strip() for column in value.split(",") if column.strip()] + if not columns: + raise argparse.ArgumentTypeError("must specify at least one spark column") + + invalid = sorted(set(columns) - set(SPARK_COLUMNS)) + if invalid: + valid = ", ".join((*SPARK_COLUMNS, "all")) + raise argparse.ArgumentTypeError( + f"invalid spark column(s): {', '.join(invalid)}; choose from {valid}" + ) + + deduped: list[str] = [] + for column in columns: + if column not in deduped: + deduped.append(column) + return deduped + + +def parse_args() -> argparse.Namespace: + parser = argparse.ArgumentParser( + description="Measure algo1 spark angular error from an algo1 output CSV." + ) + parser.add_argument("file", type=Path, help="CSV generated by src/algo1.py.") + parser.add_argument( + "--crank-gp", + choices=GP_COLUMNS, + default="gp0_falling", + help="GP falling column used as the crank pulse reference.", + ) + parser.add_argument( + "--sparks", + type=parse_spark_columns, + default="all", + help="Spark columns to analyze: all, one spark, or a comma-separated subset.", + ) + parser.add_argument( + "--spark-advance", + type=float, + default=25.0, + help="Test spark advance in crank degrees before TDC.", + ) + parser.add_argument( + "--ppr", + type=int, + default=DEFAULT_PPR, + help="Rotary encoder pulses per revolution, used for report context.", + ) + parser.add_argument("-o", "--output", type=Path, help="Distribution plot path.") + parser.add_argument("--log", type=Path, help="Text report path.") + parser.add_argument( + "--cycle-plots", + type=Path, + help="Optional root folder for per-spark engine-cycle plots grouped by error band.", + ) + parser.add_argument( + "--cycle-before-deg", + type=float, + default=540.0, + help="Crank-angle range before the analyzed spark for each cycle plot.", + ) + parser.add_argument( + "--cycle-after-deg", + type=float, + default=180.0, + help="Crank-angle range after the analyzed spark for each cycle plot.", + ) + parser.add_argument( + "--band-size", + type=float, + default=1.0, + help="Angular error band size in degrees for categorized cycle plot folders.", + ) + return parser.parse_args() + + +def main() -> int: + args = parse_args() + input_path = args.file.expanduser() + if not input_path.exists(): + print(f"Path does not exist: {input_path}", file=sys.stderr) + return 1 + if args.band_size <= 0.0: + print("--band-size must be greater than zero.", file=sys.stderr) + return 1 + if args.cycle_before_deg <= 0.0: + print("--cycle-before-deg must be greater than zero.", file=sys.stderr) + return 1 + if args.cycle_after_deg <= 0.0: + print("--cycle-after-deg must be greater than zero.", file=sys.stderr) + return 1 + + spark_columns = ( + args.sparks if isinstance(args.sparks, list) else parse_spark_columns(args.sparks) + ) + + try: + df = read_algo1(input_path) + except (OSError, ValueError) as exc: + print(f"Could not read input: {exc}", file=sys.stderr) + return 1 + + if args.crank_gp not in df.columns: + print(f"Missing crank GP column: {args.crank_gp}", file=sys.stderr) + return 1 + + errors = calculate_errors( + df, + crank_gp=args.crank_gp, + spark_columns=spark_columns, + spark_advance_deg=args.spark_advance, + ) + report = build_report( + input_path, + errors, + crank_gp=args.crank_gp, + spark_columns=spark_columns, + spark_advance_deg=args.spark_advance, + ppr=args.ppr, + band_size_deg=args.band_size, + ) + + log_path = args.log.expanduser() if args.log else default_output_path(input_path, "spark_error.txt") + plot_path = args.output.expanduser() if args.output else default_output_path(input_path, "spark_error.png") + + log_path.parent.mkdir(parents=True, exist_ok=True) + log_path.write_text(report) + plot_distribution( + errors, + plot_path, + spark_columns=spark_columns, + band_size_deg=args.band_size, + ) + + print(report, end="") + print(f"Wrote {log_path}") + print(f"Wrote {plot_path}") + if args.cycle_plots is not None: + cycle_root = args.cycle_plots.expanduser() + written = write_cycle_plots( + df, + errors, + cycle_root, + before_deg=args.cycle_before_deg, + after_deg=args.cycle_after_deg, + band_size_deg=args.band_size, + ) + print(f"Wrote {written} categorized cycle plots under {cycle_root}") + return 0 + + +if __name__ == "__main__": + sys.exit(main())