# -*- coding: utf-8 -*-
"""
calculate_accessibility.py

This script reads Summary_Table.csv from the current directory,
computes accessibility scores for models M0-M3 under multiple time
thresholds, and saves the final accessibility matrix to
results/accessibility_result.csv.
"""

import os
import warnings

import numpy as np
import pandas as pd

warnings.filterwarnings("ignore")


# Define file locations relative to the current script.
BASE_DIR = os.path.dirname(os.path.abspath(__file__))
DATA_PATH = os.path.join(BASE_DIR, "Summary_Table.csv")
RESULT_DIR = os.path.join(BASE_DIR, "results")
RESULT_PATH = os.path.join(RESULT_DIR, "accessibility_result.csv")

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

OBJECTIVE_COLUMNS = ["area", "transportation", "infrastructure", "ndvi"]
ALL_SUPPLY_COLUMNS = OBJECTIVE_COLUMNS + ["emotion"]
THRESHOLDS = [15, 30, 45, 60]
SCALE = 1e6


def validate_columns(df):
    """Check whether the input table contains all required fields."""
    missing_columns = [col for col in REQUIRED_COLUMNS if col not in df.columns]
    if missing_columns:
        raise ValueError(
            "Summary_Table.csv is missing required columns: "
            + ", ".join(missing_columns)
        )


def critic_weight(df, columns):
    """Compute CRITIC weights from indicator variability and conflict."""
    data = df[columns].values
    n_features = 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)
    conflict = np.sum(1 - np.abs(corr), axis=0)
    c_values = std * conflict

    if np.sum(c_values) == 0:
        weights = np.ones(n_features) / n_features
    else:
        weights = c_values / np.sum(c_values)

    return dict(zip(columns, weights))


def compute_supply(df_parks, weights, columns):
    """Aggregate normalized park indicators into a single supply score."""
    return sum(weights[col] * df_parks[col] for col in columns)


def vec_gaussian_decay(distance_series, threshold):
    """Apply the Gaussian distance-decay function within a time threshold."""
    result = np.zeros_like(distance_series, dtype=float)
    mask = distance_series <= threshold
    valid_distance = distance_series[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_accessibility_with_rj(
    df_od,
    supply_dict,
    park_ids,
    community_ids,
    df_parks,
    review_count_mean,
    demand_type="static",
    alpha=1.0,
    threshold=15,
):
    """
    Compute community accessibility using the two-step floating catchment area method.

    Step 1 calculates the park-level supply-demand ratio Rj.
    Step 2 sums the reachable park ratios for each community with distance decay.
    """
    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

        # Dynamic demand scales park review volume by the baseline demand level.
        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 = alpha * beta * review_count

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

    accessibility_values = []

    for community_id in community_ids:
        od_community = df_od[
            (df_od["OriginID"] == community_id) & (df_od["cost"] <= threshold)
        ].copy()
        if od_community.empty:
            accessibility_values.append(0.0)
            continue

        od_community["G"] = vec_gaussian_decay(
            od_community["cost"].to_numpy(), threshold
        )
        accessibility = 0.0

        for _, row in od_community.iterrows():
            park_id = row["DestinationID"]
            if park_id in rj_dict:
                accessibility += rj_dict[park_id] * row["G"]

        accessibility_values.append(accessibility)

    return pd.Series(accessibility_values, index=community_ids)


def main():
    # Load the OD summary table and validate the required schema.
    df = pd.read_csv(DATA_PATH)
    validate_columns(df)

    df["OriginID"] = df["OriginID"].astype(int)
    df["DestinationID"] = df["DestinationID"].astype(int)

    park_columns = [
        "DestinationID",
        "area",
        "review_count",
        "transportation",
        "infrastructure",
        "ndvi",
        "pos_ratio",
        "neg_ratio",
    ]
    df_parks = (
        df[park_columns]
        .drop_duplicates(subset="DestinationID")
        .set_index("DestinationID")
    )
    df_parks["review_count"] = df_parks["review_count"].fillna(0)
    # Emotion is defined as the difference between positive and negative ratios.
    df_parks["emotion"] = df_parks["pos_ratio"] - df_parks["neg_ratio"]

    community_population = (
        df[["OriginID", "population"]]
        .drop_duplicates(subset="OriginID")
        .set_index("OriginID")["population"]
    )
    community_ids = sorted(df["OriginID"].unique())
    park_ids = sorted(df["DestinationID"].unique())

    # Normalize objective supply indicators to keep them on a comparable scale.
    for column in OBJECTIVE_COLUMNS:
        min_value = df_parks[column].min()
        max_value = df_parks[column].max()
        df_parks[f"{column}_norm"] = (df_parks[column] - min_value) / (
            max_value - min_value
        )

    # Normalize the emotion indicator separately before weighting.
    min_emotion = df_parks["emotion"].min()
    max_emotion = df_parks["emotion"].max()
    df_parks["emotion_norm"] = (df_parks["emotion"] - min_emotion) / (
        max_emotion - min_emotion
    )

    objective_columns_norm = [f"{col}_norm" for col in OBJECTIVE_COLUMNS]
    all_supply_columns_norm = objective_columns_norm + ["emotion_norm"]

    weights_without_emotion = critic_weight(df_parks, objective_columns_norm)
    weights_with_emotion = critic_weight(df_parks, all_supply_columns_norm)

    # M0/M2 use objective supply only, while M1/M3 include emotion-adjusted supply.
    supply_m0 = compute_supply(
        df_parks, weights_without_emotion, objective_columns_norm
    )
    supply_m1 = compute_supply(df_parks, weights_with_emotion, all_supply_columns_norm)
    supply_m2 = supply_m0
    supply_m3 = supply_m1

    review_count_mean = df_parks["review_count"].mean()

    # Build the OD table used in the 2SFCA computation.
    df_od = df[["OriginID", "DestinationID", "cost"]].copy()
    df_od = df_od.merge(
        community_population.rename("population"),
        left_on="OriginID",
        right_index=True,
        how="left",
    )

    models = {
        "M0": {"supply": supply_m0.to_dict(), "demand": "static"},
        "M1": {"supply": supply_m1.to_dict(), "demand": "static"},
        "M2": {"supply": supply_m2.to_dict(), "demand": "dynamic"},
        "M3": {"supply": supply_m3.to_dict(), "demand": "dynamic"},
    }

    # Compute accessibility for each model under each travel-time threshold.
    accessibility_results = {}
    for model_name, params in models.items():
        for threshold in THRESHOLDS:
            column_name = f"A_{model_name}_{threshold}min"
            accessibility_results[column_name] = compute_accessibility_with_rj(
                df_od=df_od,
                supply_dict=params["supply"],
                park_ids=park_ids,
                community_ids=community_ids,
                df_parks=df_parks,
                review_count_mean=review_count_mean,
                demand_type=params["demand"],
                alpha=1.0,
                threshold=threshold,
            )

    result_df = pd.DataFrame(accessibility_results, index=community_ids)
    # Keep the same output scale as the original script.
    result_df *= SCALE
    result_df.columns = [col.replace("min", "") for col in result_df.columns]
    result_df.index.name = "OriginID"
    os.makedirs(RESULT_DIR, exist_ok=True)
    result_df.reset_index().to_csv(RESULT_PATH, index=False, encoding="utf-8-sig")

    print(f"Result file saved to: {RESULT_PATH}")


if __name__ == "__main__":
    main()
