preprocess.py

#global variables 

import copy
import numpy as np
import pandas as pd

import matplotlib.pyplot as plt

comorbid_map = \
{ \
    'Heart failure': ['389.91', '402.11', '402.91', '404.11', '404.13', '404.91', '404.93', 'I09.0', 'I11.0', 'I13.0', 'I13.2', \
        'I42.5', 'I42.6', 'I42.7', 'I42.8', 'I42.9', 'P29.0', 'I43', 'I43.0', 'I43.1', 'I43.2', 'I43.8', 'I50', 'I50.1', 'I50.2', \
        'I50.20', 'I50.21', 'I50.22', 'I50.23', 'I50.3', 'I50.30', 'I50.31', 'I50.32', 'I50.33', 'I50.4', 'I50.40', 'I50.41', 'I50.42', \
        'I50.43', 'I50.8', 'I50.81', 'I50.82', 'I50.83', 'I50.84', 'I50.89', 'I50.9', 'I50.810', 'I50.811', 'I50.812', 'I50.813', \
        'I50.814', 'I25.5', 'I42.0', '428', '428.0', '428.1', '428.2', '428.3', '428.4', '428.9', '428.20', '428.21', '428.22', \
        '428.23', '428.30', '428.31', '428.32', '428.33', '428.40', '428.41', '428.42', '428.43'], \
    'Diabetes': ['250', '250.0', '250.00', '250.01', '250.02', '250.03', '250.1', '250.10', '250.11', '250.12', '250.13', \
        '250.2', '250.20', '250.21', '250.22', '250.23', '250.3', '250.30', '250.31', '250.32', '250.33', '250.4', '250.40', '250.41', \
        '250.42', '250.43', '250.5', '250.50', '250.51', '250.52', '250.53', '250.6', '250.60', '250.61', '250.62', '250.63', \
        '250.7', '250.70', '250.71', '250.72', '250.73', '250.9', '250.90', '250.91', '250.92', '250.93', 'E10.0', 'E10.1', 'E10.9', 
        'E11.0', 'E11.1', 'E11.9', 'E12.0', 'E12.1', 'E12.9', 'E13.0', 'E13.1', 'E13.9', 'E14.0', 'E14.1', 'E14.9', \
        'E10.2', 'E10.3', 'E10.4', 'E10.5', 'E10.6', 'E10.7', 'E10.8', 'E11.2', 'E11.3', 'E11.4', 'E11.5', 'E11.6', 'E11.7', 'E11.8', \
        'E12.2', 'E12.3', 'E12.4', 'E12.5', 'E12.6', 'E12.7', 'E12.8', 'E13.2', 'E13.3', 'E13.4', 'E13.5', 'E13.6', 'E13.7', 'E13.8', \
        'E14.2', 'E14.3', 'E14.4', 'E14.5', 'E14.6', 'E14.7', 'E14.8'], \
    'Hypertension': ['401', '642.0', 'I10', '402', '403', '404', '405', '642.1', '642.2', '642.7', '642.9', 'I11', 'I12', 'I13', 'I15'], \
    'Peripheral vascular\ndisorders': ['440', '441', '442', '443.1', '443.2', '443.3', '443.4', '443.5', '443.6', '443.7', '443.8', '443.9', \
        '447.1', '557.1', '557.9', 'V43.4', 'I70', 'I71', 'I73.1', 'I73.8', 'I73.9', 'I77.1', 'I79.0', 'I79.2', 'K55.1', 'K55.8', 'K55.9', \
        'Z95.8', 'Z95.9'], \
    'Pulmonary circularion disorders': ['416', '417.9', 'I26', 'I27', 'I28.0', 'I28.8', 'I28.9', '415.0', '415.1', '417.0', '417.8'], \
    'Coagulopathy': ['286', '287.1', '287.3', '287.5', 'D65', 'D66', 'D77', 'D68', 'D69.1', 'D69.3', 'D69.4', 'D69.5', 'D69.6'], \
    'Obesity': ['278.0', 'E66'], \
    'Weight loss': ['260', '261', '262', '263', '783.2', 'E40', 'E41', 'E42', 'E43', 'E44', 'E45', 'E46', 'R63.4', 'R64', '799.4'], \
    'Drug abuse': ['202.0', '292.82', '292.83', '292.84', '292.85', '292.86', '292.87', '292.88', '292.89', '292.9', '304', '305.2', \
        '305.3', '305.4', '305.5', '305.6', '305.7', '305.8', '305.9', '648.3', 'F11', 'F12', 'F13', 'F14', 'F15', 'F16', 'F18', 'F19',
        'Z71.5', 'Z72.2', '292', 'V65.42'], \
    'Hearing loss': ['389', 'H90', 'H91'], \
    'Renal failure': ['403.11', '403.91', '404.12', '404.92', 'V42.0', 'V45.1', 'V56.0', 'V56.8', 'I12.0', 'I13.1', 'N25.0', 'Z94.0' \
        'Z99.2', '585.1', '585.2', '585.3', '585.4', '585.5', '585.6', '585.9', '586', 'N18.1', 'N18.2', 'N18.3', 'N18.30', 'N18.31', \
        'N18.32', 'N18.4', 'N18.5', 'N18.6', 'N18.9', 'N19', 'Z49.0', 'Z49.01', 'Z49.02'], \
    'Liver disease': ['070.32', '070.33', '070.54', '456.0', '456.1', '456.2', '571.0', '572.3', '572.8', 'V42.7', 'I86.4', 'I98.2', \
        'K71.1', 'K71.1', 'K76.0', 'Z94.4', '571.2', '571.3', '571.4', '571.40', '571.41', '571.42', '571.49', '571.5', '571.6', \
        '571.8', '571.9', 'B18.0', 'B18.1', 'B18.2', 'B18.8', 'B18.9', 'I85.0', 'I85.00', 'I85.01', 'I85.1', 'I85.10', 'I85.11', \
        'K70.0', 'K70.1', 'K70.10', 'K70.11', 'K70.2', 'K70.3', 'K70.30', 'K70.31', 'K70.4', 'K70.40', 'K70.41', 'K70.9', 'K71.3', \
        'K71.4', 'K71.5', 'K71.50', 'K71.51', 'K72.0', 'K72.00', 'K72.01', 'K72.1', 'K72.10', 'K72.11', 'K72.9', 'K72.90', 'K72.91', \
        'K73.0', 'K73.1', 'K73.2', 'K73.8', 'K73.9', 'K74.0', 'K74.00', 'K74.01', 'K74.02', 'K74.1', 'K74.2', 'K74.3', 'K74.4', \
        'K74.5',  'K74.6', 'K74.60', 'K74.69', 'K76.2', 'K76.3', 'K76.4', 'K76.5', 'K76.6', 'K76.7', 'K76.8', 'K76.81', 'K76.89', \
        'K76.9'], \
    'Depression': ['300.4', '301.12', '309.0', '309.1', '311', 'F20.4', 'F34.1', 'F41.2', 'F43.2', 'F31.3', 'F31.4', 'F31.5', 'F32', \
        'F32.0', 'F32.1', 'F32.2', 'F32.3', 'F32.8', 'F32.9', 'F33', 'F33.0', 'F33.1', 'F33.2', 'F33.3', 'F33.4', 'F33.8', 'F33.9'], \
    'Alcohol abuse': ['291.1', '291.2', '303.9', '305.0', 'V11.3', '291.5', '291.8', '291.81', '291.82', '291.89', '291.9', 'F10', \
        'E52', 'G62.1', 'I42.6', 'K29.2', 'K70.0', 'K70.3', 'K70.9', 'Z50.2', 'Z71.4', 'Z72.1', 'T51', 'T51.0', 'T51.1', 'T51.2', \
        'T51.3', 'T51.8', 'T51.9'], \
    'PTSD': ['F43.1', '309.81'], \
    'TBI': ['S02.0', 'S02.1', 'S06.2', 'S06.3', 'S06.8', 'S06.A', \
        'S06.0', 'S06.1', 'S09', '850', '851', '852', '853', '854', \
        'V15.52', 'Z87.8', 'S02', 'S04', 'S06', 'S07', 'S09', \
        '800', '801', '803', '804', '907', 'S06.4', 'S06.5', \
        'S06.6', 'S06.9'], \
    'Cerebrovascular\ndisease': ['I60', 'I61', 'I62', 'I63', 'I64', \
        'I65', 'I66', 'I67', 'I68', 'I69', '430', '431', '432', '433', \
        '434', '435', '436', '437', '438'], \
    '(post-operative) delerium': ['293.0', 'T81.89', 'F05'], \
    'Hearing aids': ['Z97.4', 'V53.2'], \
    'Smoking': ['305.1', 'V15.82', 'F17', 'Z87.891'], \
    'Dyslipidemia': ['272', 'E78'], \
    'Sleep apnea': ['G47.3', '327.2'], \
    'Anxiety/related\ndisorder': ['300.0', '300.2', '300.3', '308.0', '309.2', '309.8', \
                                 'F43'], \
    'Myocardial infarction': ['410', '412', 'I21', 'I22', 'I25.2'], \
    #presumed negative controls
    'Myopia': ['H52.1', '367.1'], \
    'Urine stones': ['N20', '592'], \
    'Gout': ['M10', '274'], \
}
    
medication_map = \
{ \
    'cardiovascular': ['lisinopril', 'ramipril', 'losartan', 'amiodarone', \
                    'warfarin', 'aspirin', 'bisoprolol', 'amlodipine', \
                    'simvastatin', 'digoxin', 'bendroflumethiazide', \
                     'atorvastatin', 'fluvastatin', 'rosuvastatin', '\
                     dabigatran', 'rivaroxaban', 'apixaban', 'edoxaban', \
                     'betrixaban', 'valsartan', 'nitroglycerin', 'nifedipine', \
                     'niacin', 'metoprolol', 'lovastatin', 'hydrochlorothiazide', \
                     'furosemide', 'enalapril', 'diltiazem', 'clopidogrel', \
                     'carvedilol', 'benazepril', 'atenolol'], \
    'dyslipidemia': ['fenofibrate', 'gemfibrozil'], \
    'diabetes': ['insulin', 'metformin', 'acarbose', 'miglitol', \
                 'bromocriptine', 'gliclazide', 'glipizide', 'glimepiride', 'tolbutamide'], \
    'depression': ['fluoxetine', 'paroxetine', 'fluvoxamine', 'citalopram', \
                   'escitalopram', 'sertraline', 'desvenlafaxine', 'duloxetine', \
                   'levomilnacipran', 'milnacipran', 'venlafaxine', 'nefazodone', \
                   'trazodone', 'vilazodone', 'vortioxetine', 'esketamine', \
                   'moclobemide', 'isocarboxazid', 'phenelzine', 'tranylcypromine', \
                   'quetiapine', 'mirtazapine', 'bupropion'], \
    'anxiety': ['clonazepam', 'alprazolam', 'lorazepam', 'bromazepam', \
                'oxazepam', 'chlordiazepoxide', 'diazepam', 'clorazepate'], \
}
    
zips = ['330', '331', '212', '104', '206', '207', '390', '391', '392', \
        '701', '317', '291', '922', '798', '799', '885']
zips2 = ['900', '901', '902', '903', '904', '905', '906', '907', '908', \
         '910', '911', '912', '913', '914', '915', '916', '917', '918', \
         '935', '600', '601', '602', '603', '604', '605', '606', '607', \
         '608', '850', '851', '852', '853', '770', '772', '773', '774', \
         '775', '919', '920', '921', '906', '907', '928', '927', '928', \
         '112', '111', '334']

def load_file(file_name):
    contents = pd.read_csv(file_name)
    return contents


def normalize(column):
    return (column - np.nanmin(column)) / (np.nanmax(column - np.nanmin(column)))

###############################################
#prepare population

#find how many patients are in each conversion window
def get_converters(pop, horizon, strata='convert_time'):
    tte = pop[strata]
    tte = tte[~np.isnan(tte)]
    obs_ttes, counts = np.unique(tte, return_counts=True)
    print(strata + ' and counts', obs_ttes, counts)
    
    converters = []
    for i in range(horizon):
        conv_i = pop[pop[strata] == i]
        converters.append(conv_i)
    
    smallest_pop = obs_ttes[np.argmin(counts[obs_ttes<horizon])]
    
    return converters, int(smallest_pop)


#match on age of mci onset and number of bmi measurements
def match(match_on, pops, mult=1, intersection=None, plot_name='conv'):
    matched_pop = []
    
    #find the age/bmi distribution common to all populations
    #example: 10 patients with mci onset at 70 with 3 bmi measurements at all horizons
    if intersection is None:
        intersection = []
        for i in range(len(pops)-1, -1, -1):
            current_ids = pops[i][match_on[0]].map(str).to_numpy()
            for j in range(1, len(match_on)):
                current_ids = current_ids + '-' + pops[i][match_on[j]].map(str).to_numpy()
            if len(intersection) == 0:
                intersection = current_ids
                continue
            used = []
            for j in range(current_ids.shape[0]):
                matches = np.where(intersection == current_ids[j])[0]
                for match in matches:
                    if not match in used:
                        used.append(match)
                        break
            intersection = intersection[used]
    
    print('how many all pops have in common', intersection.shape)
    
    #for each horizon, match some constant multiple patients to intersection
    for i in range(len(pops)):
        current_ids = pops[i][match_on[0]].map(str).to_numpy()
        for j in range(1, len(match_on)):
            current_ids = current_ids + '-' + pops[i][match_on[j]].map(str)
        matched_i = []
        for j in range(len(intersection)):
            matches = np.where(current_ids == intersection[j])[0]
            matches = np.setdiff1d(matches, matched_i)
            #perm = np.random.permutation(matches.shape[0])
            #matches = matches[perm]
            if matches.shape[0] > 0:
                mult_i = mult
                if isinstance(mult, list):
                    mult_i = mult[i]
                matched_i += [int(matches[k]) for k in range(min(mult_i, len(matches)))]
        print(i, len(matched_i))
        matched_pop.append(pops[i].iloc[matched_i])
    
    matched_pop = pd.concat(matched_pop)
    print('num matched', matched_pop.shape[0])
    return matched_pop, intersection


#find who is censored (lost to followup without acd onset before 10 years)
def get_censored(pop, horizon):
    tte = pop['convert_time']
    followup = pop['followup']
    
    censored = pd.concat([pop.loc[np.isnan(tte)], pop.loc[followup>=horizon]])
    censored.loc[censored['followup'] >= horizon, 'followup'] = horizon-1
    censored.loc[censored['convert_time'] >=horizon, 'followup'] = horizon-1
    print('num censored', censored.shape[0])
    return censored
    
        
#get the final population
def get_population(file_root, pop_file, horizon, match_on, vis_bin, age_bin):
    raw_pop = load_file(pop_file)
    print('raw pop size', raw_pop.shape)
    vit = 'bmi'
    bmi = pd.read_csv(file_root + vit + '.csv')
    
    raw_pop['num_outpat_orig'] = raw_pop['num_outpat']
    raw_pop['mci_age_orig'] = raw_pop['mci_age']
    raw_pop['history_len_days_orig'] = raw_pop['history_len_days']
    print('history length', np.percentile(raw_pop['history_len_days_orig'], [0, 25, 50, 75, 100]))
    
    #get rid of anyone with less than 1 years of history
    raw_pop = raw_pop[raw_pop['history_len_days_orig'] >= 365.25]
    
    #bin the bmi counts and ages for matching
    raw_pop = raw_pop.merge(bmi[['patientsid', vit + '_count']], how='left', left_on='patientID', right_on='patientsid')
    raw_pop[vit + '_count_orig'] = raw_pop[vit + '_count']
    if 'acd' in pop_file:
        raw_pop[vit + '_count'] = raw_pop[vit + '_count'] // 10 #bin the counts by 10 measurements per bin
    
    raw_pop['num_outpat'] = raw_pop['num_outpat'] // vis_bin #bins
    raw_pop['followup'] = raw_pop['followup'] // 365.25 #followup in  years
    raw_pop['mci_age'] = raw_pop['mci_age'] // (age_bin*2) #bin the ages by 2 years per bin
    
    #match over a 10 year horizon, match both cenosred and uncensored patients
    match_horizon = 10
    converters, smallest_pop = get_converters(raw_pop, match_horizon)
    pop_mult = [3, 6, 4, 3, 2, 2, 2, 2, 2, 1]
    matched_pop, intersect = match(match_on, converters, pop_mult, None, 'conv')
    
    censored_pop = get_censored(raw_pop, match_horizon)
    cen_groups, smallest_cen = get_converters(censored_pop, match_horizon, 'followup')
    cen_mult = 5
    censored_pop, _ = match(match_on, cen_groups, cen_mult, intersect, 'cen')
    
    #put matched censored and uncensored populations together
    print('match x', cen_mult)
    final_pop = pd.concat([matched_pop, censored_pop])
    print(final_pop.columns)
    final_pop = final_pop.drop(labels=vit + '_count', axis=1)
    final_pop = final_pop.drop_duplicates(subset='patientID')
    return final_pop


###############################################
#format features

#1-hot encodings for demographic features
def format_dem(file_root, pop):
    data = pd.read_csv(file_root + 'demographics.csv').fillna('not registered')
    feats = ['GenderName', 'RaceName', 'EthnicityName'] 

    ids = pop['patientID']
    dem_feats = {'patientID': ids}
    
    for i in range(len(feats)):
        feat = feats[i]
        feat_vals = data[[feat]].astype(str).to_numpy()
        feat_vals = np.unique(feat_vals)
        val_assigned = []
        for j in range(len(feat_vals)):
            if 'not registered' in feat_vals[j]:
                continue
            data_ij = data.loc[data[feat] == feat_vals[j]]['patientID']
            data_ij = np.setdiff1d(data_ij, val_assigned)
            feat_col = np.zeros(ids.shape)
            feat_col1 = np.where(np.isin(ids, data_ij))[0]
            feat_col[feat_col1] = 1
            feat_name = feat + '_' + feat_vals[j]
            dem_feats['dem_' + feat_name] = feat_col
            val_assigned = np.union1d(val_assigned, data_ij)
            
    dem_feats = pd.DataFrame(dem_feats) 
    print('demographics', dem_feats.shape[0])
    columns = list(dem_feats.columns)
    for feat in columns:
        if feat == 'patientID':
            continue
        count = np.sum(dem_feats[feat])
        print(feat, count, count/dem_feats.shape[0])
        
    return dem_feats


#1-hot encodings for diagnises and medications
def format_diag_meds(file_root, pop, category='diag'):
    if category == 'diag':
        data = pd.read_csv(file_root + 'diagnoses.csv').fillna('not registered')
        feats = list(comorbid_map.keys())
        map_dict = comorbid_map
    elif category == 'meds':
        data = pd.read_csv(file_root + 'medications.csv').fillna('not registered')
        feats = list(medication_map.keys())
        map_dict = medication_map

    ids = pop['patientID']
    diag_feats = {'patientID': ids}
    
    #consider feats up to 5 years before mci
    data = data[data['time_before'] <= 1827]
    
    for i in range(len(feats)):
        feat = feats[i]
        data_i = data.loc[data['variable_name'].str.contains('|'.join(map_dict[feat]), case=False)==True]['patientID']
        feat_col = np.zeros(ids.shape)
        feat_col1 = np.where(np.isin(ids, data_i))[0]
        feat_col[feat_col1] = 1
        diag_feats[category + '_' + feat] = feat_col
            
    diag_feats = pd.DataFrame(diag_feats) 
    print(category, diag_feats.shape[0])
    columns = list(diag_feats.columns)
    for feat in columns:
        if feat == 'patientID':
            continue
        count = np.sum(diag_feats[feat])
        print(feat, count, count/diag_feats.shape[0])
        
    return diag_feats


#healthcare utlization - make 1-hot encodings based on values binned into quintiles
def format_util(file_root, pop, vis_bin, age_bin):
    raw_pop = pd.read_csv(file_root + 'population_acd.csv')
    feats = ['patientID', 'mci_age', 'num_inpat', 'num_outpat']#, 'history_len_days']
    util_feats = raw_pop[feats]
    util_feats = util_feats.loc[util_feats['patientID'].isin(pop['patientID'])]
    util_feats = util_feats.drop_duplicates(subset=['patientID'])
    
    print('utilization', util_feats.shape, pop.shape)
    for feat in feats:
        if feat == 'patientID':
            continue
        print('quartiles', feat, np.percentile(util_feats[feat], [0, 25, 50, 75, 100]))
        quintiles = np.percentile(util_feats[feat], [0, 20, 40, 60, 80, 100])
        quintiles[-1] += 1
        print('n-tiles', feat, quintiles)
        for i in range(len(quintiles)-1):
            bin_name = 'util_' + feat + str(quintiles[i]) + '_' + str(quintiles[i+1])
            if bin_name in list(util_feats.columns) or quintiles[i] == quintiles[i+1]:
                continue
            quintiled_vals = np.zeros((util_feats.shape[0],))
            quintiled_vals[(util_feats[feat] >= quintiles[i]) & (util_feats[feat] < quintiles[i+1])] = 1
            quintiled_vals[util_feats[feat] >= quintiles[i+1]] = 0
            quintiled_vals[util_feats[feat] < quintiles[i]] = 0
            util_feats[bin_name] = quintiled_vals
        
    return util_feats


#vitals and labs - make 1-hot encodings based on values binned into quintiles
#also include a missingness indicator
def format_vitlabs(file_root, pop):
    feats = ['bmi', 'bpsys', 'bpdia', 'ast', 'alt', 'hemo', \
             'creatinine', 'ldl', 'hdl']
    stats = ['avg', 'recent', 'range', 'std', 'count']

    vitlab_feats = {'patientID': pop['patientID']}
    
    for i in range(len(feats)):
        feat = feats[i]
        feat_name = 'vitl' + feat
        raw_data = pd.read_csv(file_root + feat + '.csv')#.fillna(-1)
        
        for j in range(len(stats)):
            values = pop.merge(raw_data, how='left', left_on='patientID', right_on='patientsid')
            values = values[feat + '_' + stats[j]]
            print('quartiles', feat, stats[j], np.nanpercentile(values, [0, 25, 50, 75, 100]))
            quintiles = np.nanpercentile(values, [0, 20, 40, 60, 80, 100])
            quintiles[-1] += 1
            values = values.fillna(-1)
            print('n-tiles', feat, stats[j], quintiles)#, values[:10])
            for k in range(len(quintiles)-1):
                bin_name = feat_name + '_' + stats[j] + str(quintiles[k]) + '_' + str(quintiles[k+1])
                if quintiles[k] == quintiles[k+1]:
                    continue
                vitlab_feats[bin_name] = copy.deepcopy(values)
                vitlab_feats[bin_name][(values >= quintiles[k]) & (values < quintiles[k+1])] = 1
                vitlab_feats[bin_name][values >= quintiles[k+1]] = 0
                vitlab_feats[bin_name][values < quintiles[k]] = 0
               
            if stats[j] != 'recent':
                continue
            missing = feat_name + '_' + stats[j] + '_missing'
            vitlab_feats[missing] = copy.deepcopy(values)
            vitlab_feats[missing][values == -1] = 1
            vitlab_feats[missing][values != -1] = 0
            num_missing = np.sum(vitlab_feats[missing])
            prop_missing = num_missing / vitlab_feats[missing].shape[0]
            print(feat_name + '_missing', num_missing, prop_missing) 
            
    for col_name in list(vitlab_feats.keys()):
        vitlab_feats[col_name] = list(vitlab_feats[col_name].values)        
    vitlab_feats2 = pd.DataFrame(data=vitlab_feats)
    vitlab_feats2 = vitlab_feats2.fillna(0)
    
    return vitlab_feats2


#zip codes - bin according the previous work (see mapping above)
def format_zip(file_root, pop):
    data = pd.read_csv(file_root + 'zip.csv').fillna('not registered')
    feats = ['zip'] 

    ids = pop['patientID']
    dem_feats = {'patientID': ids}
    
    for i in range(len(feats)):
        feat = feats[i]
        feat_vals_all = data[feat].astype(str).str[:3]
        for j in range(len([1])):#feat_vals)):
            data_ij = data.loc[feat_vals_all.str.contains('|'.join(zips+zips2), case=False)==True]['patientID']
            feat_col = np.zeros(ids.shape)
            feat_col1 = np.where(np.isin(ids, data_ij))[0]
            feat_col[feat_col1] = 1
            feat_name = feat 
            dem_feats[feat_name] = feat_col
            
    dem_feats = pd.DataFrame(dem_feats) 
    print('zip codes', dem_feats.shape)
    columns = list(dem_feats.columns)
    for feat in columns:
        if feat == 'patientID':
            continue
        count = np.sum(dem_feats[feat])
        print(feat, count, count/dem_feats.shape[0])
        
    return dem_feats


#put together all of the features (call functions above)
def get_features(pop, file_root, vis_bin, age_bin, util_only=False, preload=False):
    if preload and not util_only: #if loading saved features, don't have to preprocess everything
        file_root = file_root[:-1] + '_45years/'
        feats = pd.read_csv(file_root + 'processed_feats_240528.csv')
        feats = feats.merge(pop[['patientID']])
        feat_names = list(feats.columns)
        feat_names_filt = []
        for feat_name in feat_names:
            #['avg', 'recent', 'range', 'std', 'count']
            if not('avg' in feat_name or 'range' in feat_name or \
                   'std' in feat_name or 'count' in feat_name):
                feat_names_filt.append(feat_name)
        feats = feats[feat_names_filt]
        return feats
    
    #preprocess features and combine
    dem_feats = format_dem(file_root, pop)
    if not util_only:
        diag_feats = format_diag_meds(file_root, pop)
        med_feats = format_diag_meds(file_root, pop, category='meds')
    file_root = file_root[:-1] + '_45years/'
    util_feats = format_util(file_root, pop, vis_bin, age_bin)
    vitlab_feats = format_vitlabs(file_root, pop)
    zip_feats = format_zip(file_root, pop)
    print('dem', dem_feats.shape[1]-1)
    if not util_only:
        print('diag', diag_feats.shape[1]-1)
        print('med', med_feats.shape[1]-1)
    print('vitlab', vitlab_feats.shape[1]-1)
    print('util', util_feats.shape[1]-1)
    print('zip', zip_feats.shape[1]-1)
    
    if not util_only: 
        feats = dem_feats.merge(diag_feats, on='patientID')
        feats = feats.merge(med_feats, on='patientID')
        feats = feats.merge(vitlab_feats, on='patientID')
        feats = feats.merge(zip_feats, on='patientID')
        feats = feats.merge(util_feats, on='patientID')
        feats = feats.drop_duplicates(subset='patientID')
    
    else:
        feats = dem_feats.merge(util_feats, on='patientID', how='outer')
        feats = feats.merge(zip_feats, on='patientID', how='outer')
    
    feat_names = np.array(feats.columns.values.tolist())[1:]
    perm = np.random.permutation(feat_names)
    feat_names = np.concatenate([['patientID'], perm])
    feats = feats[feat_names]
    print('feature shape', feats.shape)
    #comment line below in if want to save processed features
    #save processed features for quick loading so don't have to preprocess on fly all the time
    #feats.to_csv(file_root + 'processed_feats_240528.csv') 
    
    return feats


#labels - time to event (acd onset) and censoring indicator
def get_labels(pop_file, horizon, feats):
    raw_pop = load_file(pop_file)
    
    labels = []
    
    ids = feats['patientID']
    
    for i in range(ids.shape[0]):
        pat_row = raw_pop[raw_pop['patientID']==ids.iloc[i]]
        tte = pat_row['convert_time'].iloc[0] 
        followup = pat_row['followup'].iloc[0] // 365.25 
        had, time = False, followup
        if not np.isnan(tte) and tte<horizon:
            had, time = True, tte
        else:
            time = min(time, horizon-1)
        lab_i = (had, time)
        labels.append(lab_i)
    
    labels = np.array(labels, dtype=[('Status', '?'), ('TTE', '<f8')])
    
    print('done getting labels', np.unique(labels, return_counts=True))
    return labels


###############################################
#put togehter the dataset

def get_dataset(horizon, util_only, preload):
    file_root = 'tables_cipher/' #change this where the tables are stored
    
    acd = True
    print('acd', acd)
    pop_file = file_root + 'population.csv'
    if acd:
        pop_file = file_root + 'population_acd.csv'
    
    vis_bin = 2
    age_bin = 1
    match_on = ['mci_age', 'bmi_count']
    print('horizon', horizon)
    
    pop = get_population(file_root, pop_file, horizon, match_on, vis_bin, age_bin)
    feats = get_features(pop, file_root, vis_bin, age_bin, util_only, preload)
    
    feats = feats.fillna(0)
    labs = get_labels(pop_file, horizon, feats)
    
    #don't want these features to train the models (don't want the raw values, used the binned values instead)
    to_drop = ['patientID', 'mci_age', 'num_outpat', 'num_inpat']#, 'history_len_days']
    for feat in list(feats.columns):
        if 'Unnamed' in feat or '331.83' in feat or 'Unknown' in feat:
            to_drop.append(feat)
    feats = feats.drop(labels=to_drop, axis=1)
    
    print('feature dimnensions', feats.shape, feats.columns)
    print('label dimensions', labs.shape)
    
    return feats, labs

###############################################
if __name__ == '__main__':
    horizon, util_only = 10, False
    feats, labs = get_dataset(horizon, util_only)
    print('smile yay')