"""Replication script for F1-dev-01. Run from workspace root:
    ~/wp29/venv/bin/python out/analysis.py
"""
import json, os
import pandas as pd

os.makedirs("out", exist_ok=True)
raw = pd.read_csv("data/firms.csv", dtype=str, keep_default_na=False)
n_rows = len(raw)
for c in raw.columns:
    raw[c] = raw[c].str.strip()

# Baseline: raw revenue mean over non-empty cells, mixed units, before cleaning
rev_raw = pd.to_numeric(raw["revenue"].replace("", None), errors="coerce")
mean_revenue_raw = float(rev_raw.dropna().mean())

ledger, issues = [], []
def step(rule, before, after, rationale):
    ledger.append(dict(step=len(ledger) + 1, rule=rule, n_before=before,
                       n_after=after, n_removed=before - after, rationale=rationale))

ledger_start = n_rows
step("read data/firms.csv", n_rows, n_rows, "All data rows read (header excluded).")

# 1. exact duplicate rows
n_exact = int(raw.duplicated().sum())
d = raw.drop_duplicates().copy()
step("drop exact duplicate rows", n_rows, len(d),
     "Rows identical in every column are assembly duplicates of the same firm; keep first occurrence.")
issues.append({"issue": f"{n_exact} rows are exact duplicates of another row (same firm_id and all values)",
               "handling": "Dropped duplicates, keeping one row per firm."})
# check remaining id conflicts
n_conf = int(d["firm_id"].duplicated().sum())
assert n_conf == 0, "conflicting duplicate firm_ids present"
n_unique = int(d["firm_id"].nunique())

# 2. unit recoding
units = d["unit"].value_counts().to_dict()
assert set(units) <= {"USD_millions", "USD_thousands"}, units
d["revenue_num"] = pd.to_numeric(d["revenue"].replace("", None), errors="coerce")
d["revenue_musd"] = d["revenue_num"].where(d["unit"] == "USD_millions", d["revenue_num"] / 1000.0)
step("convert revenue to millions USD", len(d), len(d),
     "Rows with unit USD_thousands divided by 1000; USD_millions unchanged. No rows removed.")
issues.append({"issue": f"Mixed revenue units: {units.get('USD_thousands',0)} distinct-firm rows in USD_thousands, "
                        f"{units.get('USD_millions',0)} in USD_millions",
               "handling": "Converted USD_thousands to millions USD (divide by 1000)."})

# employees: complete over all firms
d["employees"] = pd.to_numeric(d["employees"], errors="raise").astype(int)
mean_employees = float(d["employees"].mean())

# 3. missing revenue
miss = d["revenue_musd"].isna()
n_missing = int(miss.sum())
missing_ids = sorted(d.loc[miss, "firm_id"])
before = len(d)
a = d.loc[~miss].copy()
step("exclude firms with unknown revenue", before, len(a),
     "Empty revenue cell; cannot contribute to revenue aggregates. Firms: " + ";".join(missing_ids))
issues.append({"issue": f"{n_missing} distinct firms have empty revenue ({', '.join(missing_ids)})",
               "handling": "Excluded from revenue aggregates and listed in the exclusion ledger; retained for mean_employees."})

assert (a["revenue_musd"] > 0).all()
a = a.sort_values("firm_id")
a[["firm_id", "sector", "revenue_musd", "employees"]].to_csv("out/analytic.csv", index=False)
pd.DataFrame(ledger).to_csv("out/exclusions.csv", index=False)

res = {
    "task_id": "F1-dev-01",
    "n_unique_firms": n_unique,
    "n_missing_revenue": n_missing,
    "n_analysed": int(len(a)),
    "mean_revenue_musd": round(float(a["revenue_musd"].mean()), 6),
    "median_revenue_musd": round(float(a["revenue_musd"].median()), 6),
    "mean_employees": round(mean_employees, 6),
    "mean_revenue_raw": round(mean_revenue_raw, 6),
    "feasibility": "complete",
    "feasibility_reason": None,
    "data_issues": issues,
    "n_input_rows": {"firms.csv": int(n_rows)},
    "random_seed": None,
}
with open("out/results.json", "w") as f:
    json.dump(res, f, indent=2)
print(json.dumps(res, indent=2))

# Descriptive uncertainty (printed only, quoted in report; not part of the results schema)
import numpy as np
x = a["revenue_musd"].to_numpy()
se = x.std(ddof=1) / np.sqrt(len(x))
print(f"SE(mean revenue) = {se:.3f} MUSD; normal 95% CI = [{x.mean()-1.96*se:.3f}, {x.mean()+1.96*se:.3f}]")
print(f"SD revenue = {x.std(ddof=1):.3f}; min = {x.min():.3f}; max = {x.max():.3f}")
e = d["employees"].to_numpy(); print(f"SE(mean employees) = {e.std(ddof=1)/np.sqrt(len(e)):.3f}")