Jyotisha/scripts/research/dasha_shift_per_minute.py

#!/usr/bin/env python3
"""R3a (2026-09-26): how far Vimshottari boundaries move per minute of birth time.

Fictional moments only (random instants, no people). Three measurements:

1. analytic: Moon speed (swisseph, sidereal Raman) x lord years / 13.333 deg;
2. finite difference: first-dasha balance at t and t + 1 minute (the 09-26
   audit method, reproduced here with the same seed);
3. repo helper: `scripts.rectification.event_probes._vim_start_dates` at t and
   t + 20 minutes; all MD/AD starts in [birth+5, birth+80] shift by one
   constant, which is divided by 20.

It then derives what the production 45 / 30 day probe gate means in minutes,
and how close two candidates can be for evidence of a given date precision to
separate them. It corrects the 2026-09-14 statement "1 minute ~ 1.1 days",
which used Vimshottari 1 deg ~ 122 days; the correct mean is ~365 days/deg
(120 years over 9 nakshatras = 120 deg).

Run:
    python3 scripts/research/dasha_shift_per_minute.py
"""

from __future__ import annotations

import argparse
import json
import random
import statistics
import sys
from datetime import date, timedelta
from pathlib import Path
from typing import Any

ROOT = Path(__file__).resolve().parents[2]
for extra in (ROOT, ROOT / "scripts"):
    if str(extra) not in sys.path:
        sys.path.insert(0, str(extra))

import swisseph as swe  # noqa: E402

from scripts.rectification.event_probes import (  # noqa: E402
    MIN_BOUNDARY_DAYS,
    REFRESH_MIN_BOUNDARY_DAYS,
    _vim_start_dates,
)
from scripts.research.offline_research_20260926_lib import (  # noqa: E402
    DAYS_PER_YEAR,
    NAKSHATRA_SPAN,
    VIM_LORDS,
    VIM_YEARS,
    distribution,
    minutes_for_gap,
    nakshatra_lord,
    percentile,
    shift_days_per_minute,
    zip_misaligned,
)

REPORT_JSON = ROOT / "docs" / "research" / "dasha_shift_per_minute_2026_09_26.json"
SAMPLES = 2000
SEED = 7  # same seed as the 09-26 audit sketch, so the 3.8-day median reproduces
REPO_SAMPLES = 300
WINDOWS = (2, 4, 6, 10, 20, 30, 40, 60)
MONTH_DAYS = 30.44
OLD_CLAIM_DAYS_PER_MINUTE = 1.1
OLD_CLAIM_DEG_DAYS = 122.0


def moon(jd: float) -> tuple[float, float]:
    values, _flags = swe.calc_ut(jd, swe.MOON, swe.FLG_SIDEREAL | swe.FLG_SPEED)
    return float(values[0]), float(values[3])


def balance_days(longitude: float) -> float:
    progress = (longitude % NAKSHATRA_SPAN) / NAKSHATRA_SPAN
    return (1.0 - progress) * VIM_YEARS[nakshatra_lord(longitude)] * DAYS_PER_YEAR


def sample_instants(n: int, seed: int) -> list[float]:
    rng = random.Random(seed)
    rows = []
    for _ in range(n):
        rows.append(swe.julday(1950 + rng.randrange(55), 1, 1, 0.0) + rng.random() * 365.0)
    return rows


def jd_to_date(jd: float) -> date:
    year, month, day, _hour = swe.revjul(jd)
    return date(int(year), int(month), int(day))


def repo_shift(jd: float, minutes: int) -> dict[str, Any] | None:
    birth = jd_to_date(jd)
    birth_date = birth.isoformat()
    lo, hi = birth.year + 5, birth.year + 80
    lon_a, _ = moon(jd)
    lon_b, _ = moon(jd + minutes / 1440.0)
    if int(lon_a // NAKSHATRA_SPAN) != int(lon_b // NAKSHATRA_SPAN):
        return {"crossed_nakshatra": True}
    left = _vim_start_dates(birth_date, lon_a, lo, hi)
    right = _vim_start_dates(birth_date, lon_b, lo, hi)
    misaligned = zip_misaligned(left, right)
    # Re-align on corresponding boundaries: when a boundary crosses the lower
    # or upper year edge (MD and first-AD starts coincide, so up to two
    # entries can drop), try offsets -3..+3 and keep the one whose
    # signed differences are most constant.
    best: list[int] | None = None
    for offset in (0, 1, -1, 2, -2, 3, -3):
        a = left[offset:] if offset > 0 else left
        b = right[-offset:] if offset < 0 else right
        diffs = [(one - two).days for one, two in zip(a, b)]
        if diffs and (best is None or max(diffs) - min(diffs) < max(best) - min(best)):
            best = diffs
    diffs = best or []
    if not diffs:
        return None
    return {
        "crossed_nakshatra": False,
        "per_minute": abs(statistics.median(diffs)) / minutes,
        "spread_days": max(diffs) - min(diffs),
        "raw_misaligned": misaligned,
        "boundaries": len(diffs),
    }


def main() -> int:
    parser = argparse.ArgumentParser()
    parser.add_argument("--samples", type=int, default=SAMPLES)
    parser.add_argument("--repo-samples", type=int, default=REPO_SAMPLES)
    parser.add_argument("--no-write", action="store_true")
    args = parser.parse_args()
    swe.set_sid_mode(swe.SIDM_RAMAN)
    instants = sample_instants(args.samples, SEED)
    analytic: list[float] = []
    finite: list[float] = []
    by_lord: dict[str, list[float]] = {lord: [] for lord in VIM_LORDS}
    speeds: list[float] = []
    for jd in instants:
        lon, speed = moon(jd)
        speeds.append(speed)
        value = shift_days_per_minute(lon, speed)
        analytic.append(value)
        by_lord[nakshatra_lord(lon)].append(value)
        lon_next, _ = moon(jd + 1 / 1440.0)
        if int(lon // NAKSHATRA_SPAN) == int(lon_next // NAKSHATRA_SPAN):
            # boundary = birth + balance; birth itself moves +1 minute.
            finite.append(abs(balance_days(lon_next) + 1 / 1440.0 - balance_days(lon)))
    repo_rows = []
    crossed = 0
    for jd in instants[: args.repo_samples]:
        row = repo_shift(jd, 20)
        if row is None:
            continue
        if row.get("crossed_nakshatra"):
            crossed += 1
            continue
        repo_rows.append(row)
    repo_values = [row["per_minute"] for row in repo_rows]
    p10, p50, p90 = (percentile(analytic, q) for q in (10, 50, 90))
    lo_bound, hi_bound = min(analytic), max(analytic)

    gate_minutes = {}
    for label, days in (("gate_45", MIN_BOUNDARY_DAYS), ("refresh_30", REFRESH_MIN_BOUNDARY_DAYS)):
        gate_minutes[label] = distribution([minutes_for_gap(days, s) for s in analytic], 1)

    window_rows = []
    for minutes in WINDOWS:
        gaps = [minutes * s for s in analytic]
        same_year_pass = sum(1 for g in gaps if g >= MIN_BOUNDARY_DAYS) / len(gaps)
        window_rows.append({
            "candidate_gap_minutes": minutes,
            "boundary_gap_days": distribution(gaps, 1),
            "share_ge_45_days": round(same_year_pass, 3),
            "share_ge_30_days": round(sum(1 for g in gaps if g >= REFRESH_MIN_BOUNDARY_DAYS) / len(gaps), 3),
            "old_claim_days": round(minutes * OLD_CLAIM_DAYS_PER_MINUTE, 1),
        })

    # How close can two candidates be for evidence of a given precision to
    # separate them? The event date must fall between the two boundary dates.
    # Guaranteed separation needs a boundary gap of about the precision span;
    # below that, separation happens only when the gap straddles a calendar
    # edge, with probability ~ gap / span (uniform boundary position).
    precision_rows = []
    for label, span_days in (
        ("day_exact", 1.0),
        ("day_pm3", 7.0),
        ("day_pm7", 15.0),
        ("month", MONTH_DAYS),
        ("gate_45_same_year", float(MIN_BOUNDARY_DAYS)),
        ("year", DAYS_PER_YEAR),
    ):
        precision_rows.append({
            "evidence": label,
            "span_days": round(span_days, 2),
            "min_separable_minutes_p90_shift": round(span_days / p90, 1),
            "min_separable_minutes_median_shift": round(span_days / p50, 1),
            "min_separable_minutes_p10_shift": round(span_days / p10, 1),
            "min_separable_minutes_slowest": round(span_days / lo_bound, 1),
        })
    straddle = []
    for minutes in (2, 6, 10, 20):
        gap = minutes * p50
        straddle.append({
            "candidate_gap_minutes": minutes,
            "boundary_gap_days_median": round(gap, 1),
            "p_different_month": round(min(1.0, gap / MONTH_DAYS), 3),
            "p_different_year": round(min(1.0, gap / DAYS_PER_YEAR), 3),
        })
    payload = {
        "generated_at": "2026-09-26",
        "nature": "fictional instants only; pure astronomy + repo Vimshottari helper",
        "ayanamsa": "raman (swisseph SIDM_RAMAN)",
        "dasha_year_days": DAYS_PER_YEAR,
        "samples": len(instants),
        "seed": SEED,
        "old_claim": {
            "days_per_minute": OLD_CLAIM_DAYS_PER_MINUTE,
            "vimshottari_days_per_degree": OLD_CLAIM_DEG_DAYS,
            "correct_mean_days_per_degree": round(120 * DAYS_PER_YEAR / 120, 2),
            "moon_deg_per_minute_mean": round(statistics.fmean(speeds) / 1440.0, 5),
        },
        "analytic_days_per_minute": distribution(analytic),
        "finite_difference_days_per_minute": distribution(finite),
        "repo_vim_start_dates_days_per_minute": distribution(repo_values),
        "repo_samples": len(repo_rows),
        "repo_crossed_nakshatra_in_20_min": crossed,
        "repo_max_spread_days": max((row["spread_days"] for row in repo_rows), default=None),
        "repo_raw_zip_misaligned_share": round(
            sum(1 for row in repo_rows if row["raw_misaligned"]) / max(len(repo_rows), 1), 3),
        "by_lord": {
            lord: {**distribution(values), "years": VIM_YEARS[lord]}
            for lord, values in by_lord.items()
        },
        "moon_speed_deg_per_day": distribution(speeds, 3),
        "theoretical_range_days_per_minute": {
            "min": round(lo_bound, 2), "max": round(hi_bound, 2),
        },
        "gate_minutes": gate_minutes,
        "window_table": window_rows,
        "precision_table": precision_rows,
        "straddle_table": straddle,
        "command": "python3 scripts/research/dasha_shift_per_minute.py",
    }
    if not args.no_write:
        REPORT_JSON.write_text(json.dumps(payload, ensure_ascii=False, indent=2) + "\n", encoding="utf-8")
    print(json.dumps({key: payload[key] for key in (
        "analytic_days_per_minute", "finite_difference_days_per_minute",
        "repo_vim_start_dates_days_per_minute", "repo_samples", "repo_crossed_nakshatra_in_20_min",
        "repo_max_spread_days", "repo_raw_zip_misaligned_share", "gate_minutes",
        "theoretical_range_days_per_minute", "old_claim")}, ensure_ascii=False, indent=1))
    for row in window_rows:
        print(row)
    for row in precision_rows:
        print(row)
    for row in straddle:
        print(row)
    for lord, row in payload["by_lord"].items():
        print(lord, row)
    return 0


if __name__ == "__main__":
    sys.exit(main())