research(rectification): angle-timing feature library — transits to angles, secondary progressions, solar arcs (BUG-1141 A0)

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_017eEAG8HD3mm8gsKXgk8uU8
This commit is contained in:
Jesse_Chen
2026-10-01 10:55:46 +08:00
co-authored by Claude Opus 5.5
parent 07fcba78f0
commit f507570448
2 changed files with 336 additions and 0 deletions
+253
View File
@@ -0,0 +1,253 @@
"""Angle-based timing features for birth-time rectification research (BUG-1141).
Offline only. For each public v5 case and each candidate birth minute in a
±60-minute window, measure how close three degree-level timing techniques come
to an exact contact at each dated life event:
TR transiting Saturn / Jupiter / mean Rahu (and Mars for day-precision
events) to the natal Ascendant or Midheaven, 0° / 90° / 180°;
SP secondary-progressed (day-for-a-year) Ascendant / Midheaven to natal
Sun..Saturn, 0° / 90° / 180°;
SA solar-arc-directed Ascendant / Midheaven to natal Sun..Saturn,
0° / 90° / 180°.
Aspects between points in one zodiac do not depend on the zodiac, so everything
is tropical. Event precision: day → that date; month → TR takes the closest
contact over the month (daily samples), SP / SA the mid-month point (progressed
angles move ~0.08° a month); year → SP / SA only, closest contact over the year
(quarterly samples). TR skips year events (Saturn moves ~12° a year, far wider
than any orb tested). Progressed angles are exact every 5 candidate minutes and
linearly interpolated in between (error well under 0.05°).
A candidate's feature for one event is the minimum separation (degrees) over
its targets and aspects; the research compares the truth minute's count of
events within an orb against shuffled-date controls.
"""
from __future__ import annotations
import calendar
import hashlib
import json
import sys
from dataclasses import dataclass
from datetime import date, datetime, timedelta
from pathlib import Path
from random import Random
from typing import Any, Iterable, Sequence
import numpy as np
import swisseph as swe
ROOT = Path(__file__).resolve().parents[2]
if str(ROOT) not in sys.path:
sys.path.insert(0, str(ROOT))
HOLDOUT_V5 = ROOT / "references" / "real_case_calibration" / "minute_rectification_holdout_v5.json"
EVENT_CUTOFF = date(2026, 9, 14)
RADIUS_MAX = 60
ASPECTS = (0.0, 90.0, 180.0)
NATAL_TARGETS = (swe.SUN, swe.MOON, swe.MERCURY, swe.VENUS, swe.MARS, swe.JUPITER, swe.SATURN)
TRANSIT_SLOW = (swe.SATURN, swe.JUPITER, swe.MEAN_NODE)
TRANSIT_DAY_EXTRA = (swe.MARS,)
TECHNIQUES = ("TR", "SP", "SA")
ORBS = (1.0, 2.0, 3.0)
TROPICAL_YEAR = 365.24219
FLAGS = swe.FLG_SWIEPH
def load_cases(path: Path = HOLDOUT_V5) -> list[dict[str, Any]]:
return json.loads(path.read_text(encoding="utf-8"))["cases"]
def julian_ut(moment_local: datetime, tz_hours: float) -> float:
utc = moment_local - timedelta(hours=float(tz_hours))
return swe.julday(utc.year, utc.month, utc.day, utc.hour + utc.minute / 60 + utc.second / 3600)
def birth_local(case: dict[str, Any]) -> datetime:
birth = case["birth"]
return datetime.strptime(f"{birth['date']} {str(birth['time'])[:5]}", "%Y-%m-%d %H:%M")
def angle_diff(a: np.ndarray | float, b: np.ndarray | float) -> np.ndarray | float:
"""Smallest absolute difference between longitudes in degrees (0..180)."""
d = np.abs((np.asarray(a) - np.asarray(b) + 180.0) % 360.0 - 180.0)
return d
def aspect_sep(a: np.ndarray, b: np.ndarray) -> np.ndarray:
"""Min over 0/90/180 of |sep - aspect| (degrees)."""
raw = angle_diff(a, b)
return np.minimum.reduce([np.abs(raw - aspect) for aspect in ASPECTS])
def planet_lon(jd_ut: float, body: int) -> float:
return float(swe.calc_ut(jd_ut, body, FLAGS)[0][0])
def angles(jd_ut: float, lat: float, lon: float) -> tuple[float, float]:
_cusps, ascmc = swe.houses(jd_ut, float(lat), float(lon), b"P")
return float(ascmc[0]), float(ascmc[1])
@dataclass(frozen=True)
class Event:
id: str
domain: str
precision: str # day | month | year
start: date
end: date
def parse_event(raw: dict[str, Any]) -> Event | None:
text = str(raw.get("date") or "")
precision = str(raw.get("precision") or "")
try:
if precision == "day" and len(text) >= 10:
day = date.fromisoformat(text[:10]); start = end = day
elif precision == "month" and len(text) >= 7:
y, m = int(text[:4]), int(text[5:7]); start = date(y, m, 1); end = date(y, m, calendar.monthrange(y, m)[1])
elif precision == "year" and len(text) >= 4:
y = int(text[:4]); start = date(y, 1, 1); end = date(y, 12, 31)
else:
return None
except ValueError:
return None
if start > EVENT_CUTOFF:
return None
return Event(str(raw.get("id") or text), str(raw.get("domain") or ""), precision, start, min(end, EVENT_CUTOFF))
def sample_dates(event: Event, technique: str) -> list[date]:
if event.precision == "day":
return [event.start]
if event.precision == "month":
if technique == "TR":
return [event.start + timedelta(days=k) for k in range((event.end - event.start).days + 1)]
# progressed / directed angles move ~0.08° a month: the mid-month point is enough
return [event.start + timedelta(days=14)]
# year: progressed angles move ~1°/yr, quarterly samples bound the error at ~0.125°
if technique == "TR":
return []
return [event.start + timedelta(days=k) for k in (0, 91, 182, 273, (event.end - event.start).days)]
class CaseChart:
"""Natal quantities for every candidate minute of one case (offsets -60..+60)."""
def __init__(self, case: dict[str, Any]):
self.case = case
self.case_id = str(case["case_id"])
birth = case["birth"]
self.lat, self.lon, self.tz = float(birth["latitude"]), float(birth["longitude"]), float(birth["timezone_offset"])
self.local = birth_local(case)
self.offsets = np.arange(-RADIUS_MAX, RADIUS_MAX + 1)
self.jd = np.array([julian_ut(self.local + timedelta(minutes=int(o)), self.tz) for o in self.offsets])
ang = [angles(j, self.lat, self.lon) for j in self.jd]
self.asc = np.array([a for a, _ in ang]); self.mc = np.array([m for _, m in ang])
self.natal = np.array([[planet_lon(j, body) for body in NATAL_TARGETS] for j in self.jd]) # (n, 7)
self.sun = self.natal[:, 0]
def events(self, raw_events: Iterable[dict[str, Any]] | None = None) -> list[Event]:
items = [parse_event(e) for e in (raw_events if raw_events is not None else self.case["events"])]
birth_day = self.local.date()
return [e for e in items if e is not None and e.start > birth_day]
# --- per-technique minimum separation for one event, one value per candidate ----------
def sep_TR(self, event: Event) -> np.ndarray | None:
days = sample_dates(event, "TR")
if not days:
return None
bodies = TRANSIT_SLOW + (TRANSIT_DAY_EXTRA if event.precision == "day" else ())
best = np.full(len(self.offsets), np.inf)
for day in days:
jd = swe.julday(day.year, day.month, day.day, 12.0)
for body in bodies:
lon = planet_lon(jd, body)
points = [lon, (lon + 180.0) % 360.0] if body == swe.MEAN_NODE else [lon]
for p in points:
best = np.minimum(best, aspect_sep(self.asc, p))
best = np.minimum(best, aspect_sep(self.mc, p))
return best
def _age_days(self, day: date) -> float:
noon = swe.julday(day.year, day.month, day.day, 12.0)
return (noon - self.jd[RADIUS_MAX]) / TROPICAL_YEAR
def _progressed_angles(self, age_days: float) -> tuple[np.ndarray, np.ndarray]:
"""Progressed ASC / MC for every candidate: exact every 5 minutes, linear in between (unwrapped)."""
knots = np.arange(0, len(self.offsets), 5)
if knots[-1] != len(self.offsets) - 1:
knots = np.append(knots, len(self.offsets) - 1)
vals = np.array([angles(self.jd[k] + age_days, self.lat, self.lon) for k in knots])
idx = np.arange(len(self.offsets))
asc = np.interp(idx, knots, np.unwrap(np.radians(vals[:, 0]))); mc = np.interp(idx, knots, np.unwrap(np.radians(vals[:, 1])))
return np.degrees(asc) % 360.0, np.degrees(mc) % 360.0
def sep_SP(self, event: Event) -> np.ndarray | None:
best = np.full(len(self.offsets), np.inf)
for day in sample_dates(event, "SP"):
asc, mc = self._progressed_angles(self._age_days(day)) # progressed days == age in years
for col in range(self.natal.shape[1]):
t = self.natal[:, col]
best = np.minimum(best, aspect_sep(asc, t)); best = np.minimum(best, aspect_sep(mc, t))
return best
def sep_SA(self, event: Event) -> np.ndarray | None:
best = np.full(len(self.offsets), np.inf)
for day in sample_dates(event, "SA"):
age_days = self._age_days(day)
progressed_sun = planet_lon(self.jd[RADIUS_MAX] + age_days, swe.SUN) # Sun moves ~0.04°/h: one value per sample
arc = (progressed_sun - self.sun) % 360.0
asc = (self.asc + arc) % 360.0; mc = (self.mc + arc) % 360.0
for col in range(self.natal.shape[1]):
t = self.natal[:, col]
best = np.minimum(best, aspect_sep(asc, t)); best = np.minimum(best, aspect_sep(mc, t))
return best
def separations(self, technique: str, events: Sequence[Event]) -> np.ndarray:
"""(n_events_used, n_candidates) matrix of minimum separations; events the technique skips are dropped."""
fn = {"TR": self.sep_TR, "SP": self.sep_SP, "SA": self.sep_SA}[technique]
rows = [r for r in (fn(e) for e in events) if r is not None]
return np.array(rows) if rows else np.zeros((0, len(self.offsets)))
def hit_counts(seps: np.ndarray, orb: float) -> np.ndarray:
return (seps <= orb).sum(axis=0) if seps.size else np.zeros(seps.shape[1] if seps.ndim == 2 else 0, dtype=int)
def truth_rank(scores: np.ndarray, radius: int) -> float:
"""Fraction of other window minutes that beat the truth minute (ties count half); 0 = truth best, 0.5 = chance."""
window = scores[RADIUS_MAX - radius: RADIUS_MAX + radius + 1]
truth = window[radius]
others = np.delete(window, radius)
return float(((others > truth).sum() + 0.5 * (others == truth).sum()) / len(others))
def shuffled_events(case: dict[str, Any], rng: Random) -> list[dict[str, Any]]:
"""Same events, precision and domain, with dates drawn uniformly over the observed life span."""
birth = birth_local(case).date()
parsed = [parse_event(e) for e in case["events"]]
years = [p.start.year for p in parsed if p is not None]
lo = birth.year + 1
hi = max(max(years) if years else birth.year + 20, lo + 1)
out = []
for raw in case["events"]:
p = str(raw.get("precision") or "")
y = rng.randint(lo, hi)
if p == "day":
m = rng.randint(1, 12); d = rng.randint(1, calendar.monthrange(y, m)[1]); text = f"{y:04d}-{m:02d}-{d:02d}"
elif p == "month":
m = rng.randint(1, 12); text = f"{y:04d}-{m:02d}"
else:
text = f"{y:04d}"
out.append({**raw, "date": text})
return out
def stable_json(payload: Any) -> str:
return json.dumps(payload, ensure_ascii=False, sort_keys=True, indent=2) + "\n"
def sha256_of(path: Path) -> str:
return hashlib.sha256(path.read_bytes()).hexdigest()
+83
View File
@@ -0,0 +1,83 @@
"""Angle-timing research helpers (BUG-1141). Offline; public v5 cases only."""
from __future__ import annotations
from datetime import date
from random import Random
import numpy as np
import pytest
from scripts.research import angle_timing_lib as at
def test_angle_diff_wraps_and_aspect_sep_takes_nearest_aspect():
assert at.angle_diff(359.0, 1.0) == pytest.approx(2.0)
assert at.angle_diff(10.0, 190.0) == pytest.approx(180.0)
seps = at.aspect_sep(np.array([0.0, 91.0, 178.5, 45.0]), np.array([0.0, 0.0, 0.0, 0.0]))
assert list(np.round(seps, 3)) == [0.0, 1.0, 1.5, 45.0]
@pytest.mark.parametrize("raw, precision, start, end", [
({"date": "2001-03-04", "precision": "day"}, "day", date(2001, 3, 4), date(2001, 3, 4)),
({"date": "2001-02", "precision": "month"}, "month", date(2001, 2, 1), date(2001, 2, 28)),
({"date": "2001", "precision": "year"}, "year", date(2001, 1, 1), date(2001, 12, 31)),
])
def test_parse_event_spans(raw, precision, start, end):
event = at.parse_event({**raw, "id": "x", "domain": "career"})
assert event is not None and (event.precision, event.start, event.end) == (precision, start, end)
def test_parse_event_rejects_future_and_malformed():
assert at.parse_event({"date": "2030", "precision": "year"}) is None
assert at.parse_event({"date": "2001-13", "precision": "month"}) is None
assert at.parse_event({"date": "", "precision": "day"}) is None
def test_sample_dates_by_technique_and_precision():
month = at.parse_event({"date": "2001-02", "precision": "month"})
year = at.parse_event({"date": "2001", "precision": "year"})
assert len(at.sample_dates(month, "TR")) == 28
assert len(at.sample_dates(month, "SP")) == 1
assert at.sample_dates(year, "TR") == []
assert len(at.sample_dates(year, "SA")) == 5
def test_truth_rank_semantics():
scores = np.zeros(121)
scores[at.RADIUS_MAX] = 3
assert at.truth_rank(scores, 10) == 0.0
flat = np.ones(121)
assert at.truth_rank(flat, 10) == pytest.approx(0.5)
worst = np.ones(121); worst[at.RADIUS_MAX] = 0
assert at.truth_rank(worst, 10) == 1.0
def test_shuffle_keeps_precision_count_and_is_seeded():
case = at.load_cases()[0]
a = at.shuffled_events(case, Random("seed-1"))
b = at.shuffled_events(case, Random("seed-1"))
assert a == b
assert [e["precision"] for e in a] == [e["precision"] for e in case["events"]]
assert [len(e["date"]) for e in a] == [len(str(e["date"])[:{"day": 10, "month": 7}.get(e["precision"], 4)]) for e in case["events"]]
birth_year = int(case["birth"]["date"][:4])
assert all(int(e["date"][:4]) > birth_year for e in a)
def test_natal_angles_move_about_a_quarter_degree_per_minute():
chart = at.CaseChart(at.load_cases()[0])
step_asc = at.angle_diff(chart.asc[1:], chart.asc[:-1])
step_mc = at.angle_diff(chart.mc[1:], chart.mc[:-1])
assert 0.05 < float(np.median(step_asc)) < 1.0
assert 0.2 < float(np.median(step_mc)) < 0.3
# the truth minute is the window centre
assert chart.offsets[at.RADIUS_MAX] == 0
def test_separation_matrices_have_one_row_per_usable_event():
chart = at.CaseChart(at.load_cases()[0])
events = chart.events()
tr = chart.separations("TR", events)
sp = chart.separations("SP", events)
assert tr.shape[0] == sum(1 for e in events if e.precision != "year")
assert sp.shape == (len(events), len(chart.offsets))
assert np.isfinite(sp).all() and (sp >= 0).all() and (sp <= 45).all()