--- a
+++ b/preprocessing/FIDDLE_steps_2.py
@@ -0,0 +1,649 @@
+"""
+FIDDLE Preprocessing steps
+1. Pre-filter
+2. Transform
+3. Post-filter
+"""
+try:
+    from FIDDLE_helpers import *
+except:
+    from .FIDDLE_helpers import *
+import time
+import json
+import pickle5 as pickle
+import joblib
+from datetime import datetime
+import multiprocessing
+
+# import pickle4reducer
+# import multiprocessing
+# ctx = multiprocessing.get_context()
+# ctx.reducer = pickle4reducer.Pickle4Reducer()
+
+def detect_variable_data_type(df_data, value_type_override, args):
+    data_path = args.data_path
+    print_header('*) Detecting value types', char='-')
+    
+    data_types = []
+    df = df_data
+    assert val_col in df.columns
+
+    # Collect the unique values of each variable
+    # values_by_variable: dict(variable_name -> [value1, value2, ...])
+    d = df[[var_col, val_col]].drop_duplicates().sort_values(by=[var_col, val_col])
+    values_by_variable = defaultdict(list)
+    for n,v in zip(d[var_col], d[val_col]):
+        values_by_variable[n].append(v)
+
+    # Determine type of each variable
+    for variable, values in sorted(values_by_variable.items()):
+        # Manual override type in config
+        if variable in value_type_override:
+            data_types.append((variable, value_type_override[variable]))
+            # Force categorical values to be a string
+            if value_type_override[variable] == 'Categorical' and \
+                any(is_numeric(v) for v in values if not pd.isnull(v)):
+                m_var = df[var_col] == variable
+                df.loc[m_var, val_col] = df.loc[m_var, val_col].apply(lambda s: '_' + str(s))
+        else:
+            if len(values) == 1 and pd.isnull(values[0]):
+                data_types.append((variable, 'None'))
+            elif all(is_numeric(v) for v in values if not pd.isnull(v)):
+                data_types.append((variable, 'Numeric'))
+            elif any(is_numeric(v) for v in values if not pd.isnull(v)):
+                data_types.append((variable, 'Numeric + Categorical'))
+            else:
+                data_types.append((variable, 'Categorical'))
+    
+    df_types = pd.DataFrame(data_types, columns=['variable_name', 'value_type'])
+    df_types[var_col] = df_types[var_col].astype(str)
+    df_types = df_types.set_index(var_col)
+    fpath = data_path + 'value_types.csv'
+    df_types.to_csv(fpath, quoting=1)
+    print('Saved as:', fpath)
+    return df, df_types['value_type']
+
+
+def split_by_timestamp_type(df):
+    print_header('*) Separate time-invariant and time-dependent', char='-')
+    
+    variables_inv = df[pd.isnull(df[t_col])][var_col].unique() # Invariant variables have t = NULL
+    df_time_invariant = df[df[var_col].isin(variables_inv)]
+    df_time_series = df[~df[var_col].isin(variables_inv)]
+    
+    print('Variables (time-invariant):', len(variables_inv))
+    print('Variables (time-dependent):', df[var_col].nunique() - len(variables_inv))
+    print('# rows    (time-invariant):', len(df_time_invariant))
+    print('# rows    (time-dependent):', len(df_time_series))
+    return df_time_invariant, df_time_series
+
+
+def process_time_dependent(df_data_time_series, args):
+    data_path = args.data_path
+    theta_2 = args.theta_2
+    
+    print_header('2-B) Transform time-dependent data', char='-')
+    dir_path = data_path + '/'
+    start_time = time.time()
+
+    ## Create NxLxD^ table
+    df_time_series, dtypes_time_series = transform_time_series_table(df_data_time_series, args)
+    print('Time elapsed: %f seconds' % (time.time() - start_time), flush=True)
+
+    ##############
+    joblib.dump(df_time_series, args.data_path + 'df_time_series,{}.joblib'.format(datetime.now().isoformat()))
+    joblib.dump(dtypes_time_series, args.data_path + 'dtypes_time_series,{}.joblib'.format(datetime.now().isoformat()))
+    ##############
+
+    ## Map variables to features
+    X_all, X_all_feature_names = map_time_series_features(df_time_series, dtypes_time_series, args)
+    sparse.save_npz(dir_path + 'X_all.npz', X_all)
+    with open(dir_path + 'X_all.feature_names.json', 'w') as f:
+        json.dump(list(X_all_feature_names), f, sort_keys=True)
+    print('Time elapsed: %f seconds' % (time.time() - start_time), flush=True)
+    
+    ## Filter features
+    if not args.postfilter:
+        return X_all, X_all_feature_names, {}
+    print_header('3-B) Post-filter time-dependent data', char='-')
+    print(X_all.shape, X_all.density)
+    X, X_feature_names, X_feature_aliases = post_filter_time_series(X_all, X_all_feature_names, theta_2, args)
+    print(X.shape, X.density)
+    print('Time elapsed: %f seconds' % (time.time() - start_time))
+
+    ## Save output
+    print()
+    print('Output')
+    print('X: shape={}, density={:.3f}'.format(X.shape, X.density))
+    sparse.save_npz(dir_path + 'X.npz', X)
+    with open(dir_path + 'X.feature_names.json', 'w') as f:
+        json.dump(list(X_feature_names), f, sort_keys=True)
+    with open(dir_path + 'X.feature_aliases.json', 'w') as f:
+        json.dump(X_feature_aliases, f, sort_keys=True)
+    
+    print('Total time: %f seconds' % (time.time() - start_time))
+    print('', flush=True)
+    return X, X_feature_names, X_feature_aliases
+
+######
+# Time-series routines
+######
+def func_encode_single_time_series(i, g, variables, variables_num_freq, T, dt, stats_functions, impute=True):
+    try:
+        assert g.index.nunique() == 1
+        assert g.index.unique()[0] == i
+        # non-frequent
+        variables_non = sorted(set(variables) - set(variables_num_freq))
+        if len(variables_non) > 0:
+            df_j = pivot_event_table(g).reindex(columns=variables_non).sort_index()
+            df_values_j = most_recent_values(df_j, variables, T, dt)
+            df_out = df_values_j
+        
+        if len(variables_num_freq) > 0:
+            # frequent
+            # we're only producing mask, ffill, and statistics if the data is measured frequently enough
+            df_i = pivot_event_table(g).reindex(columns=variables_num_freq).sort_index()
+            mask_i = presence_mask(df_i, variables_num_freq, T, dt)
+#             delta_t_i = get_delta_time(mask_i)
+#             df_i = impute_ffill(df_i, variables_num_freq, T, dt, mask_i)
+            df_stats_i = summary_statistics(df_i, variables_num_freq, stats_functions, T, dt)
+            df_values_i = most_recent_values(df_i, variables, T, dt)
+#             if impute:
+#                 check_imputed_output(df_values_i)
+#                 check_imputed_output(df_stats_i)
+            if len(variables_non) > 0:
+                df_out = df_out.join([mask_i, df_values_i, df_stats_i])
+            else:
+                df_out = mask_i.join([df_values_i, df_stats_i])
+#             df_out = df_out.join([mask_i, delta_t_i, df_values_i, df_stats_i])
+    except:
+        print(i)
+        print(g)
+        raise
+    return i, df_out
+
+def divide_chunks(l, n): 
+    # looping till length l 
+    for i in range(0, len(l), n):  
+        yield l[i:i + n]
+
+def form_batches_of_examples(df_in, args, batch_size=1000):
+    grouped = df_in.set_index(ID_col)
+    IDs = list(grouped.index.unique())
+    batches_IDs = list(divide_chunks(IDs, batch_size))
+    batches = [grouped.loc[chunk] for chunk in batches_IDs]
+    return batches, batches_IDs
+
+def process_batch_time_series(first_arg):
+    batch, batch_IDs, args = first_arg
+    variables, variables_num_freq = args.variables, args.variables_num_freq
+    out = dict(
+        func_encode_single_time_series(i, batch.loc[i:i], variables, variables_num_freq, args.T, args.dt, args.stats_functions) 
+        for i in batch_IDs
+    )
+    return out
+
+def transform_time_series_table(df_in, args):
+    data_path = args.data_path
+    theta_freq = args.theta_freq
+    stats_functions = args.stats_functions
+    N, L = args.N, args.L
+    df_population = args.df_population
+    parallel = args.parallel
+    
+    ## TODO: asserts shape of df_in
+
+    # Determine all unique variable names
+    variables = args.variables
+#     variables = get_unique_variables(df_in) ############
+#     assert df_in[var_col].nunique() == len(variables)
+    print('Total variables    :', len(variables), flush=True)
+    
+    # Determine frequent variables -> we'll calculate statistics, mask, and delta_time only on these
+    variables_num_freq = args.variables_num_freq
+#     variables_num_freq = get_frequent_numeric_variables(df_in, variables, theta_freq, args)
+    print('Frequent variables :', list(variables_num_freq))
+    print('{} = {}'.format('M\u2081', len(variables_num_freq)))
+    print('{} = {}'.format('M\u2082', len(variables) - len(variables_num_freq)))
+    print('{} = {} {}'.format('k ', len(stats_functions), stats_functions))
+    
+    print()
+    print('Transforming each example...', flush=True)
+    args.variables = variables
+    args.variables_num_freq = variables_num_freq
+    
+    # Encode time series table for each patient
+    batches, batches_IDs = form_batches_of_examples(df_in, args, batch_size=1000)
+    print('Batches of size 1000: ', len(batches), flush=True)
+    
+### TEST
+#     for batch, batch_IDs in zip(batches, batches_IDs):
+#         process_batch_time_series((batch, batch_IDs, args))
+#         print('done one')
+#         break
+#     exit()
+### TEST
+
+    pool = multiprocessing.Pool(multiprocessing.cpu_count())
+    out = list(tqdm(pool.imap_unordered(
+        process_batch_time_series,
+        zip(batches, batches_IDs, [args]*len(batches))), total=len(batches)
+    ))
+    pool.close()
+    pool.join()
+    
+    out = dict((key, d[key]) for d in out for key in d)
+    print()
+    print('Parallel processing done', flush=True)
+    
+    # Handle IDs not in the table
+#     df_original = list(out.values())[0]
+#     df_copy = pd.DataFrame().reindex_like(df_original)
+#     for i, j in df_original.dtypes.iteritems():
+#         if i.endswith('_mask'):
+#             assert j == bool
+#             df_copy[i] = False
+#             df_copy[i] = df_copy[i].astype(bool)
+#         if i.endswith('_delta_time'):
+#             df_copy[i] = 0
+#             df_copy[i] = df_copy[i].astype(int)
+#         if j == 'object':
+#             df_copy[i] = df_copy[i].astype('object')
+
+#     for ID in tqdm(df_population.index.values[:N]):
+#         if ID not in out:
+#             out[ID] = df_copy.copy()
+
+#     out = {ID: out[ID] for ID in df_population.index.values}
+    N = len(out)
+#     assert len(out) == N
+    D_timeseries = out
+
+    print('Filled no-data examples done', flush=True)
+#     # check each example have identical LxD table structure
+#     ID0 = sorted(D_timeseries.keys())[0]
+#     df0 = D_timeseries[ID0]
+#     for ID, df_i in D_timeseries.items():
+#         pd.testing.assert_index_equal(df_i.index, df0.index)
+#         pd.testing.assert_index_equal(df_i.columns, df0.columns)
+
+    D_timeseries = out
+    D_ = len(list(D_timeseries.values())[0].columns)
+    
+    ########
+#     joblib.dump(D_timeseries, args.data_path + 'D_timeseries,{}.joblib'.format(datetime.now().isoformat()))
+    ########
+    
+#     # (N*L)xD^ table
+#     ## Create MultiIndex of (ID, time_bin)
+#     index = sum([ 
+#         [(ID, t_) for t_ in list(df_.index)]
+#         for ID, df_ in sorted(D_timeseries.items()) 
+#     ], [])
+#     index = pd.Index(index)
+    index = [ID for ID, df_ in sorted(D_timeseries.items())]
+    index = pd.Index(index)
+# #     assert len(index) == N * L
+    
+    ## Assume all dataframes have the same columns, used after concatenation
+    columns = list(sorted(D_timeseries.items())[0][1].columns)
+    columns = np.array(columns)
+    dtypes = sorted(D_timeseries.items())[0][1].dtypes
+    
+    ## Convert each df to a numpy array
+    ## Concatenate **sorted** numpy arrays (faster than calling pd.concat)
+    feature_values = [(ID, df_.to_numpy()) for ID, df_ in sorted(D_timeseries.items())]
+    time_series = np.concatenate([feat_val[1] for feat_val in feature_values])
+#     assert time_series.shape == (len(index), len(columns))
+    
+    df_time_series = pd.DataFrame(data=time_series, index=index, columns=columns)
+    
+    print()
+    print('(N \u00D7 L \u00D7 ^D) table :\t', (N, L, len(columns)))
+    return df_time_series, dtypes
+
+def map_time_invariant_features(df, bin_numeric=True):
+    # Categorical -> binary features
+    # Numeric -> binary/float-valued features
+    if bin_numeric:
+#         df_mixed = df.apply(smart_qcut, q=5)
+#         features_mixed = pd.get_dummies(df_mixed, columns=df_mixed.columns, prefix_sep=':')
+#         time_invariant_features = features_mixed
+#         assert time_invariant_features.astype(int).dtypes.nunique() == 1
+        print('start discretization', flush=True)
+#         out = [smart_qcut_dummify(df[col], q=5) for col in tqdm(df.columns)]
+        
+#         pool = multiprocessing.Pool(multiprocessing.cpu_count())
+#         out = list(tqdm(pool.imap_unordered(
+#             smart_qcut_dummify_5,
+#             [df[col] for col in tqdm(df.columns)]), total=len(df.columns)
+#         ))
+#         pool.close()
+#         pool.join()
+
+        pool = multiprocessing.Pool(multiprocessing.cpu_count())
+        out_0 = list(tqdm(pool.imap_unordered(
+            smart_qcut,
+            [df[col] for col in df.columns]), total=len(df.columns)
+        ))
+        cols_data, discretization_bins = zip(*out_0)
+        out = list(tqdm(pool.imap_unordered(dummify, cols_data), total=len(df.columns)
+        ))
+        pool.close()
+        pool.join()
+        
+        time_invariant_features = pd.concat(out, axis=1).sort_index(axis=1)
+        feature_names_all = time_invariant_features.columns.values
+        sdf = time_invariant_features.astype(pd.SparseDtype(int, fill_value=0))
+        s_ = sparse.COO(sdf.sparse.to_coo())
+    else:
+        raise NotImplemented
+        # Split a mixed column into numeric and string columns
+        for col in df.columns:
+            col_data = df[col]
+            col_is_numeric = [is_numeric(v) for v in col_data if not pd.isnull(v)]
+            if not all(col_is_numeric) and any(col_is_numeric): # have mixed type values
+                numeric_mask = col_data.apply(is_numeric)
+                df[col+'_str'] = df[col].copy()
+                df.loc[~numeric_mask, col] = np.nan
+                df.loc[numeric_mask, col+'_str'] = np.nan
+        
+        out = [smart_dummify_impute(df[col]) for col in tqdm(df.columns)]
+        time_invariant_features = pd.concat(out, axis=1)
+        feature_names_all = time_invariant_features.columns.values
+        sdf = time_invariant_features.astype(pd.SparseDtype(float, fill_value=0))
+        s_ = sparse.COO(sdf.sparse.to_coo())
+    
+    print()
+    print('Output')
+    print('s_all, binary features    :\t', s_.shape)
+    return s_, feature_names_all, dict(discretization_bins)
+
+##########
+
+def map_time_series_features(df_time_series, dtypes, args):
+    N, L = args.N, args.L
+    
+    df_time_series = df_time_series.dropna(axis='columns', how='all').sort_index()
+
+    print('Discretizing features...', flush=True)
+    ts_mask = select_dtype(df_time_series, 'mask', dtypes)
+    ts_mixed = select_dtype(df_time_series, '~mask', dtypes)
+    assert len(ts_mixed.columns) + len(ts_mask.columns) == len(df_time_series.columns)
+    ts_feature_mask = ts_mask.astype(int)
+    ts_mixed_cols = [ts_mixed[col] for col in ts_mixed.columns]
+    
+    print()
+    if args.binarize:
+        dtype = int
+        print('Processing', len(ts_mixed_cols), 'non-boolean variable columns...')
+
+        print('    Binning numeric variables by quintile...')
+        print('    Converting variables to binary features')
+        if parallel:
+#             out = Parallel(n_jobs=n_jobs, verbose=10)( # Need to share global variables
+#                 delayed(smart_qcut_dummify)(col_data, q=5) for col_data in ts_mixed_cols
+#             )
+#             out = [smart_qcut_dummify(col_data, q=5) for col_data in tqdm(ts_mixed_cols)]
+    
+            pool = multiprocessing.Pool(multiprocessing.cpu_count())
+            out_0 = pool.starmap(smart_qcut, ts_mixed_cols)
+            cols_data, discretization_bins = zip(out_0)
+            out = pool.starmap(dummify, cols_data)
+            pool.close()
+            pool.join()
+        else:
+            out = [smart_qcut_dummify(col_data, q=5) for col_data in tqdm(ts_mixed_cols)]
+    else:
+        raise NotImplemented
+        dtype = float
+        df = ts_mixed.copy()
+        
+        # Split a mixed column into numeric and string columns
+        for col in df.columns:
+            col_data = df[col]
+            col_is_numeric = [is_numeric(v) for v in col_data if not pd.isnull(v)]
+            if not all(col_is_numeric) and any(col_is_numeric): # have mixed type values
+                numeric_mask = col_data.apply(is_numeric)
+                df[col+'_str'] = df[col].copy()
+                df.loc[~numeric_mask, col] = np.nan
+                df.loc[numeric_mask, col+'_str'] = np.nan
+        
+        ts_mixed_cols = [df[col] for col in df.columns]
+        
+        print('Discretizing categorical features...')
+        if parallel:
+#             out = Parallel(n_jobs=n_jobs, verbose=10)( # Need to share global variables?
+#                 delayed(smart_dummify_impute)(col_data) for col_data in ts_mixed_cols
+#             )
+            out = [smart_dummify_impute(col_data) for col_data in tqdm(ts_mixed_cols)]
+        else:
+            out = [smart_dummify_impute(col_data) for col_data in tqdm(ts_mixed_cols)]
+    
+    out = [ts_feature_mask, *out]
+    D_all = sum(len(df_i.columns) for df_i in out)
+    X_all_feature_names = np.asarray(sum([list(df_i.columns) for df_i in out], []))
+    X_dense = np.concatenate([df_i.values for df_i in out], axis=1).astype(dtype)
+    X_all = sparse.COO(X_dense)
+    
+    print('Finished discretizing features', flush=True)
+    assert X_all.shape[0] == N * L
+    X_all = X_all.reshape((N, L, D_all))
+    
+    print()
+    print('Output')
+    print('X_all: shape={}, density={:.3f}'.format(X_all.shape, X_all.density))
+    return X_all, X_all_feature_names, dict(discretization_bins)
+
+
+def print_metadata():
+    # Print metadata
+    print('DONE: Transforming each example...')
+    ## Freq: Count missing entries using mask
+    ts_mask = df_time_series[[col for col in df_time_series if col.endswith('_mask')]]
+    ts_mask.columns = [col.replace('_mask', '') for col in ts_mask.columns]
+    print('(freq) number of missing entries :\t', 
+          '{} out of {}={} total'.format(
+              (1-ts_mask).astype(int).sum().sum(), 
+              '\u00D7'.join(str(i) for i in [N,L,ts_mask.shape[1]]), ts_mask.size))
+    
+    ## Freq: Count imputed entries using mask and dt
+    ts_delta_time = df_time_series[[col for col in df_time_series if col.endswith('_delta_time')]]
+    ts_delta_time.columns = [col.replace('_delta_time', '') for col in ts_delta_time.columns]
+    
+    imputed = (1-ts_mask).astype(bool) & (ts_delta_time > 0)
+    print('(freq) number of imputed entries :\t', 
+          '{}'.format(imputed.sum().sum(), ts_delta_time.size))
+    imputed.sum().rename('count').to_csv(data_path + '/' + 'freq_imputed.csv')
+    
+    not_imputed = (1-ts_mask).astype(bool) & (ts_delta_time == 0)
+    print('(freq) number of not imputed entries :\t', 
+          '{}'.format(not_imputed.sum().sum(), ts_delta_time.size))
+    not_imputed.sum().rename('count').to_csv(data_path + '/' + 'freq_not_imputed.csv')
+    
+    ## Non-Freq: Count missing entries
+    non_freq_cols = sorted([c + '_value' for c in set(variables) - set(variables_num_freq)])
+    non_freqs = df_time_series[non_freq_cols]
+    print('(non-freq) number of missing entries :\t',
+          '{} out of {}={} total'.format(
+              non_freqs.isna().sum().sum(), 
+              '\u00D7'.join(str(i) for i in [N,L,non_freqs.shape[1]]), non_freqs.size))
+
+
+def post_filter_time_series(X_all, feature_names_all, threshold, args):
+    N, L = args.N, args.L
+    assert X_all.shape[0] == N
+    assert X_all.shape[1] == L
+#     assert X_all.dtype == int
+    start_time = time.time()
+    
+    X0 = X_all
+    feature_names_0 = feature_names_all
+    print('Original :', len(feature_names_0))
+    
+    ## Remove nearly-constant features (with low variance)
+    sel_const = FrequencyThreshold_temporal(threshold=threshold, L=L)
+    sel_const.fit(X0.reshape((N*L, -1)))
+    m_ts_const = sel_const.get_support()
+    assert len(m_ts_const) == X0.shape[-1]
+    X1 = X0[:, :, m_ts_const]
+    feature_names_1 = feature_names_0[m_ts_const]
+    print('Nearly-constant:', len(feature_names_0) - len(feature_names_1))
+    print('*** time: ', time.time() - start_time)
+    
+    ## Keep only first of pairwise perfectly correlated features
+    sel_ts_corr = CorrelationSelector()
+    sel_ts_corr.fit(X1.reshape((N*L, -1)))
+    m_ts_corr = sel_ts_corr.get_support()
+    assert len(m_ts_corr) == X1.shape[-1]
+    X2 = X1[:, :, m_ts_corr]
+    feature_names_2 = feature_names_1[m_ts_corr]
+    feature_aliases = sel_ts_corr.get_feature_aliases(feature_names_1)
+    print('Correlated     :', len(feature_names_1) - len(feature_names_2))
+    print('*** time: ', time.time() - start_time)
+    
+    X = sparse.COO(X2)
+    feature_names = feature_names_2
+    assert X.shape == (N, L, len(feature_names))
+    
+    ## Save output
+    print()
+    print('Output')
+    print('X: shape={}, density={:.3f}'.format(X.shape, X.density))
+    
+    return X, feature_names, feature_aliases
+
+
+def pre_filter(df, threshold, df_population, args):
+    T = int(args.T)
+    theta_1 = args.theta_1
+    df_population = args.df_population
+    
+    # Remove rows not in population
+    print('Remove rows not in population')
+    df = df[df['ID'].isin(df_population.index)]
+    
+    # Remove rows with t outside of [0, T)
+    print('Remove rows with t outside of [0, {}]'.format(T))
+    df = df[pd.isnull(df[t_col]) | ((0 <= df[t_col]) & (df[t_col] < T))]
+    
+    # Data tables should not contain duplicate rows
+    # Check for inconsistencies
+    dups = df.duplicated(subset=[ID_col, t_col, var_col], keep=False)
+    if any(dups):
+        print(df[dups].head())
+        raise Exception('Inconsistent values recorded')
+    
+    # Remove variables that occur too rarely as defined by the threshold
+    print('Remove rare variables (<= {})'.format(threshold))
+    
+    ## Calculate overall occurrence rate of each variable based on IDs
+    df_count = calculate_variable_counts(df, df_population) # (N x |var|) table of counts
+    df_bool = df_count.astype(bool) # convert counts to boolean
+    
+    ## Keep variables that are recorded for more than threshold fraction of IDs
+    variables_keep = df_bool.columns[df_bool.mean(axis=0) > threshold]
+    df_out = df[df[var_col].isin(variables_keep)]
+    assert set(variables_keep) == set(df_out[var_col].unique())
+    
+    variables = sorted(df_bool.columns)
+    variables_remove = sorted(set(variables) - set(variables_keep))
+    print('Total variables     :', len(variables))
+    print('Rare variables      :', len(variables_remove))
+    print('Remaining variables :', len(variables_keep))
+    print('# rows (original)   :', len(df))
+    print('# rows (filtered)   :', len(df_out))
+    return df_out
+
+
+
+######
+# Time-invariant routines
+######
+def transform_time_invariant_table(df_in, df_population):
+    df_in = df_in.copy()
+    
+    # Recorded Value (np.nan if not recorded)
+    df_value = pd.pivot_table(df_in, val_col, ID_col, var_col, 'last', np.nan)
+    df_value = df_value.reindex(index=df_population.index, fill_value=np.nan)
+    df_value.columns = [str(col) + '_value' for col in df_value.columns]
+    
+    print('(N \u00D7 ^d) table            :\t', df_value.shape)
+    print('number of missing entries :\t', '{} out of {} total'.format(df_value.isna().sum().sum(), df_value.size))
+    return df_value
+
+def smart_qcut_dummify_5(x):
+    return smart_qcut_dummify(x, q=5)
+
+
+def post_filter(s_, s_feature_names_all, threshold):
+    # Filter features (optional)
+    assert s_.shape[1] == len(s_feature_names_all)
+    feature_names_0 = s_feature_names_all
+    s0 = s_.to_scipy_sparse()
+    print('Original       :', len(feature_names_0))
+    
+    ## Remove nearly-constant features (with low variance)
+    ## a binary feature is removed if =0 (or =1) for >th fraction of examples
+    ## i.e., variance <= (th * (1 - th))
+    sel_rare = VarianceThreshold(threshold=(threshold * (1 - threshold)))
+    s1 = sel_rare.fit_transform(s0)
+    feature_names_1 = feature_names_0[sel_rare.get_support()]
+    print('Nearly-constant:', len(feature_names_0) - len(feature_names_1))
+    
+    ## Keep only first of pairwise perfectly correlated features
+    sel_corr = CorrelationSelector()
+    s2 = sel_corr.fit_transform(s1)
+    feature_names_2 = feature_names_1[sel_corr.get_support()]
+    feature_aliases = sel_corr.get_feature_aliases(feature_names_1)
+    print('Correlated     :', len(feature_names_1) - len(feature_names_2))
+    
+    s = sparse.COO(s2)
+    feature_names = feature_names_2
+    assert s.shape[1] == len(feature_names)
+    
+    return s, feature_names, feature_aliases
+
+
+def process_time_invariant(df_data_time_invariant, args):
+    data_path = args.data_path
+    df_population = args.df_population
+    theta_2 = args.theta_2
+    
+    print_header('2-A) Transform time-invariant data', char='-')
+    dir_path = data_path + '/'
+    start_time = time.time()
+
+    ## Create Nxd^ table
+    df_time_invariant = transform_time_invariant_table(df_data_time_invariant, df_population)
+    print('Time elapsed: %f seconds' % (time.time() - start_time))
+
+    ## Discretize
+    s_all, s_all_feature_names, s_discretization = map_time_invariant_features(df_time_invariant, args.binarize)
+    sparse.save_npz(dir_path + 's_all.npz', s_all)
+    with open(dir_path + 's_all.feature_names.json', 'w') as f:
+        json.dump(list(s_all_feature_names), f, sort_keys=True)
+    print('Time elapsed: %f seconds' % (time.time() - start_time))
+    json.dump(s_discretization, open(dir_path + 'discretization.json', 'w'))
+    
+    print_header('3-A) Post-filter time-invariant data', char='-')
+    
+    ## Filter
+    s, s_feature_names, s_feature_aliases = post_filter(s_all, s_all_feature_names, theta_2)
+    print('Time elapsed: %f seconds' % (time.time() - start_time))
+    
+    ## Save output
+    print()
+    print('Output')
+    print('s: shape={}, density={:.3f}'.format(s.shape, s.density))
+    sparse.save_npz(dir_path + 's.npz', s)
+    
+    with open(dir_path + 's.feature_names.json', 'w') as f:
+        json.dump(list(s_feature_names), f, sort_keys=True)
+    with open(dir_path + 's.feature_aliases.json', 'w') as f:
+        json.dump(s_feature_aliases, f, sort_keys=True)
+    
+    print('Total time: %f seconds' % (time.time() - start_time))
+    print('', flush=True)
+    return s, s_feature_names, s_feature_aliases
+