"""RehabTools paired measurement comparison 1.0.0. Runs locally in a worker.
ICC equations: McGraw/Wong, as tabulated by Koo & Li (2016), doi:10.1016/j.jcm.2016.02.012.
Plots distinguish correlation, absolute agreement and paired differences.
"""
import json, math
import numpy as np
from scipy import stats

def reliability(data):
    y=np.asarray(data.get('rows'),dtype=float)
    if y.ndim!=2 or y.shape[1]!=2 or not 3<=len(y)<=5000 or not np.isfinite(y).all():
        raise ValueError('Provide 3–5,000 complete numeric pairs.')
    a,b=y[:,0],y[:,1];n,k=y.shape;grand=float(y.mean());diff=b-a;means=y.mean(axis=1)
    msr=float(k*np.sum((means-grand)**2)/(n-1))
    msc=float(n*np.sum((y.mean(axis=0)-grand)**2)/(k-1))
    mse=float(np.sum((y-means[:,None]-y.mean(axis=0)[None,:]+grand)**2)/((n-1)*(k-1)))
    def divide(a,b):return float(a/b) if abs(b)>1e-14 else None
    icc_a1=divide(msr-mse,msr+(k-1)*mse+k*(msc-mse)/n)
    icc_ak=divide(msr-mse,msr+(msc-mse)/n)
    icc_c1=divide(msr-mse,msr+(k-1)*mse)
    icc_ck=divide(msr-mse,msr)
    sd_a=float(a.std(ddof=1));sd_b=float(b.std(ddof=1));sd_diff=float(diff.std(ddof=1));bias=float(diff.mean())
    pooled=math.sqrt((sd_a**2+sd_b**2)/2)
    sem=pooled*math.sqrt(1-icc_a1) if icc_a1 is not None and icc_a1<=1 else None
    mdc=sem*1.96*math.sqrt(2) if sem is not None else None
    normal=lambda result:{'statistic':float(result.statistic),'p':float(result.pvalue)}
    tests={}
    tests['paired_t']=normal(stats.ttest_rel(b,a)) if sd_diff>0 else {'statistic':None,'p':1.0 if bias==0 else None}
    tests['wilcoxon']=normal(stats.wilcoxon(diff,zero_method='wilcox',method='auto')) if np.any(diff) else {'statistic':0.0,'p':1.0}
    tests['shapiro_a']=normal(stats.shapiro(a)) if sd_a>0 else {'statistic':None,'p':None}
    tests['shapiro_b']=normal(stats.shapiro(b)) if sd_b>0 else {'statistic':None,'p':None}
    tests['shapiro_diff']=normal(stats.shapiro(diff)) if sd_diff>0 else {'statistic':None,'p':None}
    tests['levene']=normal(stats.levene(a,b,center='median')) if sd_a>0 or sd_b>0 else {'statistic':None,'p':None}
    pearson=normal(stats.pearsonr(a,b)) if sd_a>0 and sd_b>0 else {'statistic':None,'p':None}
    def fitted(x,y):
        if np.ptp(x)==0:return None
        fit=stats.linregress(x,y);grid=np.linspace(float(min(x)),float(max(x)),80)
        residual=y-(fit.intercept+fit.slope*x);s=float(np.sqrt(np.sum(residual**2)/(n-2)))
        error=s*np.sqrt(1/n+(grid-x.mean())**2/np.sum((x-x.mean())**2));critical=float(stats.t.ppf(.975,n-2))
        pred=fit.intercept+fit.slope*grid
        return {'slope':float(fit.slope),'intercept':float(fit.intercept),'p':float(fit.pvalue),'r_squared':float(fit.rvalue**2),'residual_sd':s,'line':np.column_stack([grid,pred,pred-critical*error,pred+critical*error]).tolist()}
    def histogram(x):
        counts,edges=np.histogram(x,bins=min(20,max(4,int(math.sqrt(n)))))
        grid=np.linspace(float(edges[0]),float(edges[-1]),80);sd=float(x.std(ddof=1))
        curve=stats.norm.pdf(grid,loc=x.mean(),scale=sd)*n*(edges[1]-edges[0]) if sd>0 else np.zeros(80)
        return {'bins':np.column_stack([edges[:-1],edges[1:],counts]).tolist(),'normal':np.column_stack([grid,curve]).tolist()}
    def qq(x):
        theoretical,ordered=stats.probplot(x,dist='norm',fit=False)
        return np.column_stack([theoretical,ordered]).tolist()
    out={'version':'1.0.0','n':n,'mean_a':float(a.mean()),'mean_b':float(b.mean()),'sd_a':sd_a,'sd_b':sd_b,'mean_difference_b_minus_a':bias,'sd_difference':sd_diff,'limits_of_agreement':[bias-1.96*sd_diff,bias+1.96*sd_diff],'icc':{'absolute_single_2_1':icc_a1,'absolute_average_2_2':icc_ak,'consistency_single_c1':icc_c1,'consistency_average_c2':icc_ck},'anova':{'MSR':msr,'MSC':msc,'MSE':mse},'sem_sd_pooled_icc_a1':sem,'mdc95':mdc,'sem_percent':100*sem/grand if sem is not None and grand>0 else None,'mdc95_percent':100*mdc/grand if mdc is not None and grand>0 else None,'paired_cohens_dz':divide(bias,sd_diff),'cohens_d_pooled':divide(bias,pooled),'pearson':pearson,'tests':tests,'scatter':np.column_stack([a,b]).tolist(),'scatter_fit':fitted(a,b),'bland_altman':np.column_stack([means,diff]).tolist(),'bland_altman_fit':fitted(means,diff),'qq_a':qq(a),'qq_b':qq(b),'histogram_a':histogram(a),'histogram_b':histogram(b)}
    def clean(v):
        if isinstance(v,dict):return {k:clean(x) for k,x in v.items()}
        if isinstance(v,list):return [clean(x) for x in v]
        if isinstance(v,float) and not math.isfinite(v):return None
        return v
    return clean(out)