cd ~

# Import standard packages.
import collections
import functools
import os
import pdb # Debug with pdb.
import subprocess
import sys
import time
# Import installed packages.
import numpy as np
import pandas as pd
# Import local packages.
# Insert current directory into module search path.
# Autoreload local packages after editing.
# `dsdemos` version: https://github.com/stharrold/dsdemos/releases/tag/v0.0.3
sys.path.insert(0, os.path.join(os.path.curdir, r'dsdemos'))
%reload_ext autoreload
%autoreload 2
import dsdemos as dsd

print("Timestamp:")
print(time.strftime(r'%Y-%m-%dT%H:%M:%S%Z', time.gmtime()))
print()
print("Versions:")
print("Python:", sys.version_info)
print("numpy:", np.__version__)
print("pandas:", pd.__version__)

"""
Explanation: 20160110-etl-census-with-python
Related post:
https://stharrold.github.io/20160110-etl-census-with-python.html
Data documentation:
https://www.census.gov/programs-surveys/acs/technical-documentation/pums/documentation.2013.html
Initialization
Imports
End of explanation
"""


# File paths
path_static = os.path.join(os.path.expanduser(r'~'), r'stharrold.github.io/content/static')
basename = r'20160110-etl-census-with-python'
filename = basename
path_ipynb = os.path.join(path_static, basename, filename+'.ipynb')
path_disk = os.path.abspath(r'/mnt/disk-20151227t211000z/')
path_acs = os.path.join(path_disk, r'www2-census-gov/programs-surveys/acs/')
path_pcsv = os.path.join(path_acs, r'data/pums/2013/5-Year/ss13pdc.csv') # 'pdc' = 'person DC'
path_hcsv = os.path.join(path_acs, r'data/pums/2013/5-Year/ss13hdc.csv') # 'hdc' = 'housing DC'
path_ecsv = os.path.join(path_acs, r'tech_docs/pums/estimates/pums_estimates_9_13.csv')
path_dtxt = os.path.join(path_acs, r'tech_docs/pums/data_dict/PUMS_Data_Dictionary_2009-2013.txt')

# Weights
pwt = 'PWGTP' # person weight
pwts = [pwt+str(inum) for inum in range(1, 81)]
hwt = 'WGTP' # housing weight
hwts = [hwt+str(inum) for inum in range(1, 81)]

"""
Explanation: Globals
File sources:
* 2013 5-year PUMS data dictionary: PUMS_Data_Dictionary_2009-2013.txt (<1&nbsp;MB)
* 2013 5-year PUMS person and housing records for Washington DC:
    * Person records: csv_pdc.zip (5&nbsp;MB compressed, 30&nbsp;MB decompressed)
    * Housing records: csv_hdc.zip (2&nbsp;MB compressed, 13&nbsp;MB decompressed)
* 2013 5-year PUMS estimates for user verification: pums_estimates_9_13.csv (<1&nbsp;MB)
End of explanation
"""


print("`ddict`: Load the data dictionary and display the hierarchical structure.")
# Only `ddict` is used below.
# The hierarchical data frame is only for display. 
ddict = dsd.census.parse_pumsdatadict(path=path_dtxt)
tmp = dict()
for record_type in ddict['record_types']:
    tmp[record_type] = pd.DataFrame.from_dict(ddict['record_types'][record_type], orient='index')
pd.concat(tmp, names=['record_type', 'var_name']).head()

print("`ddict`: First 10 unstructured notes from end of file.")
ddict['notes'][:10]

"""
Explanation: Extract, transform, and load
Data dictionary
End of explanation
"""


print("`dfp`, `dfh`: Load person and housing records.")
time_start = time.perf_counter()
for path in [path_pcsv, path_hcsv]:
    with open(path) as fobj:
        nlines = sum(1 for _ in fobj)
    with open(path) as fobj:
        first_line = fobj.readline()
        ncols = first_line.count(',')+1
    print("{path}:".format(path=path))
    print("    size (MB)   = {size:.1f}".format(size=os.path.getsize(path)/1e6))
    print("    num lines   = {nlines}".format(nlines=nlines))
    print("    num columns = {ncols}".format(ncols=ncols))
print()

# For ss13pdc.csv, low_memory=False since otherwise pandas raises DtypeWarning.
dfp = pd.read_csv(path_pcsv, low_memory=False)
dfh = pd.read_csv(path_hcsv, low_memory=True)
for (name, df) in [('dfp', dfp), ('dfh', dfh)]:
    print("{name} RAM usage (MB) = {mem:.1f}".format(
            name=name, mem=df.memory_usage().sum()/1e6))
time_stop = time.perf_counter()
print()
print("Time elapsed (sec) = {diff:.1f}".format(diff=time_stop-time_start))

print("`dfp`: First 5 person records.")
dfp.head()

print("`dfp`: First 5 housing records.")
dfh.head()

print(
r"""`dfp`, `dfh`, `ddict`: Describe all columns ('variables') that aren't weights or flags.
Printed format:
[PERSON, HOUSING] RECORD
COL: Column name.
    Column description.
    Multi-line optional column notes.
    1-3 line description of value meanings ('variable codes').
    Multi-line statistical description and data type.
...
num columns described = ncols""")
print()
records_dfs = collections.OrderedDict([
    ('PERSON RECORD',  {'dataframe': dfp, 'weight': pwt, 'replicate_weights': pwts}),
    ('HOUSING RECORD', {'dataframe': dfh, 'weight': hwt, 'replicate_weights': hwts})])
for record_type in records_dfs:
    print(record_type)
    df = records_dfs[record_type]['dataframe']
    ncols_desc = 0 # number of columns described
    for col in df.columns:
        if col in ddict['record_types'][record_type]:
            col_dict = ddict['record_types'][record_type][col]
            desc = col_dict['description']
        else:
            col_dict = None
            desc = 'Column not in data dictionary.'
        if not (
            (col.startswith('F') and (desc.endswith(' flag') or desc.endswith(' edit')))
            or ('WGTP' in col and "Weight replicate" in desc)):
                print("{col}: {desc}".format(col=col, desc=desc))
                ncols_desc += 1
                if col_dict is not None:
                    if 'notes' in col_dict:
                        print("    {notes}".format(notes=col_dict['notes']))
                    for (inum, var_code) in enumerate(col_dict['var_codes']):
                        var_code_desc = col_dict['var_codes'][var_code]
                        print("    {vc}: {vcd}".format(vc=var_code, vcd=var_code_desc))
                        if inum >= 2:
                            print("    ...")
                            break
                print('    '+repr(df[col].describe()).replace('\n', '\n    '))
    print("num columns described = {ncd}".format(ncd=ncols_desc))
    print()

"""
Explanation: PUMS data
End of explanation
"""


print("`dfe`: Estimates for user verification filtered for 'District of Columbia'.")
dfe = pd.read_csv(path_ecsv)
tfmask_dc = dfe['state'] == 'District of Columbia'
dfe_dc = dfe.loc[tfmask_dc]
dfe_dc

print("`dfe`: Verify characteristic estimates, direct standard errors, and margin of error.")
# Verify the estimates following
# https://www.census.gov/programs-surveys/acs/
#     technical-documentation/pums/documentation.2013.html
#     tech_docs/pums/accuracy/2009_2013AccuracyPUMS.pdf
print()
tfmask_test_strs = collections.OrderedDict([
    ('PERSON RECORD', collections.OrderedDict([
        ('Total population', "np.asarray([True]*len(dfp))"),
        ('Housing unit population (RELP=0-15)',"np.logical_and(0 <= dfp['RELP'], dfp['RELP'] <= 15)"),
        ('GQ population (RELP=16-17)', "np.logical_and(16 <= dfp['RELP'], dfp['RELP'] <= 17)"),
        ('GQ institutional population (RELP=16)', "dfp['RELP'] == 16"),
        ('GQ noninstitutional population (RELP=17)', "dfp['RELP'] == 17"),
        ('Total males (SEX=1)', "dfp['SEX'] == 1"),
        ('Total females (SEX=2)', "dfp['SEX'] == 2"),
        ('Age 0-4', "np.logical_and(0 <= dfp['AGEP'], dfp['AGEP'] <= 4)"),
        ('Age 5-9', "np.logical_and(5 <= dfp['AGEP'], dfp['AGEP'] <= 9)"),
        ('Age 10-14', "np.logical_and(10 <= dfp['AGEP'], dfp['AGEP'] <= 14)"),
        ('Age 15-19', "np.logical_and(15 <= dfp['AGEP'], dfp['AGEP'] <= 19)"),
        ('Age 20-24', "np.logical_and(20 <= dfp['AGEP'], dfp['AGEP'] <= 24)"),
        ('Age 25-34', "np.logical_and(25 <= dfp['AGEP'], dfp['AGEP'] <= 34)"),
        ('Age 35-44', "np.logical_and(35 <= dfp['AGEP'], dfp['AGEP'] <= 44)"),
        ('Age 45-54', "np.logical_and(45 <= dfp['AGEP'], dfp['AGEP'] <= 54)"),
        ('Age 55-59', "np.logical_and(55 <= dfp['AGEP'], dfp['AGEP'] <= 59)"),
        ('Age 60-64', "np.logical_and(60 <= dfp['AGEP'], dfp['AGEP'] <= 64)"),
        ('Age 65-74', "np.logical_and(65 <= dfp['AGEP'], dfp['AGEP'] <= 74)"),
        ('Age 75-84', "np.logical_and(75 <= dfp['AGEP'], dfp['AGEP'] <= 84)"),
        ('Age 85 and over', "85 <= dfp['AGEP']")])),
    ('HOUSING RECORD', collections.OrderedDict([
        ('Total housing units (TYPE=1)', "dfh['TYPE'] == 1"),
        ('Total occupied units', "dfh['TEN'].notnull()"),
        ('Owner occupied units (TEN in 1,2)', "np.logical_or(dfh['TEN'] == 1, dfh['TEN'] == 2)"),
        ('Renter occupied units (TEN in 3,4)', "np.logical_or(dfh['TEN'] == 3, dfh['TEN'] == 4)"),
        ('Owned with a mortgage (TEN=1)', "dfh['TEN'] == 1"),
        ('Owned free and clear (TEN=2)', "dfh['TEN'] == 2"),
        ('Rented for cash (TEN=3)', "dfh['TEN'] == 3"),
        ('No cash rent (TEN=4)', "dfh['TEN'] == 4"),
        ('Total vacant units', "dfh['TEN'].isnull()"),
        ('For rent (VACS=1)', "dfh['VACS'] == 1"),
        ('For sale only (VACS=3)', "dfh['VACS'] == 3"),
        ('All Other Vacant (VACS in 2,4,5,6,7)',
         "functools.reduce(np.logical_or, (dfh['VACS'] == vacs for vacs in [2,4,5,6,7]))")]))])
for record_type in records_dfs:
    print("'{rt}'".format(rt=record_type))
    df = records_dfs[record_type]['dataframe']
    wt = records_dfs[record_type]['weight']
    wts = records_dfs[record_type]['replicate_weights']
    for char in tfmask_test_strs[record_type]:
        print("    '{char}'".format(char=char))
        # Select the reference verification data
        # and the records for the characteristic.
        tfmask_ref = dfe_dc['characteristic'] == char        
        tfmask_test = eval(tfmask_test_strs[record_type][char])
        # Calculate and verify the estimate ('est') for the characteristic.
        # The estimate is the sum of the sample weights 'WGTP'.
        col = 'pums_est_09_to_13'
        print("        '{col}':".format(col=col), end=' ')
        ref_est = int(dfe_dc.loc[tfmask_ref, col].values[0].replace(',', ''))
        test_est = df.loc[tfmask_test, wt].sum()
        assert np.isclose(ref_est, test_est, rtol=0, atol=1)
        print("(ref, test) = {tup}".format(tup=(ref_est, test_est)))
        # Calculate and verify the "direct standard error" ('se') of the estimate.
        # The direct standard error is a modified root-mean-square deviation
        # using the "replicate weights" 'WGTP[1-80]'.
        col = 'pums_se_09_to_13'
        print("        '{col}' :".format(col=col), end=' ')
        ref_se = dfe_dc.loc[tfmask_ref, col].values[0]
        test_se = ((4/80)*((df.loc[tfmask_test, wts].sum() - test_est)**2).sum())**0.5
        assert np.isclose(ref_se, test_se, rtol=0, atol=1)
        print("(ref, test) = {tup}".format(tup=(ref_se, test_se)))
        # Calculate and verify the margin of error ('moe') at the
        # 90% confidence level (+/- 1.645 standard errors).
        col = 'pums_moe_09_to_13'
        print("        '{col}':".format(col=col), end=' ')
        ref_moe = dfe_dc.loc[tfmask_ref, col].values[0]
        test_moe = 1.645*test_se
        assert np.isclose(ref_moe, test_moe, rtol=0, atol=1)
        print("(ref, test) = {tup}".format(tup=(ref_moe, test_moe)))

"""
Explanation: PUMS estimates for user verification
End of explanation
"""


# Export ipynb to html
for template in ['basic', 'full']:
    path_html = os.path.splitext(path_ipynb)[0]+'-'+template+'.html'
    cmd = ['jupyter', 'nbconvert', '--to', 'html', '--template', template, path_ipynb, '--output', path_html]
    print(' '.join(cmd))
    subprocess.run(args=cmd, check=True)
    print()

"""
Explanation: Export ipynb to html
End of explanation
"""


