diff --git a/.gitignore b/.gitignore new file mode 100644 index 0000000..b885316 --- /dev/null +++ b/.gitignore @@ -0,0 +1,7 @@ +# Python cache and bytecode +__pycache__/ +*.py[cod] + +# Generated example outputs +examples/**/output/* +!examples/**/output/.gitkeep diff --git a/examples/case01/input/input.txt b/examples/case01/input/input.txt new file mode 100644 index 0000000..a17557a --- /dev/null +++ b/examples/case01/input/input.txt @@ -0,0 +1,8 @@ +m: 1 +xa: 0 +ta: 0 +xb: 20 +tb: 10 +N: 100 +xmax: 1.0 +seed: 42 diff --git a/examples/case01/output/.gitkeep b/examples/case01/output/.gitkeep new file mode 100644 index 0000000..8b13789 --- /dev/null +++ b/examples/case01/output/.gitkeep @@ -0,0 +1 @@ + diff --git a/examples/case01/run_case.py b/examples/case01/run_case.py new file mode 100644 index 0000000..9b3fd33 --- /dev/null +++ b/examples/case01/run_case.py @@ -0,0 +1,20 @@ +from pathlib import Path +import subprocess +import sys + + +def main() -> None: + case_dir = Path(__file__).resolve().parent + root_dir = case_dir.parent.parent + input_path = case_dir / "input" / "input.txt" + output_path = case_dir / "output" / "result.txt" + program_path = root_dir / "least_action.py" + + subprocess.run( + [sys.executable, str(program_path), str(input_path), str(output_path)], + check=True, + ) + + +if __name__ == "__main__": + main() diff --git a/least_action.py b/least_action.py new file mode 100644 index 0000000..30a1d2d --- /dev/null +++ b/least_action.py @@ -0,0 +1,244 @@ +from __future__ import annotations + +import argparse +import random +from pathlib import Path +from typing import Dict, List + +from PIL import Image, ImageDraw, ImageFont +import yaml + + +def load_input(input_path: Path) -> Dict[str, float]: + with input_path.open("r", encoding="utf-8") as file: + data = yaml.safe_load(file) or {} + + defaults = { + "m": 1.0, + "xa": 0.0, + "ta": 0.0, + "xb": 20.0, + "tb": 10.0, + "N": 100, + "xmax": 1.0, + "seed": None, + } + defaults.update(data) + return defaults + + +def build_straight_path(xa: float, xb: float, n_points: int) -> List[float]: + if n_points < 2: + raise ValueError("N must be at least 2.") + + return [ + xa + (xb - xa) * step / (n_points - 1) + for step in range(n_points) + ] + + +def build_random_path(straight_path: List[float], xmax: float) -> List[float]: + if len(straight_path) < 2: + raise ValueError("N must be at least 2.") + + path = [straight_path[0]] + for base_x in straight_path[1:-1]: + path.append(base_x + random.uniform(-1.0, 1.0) * xmax) + path.append(straight_path[-1]) + return path + + +def kinetic_action(path: List[float], mass: float, dt: float) -> float: + action = 0.0 + for left, right in zip(path[:-1], path[1:]): + velocity = (right - left) / dt + kinetic_energy = 0.5 * mass * velocity * velocity + action += kinetic_energy * dt + return action + + +def build_time_points(ta: float, tb: float, n_points: int) -> List[float]: + if n_points < 2: + raise ValueError("N must be at least 2.") + + dt = (tb - ta) / (n_points - 1) + return [ta + step * dt for step in range(n_points)] + + +def draw_trajectory_image( + image_path: Path, + time_points: List[float], + random_path: List[float], + straight_path: List[float], +) -> None: + width = 1200 + height = 800 + margin_left = 110 + margin_right = 60 + margin_top = 80 + margin_bottom = 110 + + image = Image.new("RGB", (width, height), "white") + draw = ImageDraw.Draw(image) + + try: + title_font = ImageFont.truetype("DejaVuSans.ttf", 34) + label_font = ImageFont.truetype("DejaVuSans.ttf", 24) + tick_font = ImageFont.truetype("DejaVuSans.ttf", 20) + legend_font = ImageFont.truetype("DejaVuSans.ttf", 24) + except OSError: + title_font = ImageFont.load_default() + label_font = ImageFont.load_default() + tick_font = ImageFont.load_default() + legend_font = ImageFont.load_default() + + all_x = random_path + straight_path + t_min = min(time_points) + t_max = max(time_points) + x_min = min(all_x) + x_max = max(all_x) + + if t_min == t_max: + t_max = t_min + 1.0 + if x_min == x_max: + x_max = x_min + 1.0 + + x_padding = max((x_max - x_min) * 0.08, 1.0) + x_min -= x_padding + x_max += x_padding + + plot_left = margin_left + plot_right = width - margin_right + plot_top = margin_top + plot_bottom = height - margin_bottom + + def map_point(time_value: float, position_value: float) -> tuple[float, float]: + x_ratio = (time_value - t_min) / (t_max - t_min) + y_ratio = (position_value - x_min) / (x_max - x_min) + px = plot_left + x_ratio * (plot_right - plot_left) + py = plot_bottom - y_ratio * (plot_bottom - plot_top) + return px, py + + draw.rectangle([plot_left, plot_top, plot_right, plot_bottom], outline=(0, 0, 0), width=2) + + tick_count = 5 + for tick in range(tick_count + 1): + x_value = t_min + (t_max - t_min) * tick / tick_count + y_value = x_min + (x_max - x_min) * tick / tick_count + + tick_x = plot_left + (plot_right - plot_left) * tick / tick_count + tick_y = plot_bottom - (plot_bottom - plot_top) * tick / tick_count + + draw.line([(tick_x, plot_top), (tick_x, plot_bottom)], fill=(225, 225, 225), width=1) + draw.line([(plot_left, tick_y), (plot_right, tick_y)], fill=(225, 225, 225), width=1) + + draw.line([(tick_x, plot_bottom), (tick_x, plot_bottom + 8)], fill=(0, 0, 0), width=2) + draw.line([(plot_left - 8, tick_y), (plot_left, tick_y)], fill=(0, 0, 0), width=2) + + draw.text((tick_x - 20, plot_bottom + 15), f"{x_value:.1f}", fill=(0, 0, 0), font=tick_font) + draw.text((20, tick_y - 10), f"{y_value:.1f}", fill=(0, 0, 0), font=tick_font) + + random_points = [map_point(t, x) for t, x in zip(time_points, random_path)] + straight_points = [map_point(t, x) for t, x in zip(time_points, straight_path)] + + draw.line(straight_points, fill=(40, 110, 220), width=4) + draw.line(random_points, fill=(220, 70, 70), width=3) + + for px, py in random_points: + draw.ellipse([(px - 3, py - 3), (px + 3, py + 3)], fill=(220, 70, 70)) + + for px, py in (straight_points[0], straight_points[-1]): + draw.ellipse([(px - 5, py - 5), (px + 5, py + 5)], fill=(40, 110, 220)) + + draw.text( + (margin_left, 20), + "Least Action Trajectory Comparison", + fill=(0, 0, 0), + font=title_font, + ) + draw.text((width // 2 - 45, height - 50), "time t", fill=(0, 0, 0), font=label_font) + draw.text((20, 20), "position x", fill=(0, 0, 0), font=label_font) + + legend_x = width - 290 + legend_y = 28 + draw.line([(legend_x, legend_y + 10), (legend_x + 40, legend_y + 10)], fill=(40, 110, 220), width=4) + draw.text((legend_x + 50, legend_y - 6), "straight path", fill=(0, 0, 0), font=legend_font) + draw.line([(legend_x, legend_y + 40), (legend_x + 40, legend_y + 40)], fill=(220, 70, 70), width=3) + draw.text((legend_x + 50, legend_y + 24), "random path", fill=(0, 0, 0), font=legend_font) + + image.save(image_path, format="JPEG", quality=95) + + +def solve_case(input_path: Path, output_path: Path) -> str: + params = load_input(input_path) + + mass = float(params["m"]) + xa = float(params["xa"]) + ta = float(params["ta"]) + xb = float(params["xb"]) + tb = float(params["tb"]) + n_points = int(params["N"]) + xmax = float(params["xmax"]) + seed = params.get("seed") + + if n_points < 2: + raise ValueError("N must be at least 2.") + if tb <= ta: + raise ValueError("tb must be greater than ta.") + + if seed is not None: + random.seed(seed) + + dt = (tb - ta) / (n_points - 1) + time_points = build_time_points(ta=ta, tb=tb, n_points=n_points) + + straight_path = build_straight_path(xa=xa, xb=xb, n_points=n_points) + random_path = build_random_path(straight_path=straight_path, xmax=xmax) + + random_action = kinetic_action(path=random_path, mass=mass, dt=dt) + straight_action = kinetic_action(path=straight_path, mass=mass, dt=dt) + ratio = random_action / straight_action if straight_action != 0 else float("inf") + image_path = output_path.parent / "trajectory.jpg" + + result = "\n".join( + [ + "Least Action Case 01: Free Particle", + f"input_file: {input_path}", + f"trajectory_image: {image_path}", + f"m: {mass}", + f"xa: {xa}", + f"ta: {ta}", + f"xb: {xb}", + f"tb: {tb}", + f"N: {n_points}", + f"xmax: {xmax}", + f"seed: {seed}", + f"dt: {dt}", + f"random_action: {random_action:.10f}", + f"straight_action: {straight_action:.10f}", + f"random_over_straight: {ratio:.10f}", + ] + ) + + output_path.parent.mkdir(parents=True, exist_ok=True) + draw_trajectory_image( + image_path=image_path, + time_points=time_points, + random_path=random_path, + straight_path=straight_path, + ) + output_path.write_text(result + "\n", encoding="utf-8") + return result + + +def main() -> None: + parser = argparse.ArgumentParser(description="Compute action for a free particle path.") + parser.add_argument("input", type=Path, help="YAML input file path.") + parser.add_argument("output", type=Path, help="Output result file path.") + args = parser.parse_args() + + solve_case(input_path=args.input, output_path=args.output) + + +if __name__ == "__main__": + main()