"""Full-cohort family ecology and common-coordinate DINO analysis; pinned input provenance."""
import os
for k in ['OMP_NUM_THREADS','OPENBLAS_NUM_THREADS','MKL_NUM_THREADS']:os.environ[k]='2'
import json, math, hashlib, collections, itertools, time, sys, csv
from pathlib import Path
import numpy as np
from scipy.spatial.distance import pdist,squareform
from scipy.linalg import eigh,pinvh
from scipy.ndimage import label
from scipy.stats import spearmanr,rankdata
from sklearn.isotonic import IsotonicRegression
from sklearn.decomposition import PCA
from PIL import Image
from threadpoolctl import threadpool_limits
import matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plt

P=Path('/home/masahiro/kics-zert2');BASE=P/'runs/moin_plot_all_20260928'
ROOT=P/'runs/moin_diversity_20261005';OUT=ROOT/'report';OUT.mkdir(parents=True,exist_ok=True)
for n in ['figures','images','data']:(OUT/n).mkdir(exist_ok=True)
NAMES={'fessbach':'Feßbach','heidesheim':'Heidesheim','rettmer':'Rettmer','ribbeck':'Ribbeck'}
COLORS={'fessbach':'#277959','heidesheim':'#b0812e','rettmer':'#4c78b6','ribbeck':'#ad5673'}
read=lambda p:json.loads(p.read_text())
def save(p,x):
    p.parent.mkdir(parents=True,exist_ok=True);tmp=p.with_suffix('.tmp')
    tmp.write_text(json.dumps(x,ensure_ascii=False,indent=2,allow_nan=False));tmp.replace(p)
def status(stage,**kw):
    save(ROOT/'analysis_status.json',{'stage':stage,'updated_utc':__import__('datetime').datetime.now(__import__('datetime').timezone.utc).isoformat(),**kw})
def plain(x):
    if isinstance(x,np.ndarray):return plain(x.tolist())
    if isinstance(x,(np.integer,np.floating)):return plain(x.item())
    if isinstance(x,float) and not math.isfinite(x):return None
    if isinstance(x,dict):return {str(k):plain(v) for k,v in x.items()}
    if isinstance(x,(list,tuple)):return [plain(v) for v in x]
    return x
def hill(p):
    p=np.asarray(p,float);p=p[p>0]
    if len(p)==0:return {'q0':0,'H':None,'q1':None,'q2':None,'evenness':None}
    p=p/p.sum();h=float(-np.sum(p*np.log(p)))
    return {'q0':len(p),'H':h,'q1':math.exp(h),'q2':float(1/np.sum(p*p)),'evenness':h/math.log(len(p)) if len(p)>1 else None}
def family_of(g):return g['path'][4] if g['rank'] in ['family','genus','species'] and len(g['path'])>=5 else None

def collect():
    cache=ROOT/'collected.json'
    if cache.exists():return read(cache),np.load(ROOT/'embedding_moments.npz')['means']
    result=read(BASE/'report_threshold_0p7/results.json');selection=read(BASE/'selection.json')
    meta=read(ROOT/'metadata_join.json');mm={x['id']:x for x in meta['images']}
    families=sorted({family_of(g) for r in result['images'] for g in r['taxonomy_supported']['groups'] if family_of(g)})
    fi={f:i for i,f in enumerate(families)};rows=[];means=[]
    for i,r in enumerate(result['images']):
        key=r['id'];status('collecting_common_features',completed=i,total=439,current=key)
        z=np.load(BASE/key/'tokens.npz');full=z['full'];v=z['raw'][full].astype(np.float64)
        norms=np.linalg.norm(v,axis=1);assert np.isfinite(v).all() and np.all(norms>0)
        v/=norms[:,None];mu=v.mean(0);means.append(mu);q=float(1-mu@mu)
        xy=z['xy'][full];lookup={tuple(c):j for j,c in enumerate(xy)};dots=[]
        for j,(x,y) in enumerate(xy):
            for neighbor in [(x+1,y),(x,y+1)]:
                if neighbor in lookup:dots.append(float(v[j]@v[lookup[neighbor]]))
        area=np.zeros(len(families))
        for g in r['taxonomy_supported']['groups']:
            f=family_of(g)
            if f:area[fi[f]]+=g['area']
        p=area/area.sum() if area.sum()>0 else area
        # Family configuration is measured only where the regional hypothesis is unambiguous.
        reg=np.asarray(Image.open(BASE/'taxonomy_regions'/key/'policy_labels.png'))
        excluded=np.load(BASE/key/'artifact_exclusion.npz')['excluded']
        mapping=collections.defaultdict(list)
        for g in r['taxonomy_supported']['groups']:
            for case in g['case_ids']:mapping[case].append(family_of(g))
        fmap=np.full(reg.shape,-1,np.int16)
        for c in r['taxonomic_case_summary']:
            fs=mapping[c['id']]
            if fs and None not in fs and len(set(fs))==1:fmap[np.isin(reg,c['source_cluster_ids'])&~excluded]=fi[fs[0]]
        h,w=fmap.shape;gh=round(h*256/max(h,w));gw=round(w*256/max(h,w))
        ys=np.minimum(((np.arange(gh)+.5)*h/gh).astype(int),h-1);xs=np.minimum(((np.arange(gw)+.5)*w/gw).astype(int),w-1)
        grid=fmap[ys[:,None],xs[None,:]];valid=grid>=0;joins=0;edges=0
        for left,right in [(grid[:,:-1],grid[:,1:]),(grid[:-1],grid[1:])]:
            both=(left>=0)&(right>=0);edges+=int(both.sum());joins+=int(((left==right)&both).sum())
        join=joins/edges if edges else None
        conc=float(np.sum((np.bincount(grid[valid],minlength=len(families))/valid.sum())**2)) if valid.any() else None
        clump=(join-conc)/(1-conc) if join is not None and conc is not None and conc<1-1e-12 else None
        components=sum(label(grid==f)[1] for f in np.unique(grid[valid])) if valid.any() else 0
        palette=np.vstack([[180,180,175],plt.get_cmap('tab20')(np.linspace(0,1,len(families)))[:,:3]*255]).astype(np.uint8)
        Image.fromarray(palette[grid+1]).save(OUT/'images'/f'{key}_family_configuration.png')
        original=Image.open(BASE/key/'original.jpg');original.thumbnail((380,380));original.save(OUT/'images'/f'{key}.jpg',quality=88)
        month=r['folder_date'][5:7];season='Spring (May)' if month=='05' else 'Summer (Jun–Aug)'
        row={'id':key,'site':r['site'],'campaign':r['folder_date'],'session':r['session'],'month':'2025-'+month,'season':season,'archive_path':r['archive_path'],'source_sha256':r['source_sha256'],'date_conflict':r['date_conflict'],'exif_datetime':r['exif_datetime'],'duplicate_aliases':r['aliases'],'families':{f:float(area[j]) for f,j in fi.items() if area[j]>0},'family_proportions':p.tolist(),'family_diversity':hill(p),'family_assignment_fraction':float(area.sum()/r['taxonomic_eligible_pixels']),'named_group_diversity':r['taxonomy_supported'],'old_named_rank':r['taxonomic_rank'],'old_visual_composite':r['visual_score'],'old_visual_rank':r['visual_rank'],'old_visual_components':r['visual_components'],'old_visual_rank_interval':r['visual_rank_interval'],'whole_structure':{'rao_Q_all_patches':q,'patches':len(v),'adjacent_cosine_dissimilarity':float(1-np.mean(dots)),'mean_embedding_norm':float(np.linalg.norm(mu))},'configuration':{'unambiguous_family_grid_fraction':float(valid.mean()),'adjacent_same_family_fraction':join,'expected_same_family_fraction':conc,'excess_clumping':clump,'connected_components':components,'components_per_10000_resolved_cells':components/max(int(valid.sum()),1)*10000 if valid.any() else None,'grid_shape':[gh,gw],'valid_adjacent_pairs':edges},'plot_metadata':mm.get(key)}
        rows.append(row)
    assert len(rows)==439 and len({r['source_sha256'] for r in rows})==439
    np.savez_compressed(ROOT/'embedding_moments.npz',means=np.array(means))
    save(cache,{'families':families,'images':rows,'input_hashes':{n:hashlib.sha256((BASE/'report_threshold_0p7'/n).read_bytes()).hexdigest() for n in ['results.json','selection.json']}})
    return read(cache),np.array(means)

def accumulation(sets,seed=42):
    n=len(sets);freq=collections.Counter(f for s in sets for f in s);S=len(freq)
    curve=[0.]+[sum(1-(math.comb(n-f,m)/math.comb(n,m) if n-f>=m else 0) for f in freq.values()) for m in range(1,n+1)]
    rng=np.random.default_rng(seed);sims=[]
    for _ in range(1000):
        seen=set();c=[0]
        for j in rng.permutation(n):seen.update(sets[j]);c.append(len(seen))
        sims.append(c)
    q1=sum(f==1 for f in freq.values());q2=sum(f==2 for f in freq.values())
    unseen=(n-1)/n*(q1*q1/(2*q2) if q2 else q1*(q1-1)/2) if n>1 else 0
    A=(n-1)*q1/((n-1)*q1+2*q2) if q2>0 else ((n-1)*(q1-1)/((n-1)*(q1-1)+2) if q1>1 else 0)
    T=sum(freq.values());coverage=1-q1/max(T,1)*A if T else None
    return {'n':n,'observed':S,'frequencies':dict(freq),'curve':curve,'order_band':np.quantile(sims,[.025,.975],axis=0).tolist(),'chao2_lower_bound':S+unseen,'incidence_sample_coverage':coverage,'singletons':q1,'doubletons':q2,'projection_to_2n':[S+unseen*(1-(1-q1/(n*unseen+q1))**k) if unseen else float(S) for k in range(n+1)]}

def pair_beta(A):
    A=(A>0).astype(int);shared=A@A.T;rich=A.sum(1);b=rich[:,None]-shared;c=rich[None,:]-shared
    sor=np.divide(b+c,2*shared+b+c,out=np.zeros_like(shared,float),where=(2*shared+b+c)>0)
    m=np.minimum(b,c);turn=np.divide(m,shared+m,out=np.zeros_like(shared,float),where=(shared+m)>0)
    return sor,turn,sor-turn

def nmds(distance,multiplicity,dimensions=2,starts=24,seed=20261005,max_iter=1200):
    """Weighted strong-tie nonmetric SMACOF. Exact duplicate profiles are collapsed."""
    n=len(distance);weights=np.outer(multiplicity,multiplicity).astype(float);np.fill_diagonal(weights,0)
    upper=np.triu_indices(n,1);d0=distance[upper];w0=weights[upper]
    assert np.all(d0>0), 'Collapse zero-distance profiles before NMDS.'
    V=np.diag(weights.sum(1))-weights;inv=pinvh(V);center=np.eye(n)-np.ones((n,n))/n
    val,vec=eigh(-.5*center@(distance**2)@center);init=vec[:,-dimensions:]*np.sqrt(np.maximum(val[-dimensions:],0))
    rng=np.random.default_rng(seed);fits=[];best=None
    for st in range(starts):
        X=init.copy() if st==0 else rng.normal(size=(n,dimensions))
        last=math.inf;converged=False
        for it in range(max_iter):
            dd=squareform(pdist(X));actual=dd[upper]
            target=IsotonicRegression(increasing=True,y_min=0,out_of_bounds='clip').fit_transform(d0,actual,sample_weight=w0)
            target*=math.sqrt(w0.sum()/max(float(np.sum(w0*target**2)),1e-30))
            disparities=np.zeros_like(distance);disparities[upper]=target;disparities+=disparities.T
            B=-weights*np.divide(disparities,dd,out=np.zeros_like(dd),where=dd>1e-15);np.fill_diagonal(B,-B.sum(1))
            new=inv@B@X;new-=np.average(new,axis=0,weights=multiplicity)
            stress=float(np.sum(w0*(pdist(new)-target)**2));relative=(last-stress)/max(last,1e-30) if math.isfinite(last) else math.inf
            X=new
            if 0<=relative<1e-8:converged=True;break
            last=stress
        dd=pdist(X);target=IsotonicRegression(increasing=True,y_min=0).fit_transform(d0,dd,sample_weight=w0)
        stress1=math.sqrt(float(np.sum(w0*(dd-target)**2)/max(np.sum(w0*dd**2),1e-30)))
        fits.append({'start':st,'stress1':stress1,'iterations':it+1,'converged':converged})
        if best is None or stress1<best['stress1']:best={'X':X,'stress1':stress1,'start':st,'converged':converged,'iterations':it+1,'disparities':target,'distances':dd,'dissimilarities':d0,'weights':w0}
    X=best['X'];_,_,vt=np.linalg.svd(X*np.sqrt(multiplicity[:,None]),full_matrices=False);best['X']=X@vt.T
    best['starts']=fits;return best

def writefig(fig,name):
    fig.savefig(OUT/'figures'/f'{name}.png',dpi=180,bbox_inches='tight');fig.savefig(OUT/'figures'/f'{name}.svg',bbox_inches='tight');plt.close(fig)

def analyze():
    data,means=collect();rows=data['images'];F=data['families'];p=np.array([r['family_proportions'] for r in rows]);A=p>0;positive=A.any(1)
    plt.rcParams.update({'font.family':'DejaVu Sans','font.size':10,'svg.fonttype':'none','axes.spines.top':False,'axes.spines.right':False})
    status('taxonomic_diversity_and_accumulation',completed=439,total=439)
    for field,values in [('family_rank',np.array([r['family_diversity']['q1'] or np.nan for r in rows])),('structure_rank',np.array([r['whole_structure']['rao_Q_all_patches'] for r in rows]))]:
        valid=np.isfinite(values);ranks=rankdata(-values[valid],method='min');k=0
        for i,r in enumerate(rows):r[field]=int(ranks[k]) if valid[i] else None;k+=int(valid[i])
    accumulation_records=[]
    for site in NAMES:
        sr=[r for r in rows if r['site']==site]
        for level in ['combined','campaign','month','season']:
            groups={'All campaigns':sr} if level=='combined' else {key:[r for r in sr if r[level]==key] for key in sorted({r[level] for r in sr})}
            for key,group in groups.items():accumulation_records.append({'site':site,'level':level,'group':key,**accumulation([set(r['families']) for r in group])})
    summaries=[];group_structures=[]
    for site in list(NAMES)+['all_sites']:
        ix=np.array([i for i,r in enumerate(rows) if site=='all_sites' or r['site']==site]);known=ix[positive[ix]];pp=p[known];AA=A[known]
        pooled=pp.mean(0);alpha0=float(AA.sum(1).mean());alpha1=float(math.exp(np.mean([-np.sum(x[x>0]*np.log(x[x>0])) for x in pp])));alpha2=float(1/np.mean(np.sum(pp*pp,axis=1)));gh=hill(pooled)
        sor,turn,nest=pair_beta(AA);ut=np.triu_indices(len(known),1)
        samecamp=np.array([rows[i]['campaign'] for i in known]);mask=samecamp[:,None]==samecamp[None,:];within=np.triu(mask,1)
        centroid=means[ix].mean(0);alpha=float(np.mean([rows[i]['whole_structure']['rao_Q_all_patches'] for i in ix]));beta=float(np.mean(np.sum((means[ix]-centroid)**2,axis=1)));gamma=float(1-centroid@centroid)
        assert abs(alpha+beta-gamma)<1e-10
        ra=accumulation([set(rows[i]['families']) for i in ix]);freq=AA.sum(0)
        sm={'site':site,'n_images':len(ix),'n_family_positive':len(known),'unresolved_family_images':len(ix)-len(known),'gamma_families':gh['q0'],'alpha_q0_conditional':alpha0,'alpha_q0_detected_all_images':float(A[ix].sum(1).mean()),'alpha_q1_conditional':alpha1,'alpha_q2_conditional':alpha2,'gamma_q1_conditional':gh['q1'],'gamma_q2_conditional':gh['q2'],'beta_q0_conditional':gh['q0']/alpha0,'beta_q1_conditional':gh['q1']/alpha1,'beta_q2_conditional':gh['q2']/alpha2,'mean_sorensen':float(sor[ut].mean()),'mean_turnover':float(turn[ut].mean()),'mean_nestedness':float(nest[ut].mean()),'within_campaign_mean_sorensen':float(sor[within].mean()) if within.any() else None,'within_campaign_mean_turnover':float(turn[within].mean()) if within.any() else None,'within_campaign_mean_nestedness':float(nest[within].mean()) if within.any() else None,'mean_bray_curtis':float(pdist(pp,'braycurtis').mean()),'rarefied_families_at_74_images':ra['curve'][74] if len(ix)>=74 else None,'incidence_sample_coverage':ra['incidence_sample_coverage'],'top_families':[{'family':F[j],'photos':int(freq[j]),'fraction_of_family_positive_images':float(freq[j]/len(known)),'mean_conditional_area_fraction':float(pooled[j])} for j in np.argsort(-freq)[:8]],'structural_alpha_Q':alpha,'structural_beta_Q':beta,'structural_gamma_Q':gamma,'structural_beta_fraction':beta/gamma,'configuration':{f:float(np.mean([rows[i]['configuration'][f] for i in ix if rows[i]['configuration'][f] is not None])) if any(rows[i]['configuration'][f] is not None for i in ix) else None for f in ['adjacent_same_family_fraction','excess_clumping','components_per_10000_resolved_cells','unambiguous_family_grid_fraction']}}
        summaries.append(sm)
        for campaign in sorted({rows[i]['campaign'] for i in ix}):
            ids=[i for i in ix if rows[i]['campaign']==campaign];c=means[ids].mean(0);al=np.mean([rows[i]['whole_structure']['rao_Q_all_patches'] for i in ids]);be=np.mean(np.sum((means[ids]-c)**2,axis=1));ga=1-c@c
            group_structures.append({'site':site,'campaign':campaign,'n':len(ids),'alpha_Q':float(al),'beta_Q':float(be),'gamma_Q':float(ga)})
    # Descriptive configuration dissimilarity; avoid undefined descriptors for pure single-family images.
    cfg=np.array([[r['configuration'][f] if r['configuration'][f] is not None else np.nan for f in ['adjacent_same_family_fraction','components_per_10000_resolved_cells','unambiguous_family_grid_fraction']] for r in rows])
    config_positive=np.all(np.isfinite(cfg),axis=1);ranges=np.nanmax(cfg,0)-np.nanmin(cfg,0);cfg_scaled=cfg/np.where(ranges>0,ranges,1)
    for sm in summaries:
        ix=[i for i,r in enumerate(rows) if config_positive[i] and (sm['site']=='all_sites' or r['site']==sm['site'])]
        sm['configuration']['n_comparable_images']=len(ix)
        sm['configuration']['mean_gower_configuration_distance']=float(pdist(cfg_scaled[ix],'cityblock').mean()/3) if len(ix)>1 else None
    # Ordinations on exact distinct profiles with photo-multiplicity weights.
    status('NMDS_multistart',family_positive=int(positive.sum()))
    ordinations={};fits={};known=np.flatnonzero(positive)
    for name,matrix,metric in [('incidence',A[positive].astype(float),'jaccard'),('area',p[positive],'braycurtis')]:
        unique,inverse,counts=np.unique(np.round(matrix,12),axis=0,return_inverse=True,return_counts=True);D=squareform(pdist(unique,metric))
        fit=nmds(D,counts,2);sens=nmds(D,counts,3,starts=12,max_iter=5000);fits[name]=fit
        full=fit['X'][inverse];full3=sens['X'][inverse];ordinations[name]={'metric':metric,'n_images':len(known),'distinct_profiles':len(unique),'stress1_2d':fit['stress1'],'stress1_3d':sens['stress1'],'starts':fit['starts'],'starts_3d':sens['starts'],'coordinates':[{'id':rows[i]['id'],'site':rows[i]['site'],'campaign':rows[i]['campaign'],'month':rows[i]['month'],'season':rows[i]['season'],'x':float(full[j,0]),'y':float(full[j,1]),'xyz_3d':full3[j].tolist()} for j,i in enumerate(known)],'converged':fit['converged'],'converged_3d':sens['converged'],'method':'Weighted strong-tie nonmetric SMACOF; identical profiles collapsed and multiplied by photo counts; no zero distances treated as missing.'}
        for j,i in enumerate(known):rows[i][f'nmds_{name}']=full[j].tolist()
    status('shared_PCA_and_structural_partition')
    pc=PCA(svd_solver='full').fit(means);scores=pc.transform(means)
    for i,r in enumerate(rows):r['pca']=scores[i,:5].tolist()
    sitecenters=[];weights=[]
    for site in NAMES:
        ix=[i for i,r in enumerate(rows) if r['site']==site];sitecenters.append(means[ix].mean(0));weights.append(len(ix)/439)
    site_beta=float(np.sum(np.array(weights)*np.sum((np.array(sitecenters)-means.mean(0))**2,axis=1)))
    within_site_between_photos=float(np.sum([len([r for r in rows if r['site']==s['site']])/439*s['structural_beta_Q'] for s in summaries if s['site']!='all_sites']))
    global_s=next(x for x in summaries if x['site']=='all_sites');assert abs(site_beta+within_site_between_photos-global_s['structural_beta_Q'])<1e-10
    # Descriptive centroid separation in common vector space, with no independent-site p-values.
    total=float(np.sum((means-means.mean(0))**2));site_ss=site_beta*439
    # Spatial (same campaign) versus temporal (same recorded plot) composition.
    pairs=[];byplot=collections.defaultdict(list)
    import re
    for i,r in enumerate(rows):
        m=r['plot_metadata'];loc=None;point=None;origin=None
        if m and m['status']=='unique_campaign_match':
            values=m['matches'][0];loc=values.get('Location');point=values.get('Photopoint');origin='metadata'
        elif r['site']=='ribbeck':
            match=re.match(r'(.+?)_(P_[IV]+)(?:_[A-Z]+)?_\d{2}-\d{2}-\d{4}_',Path(r['archive_path']).name)
            if match:loc,point=match.groups();origin='filename'
        if loc and point:
            point=re.sub(r'[\s_-]+','_',point.strip()).upper();r['physical_plot_key']=r['site']+'::'+loc.upper()+'::'+point;r['physical_plot_key_source']=origin
            byplot[r['physical_plot_key']].append(i)
    for key,ids in byplot.items():
        for i,j in itertools.combinations(ids,2):
            if rows[i]['campaign']==rows[j]['campaign'] or not positive[i] or not positive[j]:continue
            sor,turn,nest=pair_beta(A[[i,j]])
            pairs.append({'plot_key':key,'site':rows[i]['site'],'id_a':rows[i]['id'],'id_b':rows[j]['id'],'campaign_a':rows[i]['campaign'],'campaign_b':rows[j]['campaign'],'sorensen':float(sor[0,1]),'turnover':float(turn[0,1]),'nestedness':float(nest[0,1]),'structural_centroid_distance':float(np.linalg.norm(means[i]-means[j]))})
    analysis={'audit_date':'2026-10-05','images':rows,'families':F,'site_summaries':summaries,'accumulation':accumulation_records,'NMDS':ordinations,'PCA':{'input':'439 photo means of full L2-normalized 768D whole-image patch features; photo means not renormalized; common model coordinates; centered PCA without whitening','explained_variance_ratio':pc.explained_variance_ratio_.tolist(),'site_explained_fraction_of_between_photo_variance':site_ss/total,'hierarchy':{'within_photo':global_s['structural_alpha_Q'],'between_photo_within_site':within_site_between_photos,'between_site':site_beta,'total_gamma':global_s['structural_gamma_Q']},'loadings_file':'data/shared_PCA_loadings.npz'},'structural_campaigns':group_structures,'repeated_plot_comparisons':pairs,'input_hashes':data['input_hashes'],'definitions':{'taxonomic':'Family-only supported model hypotheses at0.7, collapsed from family/genus/species lineages. No-family photos are unresolved, not confirmed empty. Hill q1/q2 and ordinations condition on known family composition; incidence curves include all photos as effort.','structural':'Continuous visual feature variance in full normalized DINO coordinates. Per-image alpha_i=1-|mean(z_i)|²; group alpha=mean(alpha_i), beta=mean|mean(z_i)-group_mean|², gamma=1-|group_mean|², so gamma=alpha+beta. Equal weight per image within a group. Measures are visual structural proxies, not validated biological diversity.','configuration':'Family-region configuration from older accepted0.7 region evidence; unambiguous single-family pixels only, sampled on aspect-preserving grid longest side256. Gower distance is the mean absolute normalized difference of same-family join fraction, component density, and resolved grid fraction; incomplete coverage can drive it.','season':'Meteorological spring=May; summer=June–August in this sampled interval. Only Rettmer has May photos. Campaign date follows archive folder; EXIF conflict is retained.','NMDS':'Strong-tie weighted nonmetric SMACOF,24 starts in2D and12 starts in3D. Identical observed profiles collapse to one point with photo-multiplicity pair weights. No-family images excluded. Coordinate axes arbitrary; report normalized Stress1 and Shepard fit.','uncertainty':'Accumulation shading is2.5–97.5% variation of1000 permutations of known images, not ecological confidence. No independent ecological-site significance or causal management inference.','taxonomic_evidence_version':'Existing AnyUp-region0.7 evidence is retained and explicitly identified. The finalized local clustering requires a separate full-cohort computation; do not claim old taxonomic labels came from it.'}}
    np.savez_compressed(OUT/'data/shared_PCA_loadings.npz',components=pc.components_,mean=pc.mean_,photo_mean_embeddings=means,photo_ids=np.array([r['id'] for r in rows]))
    save(OUT/'data/analysis.json',plain(analysis));save(OUT/'data/site_summary.json',plain(summaries));save(OUT/'data/image_rankings.json',plain(rows))
    make_plots(analysis,fits)
    status('core_analysis_complete',photos=439,family_positive=int(positive.sum()),NMDS_stress={k:v['stress1_2d'] for k,v in ordinations.items()})
    print(json.dumps(plain({'sites':summaries,'NMDS':{k:{t:v[t] for t in ['distinct_profiles','stress1_2d','stress1_3d','converged']} for k,v in ordinations.items()},'PCA_first5':pc.explained_variance_ratio_[:5],'structural_partition':analysis['PCA']['hierarchy'],'repeated_plot_pairs':len(pairs)}),ensure_ascii=False),flush=True)

def make_plots(a,fits):
    rows=a['images'];records=a['accumulation'];sums=[s for s in a['site_summaries'] if s['site']!='all_sites']
    fig,ax=plt.subplots(figsize=(10,6))
    for r in records:
        if r['level']!='combined':continue
        n=r['n'];ax.plot(range(n+1),r['curve'],label=f"{NAMES[r['site']]} (N={n}; {r['observed']} families)",color=COLORS[r['site']]);ax.fill_between(range(n+1),*r['order_band'],alpha=.1,color=COLORS[r['site']])
    ax.set(xlabel='Number of distinct images',ylabel='Cumulative detected families',title='Family accumulation across all campaigns');ax.legend();ax.grid(alpha=.15);writefig(fig,'family_accumulation_all')
    for level in ['campaign','month','season']:
        fig,axs=plt.subplots(2,2,figsize=(12,9))
        for ax,site in zip(axs.flat,NAMES):
            rr=[r for r in records if r['site']==site and r['level']==level]
            for j,r in enumerate(rr):
                ax.plot(range(r['n']+1),r['curve'],label=f"{r['group']} (N={r['n']})",color=plt.get_cmap('viridis')(j/max(len(rr)-1,1)));ax.fill_between(range(r['n']+1),*r['order_band'],alpha=.06)
            combo=next(r for r in records if r['site']==site and r['level']=='combined');ax.plot(range(combo['n']+1),combo['curve'],ls='--',color='black',label='All combined')
            ax.set(title=NAMES[site],xlabel='Number of distinct images',ylabel='Detected families');ax.legend(fontsize=8);ax.grid(alpha=.15)
        fig.suptitle(f'Family accumulation by {level}; shading is image-order variation',fontsize=14);fig.tight_layout();writefig(fig,'family_accumulation_'+level)
    freq=np.array([[sum(f in r['families'] for r in rows if r['site']==s['site'])/s['n_images'] for f in a['families']] for s in sums])
    fig,ax=plt.subplots(figsize=(12,4));im=ax.imshow(freq,aspect='auto',vmin=0,vmax=1,cmap='YlGn');ax.set(yticks=range(4),yticklabels=[NAMES[s['site']] for s in sums],xticks=range(len(a['families'])),xticklabels=a['families'],title='Observed family detection frequency; denominator is all photos');plt.setp(ax.get_xticklabels(),rotation=60,ha='right');fig.colorbar(im,ax=ax,label='Fraction of photos');writefig(fig,'family_frequency_heatmap')
    fig,axs=plt.subplots(1,2,figsize=(12,5));x=np.arange(4)
    for ax,prefix,title in [(axs[0],'','All campaigns pooled'),(axs[1],'within_campaign_','Within campaigns: spatial heterogeneity')]:
        turn=[s[prefix+'mean_turnover'] for s in sums];nest=[s[prefix+'mean_nestedness'] for s in sums]
        ax.bar(x,turn,label='Turnover',color='#377c68');ax.bar(x,nest,bottom=turn,label='Nestedness component',color='#d4b871');ax.set(xticks=x,xticklabels=[NAMES[s['site']] for s in sums],ylabel='Mean pairwise Sørensen dissimilarity',ylim=(0,1),title=title);ax.legend(fontsize=9)
    fig.tight_layout();writefig(fig,'taxonomic_beta_partition')
    for name,ord in a['NMDS'].items():
        fig,axs=plt.subplots(1,2,figsize=(12,5))
        for site in NAMES:
            rr=[r for r in ord['coordinates'] if r['site']==site];axs[0].scatter([r['x'] for r in rr],[r['y'] for r in rr],s=24,alpha=.55,label=NAMES[site],color=COLORS[site])
        for month in sorted({r['month'] for r in ord['coordinates']}):
            rr=[r for r in ord['coordinates'] if r['month']==month];axs[1].scatter([r['x'] for r in rr],[r['y'] for r in rr],s=24,alpha=.55,label=month)
        for ax in axs:ax.set(xlabel='NMDS1',ylabel='NMDS2');ax.legend(fontsize=9);ax.axhline(0,color='#ccc',lw=.5);ax.axvline(0,color='#ccc',lw=.5)
        fig.suptitle(f"Family {name} NMDS · Stress1={ord['stress1_2d']:.3f} · N={ord['n_images']}");fig.tight_layout();writefig(fig,'nmds_'+name)
        fig=plt.figure(figsize=(12,5));axs=[fig.add_subplot(121,projection='3d'),fig.add_subplot(122,projection='3d')]
        for field,ax in zip(['site','month'],axs):
            for group in sorted({r[field] for r in ord['coordinates']}):
                cc=np.array([r['xyz_3d'] for r in ord['coordinates'] if r[field]==group]);ax.scatter(cc[:,0],cc[:,1],cc[:,2],s=15,alpha=.5,label=NAMES.get(group,group),color=COLORS.get(group,None))
            ax.set(xlabel='NMDS1',ylabel='NMDS2',zlabel='NMDS3');ax.legend(fontsize=8)
        fig.suptitle(f"Family {name} 3D NMDS · Stress1={ord['stress1_3d']:.3f} · N={ord['n_images']}");fig.tight_layout();writefig(fig,'nmds_'+name+'_3d')
        fig,axs=plt.subplots(1,2,figsize=(12,5));fit=fits[name]
        axs[0].scatter(fit['dissimilarities'],fit['distances'],s=5,alpha=.15);order=np.argsort(fit['dissimilarities']);axs[0].plot(fit['dissimilarities'][order],fit['disparities'][order],color='#c26437');axs[0].set(xlabel='Original dissimilarity',ylabel='Ordination distance',title='Shepard diagram: monotone fit')
        sts=[f['stress1'] for f in fit['starts']];axs[1].plot(range(len(sts)),sts,'.-');axs[1].set(xlabel='Independent start',ylabel='Stress1',title='NMDS optimization stability');fig.tight_layout();writefig(fig,'nmds_'+name+'_diagnostics')
    ev=a['PCA']['explained_variance_ratio'];fig,axs=plt.subplots(1,2,figsize=(12,5))
    for site in NAMES:
        rr=[r for r in rows if r['site']==site];axs[0].scatter([r['pca'][0] for r in rr],[r['pca'][1] for r in rr],s=24,alpha=.55,color=COLORS[site],label=NAMES[site])
    axs[0].set(xlabel=f'PC1 ({ev[0]:.1%})',ylabel=f'PC2 ({ev[1]:.1%})',title='Shared DINO photo-embedding PCA');axs[0].legend()
    axs[1].plot(np.arange(1,31),np.cumsum(ev[:30]),'.-');axs[1].set(xlabel='Number of PCs',ylabel='Cumulative variance explained',ylim=(0,1),title='Between-photo embedding variance');axs[1].grid(alpha=.15);fig.tight_layout();writefig(fig,'shared_dino_PCA')
    fig,axs=plt.subplots(2,2,figsize=(11,9))
    for ax,site in zip(axs.flat,NAMES):
        rr=[r for r in rows if r['site']==site]
        for month in sorted({r['month'] for r in rr}):
            q=[r for r in rr if r['month']==month];ax.scatter([r['pca'][0] for r in q],[r['pca'][1] for r in q],s=25,alpha=.6,label=month)
        ax.set(title=NAMES[site],xlabel='PC1',ylabel='PC2');ax.legend(fontsize=8)
    fig.suptitle('One shared PCA, campaigns separated within each site');fig.tight_layout();writefig(fig,'shared_dino_PCA_campaigns')
    fig,ax=plt.subplots(figsize=(10,5));al=[s['structural_alpha_Q'] for s in sums];be=[s['structural_beta_Q'] for s in sums]
    ax.bar(x,al,label='Alpha: within-image feature variation',color='#43866e');ax.bar(x,be,bottom=al,label='Beta: between-image centroid variation',color='#d3ad60');ax.set(xticks=x,xticklabels=[NAMES[s['site']] for s in sums],ylabel='Rao / continuous feature variance',title='Structural visual diversity: gamma = alpha + beta');ax.legend();writefig(fig,'structural_alpha_beta_gamma')
    fig,axs=plt.subplots(1,2,figsize=(12,5))
    for ax,field,title in [(axs[0],'adjacent_same_family_fraction','Family aggregation within photos'),(axs[1],'components_per_10000_resolved_cells','Family patch fragmentation')]:
        values=[[r['configuration'][field] for r in rows if r['site']==site and r['configuration'][field] is not None] for site in NAMES];ax.boxplot(values,tick_labels=list(NAMES.values()),showfliers=False);ax.set(title=title,ylabel=field.replace('_',' '))
    fig.tight_layout();writefig(fig,'community_configuration')
    # Export representative photo panels chosen on shared PC positions, not subjective appearance.
    for axis in [0,1]:
        order=sorted(rows,key=lambda r:r['pca'][axis]);chosen=[order[j] for j in np.linspace(0,438,7).round().astype(int)]
        fig,axs=plt.subplots(1,7,figsize=(15,3))
        for ax,r in zip(axs,chosen):ax.imshow(Image.open(OUT/'images'/f"{r['id']}.jpg"));ax.axis('off');ax.set_title(f"{r['id']}\n{NAMES[r['site']]}\nPC{axis+1}={r['pca'][axis]:.3f}",fontsize=9)
        fig.suptitle(f'Representative photographs along shared PC{axis+1}');fig.tight_layout();writefig(fig,f'PCA_photo_gradient_{axis+1}')

if __name__=='__main__':
    try:
        with threadpool_limits(limits=2):analyze()
    except Exception as ex:status('failed',error=str(ex));raise