# -*- coding: utf-8 -*-
"""
Generate CSV tables for accessibility analysis.

This script reads community-level accessibility outputs from
results/accessibility_result.csv and recalculates park-level supply and
supply-demand ratios from Summary_Table.csv when needed.
"""

import os
import warnings

import numpy as np
import pandas as pd
from scipy.stats import kurtosis, skew, spearmanr

warnings.filterwarnings("ignore")


BASE_DIR = os.path.dirname(os.path.abspath(__file__))
SUMMARY_PATH = os.path.join(BASE_DIR, "Summary_Table.csv")
ACCESSIBILITY_PATH = os.path.join(BASE_DIR, "results", "accessibility_result.csv")
RESULT_DIR = os.path.join(BASE_DIR, "results")

THRESHOLDS = [15, 30, 45, 60]
GRADE_BINS = [0, 1.244, 3.003, 5.556, 8.362, float("inf")]
GRADE_LABELS = [
    "No accessibility",
    "Level 1",
    "Level 2",
    "Level 3",
    "Level 4",
    "Level 5",
]

SUMMARY_REQUIRED_COLUMNS = [
    "OriginID",
    "DestinationID",
    "cost",
    "population",
    "area",
    "transportation",
    "infrastructure",
    "ndvi",
    "pos_ratio",
    "neg_ratio",
    "review_count",
]

ACCESSIBILITY_REQUIRED_COLUMNS = [
    "OriginID",
    "A_M0_15",
    "A_M0_30",
    "A_M0_45",
    "A_M0_60",
    "A_M1_15",
    "A_M1_30",
    "A_M1_45",
    "A_M1_60",
    "A_M2_15",
    "A_M2_30",
    "A_M2_45",
    "A_M2_60",
    "A_M3_15",
    "A_M3_30",
    "A_M3_45",
    "A_M3_60",
]

OBJECTIVE_COLUMNS = ["area", "transportation", "infrastructure", "ndvi"]


def validate_columns(df, required_columns, source_name):
    missing_columns = [column for column in required_columns if column not in df.columns]
    if missing_columns:
        raise ValueError(
            f"{source_name} is missing required columns: "
            + ", ".join(missing_columns)
        )


def normalize_series(series):
    min_value = series.min()
    max_value = series.max()
    if pd.isna(min_value) or pd.isna(max_value) or max_value == min_value:
        return pd.Series(np.zeros(len(series)), index=series.index, dtype=float)
    return (series - min_value) / (max_value - min_value)


def critic_weight(df, columns):
    data = df[columns].values
    feature_count = len(columns)
    std = np.std(data, axis=0, ddof=0)
    corr = pd.DataFrame(data).corr(method="spearman").values
    corr = np.nan_to_num(corr, nan=0.0)
    np.fill_diagonal(corr, 0.0)
    conflict = np.sum(1 - np.abs(corr), axis=0)
    critic_score = std * conflict

    if np.sum(critic_score) == 0:
        weights = np.ones(feature_count) / feature_count
    else:
        weights = critic_score / np.sum(critic_score)

    return dict(zip(columns, weights))


def compute_supply(df_parks, weights, columns):
    return sum(weights[column] * df_parks[column] for column in columns)


def vec_gaussian_decay(distance_series, threshold):
    distance_array = np.asarray(distance_series, dtype=float)
    result = np.zeros_like(distance_array, dtype=float)
    mask = distance_array <= threshold
    valid_distance = distance_array[mask]
    result[mask] = (
        np.exp(-0.5 * (valid_distance / threshold) ** 2) - np.exp(-0.5)
    ) / (1 - np.exp(-0.5))
    return np.clip(result, 0, 1)


def compute_rj(df_od, df_parks, park_ids, supply_dict, review_count_mean, threshold, demand_type):
    rj_dict = {}

    for park_id in park_ids:
        od_park = df_od[
            (df_od["DestinationID"] == park_id) & (df_od["cost"] <= threshold)
        ].copy()
        if od_park.empty:
            continue

        od_park["G"] = vec_gaussian_decay(od_park["cost"].to_numpy(), threshold)
        static_demand = np.sum(od_park["population"] * od_park["G"])
        dynamic_demand = 0.0

        if demand_type == "dynamic":
            review_count = df_parks.loc[park_id, "review_count"]
            if static_demand > 0 and review_count_mean > 0:
                beta = static_demand / review_count_mean
                dynamic_demand = beta * review_count

        total_demand = static_demand + dynamic_demand
        rj_dict[park_id] = supply_dict[park_id] / total_demand if total_demand > 0 else 0.0

    return rj_dict


def l1_diff(a, b):
    return np.sum(np.abs(np.asarray(a) - np.asarray(b)))


def mse(a, b):
    return np.mean((np.asarray(a) - np.asarray(b)) ** 2)


def mean_relative_change(a, b):
    a_array = np.asarray(a, dtype=float)
    b_array = np.asarray(b, dtype=float)
    return np.mean(np.abs((b_array - a_array) / (a_array + 1e-8))) * 100


def cohens_d(a, b):
    a_array = np.asarray(a, dtype=float)
    b_array = np.asarray(b, dtype=float)
    pooled_std = np.sqrt(
        (
            (len(a_array) - 1) * np.var(a_array, ddof=1)
            + (len(b_array) - 1) * np.var(b_array, ddof=1)
        )
        / (len(a_array) + len(b_array) - 2)
    )
    if pooled_std == 0 or np.isnan(pooled_std):
        return 0.0
    return (np.mean(b_array) - np.mean(a_array)) / pooled_std


def safe_change_rate(new_value, old_value):
    if pd.isna(old_value) or old_value == 0:
        return np.nan
    return (new_value - old_value) / old_value * 100


def classify_accessibility(series):
    levels = pd.Series("No accessibility", index=series.index, dtype=object)
    positive_mask = series > 0
    levels.loc[positive_mask] = pd.cut(
        series.loc[positive_mask],
        bins=GRADE_BINS,
        labels=GRADE_LABELS[1:],
        right=False,
        include_lowest=False,
    ).astype(str)
    return pd.Categorical(levels, categories=GRADE_LABELS, ordered=True)


def main():
    os.makedirs(RESULT_DIR, exist_ok=True)

    summary_df = pd.read_csv(SUMMARY_PATH)
    accessibility_df = pd.read_csv(ACCESSIBILITY_PATH)

    validate_columns(summary_df, SUMMARY_REQUIRED_COLUMNS, "Summary_Table.csv")
    validate_columns(
        accessibility_df,
        ACCESSIBILITY_REQUIRED_COLUMNS,
        "results/accessibility_result.csv",
    )

    summary_df["OriginID"] = summary_df["OriginID"].astype(int)
    summary_df["DestinationID"] = summary_df["DestinationID"].astype(int)
    accessibility_df["OriginID"] = accessibility_df["OriginID"].astype(int)
    accessibility_df = accessibility_df.set_index("OriginID").sort_index()

    park_columns = [
        "DestinationID",
        "area",
        "review_count",
        "transportation",
        "infrastructure",
        "ndvi",
        "pos_ratio",
        "neg_ratio",
    ]
    df_parks = (
        summary_df[park_columns]
        .drop_duplicates(subset="DestinationID")
        .set_index("DestinationID")
        .sort_index()
    )
    df_parks["review_count"] = df_parks["review_count"].fillna(0)
    df_parks["sentiment_score"] = df_parks["pos_ratio"] - df_parks["neg_ratio"]

    for column in OBJECTIVE_COLUMNS:
        df_parks[f"{column}_norm"] = normalize_series(df_parks[column])
    df_parks["sentiment_score_norm"] = normalize_series(df_parks["sentiment_score"])

    objective_columns_norm = [f"{column}_norm" for column in OBJECTIVE_COLUMNS]
    supply_columns_norm = objective_columns_norm + ["sentiment_score_norm"]

    weights_without_sentiment = critic_weight(df_parks, objective_columns_norm)
    weights_with_sentiment = critic_weight(df_parks, supply_columns_norm)

    supply_m0 = compute_supply(df_parks, weights_without_sentiment, objective_columns_norm)
    supply_m1 = compute_supply(df_parks, weights_with_sentiment, supply_columns_norm)

    community_population = (
        summary_df[["OriginID", "population"]]
        .drop_duplicates(subset="OriginID")
        .set_index("OriginID")["population"]
    )
    df_od = summary_df[["OriginID", "DestinationID", "cost"]].copy()
    df_od = df_od.merge(
        community_population.rename("population"),
        left_on="OriginID",
        right_index=True,
        how="left",
    )
    park_ids = sorted(df_parks.index.tolist())
    review_count_mean = df_parks["review_count"].mean()

    rj_m0 = compute_rj(
        df_od=df_od,
        df_parks=df_parks,
        park_ids=park_ids,
        supply_dict=supply_m0.to_dict(),
        review_count_mean=review_count_mean,
        threshold=15,
        demand_type="static",
    )
    rj_m2 = compute_rj(
        df_od=df_od,
        df_parks=df_parks,
        park_ids=park_ids,
        supply_dict=supply_m0.to_dict(),
        review_count_mean=review_count_mean,
        threshold=15,
        demand_type="dynamic",
    )

    stats_rows = []
    for model in ["M0", "M1", "M2", "M3"]:
        series = accessibility_df[f"A_{model}_15"]
        mean_value = series.mean()
        std_value = series.std()
        stats_rows.append(
            {
                "Model": model,
                "Mean": mean_value,
                "SD": std_value,
                "CV": std_value / mean_value if mean_value != 0 else np.nan,
                "Skewness": skew(series),
                "Kurtosis": kurtosis(series, fisher=False),
            }
        )
    stats_15_df = pd.DataFrame(stats_rows)
    stats_15_df.to_csv(
        os.path.join(
            RESULT_DIR,
            "Fig.3. Comparison of statistical characteristics of accessibility across four models _15 min.csv",
        ),
        index=False,
        encoding="utf-8-sig",
    )

    base_15 = accessibility_df["A_M0_15"].values
    difference_rows = []
    for model in ["M1", "M2", "M3"]:
        comparison_values = accessibility_df[f"A_{model}_15"].values
        difference_rows.append(
            {
                "Comparison": f"{model} vs M0",
                "L1": l1_diff(base_15, comparison_values),
                "MSE": mse(base_15, comparison_values),
                "Mean relative change": mean_relative_change(base_15, comparison_values),
            }
        )
    difference_df = pd.DataFrame(difference_rows)
    difference_df.to_csv(
        os.path.join(
            RESULT_DIR,
            "Fig.5.Comparison of quantitative difference indicators between each correction model and baseline model_15 min.csv",
        ),
        index=False,
        encoding="utf-8-sig",
    )

    supply_change_df = pd.DataFrame(
        {
            "Park ID": supply_m0.index,
            "M0 Supply": supply_m0.values,
            "M1 Supply": supply_m1.values,
            "Change Rate (%)": [
                safe_change_rate(new_value, old_value)
                for old_value, new_value in zip(supply_m0.values, supply_m1.values)
            ],
            "Sentiment Score": df_parks.loc[supply_m0.index, "sentiment_score"].values,
        }
    )
    supply_change_df = supply_change_df.sort_values(
        "Change Rate (%)",
        key=lambda series: series.abs(),
        ascending=False,
    )
    supply_change_df.to_csv(
        os.path.join(
            RESULT_DIR,
            "Table 3 Changes in comprehensive park service capacity after perception feedback correction_15 min.csv",
        ),
        index=False,
        encoding="utf-8-sig",
    )

    rj_m0_series = pd.Series(rj_m0, name="M0_Rj").sort_index()
    rj_m2_series = pd.Series(rj_m2, name="M2_Rj").sort_index()
    rj_change_df = pd.DataFrame(
        {
            "Park ID": rj_m0_series.index,
            "Review Count": df_parks.loc[rj_m0_series.index, "review_count"].values,
            "M0_Rj/10⁻⁶": rj_m0_series.values * 1e6,
            "M2_Rj/10⁻⁶": rj_m2_series.reindex(rj_m0_series.index, fill_value=0.0).values
            * 1e6,
            "Change Rate (%)": [
                safe_change_rate(
                    new_value,
                    old_value,
                )
                for old_value, new_value in zip(
                    rj_m0_series.values,
                    rj_m2_series.reindex(rj_m0_series.index, fill_value=0.0).values,
                )
            ],
        }
    )
    rj_change_df = rj_change_df.sort_values(
        "Change Rate (%)",
        key=lambda series: series.abs(),
        ascending=False,
    )
    rj_change_df.to_csv(
        os.path.join(
            RESULT_DIR,
            "Table 4 Changes in park supply–demand ratio after usage popularity correction_15 min.csv",
        ),
        index=False,
        encoding="utf-8-sig",
    )

    m1_15 = accessibility_df["A_M1_15"].values
    m2_15 = accessibility_df["A_M2_15"].values
    m3_15 = accessibility_df["A_M3_15"].values
    independent_perception = l1_diff(base_15, m1_15)
    independent_usage = l1_diff(base_15, m2_15)
    combined_effect = l1_diff(base_15, m3_15)
    residual = combined_effect - independent_perception - independent_usage
    decomp_df = pd.DataFrame(
        {
            "Change Source": [
                "Independent contribution of perception feedback",
                "Independent contribution of usage popularity",
                "Combined effect",
                "Residual",
                "Perception feedback contribution",
                "Usage popularity contribution",
                "Residual contribution",
            ],
            "Value": [
                independent_perception,
                independent_usage,
                combined_effect,
                residual,
                independent_perception / combined_effect if combined_effect != 0 else np.nan,
                independent_usage / combined_effect if combined_effect != 0 else np.nan,
                residual / combined_effect if combined_effect != 0 else np.nan,
            ],
        }
    )
    decomp_df.to_csv(
        os.path.join(
            RESULT_DIR,
            "Table 5 Decomposition of comprehensive effects of bilateral supply–demand correction_15 min.csv",
        ),
        index=False,
        encoding="utf-8-sig",
    )

    level_m0 = classify_accessibility(accessibility_df["A_M0_15"])
    level_m3 = classify_accessibility(accessibility_df["A_M3_15"])
    transition_df = pd.crosstab(
        pd.Series(level_m0, index=accessibility_df.index, name="M0 Grade"),
        pd.Series(level_m3, index=accessibility_df.index, name="M3 Grade"),
        dropna=False,
    )
    transition_df.to_csv(
        os.path.join(
            RESULT_DIR,
            "Table 6 Grade transition matrix between baseline model (M0) and behavioral perception model (M3).csv",
        ),
        encoding="utf-8-sig",
    )

    multi_time_rows = []
    for threshold in THRESHOLDS:
        base_series = accessibility_df[f"A_M0_{threshold}"]
        for model in ["M0", "M1", "M2", "M3"]:
            model_series = accessibility_df[f"A_{model}_{threshold}"]
            mean_value = model_series.mean()
            std_value = model_series.std()
            multi_time_rows.extend(
                [
                    {
                        "Threshold (min)": threshold,
                        "Model": model,
                        "Metric": "Mean",
                        "Value": mean_value,
                    },
                    {
                        "Threshold (min)": threshold,
                        "Model": model,
                        "Metric": "CV",
                        "Value": std_value / mean_value if mean_value != 0 else np.nan,
                    },
                    {
                        "Threshold (min)": threshold,
                        "Model": model,
                        "Metric": "Skewness",
                        "Value": skew(model_series),
                    },
                ]
            )
        for model in ["M1", "M2", "M3"]:
            multi_time_rows.append(
                {
                    "Threshold (min)": threshold,
                    "Model": f"{model} vs M0",
                    "Metric": "L1 Total Difference",
                    "Value": l1_diff(base_series.values, accessibility_df[f"A_{model}_{threshold}"].values),
                }
            )
    multi_time_df = pd.DataFrame(multi_time_rows)
    multi_time_df.to_csv(
        os.path.join(
            RESULT_DIR,
            "Fig.10.Trends of statistical metrics for each model under different time thresholds.csv",
        ),
        index=False,
        encoding="utf-8-sig",
    )

    effect_rows = []
    for threshold in THRESHOLDS:
        base_values = accessibility_df[f"A_M0_{threshold}"].values
        values_m1 = accessibility_df[f"A_M1_{threshold}"].values
        values_m2 = accessibility_df[f"A_M2_{threshold}"].values
        values_m3 = accessibility_df[f"A_M3_{threshold}"].values
        effect_rows.append(
            {
                "Threshold (min)": threshold,
                "d_perception": cohens_d(base_values, values_m1),
                "d_usage": cohens_d(base_values, values_m2),
                "d_integrated": cohens_d(base_values, values_m3),
                "ρ_perception": spearmanr(base_values, values_m1)[0],
                "ρ_usage": spearmanr(base_values, values_m2)[0],
                "ρ_integrated": spearmanr(base_values, values_m3)[0],
            }
        )
    effect_df = pd.DataFrame(effect_rows)
    effect_df.to_csv(
        os.path.join(
            RESULT_DIR,
            "Fig.11.Effect sizes and Spearman rank correlation coefficients for each model under different time thresholds.csv",
        ),
        index=False,
        encoding="utf-8-sig",
    )

    print("Analysis tables saved to:", RESULT_DIR)


if __name__ == "__main__":
    main()
