import sqlite3, json, re, unicodedata, csv, math, os
from collections import defaultdict, Counter
import numpy as np
from scipy.stats import spearmanr, rankdata
import statsmodels.api as sm

DB='/mnt/data/pk1_full_reanalysis/tenders.db'
NEED='/mnt/data/need_snapshots_authoritative.csv'
PANEL='/mnt/data/pk1_problem_investment_reanalysis/Paneli_Stata_Ready_61_Bashki_2010_2024_v1.csv'
OUT='/mnt/data/pk1_full_reanalysis/final'
os.makedirs(OUT,exist_ok=True)

SECTORS=['Ujë','Sanitacion','Arsim','Rrugë','Ujitje','Kullim']

# INSTAT annual-average CPI growth rates. 2015 only anchors the chain; eligible tests start in 2016.
# Sources: INSTAT CPI December releases / Albania in Figures.
INFL={2016:0.013,2017:0.020,2018:0.020,2019:0.014,2020:0.016,2021:0.020,2022:0.067,2023:0.048,2024:0.022}
CPI={2015:100.0}
for y in range(2016,2025): CPI[y]=CPI[y-1]*(1+INFL[y])
REAL24={y:CPI[2024]/CPI[y] for y in CPI}


def fold(s):
    s='' if s is None else str(s)
    s=unicodedata.normalize('NFKD',s)
    s=''.join(c for c in s if not unicodedata.combining(c))
    s=s.lower().replace('ë','e').replace('ç','c')
    s=re.sub(r'[^a-z0-9]+',' ',s)
    return re.sub(r'\s+',' ',s).strip()

def canon_key(s):
    return fold(s).replace('qeverisja vendore ','').replace('qeversija vendore ','').strip()

def parse_amount(x):
    if x is None: return None
    s=str(x).strip()
    if not s: return None
    s=s.replace('\xa0','').replace(' ','').replace("'",'')
    if ',' in s: s=s.replace('.','').replace(',','.')
    elif re.fullmatch(r'[-+]?\d{1,3}(?:\.\d{3})+',s): s=s.replace('.','')
    s=re.sub(r'[^0-9.\-+]','',s)
    try:
        v=float(s); return v if np.isfinite(v) else None
    except: return None

def get_year(d):
    s=d.get('date_published_iso') or ''
    m=re.match(r'^(20\d{2})',s)
    if m: return int(m.group(1))
    m=re.search(r'(20\d{2})',d.get('data_e_njoftimit_te_tenderit') or '')
    return int(m.group(1)) if m else None

# Canonical municipality / population / area
canonical={}; pop={}; area={}
with open(PANEL,encoding='utf-8-sig',newline='') as f:
    for r in csv.DictReader(f):
        m=r['municipality']; y=int(r['year']); canonical[canon_key(m)]=m
        try: pop[(m,y)]=float(r['pop_est'])
        except: pass
        try: area[m]=float(r['area_km2'])
        except: pass
canon_names=sorted([(fold(m),m) for m in set(canonical.values())],key=lambda z:len(z[0]),reverse=True)

con=sqlite3.connect(DB)
dbmun={mid:canonical.get(canon_key(mn),mn.replace('Qeverisja Vendore ','').replace('Qeversija Vendore ',''))
       for mid,mn in con.execute('select municipality_id,municipality_name from municipalities')}

def bashkia_mentions(text):
    t=' '+fold(text)+' '
    found=[]
    for nf,m in canon_names:
        # exact phrase after Bashkia; handles Bashkisë after folding poorly, so use bashki* stem
        if re.search(r'\bbashki\w*\s+'+re.escape(nf)+r'\b',t): found.append(m)
    return sorted(set(found))

def resolve_municipality(mid,d,title):
    page=dbmun[mid]
    inst=d.get('institucioni_prokurues') or ''
    im=bashkia_mentions(inst)
    if len(im)==1:
        return im[0],'institution_direct'
    it=fold(inst)
    tm=bashkia_mentions(title)
    # Regional utility / county council: title is a better territorial locator when unique.
    if len(tm)==1 and (('rajonal' in it and ('ujesjelles' in it or 'kanaliz' in it)) or 'keshilli i qarkut' in it):
        return tm[0],'title_direct_regional'
    return page,'database_municipality'

def nonphysical_title(title):
    t=fold(title)
    # Design/supervision-only procedures accidentally coded as public works are not physical capital.
    return bool(re.match(r'^(studim\s*projektim|studim projektim|projektim|supervizim|mbikqyrje|kolaudim|oponence|hartim projekti|ekspertize)\b',t))

def classify(title):
    t=' '+fold(title)+' '
    labs=set()
    # sanitation / wastewater
    pats=[r'\bkanalizim\w*(?: te)? ujerave? te zeza\b',r'\bkuz\b',
          r'\bimpiante?\w* (?:i )?trajtimit.*\bujerave? (?:te )?(?:ndotura|zeza)\b',
          r'\bkolektor\w*.*\bujerave? te zeza\b',r'\brrjet\w*.*\bkanalizimeve\b']
    if any(re.search(p,t) for p in pats): labs.add('Sanitacion')
    # irrigation / drainage
    if re.search(r'\b(ujitj\w*|ujites\w*|vadit\w*|skeme\s+ujit\w*|kanal\w*\s+ujites\w*)\b',t): labs.add('Ujitje')
    if re.search(r'\b(kullim\w*|kullues\w*|drenazh\w*)\b',t): labs.add('Kullim')
    # roads / bridges as object
    road_direct=[
      r'\b(?:ndertim|rikonstruksion|rehabilitim|sistemim|mirembajtje|riparim)\w*\s+(?:dhe\s+)?(?:asfaltim\w*\s+)?(?:i\s+|e\s+|te\s+)?rrug\w*',
      r'\bsistemim\w*\s+(?:dhe\s+)?asfaltim\w*',r'\basfaltim\w*\s+(?:i\s+|e\s+|te\s+)?rrug\w*',
      r'\brruga\s+pedonale\b',r'\b(?:ndertim|rikonstruksion|rehabilitim|riparim)\w*\s+(?:i\s+|e\s+|te\s+)?ur(?:e|a|es)\b']
    road_excl=[r'\bndricim\w*.*\brrug',r'\bsinjalistik',r'\bpastrim\w*.*\brrug',r'\bkriposje\b',r'\bmirembajtje dimrore\b']
    if any(re.search(p,t) for p in road_direct) and not any(re.search(p,t) for p in road_excl): labs.add('Rrugë')
    # water supply as object; exclude office works and contextual completion references
    water_direct=[
      r'\b(?:ndertim|rikonstruksion|rehabilitim|riparim|mirembajtje|zgjerim|meremetim)\w*.{0,35}\bujesjelles\w*',
      r'^\s*ujesjelles\w*',r'\bfurnizim\w* me uje\b',r'\brrjet\w* (?:i )?ujit te pijshem\b',
      r'\bdepo(?:zite)? uji\b',r'\bstacion pompimi.{0,30}\buje\b',r'\blinje.{0,20}\bfurnizim.{0,15}\buje\b']
    water_context=bool(re.search(r'\brrug\w*',t) and re.search(r'\b(?:ku|pas|prane).{0,35}\bujesjelles\w*',t))
    office_context=bool(re.search(r'\b(zyr|godin administrat)\w*.*\bujesjelles\w*',t))
    if any(re.search(p,t) for p in water_direct) and not re.search(r'\bujitj\w*|\bvadit\w*',t) and not water_context and not office_context:
        labs.add('Ujë')
    # education facility itself as target; exclude nearby park/road/water references and social repurposing of former creches
    edu_fac=r'(?:shkoll\w*|gjimnaz\w*|kopsht\w*|cerdh\w*|konvikt\w*|objekt\w*\s+arsimor\w*)'
    edu_action=r'(?:ndertim\w*|rikonstruksion\w*|rehabilitim\w*|riparim\w*|mirembajtje\w*|hidroizolim\w*|lyerj\w*|shtese\w*)'
    direct=bool(re.search(r'\b'+edu_action+r'\b(?:\s+\w+){0,3}\s+\b'+edu_fac+r'\b',t) or
                re.search(r'\b'+edu_fac+r'\b(?:\s+\w+){0,3}\s+\b'+edu_action+r'\b',t) or
                re.search(r'\binfrastrukture\s+arsimore\b',t))
    context=bool(re.search(r'\b(?:rrug|uje|ujesjelles|kanal|park|lulishte|linje)\w*\b.{0,35}\b'+edu_fac+r'\b',t))
    social_repurpose=bool(re.search(r'\bcerdh\w*.{0,40}\b(strehim|dhunuar|social)\w*',t))
    if direct and not context and not social_repurpose: labs.add('Arsim')
    return labs

# Procurement extraction. Coverage is based on source-page municipality and ALL tenders, not sector outcomes.
coverage=defaultdict(set); records=[]; loc_counts=Counter(); physical_counts=Counter()
for ref,mid,s in con.execute('select ti.ref,ti.municipality_id,td.data_json from tenders_index ti join tenders_data td using(ref)'):
    d=json.loads(s); y=get_year(d)
    if not y or y<2015 or y>2024: continue
    page=dbmun[mid]; coverage[y].add(page)
    ctype=fold(d.get('lloji_i_kontrates_publike') or '')
    if not ctype.startswith('pune publike'): continue
    title=d.get('objekti_i_tenderit') or d.get('titulli') or ''
    physical=not nonphysical_title(title); physical_counts[physical]+=1
    if not physical: continue
    m,loc=resolve_municipality(mid,d,title); loc_counts[loc]+=1
    labs=classify(title)
    status=fold(d.get('statusi_i_tenderit') or '')
    aw=parse_amount(d.get('oferta_fituese_leke_pa_tvsh'))
    lim=parse_amount(d.get('vlera_fondi_limit'))
    if 'anuluar' in status: aw=None
    award=aw if aw is not None and aw>0 else None
    award_real=award*REAL24[y] if award is not None else None
    limit_aw=lim*REAL24[y] if (award is not None and lim is not None and lim>0) else None
    records.append(dict(ref=ref,municipality=m,page_municipality=page,location_method=loc,year=y,title=title,labs=labs,
                        award_real24=award_real,award_nominal=award,award_limit_real24=limit_aw,
                        status=d.get('statusi_i_tenderit') or '',institution=d.get('institucioni_prokurues') or '',
                        authority=d.get('autoritet_prokurues') or '',url=d.get('url') or ''))

# Need data
needs=[]
with open(NEED,encoding='utf-8-sig',newline='') as f:
    for r in csv.DictReader(f):
        m=canonical.get(canon_key(r['Bashkia']),r['Bashkia'])
        try: v=float(r['Vlera']); sy=int(float(r['Viti matjes'])); first=int(float(r['Viti i parë invest.']))
        except: continue
        r.update(municipality=m,need_value=v,snapshot_year=sy,first_invest=first); needs.append(r)
sector_map={'Ujë':'Ujë','Sanitacion':'Sanitacion','Arsim':'Arsim','Rrugë':'Rrugë','Ujitje':'Ujitje','Sipërfaqe e ujitur':'Ujitje','Kullim':'Kullim'}

# Aggregation by municipality-year. Broad counts explicit multi-sector projects for every sector they directly address.
# Strict counts a project only when exactly one target sector is identified.
agg=defaultdict(lambda: {'all':0.0,'all_limit':0.0,'n_awards':0})
for r in records:
    if r['award_real24'] is None: continue
    k=(r['municipality'],r['year']); agg[k]['all']+=r['award_real24']; agg[k]['n_awards']+=1
    if r['award_limit_real24'] is not None: agg[k]['all_limit']+=r['award_limit_real24']
    for lab in r['labs']:
        agg[k]['broad_'+lab]=agg[k].get('broad_'+lab,0.0)+r['award_real24']
        if r['award_limit_real24'] is not None: agg[k]['limit_'+lab]=agg[k].get('limit_'+lab,0.0)+r['award_limit_real24']
    if len(r['labs'])==1:
        lab=next(iter(r['labs'])); agg[k]['strict_'+lab]=agg[k].get('strict_'+lab,0.0)+r['award_real24']

# Analysis windows. Old 2011 needs are a delayed procurement validation because the DB becomes usable only in 2016.
def windows(first):
    if first<2016:
        return [('proc_2016',2016,2016),('proc_2016_18',2016,2018),('proc_2016_20',2016,2020),('proc_2016_24',2016,2024)]
    a=first; out=[]
    if a<=2024: out.append(('1y',a,a))
    if a+2<=2024: out.append(('3y',a,a+2))
    if a+4<=2024: out.append(('5y',a,a+4))
    if a<=2024 and (not out or out[-1][2]!=2024): out.append(('full',a,2024))
    return out

def spear(x,y):
    x=np.asarray(x,float); y=np.asarray(y,float)
    if len(x)<3 or np.nanstd(x)==0 or np.nanstd(y)==0: return np.nan,np.nan
    z=spearmanr(x,y); return float(z.statistic),float(z.pvalue)

def fast_boot_delta(x,y1,y2,B=3000,seed=20260822):
    x=np.asarray(x,float); y1=np.asarray(y1,float); y2=np.asarray(y2,float); n=len(x)
    if n<5 or np.std(y1)==0 or np.std(y2)==0: return (np.nan,)*4
    rx=rankdata(x); r1=rankdata(y1); r2=rankdata(y2)
    rng=np.random.default_rng(seed); idx=rng.integers(0,n,size=(B,n))
    def corr_rows(a,b):
        aa=a[idx]; bb=b[idx]; aa=aa-aa.mean(1,keepdims=True); bb=bb-bb.mean(1,keepdims=True)
        den=np.sqrt((aa*aa).sum(1)*(bb*bb).sum(1)); num=(aa*bb).sum(1)
        return np.divide(num,den,out=np.full(B,np.nan),where=den>0)
    d=corr_rows(rx,r1)-corr_rows(rx,r2); d=d[np.isfinite(d)]
    if not len(d): return (np.nan,)*4
    return float(np.median(d)),float(np.quantile(d,.025)),float(np.quantile(d,.975)),float(np.mean(d<=0))

def hc3(y,burden,pv,av):
    if len(y)<10: return (np.nan,)*4
    X=sm.add_constant(np.column_stack([np.log1p(burden),np.log(pv),np.log(av)]))
    try:
        fit=sm.OLS(np.log1p(y),X).fit(cov_type='HC3')
        return float(fit.params[1]),float(fit.bse[1]),float(fit.pvalues[1]),float(fit.rsquared)
    except: return (np.nan,)*4

def bh_on_indices(rows,indices,pfield,qfield):
    valid=[i for i in indices if np.isfinite(rows[i].get(pfield,np.nan))]
    if not valid: return
    ps=np.array([rows[i][pfield] for i in valid]); order=np.argsort(ps); m=len(ps)
    qord=np.array([ps[j]*m/(rank+1) for rank,j in enumerate(order)],float)
    for j in range(m-2,-1,-1): qord[j]=min(qord[j],qord[j+1])
    q=np.empty(m); q[order]=np.minimum(qord,1)
    for i,qq in zip(valid,q): rows[i][qfield]=float(qq)

def bh_group(rows,pfield,qfield):
    groups=defaultdict(list)
    for i,r in enumerate(rows): groups[(r['family'],r['dimension'],r['role'])].append(i)
    for inds in groups.values(): bh_on_indices(rows,inds,pfield,qfield)

by=defaultdict(list)
for r in needs:
    if r['Familja'] in sector_map and r.get('Roli H1') in ('PO','SENS'):
        by[(r['Familja'],r['Treguesi'],r['Dimensioni'],r['snapshot_year'],r['first_invest'],r['Roli H1'])].append(r)

results=[]
for (fam,indicator,dim,sy,first,role),rows in sorted(by.items()):
    sec=sector_map[fam]; need={r['municipality']:r['need_value'] for r in rows}
    for wl,a,b in windows(first):
        ys=list(range(a,b+1)); obs=set.intersection(*(coverage[y] for y in ys))
        mun=sorted(set(need)&obs)
        if len(mun)<31: continue
        x=[]; broad=[]; strict=[]; total=[]; limit=[]; pv=[]; av=[]
        for m in mun:
            p=pop.get((m,sy)) or pop.get((m,a)); ar=area.get(m)
            if not p or not ar: continue
            bv=sum(agg[(m,y)].get('broad_'+sec,0.0) for y in ys)
            sv=sum(agg[(m,y)].get('strict_'+sec,0.0) for y in ys)
            tv=sum(agg[(m,y)].get('all',0.0) for y in ys)
            lv=sum(agg[(m,y)].get('limit_'+sec,0.0) for y in ys)
            x.append(need[m]); broad.append(bv); strict.append(sv); total.append(tv); limit.append(lv); pv.append(p); av.append(ar)
        x=np.array(x); broad=np.array(broad); strict=np.array(strict); total=np.array(total); limit=np.array(limit); pv=np.array(pv); av=np.array(av)
        absdim='Barra absolute' in dim
        ym=broad if absdim else broad/pv
        ys_strict=strict if absdim else strict/pv
        # Explicit target-containing projects are removed from 'other', including mixed target projects.
        yo=np.maximum(total-broad,0) if absdim else np.maximum(total-broad,0)/pv
        rho,p=spear(x,ym); rhoo,po=spear(x,yo); rhos,ps=spear(x,ys_strict)
        delta=rho-rhoo if np.isfinite(rho) and np.isfinite(rhoo) else np.nan
        bmed,blo,bhi,bp=fast_boot_delta(x,ym,yo)
        share=np.divide(broad,total,out=np.zeros_like(broad),where=total>0)
        rshare,pshare=spear(x,share)
        anym=(broad>0).astype(float); rany,pany=spear(x,anym)
        ylim=limit if absdim else limit/pv; rlim,plim=spear(x,ylim)
        # strict specificity against works with no target component
        delta_strict=rhos-rhoo if np.isfinite(rhos) and np.isfinite(rhoo) else np.nan
        smed,slo,shi,sp=fast_boot_delta(x,ys_strict,yo)
        beta=se=pscale=r2=np.nan
        if absdim: beta,se,pscale,r2=hc3(broad,x,pv,av)
        results.append(dict(family=fam,indicator=indicator,dimension=dim,role=role,snapshot_year=sy,sector=sec,
                            window=wl,start=a,end=b,N=len(x),coverage_municipalities=len(obs),positive_match=int((broad>0).sum()),
                            rho_match=rho,p_match=p,rho_other=rhoo,p_other=po,delta_rho=delta,
                            boot_delta_median=bmed,boot_delta_lo95=blo,boot_delta_hi95=bhi,boot_p_delta_le0=bp,
                            rho_share=rshare,p_share=pshare,rho_any=rany,p_any=pany,
                            rho_strict_match=rhos,p_strict_match=ps,delta_strict=delta_strict,
                            strict_boot_lo95=slo,strict_boot_hi95=shi,strict_boot_p_le0=sp,
                            rho_awarded_limit=rlim,p_awarded_limit=plim,
                            beta_burden_scale_adj=beta,se_burden_scale_adj=se,p_burden_scale_adj=pscale,r2_scale_adj=r2,
                            delayed_baseline=(first<2016)))

for pf,qf in [('p_match','q_match'),('p_share','q_share'),('p_any','q_any'),('p_strict_match','q_strict_match'),('p_burden_scale_adj','q_scale_adj')]: bh_group(results,pf,qf)

# outputs
allfields=[]
for r in results:
    for k in r:
        if k not in allfields: allfields.append(k)
with open(OUT+'/tests_all.csv','w',encoding='utf-8-sig',newline='') as f:
    w=csv.DictWriter(f,fieldnames=allfields); w.writeheader(); w.writerows(results)

# Tender-level audit table, all physical public works, preserving source ref and classification
with open(OUT+'/tender_classification_all.csv','w',encoding='utf-8-sig',newline='') as f:
    ff=['ref','year','page_municipality','municipality','location_method','sector_labels','award_nominal','award_real24','status','title','institution','authority','url']
    w=csv.DictWriter(f,fieldnames=ff); w.writeheader()
    for r in records:
        w.writerow(dict(ref=r['ref'],year=r['year'],page_municipality=r['page_municipality'],municipality=r['municipality'],location_method=r['location_method'],
                        sector_labels='|'.join(sorted(r['labs'])),award_nominal=r['award_nominal'],award_real24=r['award_real24'],status=r['status'],title=r['title'],institution=r['institution'],authority=r['authority'],url=r['url']))

# coverage / portfolio
portfolio=[]
for y in range(2015,2025):
    yearr=[r for r in records if r['year']==y]; aw=[r for r in yearr if r['award_real24'] is not None]
    for sec in SECTORS:
        ss=[r for r in aw if sec in r['labs']]
        portfolio.append(dict(year=y,municipalities_observed=len(coverage[y]),physical_public_works=len(yearr),awarded_public_works=len(aw),sector=sec,
                              sector_awards=len(ss),sector_positive_municipalities=len(set(r['municipality'] for r in ss)),sector_real24_all=sum(r['award_real24'] for r in ss)))
with open(OUT+'/portfolio_coverage.csv','w',encoding='utf-8-sig',newline='') as f:
    w=csv.DictWriter(f,fieldnames=portfolio[0].keys()); w.writeheader(); w.writerows(portfolio)

# compact results chosen for thesis: all primary plus a decision-oriented summary
summary=[]
for fam in sorted(set(r['family'] for r in results)):
    rr=[r for r in results if r['family']==fam and r['role']=='PO']
    if not rr: continue
    summary.append(dict(family=fam,tests=len(rr),positive_matched=sum(np.isfinite(r['rho_match']) and r['rho_match']>0 for r in rr),
                        matched_q05=sum(r.get('q_match',1)<.05 for r in rr),matched_q10=sum(r.get('q_match',1)<.10 for r in rr),
                        specificity_positive_ci=sum(np.isfinite(r['boot_delta_lo95']) and r['boot_delta_lo95']>0 for r in rr),
                        specificity_negative_ci=sum(np.isfinite(r['boot_delta_hi95']) and r['boot_delta_hi95']<0 for r in rr),
                        positive_share_q05=sum(r.get('q_share',1)<.05 and r['rho_share']>0 for r in rr),
                        any_positive_q05=sum(r.get('q_any',1)<.05 and r['rho_any']>0 for r in rr),
                        scale_adj_positive_p05=sum(np.isfinite(r['p_burden_scale_adj']) and r['p_burden_scale_adj']<.05 and r['beta_burden_scale_adj']>0 for r in rr)))
with open(OUT+'/family_summary.csv','w',encoding='utf-8-sig',newline='') as f:
    w=csv.DictWriter(f,fieldnames=summary[0].keys()); w.writeheader(); w.writerows(summary)

# source / method metadata
with open(OUT+'/METHOD_README.txt','w',encoding='utf-8') as f:
    f.write('''PK1/H1 granular procurement validation - FINAL\n\nUniverse: tenders.db, 61 local-government procurement pages.\nYears: 2015-2024 loaded; inferential windows require >50% of municipalities and complete municipality coverage in every year of the window.\nMissing municipality-year is never coded as zero. A zero sector award is used only when the municipality is observed in all years of the window but no awarded physical public-work procedure matches the target sector.\nMain amount: winning offer excluding VAT, restricted to non-cancelled procedures, deflated to 2024 prices using INSTAT annual-average CPI growth.\nRobustness amount: fund limit, only for procedures that reached an award and have a winning offer.\nPublic works: physical works only; design/supervision-only titles are removed.\nSector classification: conservative title-based rules; a project may have multiple labels when it explicitly contains multiple sector components. Strict robustness retains only single-label target projects.\nTerritory: database municipality except when contracting institution identifies another municipality directly, or a regional utility/county council tender title uniquely identifies the beneficiary municipality.\nNeed: authoritative Need_Snapshots_Final from PK1_H1_Analiza_PERFUNDIMTARE_AUDITUAR.xlsx.\nRelative/rank need is compared with matched award per snapshot-year population; absolute burden with absolute matched award.\nSpecificity: rho(problem, matched sector) - rho(problem, public works without the target sector), with 3,000 municipality bootstraps.\nScale: log1p(matched award) on log1p(absolute burden), ln(population), ln(area), HC3.\nBH-FDR: within family x need dimension x H1 role, separately for each outcome family.\n2011 water/sanitation procurement tests are explicitly delayed validation starting in 2016; they do not recover missing 2012-2015 procurement and are not treated as immediate t+1 tests.\n''')

print('CPI/real24', {y:round(REAL24[y],4) for y in range(2016,2025)})
print('location methods',loc_counts)
print('physical public works',len(records),'awarded',sum(r['award_real24'] is not None for r in records))
print('coverage', {y:len(coverage[y]) for y in range(2015,2025)})
print('sector portfolio')
for sec in SECTORS:
    rr=[r for r in records if r['award_real24'] is not None and sec in r['labs']]
    print(sec,len(rr),len(set(r['municipality'] for r in rr)),round(sum(r['award_real24'] for r in rr)/1e9,2),'bn real24')
print('family summary')
for r in summary: print(r)
print('selected significant/suggestive matched or specificity')
for r in results:
    if r['role']!='PO': continue
    if r.get('q_match',1)<.10 or (np.isfinite(r['boot_delta_lo95']) and r['boot_delta_lo95']>0) or (np.isfinite(r['boot_delta_hi95']) and r['boot_delta_hi95']<0):
        print(r['family'],r['snapshot_year'],r['indicator'],r['dimension'],r['window'],r['N'],
              'rhoM',round(r['rho_match'],3),'q',round(r.get('q_match',np.nan),3),'rhoO',round(r['rho_other'],3),'d',round(r['delta_rho'],3),
              'CI',round(r['boot_delta_lo95'],3),round(r['boot_delta_hi95'],3),'share',round(r['rho_share'],3),'qS',round(r.get('q_share',np.nan),3),
              'any q',round(r.get('q_any',np.nan),3),'scale',round(r['beta_burden_scale_adj'],3) if np.isfinite(r['beta_burden_scale_adj']) else None,round(r.get('q_scale_adj',np.nan),3) if np.isfinite(r.get('q_scale_adj',np.nan)) else None)