Add least action example and output ignore rules
This commit is contained in:
+244
@@ -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()
|
||||
Reference in New Issue
Block a user