"""
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