diff --git a/scripts/research/angle_timing_lib.py b/scripts/research/angle_timing_lib.py new file mode 100644 index 00000000..799f8d33 --- /dev/null +++ b/scripts/research/angle_timing_lib.py @@ -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() diff --git a/tests/test_angle_timing_research.py b/tests/test_angle_timing_research.py new file mode 100644 index 00000000..208552a3 --- /dev/null +++ b/tests/test_angle_timing_research.py @@ -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()