#!/usr/bin/env python3
"""Reproduce Earth-likeness rankings from an unmodified NASA TAP CSV."""
from pathlib import Path
import json
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt

ROOT = Path(__file__).resolve().parents[1]
RAW = ROOT / "data" / "raw" / "nasa_exoplanet_archive_pscomppars.csv"
DERIVED = ROOT / "data" / "derived"
FIGURES = ROOT / "figures"
DOCS = ROOT / "docs"
DERIVED.mkdir(parents=True, exist_ok=True)
FIGURES.mkdir(parents=True, exist_ok=True)
DOCS.mkdir(parents=True, exist_ok=True)

EARTH_RADIUS = 1.0
EARTH_EQT_K = 255.0
PROFILES = {
    "balanced": (1/3, 1/3, 1/3),
    "earth_analog": (0.45, 0.45, 0.10),
    "exploration_distance": (0.15, 0.15, 0.70),
    "size_first": (0.70, 0.15, 0.15),
    "temperature_first": (0.15, 0.70, 0.15),
}

df = pd.read_csv(RAW)
for c in ["pl_rade", "pl_eqt", "sy_dist", "pl_radelim", "pl_eqtlim"]:
    df[c] = pd.to_numeric(df[c], errors="coerce")

valid_radius = df["pl_rade"].gt(0) & df["pl_radelim"].fillna(0).eq(0)
valid_temp = df["pl_eqt"].gt(0) & df["pl_eqtlim"].fillna(0).eq(0)
valid_dist = df["sy_dist"].gt(0)

# Log deviation treats equal fractional departures above/below Earth symmetrically.
df["size_deviation_log"] = np.where(
    valid_radius, np.abs(np.log(df["pl_rade"] / EARTH_RADIUS)), np.nan)
df["temperature_deviation_log"] = np.where(
    valid_temp, np.abs(np.log(df["pl_eqt"] / EARTH_EQT_K)), np.nan)
df["distance_pc"] = np.where(valid_dist, df["sy_dist"], np.nan)

# method='min' gives every exact tie the same winning rank.
df["size_rank"] = df["size_deviation_log"].rank(method="min")
df["temperature_rank"] = df["temperature_deviation_log"].rank(method="min")
df["distance_rank"] = df["distance_pc"].rank(method="min")

complete_mask = valid_radius & valid_temp & valid_dist
for metric in ["size_deviation_log", "temperature_deviation_log", "distance_pc"]:
    pct = df.loc[complete_mask, metric].rank(method="average", pct=True)
    df.loc[complete_mask, metric + "_percentile"] = pct

for name, (w_size, w_temp, w_dist) in PROFILES.items():
    df.loc[complete_mask, name + "_score"] = (
        w_size * df.loc[complete_mask, "size_deviation_log_percentile"] +
        w_temp * df.loc[complete_mask, "temperature_deviation_log_percentile"] +
        w_dist * df.loc[complete_mask, "distance_pc_percentile"])
    df.loc[complete_mask, name + "_rank"] = (
        df.loc[complete_mask, name + "_score"].rank(method="min"))

df.to_csv(DERIVED / "calculated_planet_metrics.csv", index=False)

ranking_specs = [
    ("size_only", "size_rank"),
    ("temperature_only", "temperature_rank"),
    ("distance_only", "distance_rank"),
] + [(name, name + "_rank") for name in PROFILES]
top_frames = []
tie_frames = []
for label, rank_col in ranking_specs:
    ranked = df[df[rank_col].notna()].sort_values(
        [rank_col, "pl_name"], kind="stable").copy()
    top = ranked[ranked[rank_col] <= 10].copy()
    top.insert(0, "ranking_definition", label)
    top.insert(1, "reported_rank", top[rank_col])
    top_frames.append(top)
    winners = ranked[ranked[rank_col] == ranked[rank_col].min()].copy()
    winners.insert(0, "ranking_definition", label)
    winners.insert(1, "tie_count", len(winners))
    tie_frames.append(winners)

selected = ["ranking_definition", "reported_rank", "pl_name", "hostname",
            "pl_rade", "pl_eqt", "sy_dist", "size_deviation_log",
            "temperature_deviation_log", "distance_pc"] +            [x + "_score" for x in PROFILES]
pd.concat(top_frames, ignore_index=True)[selected].to_csv(
    DERIVED / "top_rankings.csv", index=False)
tie_selected = ["ranking_definition", "tie_count", "pl_name", "hostname",
                "pl_rade", "pl_eqt", "sy_dist"]
ties = pd.concat(tie_frames, ignore_index=True)[tie_selected]
ties.to_csv(DERIVED / "winner_ties.csv", index=False)

# Sensitivity figure: show all co-winners, not an arbitrary single row.
fig, ax = plt.subplots(figsize=(10.5, 6.0))
labels, annotations = [], []
for label, rank_col in ranking_specs:
    winners = df[df[rank_col] == df[rank_col].min()].sort_values("pl_name")
    labels.append(label.replace("_", " ").title())
    names = ", ".join(winners["pl_name"].tolist())
    annotations.append(names if len(names) <= 90 else names[:87] + "…")
y = np.arange(len(labels))
ax.scatter(np.zeros(len(y)), y, s=85, color="#2b8cbe")
for yi, txt in zip(y, annotations):
    ax.text(0.03, yi, txt, va="center", fontsize=9.5)
ax.set_yticks(y, labels)
ax.set_xlim(-0.02, 1)
ax.set_xticks([])
ax.set_title("Winner sensitivity to the definition of ‘Earth-like’", pad=14)
for spine in ax.spines.values():
    spine.set_visible(False)
fig.tight_layout()
fig.savefig(FIGURES / "ranking_winners_sensitivity.png", dpi=220,
            bbox_inches="tight")
plt.close(fig)

# Top balanced candidates figure.
balanced = df[df["balanced_rank"].notna()].nsmallest(10, "balanced_score").copy()
balanced = balanced.sort_values("balanced_score", ascending=True)
fig, ax = plt.subplots(figsize=(8.5, 5.5))
ax.barh(balanced["pl_name"], balanced["balanced_score"], color="#41ab5d")
ax.invert_yaxis()
ax.set_xlabel("Balanced percentile score (lower is better)")
ax.set_title("Top ten under equal size–temperature–distance weighting")
fig.tight_layout()
fig.savefig(FIGURES / "balanced_top10.png", dpi=220, bbox_inches="tight")
plt.close(fig)

tie_lines = []
for definition, group in ties.groupby("ranking_definition", sort=False):
    names = ", ".join(group["pl_name"])
    tie_lines.append(f"- **{definition.replace('_',' ').title()}**: {len(group)} winner(s) — {names}")

winner_rows = []
for label, rank_col in ranking_specs:
    w = df[df[rank_col] == df[rank_col].min()].sort_values("pl_name")
    winner_rows.append((label, "; ".join(w.pl_name), len(w)))

report = f"""# Which planet is most like Earth?

## Result

Changing the value function changes the winner. The equal-weighted ranking favors
**{df.loc[df['balanced_rank'].idxmin(), 'pl_name']}**; emphasizing Earth-like size
and equilibrium temperature favors **{df.loc[df['earth_analog_rank'].idxmin(), 'pl_name']}**;
and emphasizing proximity favors **{df.loc[df['exploration_distance_rank'].idxmin(), 'pl_name']}**.

The archive snapshot contains **{len(df):,}** confirmed planets. Radius is usable
for **{valid_radius.sum():,}**, equilibrium temperature for **{valid_temp.sum():,}**,
distance for **{valid_dist.sum():,}**, and all three for **{complete_mask.sum():,}**.

## Winner and tie audit

{chr(10).join(tie_lines)}

The size-only result is therefore an eight-way tie, not a unique win. Several
planets also tie on distance because planets in the same stellar system share
the same system distance. Composite profiles have unique winners in this snapshot.

## Scoring

- Size deviation: `abs(ln(radius / 1 R_earth))`
- Temperature deviation: `abs(ln(T_eq / 255 K))`
- Distance: archive system distance in parsecs
- Combined criteria are percentile-ranked over complete cases.
- Balanced: 1/3 each.
- Earth analogue: 45% size, 45% temperature, 10% distance.
- Exploration distance: 15%, 15%, 70%.
- Size first: 70%, 15%, 15%.
- Temperature first: 15%, 70%, 15%.
- Lower scores are better. Exact ties use minimum rank.

## Interpretation

These are rankings of three archive fields, not habitability estimates.
Equilibrium temperature is a black-body model rather than surface temperature.
Radius can be measured or model-derived. Distance measures accessibility, not
biological suitability. Atmospheres, stellar activity, composition, water,
mass, tidal state, and parameter uncertainty are not included.

The PSCompPars table has one relatively complete row per confirmed planet, but
its composite parameters may come from different references and need not form
a self-consistent physical solution.

## Reproduction

From the package root:

```sh
python -m pip install -r requirements.txt
python code/analyze_rankings.py
```

The script reads only the preserved raw CSV and regenerates the calculated
table, rankings, tie audit, figures, and this report.
"""
(DOCS / "full_report.md").write_text(report, encoding="utf-8")

run_summary = {
    "raw_rows": int(len(df)),
    "usable_radius": int(valid_radius.sum()),
    "usable_temperature": int(valid_temp.sum()),
    "usable_distance": int(valid_dist.sum()),
    "complete_cases": int(complete_mask.sum()),
    "winner_ties": {
        label: df.loc[df[rank_col] == df[rank_col].min(), "pl_name"].sort_values().tolist()
        for label, rank_col in ranking_specs
    }
}
(DERIVED / "run_summary.json").write_text(
    json.dumps(run_summary, indent=2), encoding="utf-8")
print(json.dumps(run_summary, indent=2))
