Upload agents/statistical_agent.py with huggingface_hub

agents/statistical_agent.py

@@ -1,362 +1,308 @@
-"""
-FORENSIQ — Statistical Priors Agent
-Tests natural image statistics violations:
-  - DCT coefficient distribution (Laplacian vs Gaussian)
-  - Benford's law on first digits of DCT coefficients
-  - Gradient sparsity (kurtosis > 3 for natural images)
-"""
-
+"""FORENSIQ — Statistical Priors Agent (22 features)"""
 import numpy as np
 from PIL import Image
 from scipy.fftpack import dct
-from scipy.stats import kurtosis as scipy_kurtosis, entropy
-from scipy.ndimage import gaussian_filter
+from scipy.stats import kurtosis as sp_kurt, skew as sp_skew, entropy as sp_entropy
+from scipy.ndimage import gaussian_filter, sobel, uniform_filter
 from typing import Dict, Any
-
 from agents.optical_agent import AgentEvidence
 
-
-# ─── DCT Coefficient Distribution ───────────────────────────────────
-def analyze_dct_distribution(img: Image.Image) -> Dict[str, Any]:
-    """
-    Natural image DCT coefficients follow a Laplacian (heavy-tailed)
-    distribution. AI-generated images often follow a Gaussian.
-    """
-    gray = np.array(img.convert("L")).astype(np.float64)
-    h, w = gray.shape
-    h_crop, w_crop = (h // 8) * 8, (w // 8) * 8
-    gray = gray[:h_crop, :w_crop]
-
-    coeffs = []
-    for i in range(0, h_crop, 8):
-        for j in range(0, w_crop, 8):
-            block = gray[i:i + 8, j:j + 8]
-            dct_block = dct(dct(block.T, norm="ortho").T, norm="ortho")
-            # Skip DC coefficient
-            ac = dct_block.copy()
-            ac[0, 0] = 0
-            coeffs.extend(ac.flatten().tolist())
-
-    coeffs = np.array(coeffs)
-    coeffs = coeffs[coeffs != 0]
-
-    if len(coeffs) < 100:
-        return {"test": "DCT Distribution", "score": 0.0, "note": "Insufficient data"}
-
-    # Kurtosis: Laplacian ≈ 6, Gaussian ≈ 3
-    kurt = float(scipy_kurtosis(coeffs, fisher=True))
-
-    if kurt > 4.5:
-        score = -0.4
-        note = f"DCT kurtosis={kurt:.2f} (Laplacian-like, consistent with natural images)"
-    elif kurt < 2.0:
-        score = 0.5
-        note = f"DCT kurtosis={kurt:.2f} (Gaussian-like, inconsistent with natural images)"
-    elif kurt < 3.5:
-        score = 0.2
-        note = f"DCT kurtosis={kurt:.2f} (borderline, mildly Gaussian)"
-    else:
-        score = -0.1
-        note = f"DCT kurtosis={kurt:.2f} (near-natural)"
-
-    return {
-        "test": "DCT Distribution",
-        "kurtosis": round(kurt, 4),
-        "mean": round(float(np.mean(coeffs)), 4),
-        "std": round(float(np.std(coeffs)), 4),
-        "score": score,
-        "note": note,
-    }
-
-
-# ─── Benford's Law ──────────────────────────────────────────────────
-def analyze_benford(img: Image.Image) -> Dict[str, Any]:
-    """
-    First-digit distribution of DCT coefficients should follow
-    Benford's Law in natural images. AI images deviate.
-    """
-    gray = np.array(img.convert("L")).astype(np.float64)
-    h, w = gray.shape
-    h_crop, w_crop = (h // 8) * 8, (w // 8) * 8
-    gray = gray[:h_crop, :w_crop]
-
-    coeffs = []
-    for i in range(0, h_crop, 8):
-        for j in range(0, w_crop, 8):
-            block = gray[i:i + 8, j:j + 8]
-            dct_block = dct(dct(block.T, norm="ortho").T, norm="ortho")
-            coeffs.extend(np.abs(dct_block.flatten()).tolist())
-
-    coeffs = np.array(coeffs)
-    nonzero = coeffs[coeffs > 0]
-
-    if len(nonzero) < 100:
-        return {"test": "Benford's Law", "score": 0.0, "note": "Insufficient data"}
-
-    # Extract first digits
-    log_vals = np.floor(np.log10(nonzero + 1e-12))
-    first_digits = np.floor(nonzero / (10 ** log_vals)).astype(int)
-    first_digits = first_digits[(first_digits >= 1) & (first_digits <= 9)]
-
-    observed = np.array([np.sum(first_digits == d) for d in range(1, 10)], dtype=np.float64)
-    observed = observed / (observed.sum() + 1e-9)
-
-    # Benford's expected distribution
-    benford = np.log10(1 + 1.0 / np.arange(1, 10))
-
-    # Chi-squared statistic
-    chi2 = float(np.sum((observed - benford) ** 2 / (benford + 1e-9)))
-
-    # KL divergence
-    kl_div = float(np.sum(observed * np.log((observed + 1e-9) / (benford + 1e-9))))
-
-    if chi2 < 0.005:
-        score = -0.4
-        note = f"Excellent Benford's law fit (χ²={chi2:.5f}, natural image)"
-    elif chi2 < 0.02:
-        score = -0.1
-        note = f"Good Benford's law fit (χ²={chi2:.5f})"
-    elif chi2 < 0.05:
-        score = 0.3
-        note = f"Moderate Benford's deviation (χ²={chi2:.5f})"
-    else:
-        score = 0.6
-        note = f"Strong Benford's law violation (χ²={chi2:.5f}, AI-like)"
-
-    return {
-        "test": "Benford's Law",
-        "chi_squared": round(chi2, 6),
-        "kl_divergence": round(kl_div, 6),
-        "observed": observed.tolist(),
-        "benford_expected": benford.tolist(),
-        "score": score,
-        "note": note,
-    }
-
-
-# ─── Gradient Sparsity ──────────────────────────────────────────────
-def analyze_gradient_sparsity(img: Image.Image) -> Dict[str, Any]:
-    """
-    Natural images have sparse gradients (kurtosis > 3).
-    AI images often have smoother gradients with lower kurtosis.
-    """
-    gray = np.array(img.convert("L")).astype(np.float64)
-
-    # Compute gradients
-    gx = np.diff(gray, axis=1)
-    gy = np.diff(gray, axis=0)
-
-    # Combine
-    gx_flat = gx.ravel()
-    gy_flat = gy.ravel()
-    all_grads = np.concatenate([gx_flat, gy_flat])
-
-    kurt_val = float(scipy_kurtosis(all_grads, fisher=True))
-    
-    # Sparsity: fraction of near-zero gradients
-    threshold = np.std(all_grads) * 0.1
-    sparsity = float(np.mean(np.abs(all_grads) < threshold))
-
-    if kurt_val > 5.0 and sparsity > 0.4:
-        score = -0.4
-        note = f"Sparse gradients (kurtosis={kurt_val:.2f}, sparsity={sparsity:.2f}, natural)"
-    elif kurt_val < 2.0:
-        score = 0.5
-        note = f"Low gradient kurtosis ({kurt_val:.2f}), unnaturally smooth"
-    elif kurt_val < 3.5:
-        score = 0.2
-        note = f"Borderline gradient statistics (kurtosis={kurt_val:.2f})"
-    else:
-        score = -0.1
-        note = f"Normal gradient statistics (kurtosis={kurt_val:.2f})"
-
-    return {
-        "test": "Gradient Sparsity",
-        "kurtosis": round(kurt_val, 4),
-        "sparsity": round(sparsity, 4),
-        "gradient_mean": round(float(np.mean(np.abs(all_grads))), 4),
-        "score": score,
-        "note": note,
-    }
-
-
-# ─── Local Kurtosis Map ──────────────────────────────────────────────
-def analyze_local_kurtosis(img: Image.Image) -> Dict[str, Any]:
-    """
-    Natural images have spatially varying kurtosis (textured vs smooth).
-    AI images often have unnaturally uniform local statistics.
-    """
-    gray = np.array(img.convert("L")).astype(np.float64)
-    h, w = gray.shape
-    block_size = 32
-    h_crop, w_crop = (h // block_size) * block_size, (w // block_size) * block_size
-    gray = gray[:h_crop, :w_crop]
-
-    local_kurts = []
-    for i in range(0, h_crop, block_size):
-        for j in range(0, w_crop, block_size):
-            block = gray[i:i + block_size, j:j + block_size].ravel()
-            if np.std(block) > 1:
-                local_kurts.append(float(scipy_kurtosis(block, fisher=True)))
-
-    if len(local_kurts) < 10:
-        return {"test": "Local Kurtosis Map", "score": 0.0, "note": "Insufficient blocks"}
-
-    local_kurts = np.array(local_kurts)
-    kurt_std = float(np.std(local_kurts))
-    kurt_mean = float(np.mean(local_kurts))
-
-    # Natural images: high variation in local kurtosis
-    if kurt_std > 3.0:
-        score = -0.3
-        note = f"High local kurtosis variation (σ={kurt_std:.2f}, natural spatial statistics)"
-    elif kurt_std < 1.0:
-        score = 0.4
-        note = f"Unnaturally uniform local statistics (σ={kurt_std:.2f}, AI-like)"
-    else:
-        score = 0.0
-        note = f"Moderate local kurtosis variation (σ={kurt_std:.2f})"
-
-    return {
-        "test": "Local Kurtosis Map",
-        "kurtosis_std": round(kurt_std, 4),
-        "kurtosis_mean": round(kurt_mean, 4),
-        "score": score,
-        "note": note,
-    }
-
-
-# ─── Color Histogram Analysis ───────────────────────────────────────
-def analyze_color_histogram(img: Image.Image) -> Dict[str, Any]:
-    """
-    Natural images have smooth, continuous color histograms.
-    AI/GAN images can show comb-like gaps or unusual peaks in histograms.
-    """
-    rgb = np.array(img.convert("RGB"))
-
-    anomaly_scores = []
-    for c, name in enumerate(["Red", "Green", "Blue"]):
-        hist, _ = np.histogram(rgb[:, :, c].ravel(), bins=256, range=(0, 256))
-        hist = hist.astype(np.float64)
-
-        # Check for gaps (zero bins surrounded by non-zero)
-        zero_bins = np.sum(hist == 0)
-
-        # Check for comb pattern (alternating zero/nonzero)
-        diffs = np.diff((hist > 0).astype(int))
-        transitions = int(np.sum(np.abs(diffs)))
-
-        # Smoothness: ratio of histogram derivative to histogram
-        hist_smooth = gaussian_filter(hist.astype(np.float64), sigma=2)
-        smoothness = float(np.mean(np.abs(hist - hist_smooth)) / (np.mean(hist) + 1e-9))
-
-        anomaly_scores.append(smoothness)
-
-    avg_smoothness = float(np.mean(anomaly_scores))
-
-    if avg_smoothness < 0.3:
-        score = -0.2
-        note = f"Smooth color histograms (smoothness={avg_smoothness:.3f}, natural)"
-    elif avg_smoothness > 0.8:
-        score = 0.4
-        note = f"Irregular color histograms (smoothness={avg_smoothness:.3f}, manipulation artifact)"
-    else:
-        score = 0.0
-        note = f"Normal histogram smoothness ({avg_smoothness:.3f})"
-
-    return {
-        "test": "Color Histogram Analysis",
-        "histogram_smoothness": round(avg_smoothness, 4),
-        "score": score,
-        "note": note,
-    }
-
-
-# ─── Wavelet Coefficient Distribution ───────────────────────────────
-def analyze_wavelet_coefficients(img: Image.Image) -> Dict[str, Any]:
-    """
-    Natural image wavelet coefficients follow generalized Gaussian.
-    AI images show deviations, especially at high-frequency subbands.
-    Uses Haar wavelet (simple, no pywt dependency needed).
-    """
-    gray = np.array(img.convert("L")).astype(np.float64)
-    h, w = gray.shape
-    h2, w2 = h // 2 * 2, w // 2 * 2
-    gray = gray[:h2, :w2]
-
-    # Simple Haar wavelet decomposition
-    # Level 1
-    ll = (gray[0::2, 0::2] + gray[0::2, 1::2] + gray[1::2, 0::2] + gray[1::2, 1::2]) / 4
-    lh = (gray[0::2, 0::2] + gray[0::2, 1::2] - gray[1::2, 0::2] - gray[1::2, 1::2]) / 4
-    hl = (gray[0::2, 0::2] - gray[0::2, 1::2] + gray[1::2, 0::2] - gray[1::2, 1::2]) / 4
-    hh = (gray[0::2, 0::2] - gray[0::2, 1::2] - gray[1::2, 0::2] + gray[1::2, 1::2]) / 4
-
-    # Analyze high-frequency subbands
-    hf_coeffs = np.concatenate([lh.ravel(), hl.ravel(), hh.ravel()])
-    hf_coeffs = hf_coeffs[hf_coeffs != 0]
-
-    if len(hf_coeffs) < 100:
-        return {"test": "Wavelet Coefficients", "score": 0.0, "note": "Insufficient data"}
-
-    kurt = float(scipy_kurtosis(hf_coeffs, fisher=True))
-    # Generalized Gaussian: natural images have kurtosis > 3
-    # AI images: often lower kurtosis (more Gaussian-like)
-
-    if kurt > 5.0:
-        score = -0.3
-        note = f"Heavy-tailed wavelet coefficients (kurtosis={kurt:.2f}, natural)"
-    elif kurt < 1.5:
-        score = 0.4
-        note = f"Gaussian-like wavelet coefficients (kurtosis={kurt:.2f}, AI-like)"
-    else:
-        score = 0.0
-        note = f"Wavelet kurtosis={kurt:.2f}"
-
-    return {
-        "test": "Wavelet Coefficients",
-        "hf_kurtosis": round(kurt, 4),
-        "score": score,
-        "note": note,
-    }
-
-
-# ─── Main Agent Entry Point ─────────────────────────────────────────
-def run_statistical_agent(img: Image.Image) -> AgentEvidence:
-    """Run all statistical priors tests."""
-    findings = []
-    scores = []
-
-    for fn in [analyze_dct_distribution, analyze_benford, analyze_gradient_sparsity,
-               analyze_local_kurtosis, analyze_color_histogram, analyze_wavelet_coefficients]:
-        try:
-            result = fn(img)
-            findings.append(result)
-            scores.append(result["score"])
-        except Exception as e:
-            findings.append({"test": fn.__name__, "error": str(e), "score": 0})
-
-    avg_score = float(np.mean(scores)) if scores else 0.0
-    confidence = min(1.0, 0.5 + 0.5 * abs(avg_score))
-
-    violations = [f["test"] for f in findings if f.get("score", 0) > 0.2]
-    compliant = [f["test"] for f in findings if f.get("score", 0) < -0.1]
-
-    if violations:
-        rationale = f"Statistical violations: {', '.join(violations)}."
-    elif compliant:
-        rationale = f"Natural statistics confirmed: {', '.join(compliant)}."
-    else:
-        rationale = "Statistical analysis inconclusive."
-
+def _g(img): return np.array(img.convert("L")).astype(np.float64)
+def _rgb(img): return np.array(img.convert("RGB")).astype(np.float64)
+
+def t01_dct_kurtosis(img):
+    gray=_g(img); h,w=gray.shape; hc,wc=(h//8)*8,(w//8)*8; gray=gray[:hc,:wc]
+    coeffs=[]
+    for i in range(0,hc,8):
+        for j in range(0,wc,8):
+            b=gray[i:i+8,j:j+8]; d=dct(dct(b.T,norm="ortho").T,norm="ortho"); ac=d.copy(); ac[0,0]=0
+            coeffs.extend(ac.ravel().tolist())
+    c=np.array(coeffs); c=c[c!=0]
+    if len(c)<100: return {"test":"DCT Kurtosis","score":0.0,"note":"Insufficient data"}
+    k=float(sp_kurt(c,fisher=True))
+    if k>4.5: s,n=-0.4,f"Laplacian DCT (κ={k:.2f})"
+    elif k<2.0: s,n=0.5,f"Gaussian DCT (κ={k:.2f})"
+    elif k<3.5: s,n=0.2,f"Borderline (κ={k:.2f})"
+    else: s,n=-0.1,f"Near-natural (κ={k:.2f})"
+    return {"test":"DCT Kurtosis","kurtosis":round(k,4),"score":s,"note":n}
+
+def t02_benford(img):
+    gray=_g(img); h,w=gray.shape; hc,wc=(h//8)*8,(w//8)*8; gray=gray[:hc,:wc]
+    coeffs=[]
+    for i in range(0,hc,8):
+        for j in range(0,wc,8):
+            coeffs.extend(np.abs(dct(dct(gray[i:i+8,j:j+8].T,norm="ortho").T,norm="ortho").ravel()).tolist())
+    c=np.array(coeffs); nz=c[c>0]
+    if len(nz)<100: return {"test":"Benford's Law","score":0.0,"note":"Insufficient"}
+    lv=np.floor(np.log10(nz+1e-12)); fd=np.floor(nz/(10**lv)).astype(int); fd=fd[(fd>=1)&(fd<=9)]
+    obs=np.array([np.sum(fd==d) for d in range(1,10)],dtype=float); obs/=(obs.sum()+1e-9)
+    ben=np.log10(1+1.0/np.arange(1,10))
+    chi2=float(np.sum((obs-ben)**2/(ben+1e-9)))
+    if chi2<0.005: s,n=-0.4,f"Excellent Benford fit (χ²={chi2:.5f})"
+    elif chi2<0.02: s,n=-0.1,f"Good fit (χ²={chi2:.5f})"
+    elif chi2<0.05: s,n=0.3,f"Moderate deviation (χ²={chi2:.5f})"
+    else: s,n=0.6,f"Strong violation (χ²={chi2:.5f})"
+    return {"test":"Benford's Law","chi2":round(chi2,6),"observed":obs.tolist(),"benford_expected":ben.tolist(),"score":s,"note":n}
+
+def t03_gradient_sparsity(img):
+    gray=_g(img); gx=np.diff(gray,axis=1).ravel(); gy=np.diff(gray,axis=0).ravel()
+    ag=np.concatenate([gx,gy]); k=float(sp_kurt(ag,fisher=True))
+    thr=np.std(ag)*0.1; sp=float(np.mean(np.abs(ag)<thr))
+    if k>5 and sp>0.4: s,n=-0.4,f"Sparse gradients (κ={k:.2f}, sp={sp:.2f})"
+    elif k<2: s,n=0.5,f"Low kurtosis ({k:.2f})"
+    elif k<3.5: s,n=0.2,f"Borderline ({k:.2f})"
+    else: s,n=-0.1,f"Normal (κ={k:.2f})"
+    return {"test":"Gradient Sparsity","kurtosis":round(k,4),"sparsity":round(sp,4),"score":s,"note":n}
+
+def t04_local_kurtosis(img):
+    gray=_g(img); h,w=gray.shape; bs=32; hc,wc=(h//bs)*bs,(w//bs)*bs; gray=gray[:hc,:wc]
+    lk=[]
+    for i in range(0,hc,bs):
+        for j in range(0,wc,bs):
+            b=gray[i:i+bs,j:j+bs].ravel()
+            if np.std(b)>1: lk.append(float(sp_kurt(b,fisher=True)))
+    if len(lk)<10: return {"test":"Local Kurtosis Map","score":0.0,"note":"Insufficient"}
+    std=float(np.std(lk))
+    if std>3: s,n=-0.3,f"High kurtosis variation (σ={std:.2f})"
+    elif std<1: s,n=0.4,f"Uniform statistics (σ={std:.2f})"
+    else: s,n=0.0,f"Moderate (σ={std:.2f})"
+    return {"test":"Local Kurtosis Map","kurtosis_std":round(std,4),"score":s,"note":n}
+
+def t05_color_histogram(img):
+    rgb=np.array(img.convert("RGB")); scores=[]
+    for c in range(3):
+        h,_=np.histogram(rgb[:,:,c].ravel(),bins=256,range=(0,256))
+        sm=gaussian_filter(h.astype(float),2)
+        scores.append(float(np.mean(np.abs(h-sm))/(np.mean(h)+1e-9)))
+    avg=float(np.mean(scores))
+    if avg<0.3: s,n=-0.2,f"Smooth histograms ({avg:.3f})"
+    elif avg>0.8: s,n=0.4,f"Irregular histograms ({avg:.3f})"
+    else: s,n=0.0,f"Histogram smoothness={avg:.3f}"
+    return {"test":"Color Histogram","smoothness":round(avg,4),"score":s,"note":n}
+
+def t06_wavelet_kurtosis(img):
+    gray=_g(img); h,w=gray.shape; h2,w2=h//2*2,w//2*2; gray=gray[:h2,:w2]
+    lh=(gray[0::2,0::2]+gray[0::2,1::2]-gray[1::2,0::2]-gray[1::2,1::2])/4
+    hl=(gray[0::2,0::2]-gray[0::2,1::2]+gray[1::2,0::2]-gray[1::2,1::2])/4
+    hh=(gray[0::2,0::2]-gray[0::2,1::2]-gray[1::2,0::2]+gray[1::2,1::2])/4
+    hf=np.concatenate([lh.ravel(),hl.ravel(),hh.ravel()]); hf=hf[hf!=0]
+    if len(hf)<100: return {"test":"Wavelet Kurtosis","score":0.0,"note":"Insufficient"}
+    k=float(sp_kurt(hf,fisher=True))
+    if k>5: s,n=-0.3,f"Heavy-tailed wavelets (κ={k:.2f})"
+    elif k<1.5: s,n=0.4,f"Gaussian wavelets (κ={k:.2f})"
+    else: s,n=0.0,f"Wavelet κ={k:.2f}"
+    return {"test":"Wavelet Kurtosis","kurtosis":round(k,4),"score":s,"note":n}
+
+def t07_entropy_map(img):
+    gray=_g(img); h,w=gray.shape; bs=32; ents=[]
+    for i in range(0,h-bs,bs):
+        for j in range(0,w-bs,bs):
+            b=gray[i:i+bs,j:j+bs].ravel().astype(int)
+            h_,_=np.histogram(b,bins=64,range=(0,256)); h_=h_.astype(float); h_/=(h_.sum()+1e-9)
+            ents.append(-float(np.sum(h_*np.log2(h_+1e-12))))
+    if len(ents)<4: return {"test":"Entropy Map","score":0.0,"note":"Too small"}
+    std=float(np.std(ents)); mn=float(np.mean(ents))
+    if std>0.5: s,n=-0.2,f"Varied local entropy (σ={std:.2f})"
+    elif std<0.15: s,n=0.3,f"Uniform entropy (σ={std:.2f})"
+    else: s,n=0.0,f"Entropy σ={std:.2f}"
+    return {"test":"Entropy Map","entropy_std":round(std,4),"mean":round(mn,4),"score":s,"note":n}
+
+def t08_edge_orientation(img):
+    gray=_g(img); gx=sobel(gray,1); gy=sobel(gray,0); mag=np.hypot(gx,gy)
+    strong=mag>np.percentile(mag,80); angles=np.arctan2(gy[strong],gx[strong])
+    hist,_=np.histogram(angles,bins=36,range=(-np.pi,np.pi)); hist=hist.astype(float); hist/=(hist.sum()+1e-9)
+    ent=-float(np.sum(hist*np.log(hist+1e-9)))
+    max_ent=np.log(36)
+    norm_ent=ent/max_ent
+    if norm_ent<0.85: s,n=-0.2,f"Directional edges (entropy={norm_ent:.3f})"
+    elif norm_ent>0.95: s,n=0.2,f"Isotropic edges ({norm_ent:.3f})"
+    else: s,n=0.0,f"Edge entropy={norm_ent:.3f}"
+    return {"test":"Edge Orientation","entropy":round(norm_ent,4),"score":s,"note":n}
+
+def t09_lbp_distribution(img):
+    gray=np.array(img.convert("L")); h,w=gray.shape
+    # Simplified LBP
+    lbp=np.zeros((h-2,w-2),dtype=int)
+    for dy,dx,bit in [(-1,-1,0),(-1,0,1),(-1,1,2),(0,1,3),(1,1,4),(1,0,5),(1,-1,6),(0,-1,7)]:
+        lbp|=((gray[1+dy:h-1+dy,1+dx:w-1+dx]>=gray[1:h-1,1:w-1]).astype(int)<<bit)
+    hist,_=np.histogram(lbp.ravel(),bins=256,range=(0,256)); hist=hist.astype(float); hist/=(hist.sum()+1e-9)
+    # Uniform LBP patterns (≤2 transitions) dominate in natural images
+    uniform=0
+    for v in range(256):
+        b=format(v,'08b'); t=sum(1 for i in range(7) if b[i]!=b[i+1])+int(b[0]!=b[7])
+        if t<=2: uniform+=hist[v]
+    if uniform>0.6: s,n=-0.2,f"Natural LBP (uniform={uniform:.2%})"
+    elif uniform<0.3: s,n=0.3,f"Non-uniform LBP ({uniform:.2%})"
+    else: s,n=0.0,f"LBP uniform={uniform:.2%}"
+    return {"test":"LBP Distribution","uniform_ratio":round(uniform,4),"score":s,"note":n}
+
+def t10_cooccurrence(img):
+    gray=(np.array(img.convert("L"))//16).astype(int); h,w=gray.shape  # Quantize to 16 levels
+    glcm=np.zeros((16,16))
+    for i in range(h):
+        for j in range(w-1):
+            glcm[gray[i,j],gray[i,j+1]]+=1
+    glcm/=(glcm.sum()+1e-9)
+    # Energy (angular second moment)
+    energy=float(np.sum(glcm**2))
+    # Homogeneity
+    I,J=np.mgrid[0:16,0:16]; homog=float(np.sum(glcm/(1+np.abs(I-J))))
+    if energy<0.05 and homog>0.5: s,n=-0.2,f"Natural texture (E={energy:.4f}, H={homog:.3f})"
+    elif energy>0.2: s,n=0.3,f"Flat/repetitive (E={energy:.4f})"
+    else: s,n=0.0,f"GLCM E={energy:.4f}, H={homog:.3f}"
+    return {"test":"Co-occurrence Matrix","energy":round(energy,4),"homogeneity":round(homog,4),"score":s,"note":n}
+
+def t11_block_variance(img):
+    gray=_g(img); h,w=gray.shape; bs=8; hc,wc=(h//bs)*bs,(w//bs)*bs
+    gray=gray[:hc,:wc]; bvars=[]
+    for i in range(0,hc,bs):
+        for j in range(0,wc,bs):
+            bvars.append(float(np.var(gray[i:i+bs,j:j+bs])))
+    bv=np.array(bvars)
+    # ANOVA-like test: variance of variances
+    vov=float(np.std(bv))/(float(np.mean(bv))+1e-9)
+    if vov>1: s,n=-0.2,f"Varied block variance (VoV={vov:.3f})"
+    elif vov<0.3: s,n=0.3,f"Uniform block variance ({vov:.3f})"
+    else: s,n=0.0,f"VoV={vov:.3f}"
+    return {"test":"Block Variance ANOVA","vov":round(vov,4),"score":s,"note":n}
+
+def t12_gradient_magnitude(img):
+    gray=_g(img); gm=np.hypot(sobel(gray,0),sobel(gray,1))
+    k=float(sp_kurt(gm.ravel(),fisher=True)); sk=float(sp_skew(gm.ravel()))
+    if k>5: s,n=-0.2,f"Heavy-tailed gradients (κ={k:.2f})"
+    elif k<2: s,n=0.3,f"Light-tailed ({k:.2f})"
+    else: s,n=0.0,f"Gradient κ={k:.2f}"
+    return {"test":"Gradient Magnitude Dist","kurtosis":round(k,3),"skewness":round(sk,3),"score":s,"note":n}
+
+def t13_spatial_correlation(img):
+    gray=_g(img); h,w=gray.shape; step=max(1,h*w//200000)
+    ac1=float(np.corrcoef(gray[:,:-1].ravel()[::step],gray[:,1:].ravel()[::step])[0,1])
+    ac5=float(np.corrcoef(gray[:,:-5].ravel()[::step],gray[:,5:].ravel()[::step])[0,1])
+    decay=ac1-ac5
+    if 0.05<decay<0.3: s,n=-0.2,f"Natural correlation decay ({decay:.3f})"
+    elif decay<0.01: s,n=0.3,f"Flat correlation ({decay:.3f})"
+    else: s,n=0.0,f"Decay={decay:.3f}"
+    return {"test":"Spatial Correlation Decay","decay":round(decay,4),"score":s,"note":n}
+
+def t14_dct_skewness(img):
+    gray=_g(img); h,w=gray.shape; hc,wc=(h//8)*8,(w//8)*8; gray=gray[:hc,:wc]
+    coeffs=[]
+    for i in range(0,hc,8):
+        for j in range(0,wc,8):
+            d=dct(dct(gray[i:i+8,j:j+8].T,norm="ortho").T,norm="ortho"); ac=d.copy(); ac[0,0]=0
+            coeffs.extend(ac.ravel().tolist())
+    c=np.array(coeffs); c=c[c!=0]
+    if len(c)<100: return {"test":"DCT Skewness","score":0.0,"note":"Insufficient"}
+    sk=float(sp_skew(c))
+    if abs(sk)<0.1: s,n=-0.2,f"Symmetric DCT (skew={sk:.3f})"
+    elif abs(sk)>0.5: s,n=0.3,f"Skewed DCT ({sk:.3f})"
+    else: s,n=0.0,f"DCT skew={sk:.3f}"
+    return {"test":"DCT Skewness","skewness":round(sk,4),"score":s,"note":n}
+
+def t15_saturation_distribution(img):
+    rgb=np.array(img.convert("RGB")).astype(float)
+    mx=np.max(rgb,axis=-1); mn=np.min(rgb,axis=-1)
+    sat=(mx-mn)/(mx+1e-9); sat_flat=sat.ravel()
+    k=float(sp_kurt(sat_flat,fisher=True))
+    if k>3: s,n=-0.2,f"Natural saturation (κ={k:.2f})"
+    elif k<1: s,n=0.3,f"Unusual saturation ({k:.2f})"
+    else: s,n=0.0,f"Saturation κ={k:.2f}"
+    return {"test":"Saturation Distribution","kurtosis":round(k,3),"score":s,"note":n}
+
+def t16_luminance_gradient_ratio(img):
+    gray=_g(img); gx=np.abs(np.diff(gray,axis=1)); gy=np.abs(np.diff(gray,axis=0))
+    hg=float(np.mean(gx)); vg=float(np.mean(gy))
+    ratio=hg/(vg+1e-9)
+    if 0.7<ratio<1.4: s,n=-0.1,f"Balanced H/V gradients ({ratio:.3f})"
+    elif ratio>2 or ratio<0.5: s,n=0.2,f"Extreme H/V bias ({ratio:.3f})"
+    else: s,n=0.0,f"H/V ratio={ratio:.3f}"
+    return {"test":"H/V Gradient Ratio","ratio":round(ratio,3),"score":s,"note":n}
+
+def t17_pixel_uniqueness(img):
+    gray=np.array(img.convert("L")); total=gray.size; unique=len(np.unique(gray))
+    ratio=unique/256
+    if ratio>0.9: s,n=-0.1,f"Full tonal range ({unique} levels)"
+    elif ratio<0.5: s,n=0.2,f"Limited range ({unique} levels)"
+    else: s,n=0.0,f"{unique} levels"
+    return {"test":"Pixel Uniqueness","levels":unique,"score":s,"note":n}
+
+def t18_global_entropy(img):
+    gray=np.array(img.convert("L")); hist,_=np.histogram(gray,bins=256,range=(0,256))
+    hist=hist.astype(float); hist/=(hist.sum()+1e-9)
+    ent=-float(np.sum(hist*np.log2(hist+1e-12)))
+    if 6<ent<7.8: s,n=-0.2,f"Natural entropy ({ent:.3f})"
+    elif ent<5: s,n=0.3,f"Low entropy ({ent:.3f})"
+    else: s,n=0.0,f"Entropy={ent:.3f}"
+    return {"test":"Global Entropy","entropy":round(ent,4),"score":s,"note":n}
+
+def t19_power_law_fit(img):
+    gray=_g(img); gm=np.hypot(sobel(gray,0),sobel(gray,1)).ravel()
+    gm=gm[gm>1]; hist,edges=np.histogram(gm,bins=50); hist=hist.astype(float)+1
+    centers=(edges[:-1]+edges[1:])/2; valid=hist>1
+    if np.sum(valid)<5: return {"test":"Power Law Gradient","score":0.0,"note":"Insufficient"}
+    try:
+        c=np.polyfit(np.log(centers[valid]),np.log(hist[valid]),1); slope=float(c[0])
+    except: slope=0
+    if -3<slope<-1: s,n=-0.2,f"Power-law gradients (α={slope:.2f})"
+    elif slope>-0.5: s,n=0.3,f"Non-power-law ({slope:.2f})"
+    else: s,n=0.0,f"Slope={slope:.2f}"
+    return {"test":"Power Law Gradient","slope":round(slope,3),"score":s,"note":n}
+
+def t20_contrast_distribution(img):
+    gray=_g(img); h,w=gray.shape; bs=16
+    contrasts=[]
+    for i in range(0,h-bs,bs):
+        for j in range(0,w-bs,bs):
+            b=gray[i:i+bs,j:j+bs]; contrasts.append(float(np.max(b)-np.min(b)))
+    c=np.array(contrasts)
+    if len(c)<10: return {"test":"Contrast Distribution","score":0.0,"note":"Insufficient"}
+    k=float(sp_kurt(c,fisher=True))
+    if k>2: s,n=-0.2,f"Natural contrast variation (κ={k:.2f})"
+    elif k<0.5: s,n=0.2,f"Uniform contrast ({k:.2f})"
+    else: s,n=0.0,f"Contrast κ={k:.2f}"
+    return {"test":"Contrast Distribution","kurtosis":round(k,3),"score":s,"note":n}
+
+def t21_joint_histogram(img):
+    rgb=np.array(img.convert("RGB")); r,g=rgb[:,:,0].ravel(),rgb[:,:,1].ravel()
+    step=max(1,len(r)//100000)
+    h2d,_,_=np.histogram2d(r[::step],g[::step],bins=32,range=[[0,256],[0,256]])
+    h2d/=(h2d.sum()+1e-9)
+    # Mutual information
+    hr=np.sum(h2d,axis=1); hg=np.sum(h2d,axis=0)
+    mi=float(np.sum(h2d*np.log2(h2d/(np.outer(hr,hg)+1e-12)+1e-12)))
+    if mi>0.5: s,n=-0.2,f"Natural color correlation (MI={mi:.3f})"
+    elif mi<0.1: s,n=0.2,f"Weak color correlation ({mi:.3f})"
+    else: s,n=0.0,f"MI={mi:.3f}"
+    return {"test":"Joint Color Histogram","mi":round(mi,4),"score":s,"note":n}
+
+def t22_run_length(img):
+    gray=np.array(img.convert("L")); row=gray[gray.shape[0]//2,:]
+    runs=[]; cur=1
+    for i in range(1,len(row)):
+        if row[i]==row[i-1]: cur+=1
+        else: runs.append(cur); cur=1
+    runs.append(cur); runs=np.array(runs)
+    avg_run=float(np.mean(runs)); max_run=int(np.max(runs))
+    if 1<avg_run<5: s,n=-0.2,f"Natural run lengths (avg={avg_run:.2f})"
+    elif avg_run>10: s,n=0.3,f"Long runs ({avg_run:.2f}) — flat patches"
+    else: s,n=0.0,f"Run avg={avg_run:.2f}"
+    return {"test":"Run Length Analysis","avg_run":round(avg_run,3),"score":s,"note":n}
+
+ALL_TESTS=[t01_dct_kurtosis,t02_benford,t03_gradient_sparsity,t04_local_kurtosis,t05_color_histogram,
+           t06_wavelet_kurtosis,t07_entropy_map,t08_edge_orientation,t09_lbp_distribution,t10_cooccurrence,
+           t11_block_variance,t12_gradient_magnitude,t13_spatial_correlation,t14_dct_skewness,
+           t15_saturation_distribution,t16_luminance_gradient_ratio,t17_pixel_uniqueness,t18_global_entropy,
+           t19_power_law_fit,t20_contrast_distribution,t21_joint_histogram,t22_run_length]
+
+def run_statistical_agent(img):
+    findings,scores=[],[]
+    for fn in ALL_TESTS:
+        try: r=fn(img); findings.append(r); scores.append(r["score"])
+        except Exception as e: findings.append({"test":fn.__name__,"error":str(e),"score":0})
+    avg=float(np.mean(scores)) if scores else 0.0; conf=min(1.0,0.5+0.5*abs(avg))
+    viol=[f["test"] for f in findings if f.get("score",0)>0.2]
+    comp=[f["test"] for f in findings if f.get("score",0)<-0.1]
+    rat=f"Statistical violations: {', '.join(viol)}." if viol else f"Natural statistics: {', '.join(comp)}." if comp else "Statistical inconclusive."
     for f in findings:
-        if f.get("note"):
-            rationale += f" [{f['test']}]: {f['note']}."
-
-    return AgentEvidence(
-        agent_name="Statistical Priors Agent",
-        violation_score=np.clip(avg_score, -1, 1),
-        confidence=confidence,
-        failure_prob=max(0.0, 1.0 - len(scores) / 6),
-        rationale=rationale,
-        sub_findings=findings,
-    )
+        if f.get("note"): rat+=f" [{f['test']}]: {f['note']}."
+    return AgentEvidence("Statistical Priors Agent",np.clip(avg,-1,1),conf,max(0,1-len(scores)/len(ALL_TESTS)),rat,findings)