Files
least_action/least_action.py
T
2026-06-03 16:36:23 +08:00

732 lines
23 KiB
Python
Raw Blame History

This file contains ambiguous Unicode characters
This file contains Unicode characters that might be confused with other characters. If you think that this is intentional, you can safely ignore this warning. Use the Escape button to reveal them.
from __future__ import annotations
import argparse
import html
import random
from pathlib import Path
from typing import Dict, List, Sequence
from PIL import Image, ImageDraw, ImageFont
import yaml
def load_input(input_path: Path) -> Dict[str, object]:
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 parse_xmax_values(raw_value: object) -> List[float]:
if isinstance(raw_value, (int, float)):
return [float(raw_value)]
if isinstance(raw_value, str):
normalized = raw_value.replace(",", " ")
values = [float(part) for part in normalized.split() if part]
if values:
return values
if isinstance(raw_value, list):
values = [float(item) for item in raw_value]
if values:
return values
raise ValueError("xmax must be a number, a YAML list, or a space-separated string of numbers.")
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_noise_profile(n_points: int, seed: object) -> List[float]:
if n_points < 2:
raise ValueError("N must be at least 2.")
rng = random.Random(seed)
return [rng.uniform(-1.0, 1.0) for _ in range(n_points - 2)]
def build_random_path(straight_path: Sequence[float], xmax: float, noise_profile: Sequence[float]) -> List[float]:
if len(straight_path) < 2:
raise ValueError("N must be at least 2.")
if len(noise_profile) != len(straight_path) - 2:
raise ValueError("Noise profile length must match the number of interior points.")
path = [straight_path[0]]
for base_x, noise in zip(straight_path[1:-1], noise_profile):
path.append(base_x + noise * xmax)
path.append(straight_path[-1])
return path
def kinetic_action(path: Sequence[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_circle(draw: ImageDraw.ImageDraw, center: tuple[float, float], size: int, color: tuple[int, int, int]) -> None:
px, py = center
draw.ellipse([(px - size, py - size), (px + size, py + size)], fill=color, outline=color)
def draw_square(draw: ImageDraw.ImageDraw, center: tuple[float, float], size: int, color: tuple[int, int, int]) -> None:
px, py = center
draw.rectangle([(px - size, py - size), (px + size, py + size)], fill=color, outline=color)
def draw_triangle(draw: ImageDraw.ImageDraw, center: tuple[float, float], size: int, color: tuple[int, int, int]) -> None:
px, py = center
points = [(px, py - size), (px - size, py + size), (px + size, py + size)]
draw.polygon(points, fill=color, outline=color)
def draw_diamond(draw: ImageDraw.ImageDraw, center: tuple[float, float], size: int, color: tuple[int, int, int]) -> None:
px, py = center
points = [(px, py - size), (px - size, py), (px, py + size), (px + size, py)]
draw.polygon(points, fill=color, outline=color)
def draw_cross(draw: ImageDraw.ImageDraw, center: tuple[float, float], size: int, color: tuple[int, int, int]) -> None:
px, py = center
draw.line([(px - size, py - size), (px + size, py + size)], fill=color, width=2)
draw.line([(px - size, py + size), (px + size, py - size)], fill=color, width=2)
def draw_plus(draw: ImageDraw.ImageDraw, center: tuple[float, float], size: int, color: tuple[int, int, int]) -> None:
px, py = center
draw.line([(px - size, py), (px + size, py)], fill=color, width=2)
draw.line([(px, py - size), (px, py + size)], fill=color, width=2)
def draw_marker(
draw: ImageDraw.ImageDraw,
center: tuple[float, float],
marker_name: str,
size: int,
color: tuple[int, int, int],
) -> None:
if marker_name == "circle":
draw_circle(draw, center, size, color)
elif marker_name == "square":
draw_square(draw, center, size, color)
elif marker_name == "triangle":
draw_triangle(draw, center, size, color)
elif marker_name == "diamond":
draw_diamond(draw, center, size, color)
elif marker_name == "cross":
draw_cross(draw, center, size, color)
else:
draw_plus(draw, center, size, color)
def draw_trajectory_image(
image_path: Path,
time_points: Sequence[float],
straight_path: Sequence[float],
random_cases: Sequence[Dict[str, object]],
) -> None:
width = 1800
height = 1200
margin_left = 150
margin_right = 460
margin_top = 110
margin_bottom = 150
image = Image.new("RGB", (width, height), "white")
draw = ImageDraw.Draw(image)
try:
title_font = ImageFont.truetype("DejaVuSans.ttf", 50)
label_font = ImageFont.truetype("DejaVuSans.ttf", 34)
tick_font = ImageFont.truetype("DejaVuSans.ttf", 28)
legend_title_font = ImageFont.truetype("DejaVuSans.ttf", 34)
legend_font = ImageFont.truetype("DejaVuSans.ttf", 28)
except OSError:
title_font = ImageFont.load_default()
label_font = ImageFont.load_default()
tick_font = ImageFont.load_default()
legend_title_font = ImageFont.load_default()
legend_font = ImageFont.load_default()
all_positions = list(straight_path)
for case in random_cases:
all_positions.extend(case["path"])
t_min = min(time_points)
t_max = max(time_points)
x_min = min(all_positions)
x_max = max(all_positions)
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 + 12)], fill=(0, 0, 0), width=3)
draw.line([(plot_left - 12, tick_y), (plot_left, tick_y)], fill=(0, 0, 0), width=3)
draw.text((tick_x - 28, plot_bottom + 24), f"{x_value:.1f}", fill=(0, 0, 0), font=tick_font)
draw.text((28, tick_y - 14), f"{y_value:.1f}", fill=(0, 0, 0), font=tick_font)
straight_points = [map_point(t, x) for t, x in zip(time_points, straight_path)]
draw.line(straight_points, fill=(35, 95, 210), width=6)
for point in (straight_points[0], straight_points[-1]):
draw_circle(draw, point, 8, (35, 95, 210))
marker_step = max(1, len(time_points) // 12)
for case in random_cases:
points = [map_point(t, x) for t, x in zip(time_points, case["path"])]
draw.line(points, fill=case["color"], width=4)
for point in points[::marker_step]:
draw_marker(draw, point, case["marker"], 6, case["color"])
draw.text((margin_left, 28), "Least Action Trajectory Comparison", fill=(0, 0, 0), font=title_font)
draw.text((width // 2 - 70, height - 72), "time t", fill=(0, 0, 0), font=label_font)
draw.text((24, 28), "position x", fill=(0, 0, 0), font=label_font)
legend_left = plot_right + 24
legend_top = margin_top
legend_bottom = plot_bottom
draw.rectangle(
[legend_left - 10, legend_top, width - 30, legend_bottom],
fill=(250, 248, 242),
outline=(210, 205, 195),
width=2,
)
draw.text((legend_left, legend_top + 22), "Legend", fill=(0, 0, 0), font=legend_title_font)
item_y = legend_top + 102
draw.line([(legend_left, item_y + 14), (legend_left + 66, item_y + 14)], fill=(35, 95, 210), width=6)
draw.text((legend_left + 86, item_y - 4), "straight path", fill=(0, 0, 0), font=legend_font)
item_y += 64
for case in random_cases:
draw.line([(legend_left, item_y + 14), (legend_left + 66, item_y + 14)], fill=case["color"], width=4)
draw_marker(draw, (legend_left + 33, item_y + 14), case["marker"], 6, case["color"])
label = f"xmax = {case['xmax']:.6g}"
draw.text((legend_left + 86, item_y - 4), label, fill=(0, 0, 0), font=legend_font)
item_y += 60
image.save(image_path, format="JPEG", quality=95)
def render_stats_cards(case_count: int, straight_action: float, min_xmax: float, best_ratio: float) -> str:
return "\n".join(
[
'<div class="stats">',
' <div class="stat"><span class="stat-label">Straight Action</span>'
f'<span class="stat-value">{straight_action:.4f}</span></div>',
' <div class="stat"><span class="stat-label">Random Cases</span>'
f'<span class="stat-value">{case_count}</span></div>',
' <div class="stat"><span class="stat-label">Smallest xmax</span>'
f'<span class="stat-value">{min_xmax:.6g}</span></div>',
' <div class="stat"><span class="stat-label">Best Ratio</span>'
f'<span class="stat-value">{best_ratio:.4f}</span></div>',
'</div>',
]
)
def build_summary_rows(straight_action: float, random_cases: Sequence[Dict[str, object]]) -> str:
rows: List[str] = []
for case in random_cases:
rows.append(
"<tr>"
f"<td>{case['xmax']:.6g}</td>"
f"<td>{case['random_action']:.6f}</td>"
f"<td>{straight_action:.6f}</td>"
f"<td>{case['ratio']:.6f}</td>"
f"<td>{case['difference']:.6f}</td>"
"</tr>"
)
return "\n".join(rows)
def build_observation(random_cases: Sequence[Dict[str, object]]) -> str:
first_case = random_cases[0]
last_case = random_cases[-1]
return (
f"在这组同源扰动下,当 xmax 从 {first_case['xmax']:.6g} 逐步减小到 {last_case['xmax']:.6g} 时,"
f"随机路径作用量从 {first_case['random_action']:.4f} 收敛到 {last_case['random_action']:.4f}"
f"并逐渐靠近匀速直线轨迹的作用量。"
)
def write_index_html(
output_dir: Path,
mass: float,
n_points: int,
seed: object,
straight_action: float,
random_cases: Sequence[Dict[str, object]],
) -> None:
min_xmax = min(case["xmax"] for case in random_cases)
best_ratio = min(case["ratio"] for case in random_cases)
stats_cards = render_stats_cards(len(random_cases), straight_action, min_xmax, best_ratio)
summary_rows = build_summary_rows(straight_action, random_cases)
observation = html.escape(build_observation(random_cases))
html_text = f"""<!DOCTYPE html>
<html lang="zh-CN">
<head>
<meta charset="utf-8" />
<meta name="viewport" content="width=device-width, initial-scale=1" />
<title>Least Action Case 01</title>
<style>
:root {{
--bg: #f4f1ea;
--panel: #fffdf8;
--ink: #1c1b1a;
--muted: #5e5a54;
--line: #d8d1c4;
--accent: #245eb8;
--accent-2: #c44747;
}}
* {{
box-sizing: border-box;
}}
body {{
margin: 0;
font-family: "Georgia", "Noto Serif SC", serif;
background:
radial-gradient(circle at top right, rgba(36, 94, 184, 0.10), transparent 24rem),
radial-gradient(circle at left 20%, rgba(196, 71, 71, 0.08), transparent 20rem),
var(--bg);
color: var(--ink);
}}
main {{
width: min(1180px, calc(100vw - 32px));
margin: 32px auto 56px;
}}
.hero,
.panel {{
background: var(--panel);
border: 1px solid var(--line);
border-radius: 24px;
box-shadow: 0 16px 48px rgba(45, 39, 29, 0.08);
}}
.hero {{
padding: 28px 30px;
margin-bottom: 24px;
}}
.panel {{
padding: 24px;
margin-bottom: 24px;
}}
.eyebrow {{
margin: 0 0 8px;
color: var(--muted);
font-size: 14px;
letter-spacing: 0.08em;
text-transform: uppercase;
}}
h1, h2 {{
margin: 0;
font-weight: 700;
}}
h1 {{
font-size: clamp(30px, 4vw, 44px);
line-height: 1.1;
}}
h2 {{
font-size: 28px;
margin-bottom: 16px;
}}
.subtitle,
.note,
.observation {{
color: var(--muted);
line-height: 1.7;
font-size: 17px;
}}
.stats {{
display: grid;
grid-template-columns: repeat(auto-fit, minmax(180px, 1fr));
gap: 14px;
margin-top: 26px;
}}
.stat {{
padding: 16px 18px;
background: #faf7ef;
border: 1px solid var(--line);
border-radius: 18px;
}}
.stat-label {{
display: block;
color: var(--muted);
font-size: 13px;
margin-bottom: 8px;
text-transform: uppercase;
letter-spacing: 0.05em;
}}
.stat-value {{
font-size: 24px;
font-weight: 700;
}}
.trajectory-image {{
width: 100%;
display: block;
border-radius: 18px;
border: 1px solid var(--line);
background: white;
}}
table {{
width: 100%;
border-collapse: collapse;
overflow: hidden;
border-radius: 18px;
border: 1px solid var(--line);
background: white;
}}
th, td {{
padding: 14px 16px;
border-bottom: 1px solid var(--line);
text-align: left;
}}
th {{
background: #f7f1e7;
}}
.formula-grid {{
display: grid;
grid-template-columns: repeat(auto-fit, minmax(320px, 1fr));
gap: 18px;
}}
.formula-card {{
padding: 18px 20px;
border: 1px solid var(--line);
border-radius: 18px;
background: #fffdfa;
}}
.formula-card h3 {{
margin: 0 0 12px;
font-size: 22px;
}}
.formula-card p {{
margin: 0 0 12px;
color: var(--muted);
line-height: 1.7;
}}
.formula-block {{
margin: 14px 0;
padding: 14px 16px;
border-radius: 14px;
background: #faf7ef;
overflow-x: auto;
}}
.formula-block.blue {{
border-left: 4px solid var(--accent);
}}
.formula-block.red {{
border-left: 4px solid var(--accent-2);
}}
code {{
font-family: "Consolas", "SFMono-Regular", monospace;
font-size: 0.95em;
}}
</style>
<script>
window.MathJax = {{
tex: {{
inlineMath: [['\\\\(', '\\\\)']],
displayMath: [['\\\\[', '\\\\]']],
packages: {{'[+]': ['ams', 'textmacros']}}
}},
options: {{
skipHtmlTags: ['script', 'noscript', 'style', 'textarea', 'pre', 'code']
}}
}};
</script>
<script defer src="https://cdn.jsdelivr.net/npm/mathjax@3/es5/tex-mml-chtml.js"></script>
</head>
<body>
<main>
<section class="hero">
<p class="eyebrow">Least Action / Case 01</p>
<h1>自由粒子最小作用量收敛示例</h1>
<p class="subtitle">
页面使用同一组基础扰动,分别按不同的 <code>xmax</code> 缩放随机路径,并比较它们的离散作用量如何逐步收敛到匀速直线轨迹的作用量。
</p>
{stats_cards}
<p class="note">
参数:<code>m={mass:.6g}</code>
<code>N={n_points}</code>
<code>seed={seed}</code>
<code>xmax={", ".join(f"{case['xmax']:.6g}" for case in random_cases)}</code>
</p>
</section>
<section class="panel">
<h2>轨迹对比图</h2>
<img class="trajectory-image" src="trajectory.jpg" alt="不同 xmax 下的随机轨迹与匀速直线轨迹对比图" />
</section>
<section class="panel">
<h2>收敛汇总</h2>
<p class="observation">{observation}</p>
<table>
<thead>
<tr>
<th>xmax</th>
<th>random action</th>
<th>straight action</th>
<th>random / straight</th>
<th>difference</th>
</tr>
</thead>
<tbody>
{summary_rows}
</tbody>
</table>
</section>
<section class="panel">
<h2>作用量公式</h2>
<div class="formula-grid">
<article class="formula-card">
<h3>匀速直线轨迹</h3>
<p>先按起点和终点做等时间间隔的线性插值,再基于离散速度计算作用量。</p>
<div class="formula-block blue">
\\[
x_i = x_a + \\frac{{i}}{{N-1}}(x_b - x_a), \\qquad i = 0,1,\\dots,N-1
\\]
</div>
<div class="formula-block blue">
\\[
\\Delta t = \\frac{{t_b - t_a}}{{N-1}}, \\qquad
v_i = \\frac{{x_{{i+1}} - x_i}}{{\\Delta t}}
\\]
</div>
<div class="formula-block blue">
\\[
S_{{\\mathrm{{straight}}}} = \\sum_{{i=0}}^{{N-2}}\\frac{{1}}{{2}}m v_i^2 \\Delta t
\\]
</div>
</article>
<article class="formula-card">
<h3>随机扰动轨迹</h3>
<p>以直线路径为基线,在每个中间采样点加入均匀分布扰动,并保持起点终点不变。</p>
<div class="formula-block red">
\\[
x_i^{{(\\mathrm{{rand}})}} = x_i + \\epsilon_i, \\qquad
\\epsilon_i \\sim U(-x_{{\\max}}, x_{{\\max}})
\\]
</div>
<div class="formula-block red">
\\[
x_0^{{(\\mathrm{{rand}})}} = x_a, \\qquad
x_{{N-1}}^{{(\\mathrm{{rand}})}} = x_b
\\]
</div>
<div class="formula-block red">
\\[
v_i^{{(\\mathrm{{rand}})}} =
\\frac{{x_{{i+1}}^{{(\\mathrm{{rand}})}} - x_i^{{(\\mathrm{{rand}})}}}}{{\\Delta t}}
\\]
</div>
<div class="formula-block red">
\\[
S_{{\\mathrm{{rand}}}} =
\\sum_{{i=0}}^{{N-2}}\\frac{{1}}{{2}}m\\left(v_i^{{(\\mathrm{{rand}})}}\\right)^2 \\Delta t
\\]
</div>
</article>
</div>
<p class="note">
页面中的公式由 MathJax 在浏览器中渲染。若本地打开后暂时没有看到公式样式,请确认浏览器能够访问
<code>cdn.jsdelivr.net</code>。
</p>
</section>
</main>
</body>
</html>
"""
(output_dir / "index.html").write_text(html_text, encoding="utf-8")
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_values = parse_xmax_values(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.")
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)
straight_action = kinetic_action(path=straight_path, mass=mass, dt=dt)
noise_profile = build_noise_profile(n_points=n_points, seed=seed)
palette = [
(220, 70, 70),
(28, 150, 110),
(220, 145, 40),
(125, 85, 220),
(40, 135, 205),
(180, 95, 145),
]
markers = ["circle", "square", "triangle", "diamond", "cross", "plus"]
random_cases: List[Dict[str, object]] = []
for index, xmax in enumerate(xmax_values):
random_path = build_random_path(straight_path=straight_path, xmax=xmax, noise_profile=noise_profile)
random_action = kinetic_action(path=random_path, mass=mass, dt=dt)
ratio = random_action / straight_action if straight_action != 0 else float("inf")
random_cases.append(
{
"xmax": xmax,
"path": random_path,
"random_action": random_action,
"ratio": ratio,
"difference": random_action - straight_action,
"color": palette[index % len(palette)],
"marker": markers[index % len(markers)],
}
)
image_path = output_path.parent / "trajectory.jpg"
result_lines = [
"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_values: {', '.join(f'{value:.6g}' for value in xmax_values)}",
f"seed: {seed}",
f"dt: {dt}",
f"straight_action: {straight_action:.10f}",
]
for case in random_cases:
result_lines.append(f"xmax={case['xmax']:.6g} random_action={case['random_action']:.10f}")
result_lines.append(f"xmax={case['xmax']:.6g} random_over_straight={case['ratio']:.10f}")
output_path.parent.mkdir(parents=True, exist_ok=True)
draw_trajectory_image(
image_path=image_path,
time_points=time_points,
straight_path=straight_path,
random_cases=random_cases,
)
write_index_html(
output_dir=output_path.parent,
mass=mass,
n_points=n_points,
seed=seed,
straight_action=straight_action,
random_cases=random_cases,
)
result_text = "\n".join(result_lines)
output_path.write_text(result_text + "\n", encoding="utf-8")
return result_text
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()