"""
Reproducibility package for:
Maladaptive Personality Traits and Problematic Pornography Use:
Serial Indirect Associations Through Personality Functioning and Emotional Avoidance

Expected input file in the SAME FOLDER:
    Personality_Pornography_Project_Data.csv

Outputs created automatically:
    Personality_Pornography_Project_Output.txt
    Personality_Pornography_Project_TableS2.csv
    Personality_Pornography_Project_TableS2_Descriptives.csv
    Personality_Pornography_Project_TableS3.csv
    Personality_Pornography_Project_TableS4.csv
    Personality_Pornography_Project_TableS5.csv
    Personality_Pornography_Project_TableS6.csv
    Personality_Pornography_Project_TableS7.csv

Historical screening note:
The retained dataset contains N=591. Project records document 625 submitted
questionnaires and 34 exclusions for same-option/highly uniform responding.
Granular historical response-speed and long-string metrics were not archived,
so those exclusions cannot be algorithmically recreated from the retained file.
No post-hoc screening threshold is introduced in this script.
"""

from pathlib import Path
import sys, platform
import numpy as np
import pandas as pd
import scipy
from scipy.stats import norm, skew, kurtosis, pearsonr
import statsmodels
import statsmodels.api as sm
from statsmodels.stats.outliers_influence import variance_inflation_factor

BASE = Path(__file__).resolve().parent
DATA_FILE = BASE / "Personality_Pornography_Project_Data.csv"

B = 5000
SEED_PPUS12 = 20260825
SEED_PPUS9 = 20260826
SEED_ALT = 20260827
SEED_ADJ12 = 20260828
SEED_ADJ9 = 20260829

assert DATA_FILE.name == "Personality_Pornography_Project_Data.csv"

if not DATA_FILE.exists():
    raise FileNotFoundError(
        f"Expected data file not found: {DATA_FILE.name}\n"
        f"Place the CSV in the same folder as this script and rerun."
    )

df = pd.read_csv(DATA_FILE)
# Demographic coding in the final analytic file:
# Age: 1=18-23, 2=24-29, 3=30-35, 4=36-40 in this analytic sample.
# The original survey coding was open-ended at the tails: respondents <18 also
# mapped to code 1 and respondents >40 to code 4, but this study used ages 18-40.
# Gender: 1=female (reference), 2=male.
# Education: 1=below diploma, 2=diploma, 3=bachelor's, 4=master's, 5=doctorate.
# Marital_Status: 1=single, 2=married, 3=divorced.
# Employment_Status: 1=unemployed, 2=full-time, 3=part-time, 4=homemaker, 5=student.


X = "PID5_Total"
M1 = "LPFS_Total"
M2 = "PCI_Emotional_Avoidance"
Y = "PPUS_Total_12"
Y9 = "PPUS9_Nonavoid"
YA = "PPUS_Avoidance_3"

required = [X, M1, M2, Y, Y9, YA]
required += [f"PID{i}" for i in range(1,26)]
required += [f"LPFS{i}" for i in range(1,13)]
required += [f"PCI{i}" for i in range(1,16)]
required += [f"PPU{i}" for i in range(1,13)]
required += ["Age","Gender","Marital_Status","Education","Employment_Status"]

missing = [c for c in required if c not in df.columns]
if missing:
    raise ValueError(f"Missing required columns: {missing}")

assert len(df) == 591, f"Expected N=591, found N={len(df)}"
assert df[[X,M1,M2,Y,Y9]].isna().sum().sum() == 0
assert np.allclose(df[X], df[[f"PID{i}" for i in range(1,26)]].sum(axis=1))
assert np.allclose(df[M1], df[[f"LPFS{i}" for i in range(1,13)]].sum(axis=1))
assert np.allclose(df[M2], df[["PCI2","PCI3","PCI10","PCI12","PCI15"]].sum(axis=1))
assert np.allclose(df[Y], df[[f"PPU{i}" for i in range(1,13)]].sum(axis=1))
assert np.allclose(df[YA], df[["PPU10","PPU11","PPU12"]].sum(axis=1))
assert np.allclose(df[Y9], df[[f"PPU{i}" for i in range(1,10)]].sum(axis=1))

def alpha(cols):
    z = df[cols].astype(float)
    k = z.shape[1]
    return k/(k-1) * (1 - z.var(ddof=1).sum()/z.sum(axis=1).var(ddof=1))

scale_defs = {
    "PID-5": (X, [f"PID{i}" for i in range(1,26)]),
    "LPFS": (M1, [f"LPFS{i}" for i in range(1,13)]),
    "PCI_EA": (M2, ["PCI2","PCI3","PCI10","PCI12","PCI15"]),
    "PPUS12": (Y, [f"PPU{i}" for i in range(1,13)]),
    "PPUS_Avoidance": (YA, ["PPU10","PPU11","PPU12"]),
    "PPUS9": (Y9, [f"PPU{i}" for i in range(1,10)]),
}

# ---------- Descriptives and correlations ----------
desc_rows = []
for label,(v,cols) in scale_defs.items():
    desc_rows.append({
        "Variable": label,
        "N": int(df[v].count()),
        "Mean": df[v].mean(),
        "SD": df[v].std(ddof=1),
        "Min": df[v].min(),
        "Max": df[v].max(),
        "Skewness": skew(df[v], bias=False),
        "Excess_Kurtosis": kurtosis(df[v], fisher=True, bias=False),
        "Cronbach_alpha": alpha(cols)
    })
desc = pd.DataFrame(desc_rows)

corr_vars = [X,M1,M2,Y,YA,Y9]
corr = df[corr_vars].corr()

# Reproducible Pearson r and exact two-sided p-values for every Table S2 pair.
corr_rows = []
for i, v1 in enumerate(corr_vars):
    for j, v2 in enumerate(corr_vars):
        if i == j:
            r, p = 1.0, np.nan
        else:
            r, p = pearsonr(df[v1], df[v2])
        corr_rows.append({"Variable1": v1, "Variable2": v2, "Pearson_r": r, "p_two_sided": p})
corr_long = pd.DataFrame(corr_rows)

def format_p_for_report(p):
    # Floating-point underflow can yield 0.0 for an extremely small probability.
    # For publication, such values are never interpreted as exact zero.
    if pd.isna(p):
        return "—"
    if p == 0.0 or p < .001:
        return "< .001"
    return f"{p:.3f}"

corr_long["p_report"] = corr_long["p_two_sided"].apply(format_p_for_report)
desc.to_csv(BASE/"Personality_Pornography_Project_TableS2_Descriptives.csv", index=False)
corr_long.to_csv(BASE/"Personality_Pornography_Project_TableS2.csv", index=False)

def fit_std(data, outcome, predictors, covariates=None):
    covariates = covariates or []
    z = data.copy()
    for c in set([outcome] + [p for p in predictors if p in [X,M1,M2,Y,Y9]]):
        z[c] = (z[c]-z[c].mean())/z[c].std(ddof=1)
    return sm.OLS(z[outcome], sm.add_constant(z[predictors+covariates])).fit(cov_type="HC3")

def fit_unstd(data, outcome, predictors, covariates=None):
    covariates = covariates or []
    return sm.OLS(data[outcome], sm.add_constant(data[predictors+covariates])).fit(cov_type="HC3")

# ---------- Main and PPUS-9 component paths ----------
models = []
for outcome_label, outcome in [("PPUS-12",Y),("PPUS-9",Y9)]:
    m1s = fit_std(df,M1,[X]); m1u = fit_unstd(df,M1,[X])
    m2s = fit_std(df,M2,[X,M1]); m2u = fit_unstd(df,M2,[X,M1])
    m3s = fit_std(df,outcome,[X,M1,M2]); m3u = fit_unstd(df,outcome,[X,M1,M2])
    for eq,ms,mu,preds in [
        ("Eq1_LPFS",m1s,m1u,[X]),
        ("Eq2_EmotionalAvoidance",m2s,m2u,[X,M1]),
        ("Eq3_Outcome",m3s,m3u,[X,M1,M2]),
    ]:
        for p in preds:
            ci = ms.conf_int().loc[p].tolist()
            models.append({
                "Outcome": outcome_label,
                "Equation": eq,
                "Predictor": p,
                "b": mu.params[p],
                "HC3_SE_b": mu.bse[p],
                "beta": ms.params[p],
                "HC3_SE_beta": ms.bse[p],
                "p": ms.pvalues[p],
                "CI95_beta_low": ci[0],
                "CI95_beta_high": ci[1],
                "R2": ms.rsquared
            })
pd.DataFrame(models).to_csv(BASE/"Personality_Pornography_Project_TableS3.csv", index=False)

# ---------- VIF ----------
vif_rows = []
for eq,preds in [
    ("Eq1_LPFS",[X]),
    ("Eq2_EmotionalAvoidance",[X,M1]),
    ("Eq3_Outcome",[X,M1,M2]),
]:
    Z = sm.add_constant(df[preds])
    for i,p in enumerate(Z.columns):
        if p == "const":
            continue
        vif = variance_inflation_factor(Z.values,i)
        vif_rows.append({
            "Equation":eq,
            "Predictor":p,
            "VIF":vif,
            "Tolerance":1/vif
        })
pd.DataFrame(vif_rows).to_csv(BASE/"Personality_Pornography_Project_TableS7.csv", index=False)

# ---------- Bootstrap helpers ----------
def zmat(sample, cols):
    a = sample[cols].to_numpy(dtype=float)
    return (a-a.mean(0))/a.std(0,ddof=1)

def olscoef(y,x):
    return np.linalg.lstsq(x,y,rcond=None)[0]

def effects(sample,y,order="hyp"):
    z = zmat(sample,[X,M1,M2,y])
    xx,m1v,m2v,yy = z.T
    one = np.ones(len(z))
    if order == "hyp":
        q1 = olscoef(m1v,np.c_[one,xx])
        q2 = olscoef(m2v,np.c_[one,xx,m1v])
        q3 = olscoef(yy,np.c_[one,xx,m1v,m2v])
        a1=q1[1]; a2=q2[1]; d=q2[2]; cp=q3[1]; b1=q3[2]; b2=q3[3]
        i1=a1*b1; i2=a2*b2; ser=a1*d*b2; ti=i1+i2+ser
        return np.array([i1,i2,ser,ti,cp+ti])
    if order == "parallel":
        q1 = olscoef(m1v,np.c_[one,xx])
        q2 = olscoef(m2v,np.c_[one,xx])
        q3 = olscoef(yy,np.c_[one,xx,m1v,m2v])
        return np.array([q1[1]*q3[2], q2[1]*q3[3]])
    if order == "reverse":
        q1 = olscoef(m2v,np.c_[one,xx])
        q2 = olscoef(m1v,np.c_[one,xx,m2v])
        q3 = olscoef(yy,np.c_[one,xx,m2v,m1v])
        return np.array([q1[1]*q2[2]*q3[3]])
    raise ValueError(order)

def bootstrap(sample,y,order,seed):
    rng = np.random.default_rng(seed)
    theta = effects(sample,y,order)
    n = len(sample)
    draws = np.empty((B,len(theta)))
    for i in range(B):
        draws[i] = effects(sample.iloc[rng.integers(0,n,n)],y,order)
    return theta,draws

def bca(sample,y,order,theta,draws):
    n = len(sample)
    jack = np.empty((n,len(theta)))
    for i in range(n):
        jack[i] = effects(sample.drop(sample.index[i]),y,order)
    jm = jack.mean(0)
    num = np.sum((jm-jack)**3,0)
    den = 6*np.power(np.sum((jm-jack)**2,0),1.5)
    acc = np.divide(num,den,out=np.zeros_like(num),where=den!=0)
    z0 = norm.ppf(np.clip((draws<theta).mean(0),1e-8,1-1e-8))
    bounds=[]
    for za in [norm.ppf(.025),norm.ppf(.975)]:
        q = norm.cdf(z0+(z0+za)/(1-acc*(z0+za)))
        bounds.append(np.array([np.quantile(draws[:,j],q[j]) for j in range(len(theta))]))
    return np.vstack(bounds)

boot_rows=[]
for label,y,seed in [("PPUS-12",Y,SEED_PPUS12),("PPUS-9",Y9,SEED_PPUS9)]:
    theta,draws = bootstrap(df,y,"hyp",seed)
    pct = np.quantile(draws,[.025,.975],axis=0)
    bc = bca(df,y,"hyp",theta,draws)
    names=[
        "Specific_indirect_PID5_LPFS_outcome",
        "Specific_indirect_PID5_EA_outcome",
        "Serial_indirect_PID5_LPFS_EA_outcome",
        "Total_indirect",
        "Total_association",
    ]
    for i,nm in enumerate(names):
        boot_rows.append({
            "Outcome":label,"Effect":nm,"beta":theta[i],
            "Percentile95_low":pct[0,i],"Percentile95_high":pct[1,i],
            "BCa95_low":bc[0,i],"BCa95_high":bc[1,i]
        })
pd.DataFrame(boot_rows).to_csv(BASE/"Personality_Pornography_Project_TableS4.csv", index=False)

# ---------- Alternative models ----------
alt_rows=[]
# Parallel
par_m1 = fit_std(df,M1,[X])
par_m2 = fit_std(df,M2,[X])
par_y = fit_std(df,Y,[X,M1,M2])
for eq,m,preds in [
    ("Parallel_LPFS_on_PID5",par_m1,[X]),
    ("Parallel_EA_on_PID5",par_m2,[X]),
    ("Parallel_PPUS12_outcome",par_y,[X,M1,M2]),
]:
    for p in preds:
        alt_rows.append({
            "Configuration":"Parallel","Equation":eq,"Predictor":p,
            "beta":m.params[p],"HC3_SE":m.bse[p],"p":m.pvalues[p],"R2":m.rsquared
        })
thpar,drpar = bootstrap(df,Y,"parallel",SEED_ALT)
pctpar = np.quantile(drpar,[.025,.975],axis=0)
for nm,i in [("Indirect_via_LPFS",0),("Indirect_via_EA",1)]:
    alt_rows.append({
        "Configuration":"Parallel","Equation":"Indirect_effect","Predictor":nm,
        "beta":thpar[i],"HC3_SE":np.nan,"p":np.nan,"R2":np.nan,
        "Bootstrap95_low":pctpar[0,i],"Bootstrap95_high":pctpar[1,i]
    })

# Reversed
rev_m2 = fit_std(df,M2,[X])
rev_m1 = fit_std(df,M1,[X,M2])
rev_y = fit_std(df,Y,[X,M2,M1])
for eq,m,preds in [
    ("Reverse_EA_on_PID5",rev_m2,[X]),
    ("Reverse_LPFS_on_PID5_EA",rev_m1,[X,M2]),
    ("Reverse_PPUS12_outcome",rev_y,[X,M2,M1]),
]:
    for p in preds:
        alt_rows.append({
            "Configuration":"Reversed","Equation":eq,"Predictor":p,
            "beta":m.params[p],"HC3_SE":m.bse[p],"p":m.pvalues[p],"R2":m.rsquared
        })
threv,drrev = bootstrap(df,Y,"reverse",SEED_ALT)
pctrev = np.quantile(drrev,[.025,.975],axis=0)
alt_rows.append({
    "Configuration":"Reversed","Equation":"Indirect_effect",
    "Predictor":"Serial_PID5_EA_LPFS_PPUS12",
    "beta":threv[0],"HC3_SE":np.nan,"p":np.nan,"R2":np.nan,
    "Bootstrap95_low":pctrev[0,0],"Bootstrap95_high":pctrev[1,0]
})
pd.DataFrame(alt_rows).to_csv(BASE/"Personality_Pornography_Project_TableS5.csv", index=False)

# ---------- Demographic-adjusted ----------
demo = ["Age","Gender","Marital_Status","Education","Employment_Status"]
D = pd.get_dummies(df[demo].astype("category"),drop_first=True,dtype=float)
adj = pd.concat([df,D],axis=1)
covs = list(D.columns)

adj_rows=[]
for outcome_label,outcome in [("PPUS-12",Y),("PPUS-9",Y9)]:
    am1 = fit_std(adj,M1,[X],covs)
    am2 = fit_std(adj,M2,[X,M1],covs)
    amy = fit_std(adj,outcome,[X,M1,M2],covs)
    for eq,m,preds in [
        ("Adjusted_LPFS",am1,[X]+covs),
        ("Adjusted_EA",am2,[X,M1]+covs),
        ("Adjusted_Outcome",amy,[X,M1,M2]+covs),
    ]:
        for p in preds:
            adj_rows.append({
                "Outcome":outcome_label,"Equation":eq,"Predictor":p,
                "beta_scale_coefficient":m.params[p],
                "HC3_SE":m.bse[p],"p":m.pvalues[p],"R2":m.rsquared
            })

levels={c:np.sort(df[c].dropna().unique()) for c in demo}
def dummy_matrix(sample):
    mats=[]
    for c in demo:
        vals=sample[c].to_numpy()
        for lev in levels[c][1:]:
            mats.append((vals==lev).astype(float))
    return np.column_stack(mats)

def adjusted_serial(sample,y):
    z = zmat(sample,[X,M1,M2,y])
    xx,m1v,m2v,yy = z.T
    one=np.ones(len(z)); C=dummy_matrix(sample)
    q1=olscoef(m1v,np.c_[one,xx,C])
    q2=olscoef(m2v,np.c_[one,xx,m1v,C])
    q3=olscoef(yy,np.c_[one,xx,m1v,m2v,C])
    return q1[1]*q2[2]*q3[3]

def boot_adjusted(y,seed):
    rng=np.random.default_rng(seed)
    n=len(df)
    theta=adjusted_serial(df,y)
    d=np.empty(B)
    for i in range(B):
        d[i]=adjusted_serial(df.iloc[rng.integers(0,n,n)],y)
    return theta,np.quantile(d,[.025,.975])

for label,y,seed in [("PPUS-12",Y,SEED_ADJ12),("PPUS-9",Y9,SEED_ADJ9)]:
    th,ci=boot_adjusted(y,seed)
    adj_rows.append({
        "Outcome":label,"Equation":"Adjusted_serial_indirect",
        "Predictor":"PID5_LPFS_EA_outcome",
        "beta_scale_coefficient":th,
        "HC3_SE":np.nan,"p":np.nan,"R2":np.nan,
        "Bootstrap95_low":ci[0],"Bootstrap95_high":ci[1]
    })

adj_df = pd.DataFrame(adj_rows)
adj_df.to_csv(BASE/"Personality_Pornography_Project_TableS6.csv", index=False)

# ---------- Human-readable log ----------
with open(BASE/"Personality_Pornography_Project_Output.txt","w",encoding="utf-8") as f:
    def log(*args):
        print(*args,file=f)
    log("REPRODUCIBILITY RUN COMPLETED SUCCESSFULLY")
    log("="*68)
    log("Data file:",DATA_FILE.name)
    log("N =",len(df))
    log("Python:",sys.version.replace("\n"," "))
    log("Platform:",platform.platform())
    log("numpy:",np.__version__)
    log("pandas:",pd.__version__)
    log("scipy:",scipy.__version__)
    log("statsmodels:",statsmodels.__version__)
    log("Bootstrap resamples:",B)
    log("CI methods: BCa + percentile for primary PPUS-12 and overlap-reduced PPUS-9 serial models; percentile only for parallel, reverse-order, and sociodemographically adjusted sensitivity models.")
    log("Seeds:",SEED_PPUS12,SEED_PPUS9,SEED_ALT,SEED_ADJ12,SEED_ADJ9)
    log()
    log("SCALE RECONSTRUCTION: PASS")
    log()
    log("DESCRIPTIVES / RELIABILITY")
    log(desc.to_string(index=False))
    log()
    log("CORRELATION MATRIX")
    log(corr.to_string())
    log()
    log("CORRELATION P-VALUES")
    log("Numerical p=0.0 from floating-point underflow is reported as p < .001, never as exact zero.")
    log(corr_long.to_string(index=False))
    log()
    log("CORRELATION P-VALUES (publication-safe representation)")
    log("Note: any floating-point p=0.0 caused by numerical underflow is reported as p < .001, never as exact zero.")
    log(corr_long.to_string(index=False))
    log()
    log("PRIMARY / PPUS-9 PATH COEFFICIENTS")
    log(pd.DataFrame(models).to_string(index=False))
    log()
    log("VIF / TOLERANCE")
    log(pd.DataFrame(vif_rows).to_string(index=False))
    log()
    log("BOOTSTRAP EFFECTS")
    log(pd.DataFrame(boot_rows).to_string(index=False))
    log()
    log("ALTERNATIVE MODELS")
    log(pd.DataFrame(alt_rows).to_string(index=False))
    log()
    log("DEMOGRAPHIC-ADJUSTED MODELS")
    log(adj_df.to_string(index=False))
    log()
    log("END OF REPRODUCIBILITY OUTPUT")

print("Reproducibility run completed successfully.")
print("Input:", DATA_FILE.name)
print("N =", len(df))
print("Outputs written to:", BASE)