"""Build the same census profile for ANY set of geographies (a custom region).

Usage (from Python):
    from aggregate import region_profile
    p = region_profile(['2300587845','2300555225'], level='town')     # Yarmouth + North Yarmouth
Or from the command line:
    python3 aggregate.py "My Region" town 2300587845 2300555225

Method (standard Census Bureau practice):
* counts are summed; their margins of error are combined as the square root of the sum of squares;
* percentages are recomputed from the summed numerator and denominator (MOE by the Bureau's ratio formula);
* medians (household income, home value, gross rent, age) are re-estimated from the summed distribution tables
  (B19001, B25075, B25063, B01001) by linear interpolation within the bin that contains the midpoint — the same way
  the Bureau derives medians from grouped data. They are marked 'interpolated' and can differ by a few percent from
  what the Bureau would publish for the same area;
* per-capita and mean income come from aggregate income (B19025) divided by summed population / households;
* population estimates and decennial counts are summed.
Requires values_long.pkl (built by build_long.py) in the same folder or OUT.
"""
import pandas as pd, numpy as np, sys, os
HERE=os.path.dirname(os.path.abspath(__file__))
OUT=os.environ.get('CENSUS_OUT', os.path.join(HERE,'out') if os.path.exists(os.path.join(HERE,'out')) else HERE)
_P=None
def _load():
    global _P
    if _P is None:
        p=os.path.join(OUT,'values_long.pkl')
        if not os.path.exists(p): p=os.path.join(OUT,'census_values_long.csv')
        _P=pd.read_pickle(p) if p.endswith('.pkl') else pd.read_csv(p,dtype={'geoid':str,'year':str})
    return _P
# numerator variables, denominator variable, for every percentage indicator (ACS 2020-2024 data profile / detailed tables)
PCT={
 'pop_under18_pct':(['DP05_0019'],'DP05_0001'),'pop_65plus_pct':(['DP05_0024'],'DP05_0001'),'female_pct':(['DP05_0003'],'DP05_0001'),'male_pct':(['DP05_0002'],'DP05_0001'),
 'veterans_pct':(['DP02_0070'],'DP02_0069'),'disability_pct':(['DP02_0072'],'DP02_0071'),
 'white_nh_pct':(['DP05_0096'],'DP05_0089'),'hispanic_pct':(['DP05_0090'],'DP05_0089'),'black_nh_pct':(['DP05_0097'],'DP05_0089'),'aian_nh_pct':(['DP05_0098'],'DP05_0089'),'asian_nh_pct':(['DP05_0099'],'DP05_0089'),'nhpi_nh_pct':(['DP05_0100'],'DP05_0089'),'other_race_nh_pct':(['DP05_0101'],'DP05_0089'),'two_or_more_nh_pct':(['DP05_0102'],'DP05_0089'),
 'foreign_born_pct':(['DP02_0094'],'DP02_0088'),'born_in_maine_pct':(['DP02_0091'],'DP02_0088'),'language_not_english_pct':(['DP02_0114'],'DP02_0112'),
 'married_couple_hh_pct':(['DP02_0002'],'DP02_0001'),'hh_with_under18_pct':(['DP02_0014'],'DP02_0001'),'hh_with_65plus_pct':(['DP02_0015'],'DP02_0001'),'living_alone_pct':(['DP02_0008','DP02_0012'],'DP02_0001'),
 'broadband_pct':(['DP02_0154'],'DP02_0152'),'computer_pct':(['DP02_0153'],'DP02_0152'),'no_vehicle_pct':(['DP04_0058'],'DP04_0057'),
 'vacant_units_pct':(['DP04_0003'],'DP04_0001'),'owner_occupied_pct':(['DP04_0046'],'DP04_0045'),'renter_occupied_pct':(['DP04_0047'],'DP04_0045'),'single_family_detached_pct':(['DP04_0007'],'DP04_0001'),'mobile_home_pct':(['DP04_0014'],'DP04_0001'),'built_before_1940_pct':(['DP04_0026'],'DP04_0001'),
 'home_value_1m_plus_pct':(['DP04_0088'],'DP04_0080'),'home_value_500k_plus_pct':(['DP04_0087','DP04_0088'],'DP04_0080'),
 'rent_burden_35pct_plus_pct':(['DP04_0142'],'DP04_0136'),'owner_burden_35pct_plus_pct':(['DP04_0115'],'DP04_0110'),
 'heat_fuel_oil_pct':(['DP04_0066'],'DP04_0062'),'heat_wood_pct':(['DP04_0068'],'DP04_0062'),'heat_electricity_pct':(['DP04_0065'],'DP04_0062'),'heat_utility_gas_pct':(['DP04_0063'],'DP04_0062'),
 'moved_in_2020_or_later_pct':(['DP04_0051','DP04_0052'],'DP04_0050'),'moved_last_year_pct':(['DP02_0081'],'DP02_0079'),
 'poverty_rate_pct':(['B17001_002'],'B17001_001'),'child_poverty_rate_pct':([f'B17001_{i:03d}' for i in list(range(4,10))+list(range(18,24))],[f'B17001_{i:03d}' for i in list(range(4,10))+list(range(18,24))+list(range(33,39))+list(range(47,53))]),
 'hh_social_security_pct':(['DP03_0066'],'DP03_0051'),'hh_retirement_income_pct':(['DP03_0068'],'DP03_0051'),'hh_snap_pct':(['DP03_0074'],'DP03_0051'),'hh_cash_assistance_pct':(['DP03_0072'],'DP03_0051'),
 'uninsured_pct':(['DP03_0099'],'DP03_0095'),'labor_force_participation_pct':(['DP03_0002'],'DP03_0001'),'unemployment_rate_pct':(['DP03_0005'],'DP03_0003'),
 'mgmt_professional_pct':(['DP03_0027'],'DP03_0026'),'government_workers_pct':(['DP03_0048'],'DP03_0046'),'self_employed_pct':(['DP03_0049'],'DP03_0046'),'work_from_home_pct':(['DP03_0024'],'DP03_0018'),'drove_alone_pct':(['DP03_0019'],'DP03_0018'),
 'seasonal_homes_pct':(['B25004_006'],'B25001_001'),'hs_or_higher_pct':(['DP02_0067'],'DP02_0059'),'bachelors_or_higher_pct':(['DP02_0068'],'DP02_0059'),'graduate_degree_pct':(['DP02_0066'],'DP02_0059'),'some_college_or_associates_pct':(['DP02_0063','DP02_0064'],'DP02_0059'),
}
COUNTS={'pop_total':'DP05_0001','pop_18plus':'DP05_0021','citizen_voting_age_pop':'DP05_0106','households':'DP02_0001','housing_units':'DP04_0001','occupied_units':'DP04_0002',
        'pop_2020_census':('Decennial Census 2020 (DHC)','P1_001'),'pop_2010_census':('Decennial Census 2010 (SF1)','P001001'),'pop_2000_census':('Decennial Census 2000 (SF1)','P001001'),
        'pop_18plus_2020_census':('Decennial Census 2020 (PL 94-171)','P3_001'),'housing_units_2020_census':('Decennial Census 2020 (DHC)','H1_001')}
BINS={ # (table stem, [(var, lower, upper), ...]) upper=None means open-ended
 'median_household_income':('B19001',[(2,0,10000),(3,10000,15000),(4,15000,20000),(5,20000,25000),(6,25000,30000),(7,30000,35000),(8,35000,40000),(9,40000,45000),(10,45000,50000),(11,50000,60000),(12,60000,75000),(13,75000,100000),(14,100000,125000),(15,125000,150000),(16,150000,200000),(17,200000,None)]),
 'median_home_value':('B25075',[(2,0,10000),(3,10000,15000),(4,15000,20000),(5,20000,25000),(6,25000,30000),(7,30000,35000),(8,35000,40000),(9,40000,50000),(10,50000,60000),(11,60000,70000),(12,70000,80000),(13,80000,90000),(14,90000,100000),(15,100000,125000),(16,125000,150000),(17,150000,175000),(18,175000,200000),(19,200000,250000),(20,250000,300000),(21,300000,400000),(22,400000,500000),(23,500000,750000),(24,750000,1000000),(25,1000000,1500000),(26,1500000,2000000),(27,2000000,None)]),
 'median_gross_rent':('B25063',[(3,0,100),(4,100,150),(5,150,200),(6,200,250),(7,250,300),(8,300,350),(9,350,400),(10,400,450),(11,450,500),(12,500,550),(13,550,600),(14,600,650),(15,650,700),(16,700,750),(17,750,800),(18,800,900),(19,900,1000),(20,1000,1250),(21,1250,1500),(22,1500,2000),(23,2000,2500),(24,2500,3000),(25,3000,3500),(26,3500,None)]),
}
AGE_BINS=[(0,5),(5,10),(10,15),(15,18),(18,20),(20,21),(21,22),(22,25),(25,30),(30,35),(35,40),(40,45),(45,50),(50,55),(55,60),(60,62),(62,65),(65,67),(67,70),(70,75),(75,80),(80,85),(85,None)]
def _median_from_bins(counts,bins):
    total=sum(counts);
    if total<=0: return np.nan
    half=total/2; cum=0
    for c,(lo,hi) in zip(counts,bins):
        if cum+c>=half:
            if hi is None: return float(lo)
            return lo+(hi-lo)*((half-cum)/c if c else 0)
        cum+=c
    return np.nan
def region_profile(geoids, level='town', name='Custom region'):
    P=_load()
    S=P[(P.geo_level==level)&(P.geoid.isin(list(geoids)))]
    if S.geoid.nunique()!=len(set(geoids)):
        missing=set(geoids)-set(S.geoid); raise ValueError(f'not found at level {level}: {sorted(missing)}')
    acs=S[S.dataset.str.startswith('ACS 5-year 2020-2024')]
    E=acs.groupby('variable').estimate.sum(min_count=1); M=np.sqrt((acs.groupby('variable').moe.apply(lambda s:(s**2).sum(min_count=1))))
    out={'name':name,'level':level,'n_units':len(set(geoids)),'geoids':';'.join(sorted(geoids))}
    def cnt(v): return E.get(v,np.nan), M.get(v,np.nan)
    for k,v in COUNTS.items():
        if isinstance(v,tuple):
            d=S[(S.dataset==v[0])&(S.variable==v[1])]; out[k]=d.estimate.sum(min_count=1)
        else:
            out[k],out[k+'_moe']=cnt(v)
    for k,(nums,den) in PCT.items():
        n=sum(E.get(v,np.nan) for v in nums); nm=np.sqrt(sum(M.get(v,np.nan)**2 for v in nums))
        if isinstance(den,list): d=sum(E.get(v,np.nan) for v in den); dm=np.sqrt(sum(M.get(v,np.nan)**2 for v in den))
        else: d,dm=E.get(den,np.nan),M.get(den,np.nan)
        if not d or np.isnan(d): out[k]=np.nan; out[k+'_moe']=np.nan; continue
        p=n/d; inner=nm**2-(p**2)*(dm**2)
        pm=np.sqrt(inner)/d if inner>0 else np.sqrt(nm**2+(p**2)*(dm**2))/d   # Bureau proportion formula, ratio formula if negative
        out[k]=round(100*p,1); out[k+'_moe']=round(100*pm,1)
    for k,(tab,bins) in BINS.items():
        counts=[E.get(f'{tab}_{i:03d}',0) for i,lo,hi in bins]
        out[k]=round(_median_from_bins(counts,[(lo,hi) for i,lo,hi in bins])); out[k+'_method']='interpolated from bins'
    # median age from B01001 (male 003-025, female 027-049 share the same 23 bins)
    counts=[E.get(f'B01001_{i:03d}',0)+E.get(f'B01001_{i+24:03d}',0) for i in range(3,26)]
    out['median_age']=round(_median_from_bins(counts,AGE_BINS),1); out['median_age_method']='interpolated from 5-year bins'
    out['per_capita_income']=round(E.get('B19313_001',np.nan)/E.get('B01003_001',np.nan)) if E.get('B01003_001',0) else np.nan
    out['mean_household_income']=round(E.get('B19025_001',np.nan)/E.get('B11001_001',np.nan)) if E.get('B11001_001',0) else np.nan
    out['avg_household_size']=round(E.get('DP02_0018',np.nan)/E.get('DP02_0001',np.nan),2) if E.get('DP02_0001',0) else np.nan
    return out
if __name__=='__main__':
    name,level,*geos=sys.argv[1:]
    r=region_profile(geos,level,name)
    for k,v in r.items(): print(f'{k:36s} {v}')