/* phir.h — Phi_r (PhiID causal emergence) estimator, shared by the probe programs.
*
* Phi_r = I(X_t ; X_t+1) - sum_i I(X_t^i ; X_t+1)
* i.e. predictive information the WHOLE holds about its own future, beyond what its parts hold
* about that future. Following Pigozzi & Levin (arXiv 2607.28250) Materials and Methods, the
* intractable multivariate form is reduced via the MINIMUM-INFORMATION BIPARTITION to
* Phi_r = I(X_t ; X_t+1) - I(A_t ; X_t+1) - I(B_t ; X_t+1)
* with (A,B) minimising the normalised cross-cut lagged MI [I(A_t;B_t+1)+I(B_t;A_t+1)]/min(|A|,|B|).
*
* MI is estimated Gaussian (log-determinant form) with Ledoit-Wolf shrinkage toward the diagonal,
* so the regularisation is data-driven rather than a tuned knob.
*
* KNOWN LIMITATION, measured: this estimator is far more sensitive to redundancy than to synergy
* (|z| ~ 165 vs ~3.4 on matched synthetic controls). Carry that caveat into any result.
*
* The permutation null uses its own RNG so that null draws never perturb a caller's model stream.
*/
#ifndef PHIR_H
#define PHIR_H
#include <math.h>
#include <stdlib.h>
#include <string.h>
#define PHIR_MAXP 24
static unsigned long long PHIR_RNG = 0x243f6a8885a308d3ULL;
static unsigned long long phir_xs(void){ PHIR_RNG^=PHIR_RNG<<13; PHIR_RNG^=PHIR_RNG>>7; PHIR_RNG^=PHIR_RNG<<17; return PHIR_RNG; }
static void phir_seed(unsigned long long s){ PHIR_RNG = s? s : 0x243f6a8885a308d3ULL; }
/* ---------- generic stats ---------- */
static int phir_cmpd(const void*a,const void*b){ double x=*(const double*)a,y=*(const double*)b; return (x>y)-(x<y); }
static double meanv(const double *v,int n){ double s=0; for(int i=0;i<n;i++) s+=v[i]; return n? s/n : NAN; }
static double sdv(const double *v,int n){ double m=meanv(v,n),s=0; for(int i=0;i<n;i++) s+=(v[i]-m)*(v[i]-m); return sqrt(s/(n>1?n-1:1)); }
static double median(double *v,int n){ double *c=malloc((size_t)n*sizeof(double)); memcpy(c,v,(size_t)n*sizeof(double)); qsort(c,n,sizeof(double),phir_cmpd); double m=(n%2)?c[n/2]:0.5*(c[n/2-1]+c[n/2]); free(c); return m; }
/* Mann-Whitney U -> two-sided normal-approximation p (no ties correction) */
static double mannwhitney_p(const double *x,int nx,const double *y,int ny)
{
double U=0;
for(int i=0;i<nx;i++) for(int j=0;j<ny;j++){ if(x[i]>y[j]) U+=1.0; else if(x[i]==y[j]) U+=0.5; }
double mu=0.5*nx*ny, sd=sqrt((double)nx*ny*(nx+ny+1)/12.0);
if(sd<=0) return 1.0;
return erfc(fabs(U-mu)/sd/sqrt(2.0));
}
/* Spearman rank correlation (average ranks for ties) */
static void phir_ranks(const double *v,int n,double *r)
{
int *idx=malloc((size_t)n*sizeof(int));
for(int i=0;i<n;i++) idx[i]=i;
for(int i=1;i<n;i++){ int k=idx[i],j=i-1; while(j>=0 && v[idx[j]]>v[k]){ idx[j+1]=idx[j]; j--; } idx[j+1]=k; }
int i=0;
while(i<n){ int j=i; while(j+1<n && v[idx[j+1]]==v[idx[i]]) j++;
double avg=0.5*((double)i+(double)j)+1.0;
for(int k=i;k<=j;k++) r[idx[k]]=avg;
i=j+1; }
free(idx);
}
static double spearman(const double *x,const double *y,int n)
{
double *rx=malloc((size_t)n*sizeof(double)), *ry=malloc((size_t)n*sizeof(double));
phir_ranks(x,n,rx); phir_ranks(y,n,ry);
double mx=meanv(rx,n),my=meanv(ry,n),sxy=0,sxx=0,syy=0;
for(int i=0;i<n;i++){ double a=rx[i]-mx,b=ry[i]-my; sxy+=a*b; sxx+=a*a; syy+=b*b; }
free(rx);free(ry);
return (sxx>0&&syy>0)? sxy/sqrt(sxx*syy) : NAN;
}
static double spearman_p(double rho,int n) /* two-sided, t approximation */
{
if(n<4||isnan(rho)) return NAN;
double r=rho; if(r>0.999999)r=0.999999; if(r<-0.999999)r=-0.999999;
double t=fabs(r)*sqrt((n-2)/(1-r*r));
return erfc(t/sqrt(2.0)); /* normal approx; conservative enough for n>=20 */
}
/* ---------- compositional preprocessing: CLR with pseudocount, drop last component ---------- */
static void phir_clr(const double *cnt,int nb,double *out)
{
double tot=0; for(int i=0;i<nb;i++) tot+=cnt[i]+0.5;
double *lg=malloc((size_t)nb*sizeof(double)); double mean=0;
for(int i=0;i<nb;i++){ lg[i]=log((cnt[i]+0.5)/tot); mean+=lg[i]; }
mean/=nb;
for(int i=0;i<nb-1;i++) out[i]=lg[i]-mean;
free(lg);
}
/* ---------- Gaussian MI via log-determinants, Ledoit-Wolf shrinkage ---------- */
static void phir_cov_lw(const double *Z,int n,int d,double *S)
{
double *mu=calloc((size_t)d,sizeof(double));
for(int t=0;t<n;t++) for(int j=0;j<d;j++) mu[j]+=Z[(size_t)t*d+j];
for(int j=0;j<d;j++) mu[j]/=n;
double *Y=malloc((size_t)n*d*sizeof(double));
for(int t=0;t<n;t++) for(int j=0;j<d;j++) Y[(size_t)t*d+j]=Z[(size_t)t*d+j]-mu[j];
for(int i=0;i<d;i++) for(int j=0;j<d;j++){
double s=0; for(int t=0;t<n;t++) s+=Y[(size_t)t*d+i]*Y[(size_t)t*d+j];
S[i*d+j]=s/n;
}
double num=0, den=0;
for(int i=0;i<d;i++) for(int j=0;j<d;j++){
if(i==j) continue;
double sij=S[i*d+j], v=0;
for(int t=0;t<n;t++){ double e=Y[(size_t)t*d+i]*Y[(size_t)t*d+j]-sij; v+=e*e; }
v *= (double)n/((double)(n-1)*(n-1)*(n-1));
num+=v; den+=sij*sij;
}
double lam = (den>0)? num/den : 1.0;
if(lam<0) lam=0;
if(lam>1) lam=1;
for(int i=0;i<d;i++) for(int j=0;j<d;j++) if(i!=j) S[i*d+j]*=(1.0-lam);
double sc=(double)n/(double)(n-1);
for(int i=0;i<d;i++) for(int j=0;j<d;j++) S[i*d+j]*=sc;
free(mu); free(Y);
}
static double phir_logdet(const double *S,int d,const int *idx,int k,int *ok)
{
double *M=malloc((size_t)k*k*sizeof(double));
for(int a=0;a<k;a++) for(int b=0;b<k;b++) M[a*k+b]=S[idx[a]*d+idx[b]];
double tr=0; for(int a=0;a<k;a++) tr+=M[a*k+a];
double jitter=1e-10*(tr/k+1e-12);
for(int attempt=0;attempt<8;attempt++){
double *C=malloc((size_t)k*k*sizeof(double));
memcpy(C,M,(size_t)k*k*sizeof(double));
for(int a=0;a<k;a++) C[a*k+a]+=jitter;
int bad=0; double ld=0;
for(int a=0;a<k && !bad;a++){
for(int b=0;b<=a;b++){
double s=C[a*k+b];
for(int m=0;m<b;m++) s-=C[a*k+m]*C[b*k+m];
if(a==b){ if(s<=0){ bad=1; break; } C[a*k+a]=sqrt(s); ld+=2.0*log(C[a*k+a]); }
else C[a*k+b]=s/C[b*k+b];
}
}
free(C);
if(!bad){ free(M); *ok=1; return ld; }
jitter = (jitter>0)? jitter*100.0 : 1e-9*(tr/k+1e-12);
}
free(M); *ok=0; return 0.0;
}
static double phir_mi(const double *S,int d,const int *A,int ka,const int *B,int kb,int *ok)
{
int *U=malloc((size_t)(ka+kb)*sizeof(int));
memcpy(U,A,(size_t)ka*sizeof(int)); memcpy(U+ka,B,(size_t)kb*sizeof(int));
int o1,o2,o3;
double la=phir_logdet(S,d,A,ka,&o1), lb=phir_logdet(S,d,B,kb,&o2), lu=phir_logdet(S,d,U,ka+kb,&o3);
free(U);
*ok=o1&&o2&&o3;
return 0.5*(la+lb-lu);
}
/* pairs: n x (2p) row-major, [past(0..p-1) | future(p..2p-1)] */
static double phi_r(const double *pairs,int n,int p,int *ok)
{
int d=2*p;
double *S=malloc((size_t)d*d*sizeof(double));
phir_cov_lw(pairs,n,d,S);
int past[PHIR_MAXP], fut[PHIR_MAXP];
for(int i=0;i<p;i++){ past[i]=i; fut[i]=p+i; }
int o;
double I_whole=phir_mi(S,d,past,p,fut,p,&o);
if(!o){ free(S); *ok=0; return 0.0; }
int bestA[PHIR_MAXP], bestB[PHIR_MAXP], bka=0, bkb=0, found=0;
double bestcut=1e300;
int Ap[PHIR_MAXP],Bp[PHIR_MAXP],Af[PHIR_MAXP],Bf[PHIR_MAXP];
for(unsigned m=1; m<(1u<<(p-1)); m++){ /* part 0 pinned to A -> unordered bipartitions */
int ka=0,kb=0;
for(int i=0;i<p;i++){
int inA=(i==0)?1:((m>>(i-1))&1);
if(inA){ Ap[ka]=i; Af[ka]=p+i; ka++; } else { Bp[kb]=i; Bf[kb]=p+i; kb++; }
}
if(!ka||!kb) continue;
int oa,ob;
double c1=phir_mi(S,d,Ap,ka,Bf,kb,&oa);
double c2=phir_mi(S,d,Bp,kb,Af,ka,&ob);
if(!oa||!ob) continue;
int mn=(ka<kb)?ka:kb;
double cut=(c1+c2)/(double)mn;
if(cut<bestcut){ bestcut=cut; bka=ka; bkb=kb;
memcpy(bestA,Ap,(size_t)ka*sizeof(int)); memcpy(bestB,Bp,(size_t)kb*sizeof(int)); found=1; }
}
if(!found){ free(S); *ok=0; return 0.0; }
int o1,o2;
double I_A=phir_mi(S,d,bestA,bka,fut,p,&o1);
double I_B=phir_mi(S,d,bestB,bkb,fut,p,&o2);
free(S);
*ok=o1&&o2;
return I_whole-I_A-I_B;
}
/* Permutation null: shuffle the future block across rows to destroy past->future pairing. */
static void phi_r_null(const double *pairs,int n,int p,int R,double *nmean,double *nsd)
{
double *buf=malloc((size_t)n*2*p*sizeof(double));
double *vals=malloc((size_t)R*sizeof(double));
int *perm=malloc((size_t)n*sizeof(int));
int nv=0;
for(int r=0;r<R;r++){
memcpy(buf,pairs,(size_t)n*2*p*sizeof(double));
for(int i=0;i<n;i++) perm[i]=i;
for(int i=n-1;i>0;i--){ int j=(int)(phir_xs()%(unsigned long long)(i+1)); int t=perm[i]; perm[i]=perm[j]; perm[j]=t; }
for(int i=0;i<n;i++) for(int k=0;k<p;k++) buf[(size_t)i*2*p+p+k]=pairs[(size_t)perm[i]*2*p+p+k];
int ok; double v=phi_r(buf,n,p,&ok);
if(ok) vals[nv++]=v;
}
*nmean = nv? meanv(vals,nv):NAN;
*nsd = nv>1? sdv(vals,nv):NAN;
free(buf);free(vals);free(perm);
}
#endif