From f5c5e350c9178abffd1a1cc12bae18e6b7ebfab8 Mon Sep 17 00:00:00 2001 From: Ying-Li Niu <64801511@qq.com> Date: Wed, 3 Jun 2026 16:36:23 +0800 Subject: [PATCH] modified: examples/case01/input/input.txt modified: least_action.py --- examples/case01/input/input.txt | 2 +- least_action.py | 657 +++++++++++++++++++++++++++----- 2 files changed, 573 insertions(+), 86 deletions(-) diff --git a/examples/case01/input/input.txt b/examples/case01/input/input.txt index a17557a..fb67813 100644 --- a/examples/case01/input/input.txt +++ b/examples/case01/input/input.txt @@ -4,5 +4,5 @@ ta: 0 xb: 20 tb: 10 N: 100 -xmax: 1.0 +xmax: 1.0 0.5 0.1 0.01 0.001 seed: 42 diff --git a/least_action.py b/least_action.py index 30a1d2d..6eb1398 100644 --- a/least_action.py +++ b/least_action.py @@ -1,15 +1,16 @@ from __future__ import annotations import argparse +import html import random from pathlib import Path -from typing import Dict, List +from typing import Dict, List, Sequence from PIL import Image, ImageDraw, ImageFont import yaml -def load_input(input_path: Path) -> Dict[str, float]: +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 {} @@ -27,28 +28,53 @@ def load_input(input_path: Path) -> Dict[str, float]: 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) - ] + 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: +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 in straight_path[1:-1]: - path.append(base_x + random.uniform(-1.0, 1.0) * xmax) + 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: List[float], mass: float, dt: float) -> float: +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 @@ -65,38 +91,98 @@ def build_time_points(ta: float, tb: float, n_points: int) -> List[float]: 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: List[float], - random_path: List[float], - straight_path: List[float], + time_points: Sequence[float], + straight_path: Sequence[float], + random_cases: Sequence[Dict[str, object]], ) -> None: - width = 1200 - height = 800 - margin_left = 110 - margin_right = 60 - margin_top = 80 - margin_bottom = 110 + 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", 34) - label_font = ImageFont.truetype("DejaVuSans.ttf", 24) - tick_font = ImageFont.truetype("DejaVuSans.ttf", 20) - legend_font = ImageFont.truetype("DejaVuSans.ttf", 24) + 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_x = random_path + straight_path + 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_x) - x_max = max(all_x) + x_min = min(all_positions) + x_max = max(all_positions) if t_min == t_max: t_max = t_min + 1.0 @@ -125,50 +211,421 @@ def draw_trajectory_image( 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.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 - 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) - 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=(35, 95, 210), width=6) + for point in (straight_points[0], straight_points[-1]): + draw_circle(draw, point, 8, (35, 95, 210)) - draw.line(straight_points, fill=(40, 110, 220), width=4) - draw.line(random_points, fill=(220, 70, 70), width=3) + 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"]) - for px, py in random_points: - draw.ellipse([(px - 3, py - 3), (px + 3, py + 3)], fill=(220, 70, 70)) + 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) - 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, + 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((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) + draw.text((legend_left, legend_top + 22), "Legend", fill=(0, 0, 0), font=legend_title_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) + 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( + [ + '
', + '
Straight Action' + f'{straight_action:.4f}
', + '
Random Cases' + f'{case_count}
', + '
Smallest xmax' + f'{min_xmax:.6g}
', + '
Best Ratio' + f'{best_ratio:.4f}
', + '
', + ] + ) + + +def build_summary_rows(straight_action: float, random_cases: Sequence[Dict[str, object]]) -> str: + rows: List[str] = [] + for case in random_cases: + rows.append( + "" + f"{case['xmax']:.6g}" + f"{case['random_action']:.6f}" + f"{straight_action:.6f}" + f"{case['ratio']:.6f}" + f"{case['difference']:.6f}" + "" + ) + 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""" + + + + + Least Action Case 01 + + + + + +
+
+

Least Action / Case 01

+

自由粒子最小作用量收敛示例

+

+ 页面使用同一组基础扰动，分别按不同的 xmax 缩放随机路径，并比较它们的离散作用量如何逐步收敛到匀速直线轨迹的作用量。 +

+ {stats_cards} +

+ 参数：m={mass:.6g}， + N={n_points}， + seed={seed}， + xmax={", ".join(f"{case['xmax']:.6g}" for case in random_cases)} +

+
+ +
+

轨迹对比图

+ 不同 xmax 下的随机轨迹与匀速直线轨迹对比图 +
+ +
+

收敛汇总

+

{observation}

+ + + + + + + + + + + + {summary_rows} + +
xmaxrandom actionstraight actionrandom / straightdifference
+
+ +
+

作用量公式

+
+
+

匀速直线轨迹

+

先按起点和终点做等时间间隔的线性插值，再基于离散速度计算作用量。

+
+ \\[ + x_i = x_a + \\frac{{i}}{{N-1}}(x_b - x_a), \\qquad i = 0,1,\\dots,N-1 + \\] +
+
+ \\[ + \\Delta t = \\frac{{t_b - t_a}}{{N-1}}, \\qquad + v_i = \\frac{{x_{{i+1}} - x_i}}{{\\Delta t}} + \\] +
+
+ \\[ + S_{{\\mathrm{{straight}}}} = \\sum_{{i=0}}^{{N-2}}\\frac{{1}}{{2}}m v_i^2 \\Delta t + \\] +
+
+ +
+

随机扰动轨迹

+

以直线路径为基线，在每个中间采样点加入均匀分布扰动，并保持起点终点不变。

+
+ \\[ + x_i^{{(\\mathrm{{rand}})}} = x_i + \\epsilon_i, \\qquad + \\epsilon_i \\sim U(-x_{{\\max}}, x_{{\\max}}) + \\] +
+
+ \\[ + x_0^{{(\\mathrm{{rand}})}} = x_a, \\qquad + x_{{N-1}}^{{(\\mathrm{{rand}})}} = x_b + \\] +
+
+ \\[ + v_i^{{(\\mathrm{{rand}})}} = + \\frac{{x_{{i+1}}^{{(\\mathrm{{rand}})}} - x_i^{{(\\mathrm{{rand}})}}}}{{\\Delta t}} + \\] +
+
+ \\[ + S_{{\\mathrm{{rand}}}} = + \\sum_{{i=0}}^{{N-2}}\\frac{{1}}{{2}}m\\left(v_i^{{(\\mathrm{{rand}})}}\\right)^2 \\Delta t + \\] +
+
+
+

+ 页面中的公式由 MathJax 在浏览器中渲染。若本地打开后暂时没有看到公式样式，请确认浏览器能够访问 + cdn.jsdelivr.net。 +

+
+
+ + +""" + + (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) @@ -178,7 +635,7 @@ def solve_case(input_path: Path, output_path: Path) -> str: xb = float(params["xb"]) tb = float(params["tb"]) n_points = int(params["N"]) - xmax = float(params["xmax"]) + xmax_values = parse_xmax_values(params["xmax"]) seed = params.get("seed") if n_points < 2: @@ -186,49 +643,79 @@ def solve_case(input_path: Path, output_path: Path) -> str: 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") + 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 = "\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}", - ] - ) + 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, - random_path=random_path, straight_path=straight_path, + random_cases=random_cases, ) - output_path.write_text(result + "\n", encoding="utf-8") - return result + 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: