data_utils.py

import numpy as np
import pandas as pd
import pdb
import re
from time import time
import json
import random

import model

from scipy.spatial.distance import pdist, squareform
from scipy.stats import multivariate_normal, invgamma, mode
from scipy.special import gamma
from scipy.misc import imresize
from functools import partial
from math import ceil

from sklearn.metrics.pairwise import rbf_kernel
from sklearn.preprocessing import MinMaxScaler

def get_batch(samples, batch_size, batch_idx, labels=None):
    start_pos = batch_idx * batch_size
    end_pos = start_pos + batch_size
    if labels is None:
        return samples[start_pos:end_pos], None
    else:
        if type(labels) == tuple: # two sets of labels
            assert len(labels) == 2
            return samples[start_pos:end_pos], labels[0][start_pos:end_pos], labels[1][start_pos:end_pos]
        else:
            assert type(labels) == np.ndarray
            return samples[start_pos:end_pos], labels[start_pos:end_pos]

def normalise_data(train, vali, test, low=-1, high=1):
    """ Apply some sort of whitening procedure
    """
    # remember, data is num_samples x seq_length x signals
    # whiten each signal - mean 0, std 1
    mean = np.mean(np.vstack([train, vali]), axis=(0, 1))
    std = np.std(np.vstack([train-mean, vali-mean]), axis=(0, 1))

    normalised_train = (train - mean)/std
    normalised_vali = (vali - mean)/std
    normalised_test = (test - mean)/std
#    normalised_data = data - np.nanmean(data, axis=(0, 1))
#    normalised_data /= np.std(data, axis=(0, 1))    

#    # normalise samples to be between -1 and +1
    # normalise just using train and vali
#    min_val = np.nanmin(np.vstack([train, vali]), axis=(0, 1))
#    max_val = np.nanmax(np.vstack([train, vali]), axis=(0, 1))
#    
#    normalised_train = (train - min_val)/(max_val - min_val)
#    normalised_train = (high - low)*normalised_train + low
#
#    normalised_vali = (vali - min_val)/(max_val - min_val)
#    normalised_vali = (high - low)*normalised_vali + low
#
#    normalised_test = (test - min_val)/(max_val - min_val)
#    normalised_test = (high - low)*normalised_test + low
    return normalised_train, normalised_vali, normalised_test

def scale_data(train, vali, test, scale_range=(-1, 1)):
    signal_length = train.shape[1]
    num_signals = train.shape[2]
    # reshape everything
    train_r = train.reshape(-1, signal_length*num_signals)
    vali_r = vali.reshape(-1, signal_length*num_signals)
    test_r = test.reshape(-1, signal_length*num_signals)
    # fit scaler using train, vali
    scaler = MinMaxScaler(feature_range=scale_range).fit(np.vstack([train_r, vali_r]))
    # scale everything
    scaled_train = scaler.transform(train_r).reshape(-1, signal_length, num_signals)
    scaled_vali = scaler.transform(vali_r).reshape(-1, signal_length, num_signals)
    scaled_test = scaler.transform(test_r).reshape(-1, signal_length, num_signals)
    return scaled_train, scaled_vali, scaled_test

def split(samples, proportions, normalise=False, scale=False, labels=None, random_seed=None):
    """
    Return train/validation/test split.
    """
    if random_seed != None:
        random.seed(random_seed)
        np.random.seed(random_seed)
    assert np.sum(proportions) == 1
    n_total = samples.shape[0]
    n_train = ceil(n_total*proportions[0])
    n_test = ceil(n_total*proportions[2])
    n_vali = n_total - (n_train + n_test)
    # permutation to shuffle the samples
    shuff = np.random.permutation(n_total)
    train_indices = shuff[:n_train]
    vali_indices = shuff[n_train:(n_train + n_vali)]
    test_indices = shuff[(n_train + n_vali):]
    # TODO when we want to scale we can just return the indices
    assert len(set(train_indices).intersection(vali_indices)) == 0
    assert len(set(train_indices).intersection(test_indices)) == 0
    assert len(set(vali_indices).intersection(test_indices)) == 0
    # split up the samples
    train = samples[train_indices]
    vali = samples[vali_indices]
    test = samples[test_indices]
    # apply the same normalisation scheme to all parts of the split
    if normalise:
        if scale: raise ValueError(normalise, scale)        # mutually exclusive
        train, vali, test = normalise_data(train, vali, test)
    elif scale:
        train, vali, test = scale_data(train, vali, test)
    if labels is None:
        return train, vali, test
    else:
        print('Splitting labels...')
        if type(labels) == np.ndarray:
            train_labels = labels[train_indices]
            vali_labels = labels[vali_indices]
            test_labels = labels[test_indices]
            labels_split = [train_labels, vali_labels, test_labels]
        elif type(labels) == dict:
            # more than one set of labels!  (weird case)
            labels_split = dict()
            for (label_name, label_set) in labels.items():
                train_labels = label_set[train_indices]
                vali_labels = label_set[vali_indices]
                test_labels = label_set[test_indices]
                labels_split[label_name] = [train_labels, vali_labels, test_labels]
        else:
            raise ValueError(type(labels))
        return train, vali, test, labels_split

def get_data(data_type, data_options=None):
    """
    Helper/wrapper function to get the requested data.
    """
    labels = None
    pdf = None
    if data_type == 'load':
        data_dict = np.load(data_options).item()
        samples = data_dict['samples']
        pdf = data_dict['pdf']
        labels = data_dict['labels']
    elif data_type == 'sine':
        samples = sine_wave(**data_options)
    elif data_type == 'mnist':
        if data_options['full_mnist']:
            samples, labels = mnist()
        else:
            samples, labels = load_resized_mnist_0_5(14)
    elif data_type == 'gp_rbf':
        print(data_options)
        samples, pdf = GP(**data_options, kernel='rbf')
    elif data_type == 'linear':
        samples, pdf = linear(**data_options)
    elif data_type == 'eICU_task':
        samples, labels = eICU_task()
    elif data_type == 'resampled_eICU':
        samples, labels = resampled_eICU(**data_options)
    else:
        raise ValueError(data_type)
    print('Generated/loaded', len(samples), 'samples from data-type', data_type)
    return samples, pdf, labels

def make_predict_labels(samples, labels):
    """ Given two dictionaries of samples, labels (already normalised, split etc)
        append the labels on as additional signals in the data
    """
    print('Appending label to samples')
    assert not labels is None
    if len(labels['train'].shape) > 1:
        num_labels = labels['train'].shape[1]
    else:
        num_labels = 1
    seq_length = samples['train'].shape[1]
    num_signals = samples['train'].shape[2]
    new_samples = dict()
    new_labels = dict()
    for (k, X) in samples.items():
        num_samples = X.shape[0]
        lab = labels[k]
        # slow code because i am sick and don't want to try to be smart
        new_X = np.zeros(shape=(num_samples, seq_length, num_signals + num_labels))
        for row in range(num_samples):
            new_X[row, :, :] = np.hstack([X[row, :, :], np.array(seq_length*[(2*lab[row]-1).reshape(num_labels)])])
        new_samples[k] = new_X
        new_labels[k] = None
    return new_samples, new_labels

# --- specific data-types --- #

def eICU_task(predict_label=False):
    """
    Load the eICU data for the extreme-value prediction task
    """
    path = 'REDACTED'
    data = np.load(path).item()
    # convert it into similar format
    labels = {'train': data['Y_train'], 'vali': data['Y_vali'], 'test': data['Y_test']}
    samples = {'train': data['X_train'], 'vali': data['X_vali'], 'test': data['X_test']}
    # reshape
    for (k, X) in samples.items():
        samples[k] = X.reshape(-1, 16, 4)
    return samples, labels

def mnist(randomize=False):
    """ Load and serialise """
    try:
        train = np.load('./data/mnist_train.npy')
        print('Loaded mnist from .npy')
    except IOError:
        print('Failed to load from .npy, loading from csv')
        # read from the csv
        train = np.loadtxt(open('./data/mnist_train.csv', 'r'), delimiter=',')
        # scale samples from 0 to 1
        train[:, 1:] /= 255
        # scale from -1 to 1
        train[:, 1:] = 2*train[:, 1:] - 1
        # save to the npy
        np.save('./data/mnist_train.npy', train)
    # the first column is labels, kill them
    labels = train[:, 0]
    samples = train[:, 1:]
    if randomize:
        # not needed for GAN experiments...
        print('Applying fixed permutation to mnist digits.')
        fixed_permutation = np.random.permutation(28*28)
        samples = train[:, fixed_permutation]
    samples = samples.reshape(-1, 28*28, 1)     # add redundant additional signals
    return samples, labels


def load_resized_mnist_0_5(new_size, randomize=False):
    """ Load resised mnist digits from 0 to 5 """

    samples, labels = mnist()
    print('Resizing...')
    samples = samples[np.in1d(labels,[0,1,2,3,4,5])]
    labels = labels[np.in1d(labels,[0,1,2,3,4,5])]
    if new_size != 28:
        resized_imgs = [imresize(img.reshape([28,28]), [new_size,new_size], interp='lanczos').ravel()[np.newaxis].T 
                        for img in samples]
        resized_imgs = np.array(resized_imgs)
        resized_imgs = resized_imgs.astype(float)
        resized_imgs /= 255.0
        resized_imgs = 2*resized_imgs - 1
        np.save('./data/resized_mnist_1_5_samples.npy', resized_imgs)
        np.save('./data/resized_mnist_1_5_labels.npy', labels)
        return resized_imgs, labels
    else:
        return samples, labels


def load_resized_mnist(new_size, from_to_digits=(0,2), randomize=False):
    """ Load resised mnist digits from 0 to 5 """

    samples, labels = mnist()
    print('Resizing...')
    samples = samples[np.in1d(labels,np.arange(from_to_digits[0], from_to_digits[1]+1))]
    labels = labels[np.in1d(labels,np.arange(from_to_digits[0], from_to_digits[1]+1))]
    if new_size != 28:
        resized_imgs = [imresize(img.reshape([28,28]), [new_size,new_size], interp='lanczos').ravel()[np.newaxis].T
                        for img in samples]
        resized_imgs = np.array(resized_imgs)
        resized_imgs = resized_imgs.astype(float)
        resized_imgs /= 255.0
        resized_imgs = 2*resized_imgs - 1
        np.save('./data/resized_mnist_'+ str(from_to_digits[0]) + '_' + str(from_to_digits[1]) + '_5_samples.npy', resized_imgs)
        np.save('./data/resized_mnist_'+ str(from_to_digits[0]) + '_' + str(from_to_digits[1]) + '_labels.npy', labels)
        return resized_imgs, labels
    else:
        return samples, labels

def resampled_eICU(seq_length=16, resample_rate_in_min=15, 
        variables=['sao2', 'heartrate', 'respiration', 'systemicmean'], **kwargs):
    """
    Note: resampling rate is 15 minutes
    """
    print('Getting resampled eICU data')
    base_dir = 'REDACTED'
    if not resample_rate_in_min == 15:
        raise NotImplementedError(resample_rate_in_min)
    try:
        data = np.load(base_dir + 'eICU_' + str(resample_rate_in_min) + '.npy').item()
        samples = data['samples']
        pids = data['pids']
        print('Loaded from file!')
        return samples, pids
    except FileNotFoundError:
        pass
    resampled_data_path = base_dir + 'resampled_pats_' + str(resample_rate_in_min) + 'min.csv'
    resampled_pids_path = base_dir + 'cohort_complete_resampled_pats_' + str(resample_rate_in_min) + 'min.csv'

    if not os.path.isfile(resampled_data_path):
        generate_eICU_resampled_patients(resample_factor_in_min=resample_rate_in_min, upto_in_minutes=None)
        get_cohort_of_complete_downsampled_patients(time_in_hours=1.5*resample_rate_in_min*seq_length, resample_factor_in_min=resample_rate_in_min)

    pids = set(np.loadtxt(resampled_pids_path, dtype=int))
    df = pd.read_csv(resampled_data_path)
    # restrict to variables
    df_restricted = df.loc[:, variables + ['offset', 'pid']]
    # restrict to patients in the "good list"
    df_restricted = df_restricted.where(df_restricted.pid.isin(pids)).dropna()
    # assert no negative offsets
    assert np.all(df_restricted.offset >= 0)
    # restrict to 1.5 time the region length
#    df_restricted = df_restricted.loc[np.all([df_restricted.offset <= 1.5*resample_rate_in_min*seq_length, df_restricted.offset >= 0], axis=0), :]
    df_restricted = df_restricted.loc[df_restricted.offset <= 1.5*resample_rate_in_min*seq_length, :]
    # for each patient, return the first seq_length observations
    patient_starts = df_restricted.groupby('pid').head(seq_length)
    n_pats_prefilter = len(set(patient_starts.pid))
    # filter out patients who have fewer than seq_length observations
    patient_starts = patient_starts.groupby('pid').filter(lambda x: x.pid.count() == seq_length)
    n_pats_postfilter = len(set(patient_starts.pid))
    print('Removed', n_pats_prefilter - n_pats_postfilter, 'patients with <', seq_length, 'observations in the first', 1.5*resample_rate_in_min*seq_length, 'minutes, leaving', n_pats_postfilter, 'patients remaining.')
    # convert to samples - shape is [n_pats, seq_length, num_signals]
    n_patients = n_pats_postfilter
    num_signals = len(variables)
    samples = np.empty(shape=(n_patients, seq_length, num_signals))
    pats_grouped = patient_starts.groupby('pid')
    pids = []
    for (i, patient) in enumerate(pats_grouped.groups):
        samples[i, :, :] = pats_grouped.get_group(patient).loc[:, variables].values
        pids.append(patient)
    assert i == n_patients - 1
    assert np.mean(np.isnan(samples) == 0)
    np.save(base_dir + 'eICU_' + str(resample_rate_in_min) + '.npy', {'samples': samples, 'pids': pids})
    print('Saved to file!')
    return samples, pids

def sine_wave(seq_length=30, num_samples=28*5*100, num_signals=1, 
        freq_low=1, freq_high=5, amplitude_low = 0.1, amplitude_high=0.9, **kwargs):
    ix = np.arange(seq_length) + 1
    samples = []
    for i in range(num_samples):
        signals = []
        for i in range(num_signals):
            f = np.random.uniform(low=freq_high, high=freq_low)     # frequency
            A = np.random.uniform(low=amplitude_high, high=amplitude_low)        # amplitude
            # offset
            offset = np.random.uniform(low=-np.pi, high=np.pi)
            signals.append(A*np.sin(2*np.pi*f*ix/float(seq_length) + offset))
        samples.append(np.array(signals).T)
    # the shape of the samples is num_samples x seq_length x num_signals
    samples = np.array(samples)
    return samples

def periodic_kernel(T, f=1.45/30, gamma=7.0, A=0.1):
    """ 
    Calculates periodic kernel between all pairs of time points (there
    should be seq_length of those), returns the Gram matrix.
    f is frequency - higher means more peaks
    gamma is a scale, smaller makes the covariance peaks shallower (smoother)
    
    Heuristic for non-singular rbf:
        periodic_kernel(np.arange(len), f=1.0/(0.79*len), A=1.0, gamma=len/4.0)
    """
    dists = squareform(pdist(T.reshape(-1, 1)))
    cov = A*np.exp(-gamma*(np.sin(2*np.pi*dists*f)**2))
    return cov


def GP(seq_length=30, num_samples=28*5*100, num_signals=1, scale=0.1, kernel='rbf', **kwargs):
    # the shape of the samples is num_samples x seq_length x num_signals
    samples = np.empty(shape=(num_samples, seq_length, num_signals))
    #T = np.arange(seq_length)/seq_length    # note, between 0 and 1
    T = np.arange(seq_length)    # note, not between 0 and 1
    if kernel == 'periodic':
        cov = periodic_kernel(T)
    elif kernel =='rbf':
        cov = rbf_kernel(T.reshape(-1, 1), gamma=scale)
    else:
        raise NotImplementedError
    # scale the covariance
    cov *= 0.2
    # define the distribution
    mu = np.zeros(seq_length)
    print(np.linalg.det(cov))
    distribution = multivariate_normal(mean=np.zeros(cov.shape[0]), cov=cov) 
    pdf = distribution.logpdf
    # now generate samples
    for i in range(num_signals):
        samples[:, :, i] = distribution.rvs(size=num_samples)
    return samples, pdf

def linear_marginal_likelihood(Y, X, a0, b0, mu0, lambda0, log=True, **kwargs):
    """
    Marginal likelihood for linear model.
    See https://en.wikipedia.org/wiki/Bayesian_linear_regression pretty much
    """
    seq_length = Y.shape[1]         # note, y is just a line (one channel) TODO
    n = seq_length
    an = a0 + 0.5*n
    XtX = np.dot(X.T, X)
    lambdan = XtX + lambda0
    prefactor = (2*np.pi)**(-0.5*n)
    dets = np.sqrt(np.linalg.det(lambda0)/np.linalg.det(lambdan))
    marginals = np.empty(Y.shape[0])
    for (i, y) in enumerate(Y):
        y_reshaped = y.reshape(seq_length)
        betahat = np.dot(np.linalg.inv(XtX), np.dot(X.T, y_reshaped))
        mun = np.dot(np.linalg.inv(lambdan), np.dot(XtX, betahat) + np.dot(lambda0, mu0))
        bn = b0 + 0.5*(np.dot(y_reshaped.T, y_reshaped) + np.dot(np.dot(mu0.T, lambda0), mu0) - np.dot(np.dot(mun.T, lambdan), mun))
        bs = (b0**a0)/(bn**an)
        gammas = gamma(an)/gamma(a0)
        marginals[i] = prefactor*dets*bs*gammas
    if log:
        marginals = np.log(marginals)
    return marginals

def linear(seq_length=30, num_samples=28*5*100, a0=10, b0=0.01, k=2, **kwargs):
    """
    Generate data from linear trend from probabilistic model.

    The invgamma function in scipy corresponds to wiki defn. of inverse gamma:
        scipy a = wiki alpha = a0
        scipy scale = wiki beta = b0 

    k is the number of regression coefficients (just 2 here, slope and intercept)
    """
    T = np.zeros(shape=(seq_length, 2))
    T[:, 0] = np.arange(seq_length)
    T[:, 1] = 1                         # equivalent to X
    lambda0 = 0.01*np.eye(k)     # diagonal covariance for beta
    y = np.zeros(shape=(num_samples, seq_length, 1))
    sigmasq = invgamma.rvs(a=a0, scale=b0, size=num_samples)
    increasing = np.random.choice([-1, 1], num_samples)     # flip slope
    for n in range(num_samples):
        sigmasq_n = sigmasq[n]
        offset = np.random.uniform(low=-0.5, high=0.5)     # todo limits
        mu0 = np.array([increasing[n]*(1.0-offset)/seq_length, offset])
        beta = multivariate_normal.rvs(mean=mu0, cov=sigmasq_n*lambda0)
        epsilon = np.random.normal(loc=0, scale=np.sqrt(sigmasq_n), size=seq_length)
        y[n, :, :] = (np.dot(T, beta) + epsilon).reshape(seq_length, 1)
    marginal = partial(linear_marginal_likelihood, X=T, a0=a0, b0=b0, mu0=mu0, lambda0=lambda0)
    samples = y
    pdf = marginal
    return samples, pdf

def changepoint_pdf(Y, cov_ms, cov_Ms):
    """
    """
    seq_length = Y.shape[0]
    logpdf = []
    for (i, m) in enumerate(range(int(seq_length/2), seq_length-1)):
        Y_m = Y[:m, 0]
        Y_M = Y[m:, 0]
        M = seq_length - m
        # generate mean function for second part
        Ymin = np.min(Y_m)
        initial_val = Y_m[-1]
        if Ymin > 1:
            final_val = (1.0 - M/seq_length)*Ymin
        else:
            final_val = (1.0 + M/seq_length)*Ymin
        mu_M = np.linspace(initial_val, final_val, M)
        # ah yeah
        logpY_m = multivariate_normal.logpdf(Y_m, mean=np.zeros(m), cov=cov_ms[i])
        logpY_M = multivariate_normal.logpdf(Y_M, mean=mu_M, cov=cov_Ms[i])
        logpdf_m = logpY_m + logpY_M
        logpdf.append(logpdf_m)
    return logpdf

def changepoint_cristobal(seq_length=30, num_samples=28*5*100):
    """
    Porting Cristobal's code for generating data with a changepoint.
    """
    raise NotImplementedError

    basal_values_signal_a = np.random.randn(n_samples) * 0.33
    trends_seed_a = np.random.randn(n_samples) * 0.005
    trends = np.array([i*trends_seed_a for i in range(51)[1:]]).T
    signal_a = (basal_values_signal_a + trends.T).T
    time_noise = np.random.randn(n_samples, n_steps) * 0.01
    signal_a = time_noise + signal_a

    basal_values_signal_b = np.random.randn(n_samples) * 0.33
    trends_seed_b = np.random.randn(n_samples) * 0.005
    trends = np.array([i*trends_seed_b for i in range(51)[1:]]).T
    signal_b = (basal_values_signal_b + trends.T).T
    time_noise = np.random.randn(n_samples, n_steps) * 0.01
    signal_b = time_noise + signal_b

    signal_a = np.clip(signal_a, -1, 1)
    signal_b = np.clip(signal_b, -1, 1)

    # the change in the trend is based on the top extreme values of each
    # signal in the first half
    time_steps_until_change = np.max(np.abs(signal_a), axis=1) + np.max(np.abs(signal_b), axis=1)*100
    # noise added to the starting point
    time_steps_until_change += np.random.randn(n_samples) * 5
    time_steps_until_change = np.round(time_steps_until_change)
    time_steps_until_change = np.clip(time_steps_until_change, 0, n_steps-1)
    time_steps_until_change = n_steps - 1 - time_steps_until_change

    trends = np.array([i*trends_seed_a for i in range(101)[51:]]).T
    signal_a_target = (basal_values_signal_a + trends.T).T
    time_noise = np.random.randn(n_samples, n_steps) * 0.01
    signal_a_target = time_noise + signal_a_target

    trends = np.array([i*trends_seed_b for i in range(101)[51:]]).T
    signal_b_target = (basal_values_signal_b + trends.T).T
    time_noise = np.random.randn(n_samples, n_steps) * 0.01
    signal_b_target = time_noise + signal_b_target
    
    signal_multipliers = []
    for ts in time_steps_until_change:
        signal_multiplier = []
        if ts > 0:
            for i in range(int(ts)):
                signal_multiplier.append(1)
            i += 1
        else:
            i = 0
        multiplier = 1.25
        while(i<n_steps):
            signal_multiplier.append(multiplier)
            multiplier += 0.25
            i+=1
        signal_multipliers.append(signal_multiplier)
    signal_multipliers = np.array(signal_multipliers)

    for s_idx, signal_choice in enumerate(basal_values_signal_b > basal_values_signal_a):
        if signal_choice == False:
            signal_a_target[s_idx] *= signal_multipliers[s_idx]
        else:
            signal_b_target[s_idx] *= signal_multipliers[s_idx]

    signal_a_target = np.clip(signal_a_target, -1, 1)
    signal_b_target = np.clip(signal_b_target, -1, 1)

    # merging signals
    signal_a = np.swapaxes(signal_a[np.newaxis].T, 0, 1)
    signal_b = np.swapaxes(signal_b[np.newaxis].T, 0, 1)
    signal_a_target = np.swapaxes(signal_a_target[np.newaxis].T, 0, 1)
    signal_b_target = np.swapaxes(signal_b_target[np.newaxis].T, 0, 1)
    input_seqs = np.dstack((signal_a,signal_b))
    target_seqs = np.dstack((signal_a_target,signal_b_target))
    return False

def changepoint(seq_length=30, num_samples=28*5*100):
    """
    Generate data from two GPs, roughly speaking.
    The first part (up to m) is as a normal GP.
    The second part (m to end) has a linear downwards trend conditioned on the
    first part.
    """
    print('Generating samples from changepoint...')
    T = np.arange(seq_length)
    # sample breakpoint from latter half of sequence
    m_s = np.random.choice(np.arange(int(seq_length/2), seq_length-1), size=num_samples)
    samples = np.zeros(shape=(num_samples, seq_length, 1))
    # kernel parameters and stuff
    gamma=5.0/seq_length
    A = 0.01
    sigmasq = 0.8*A
    lamb = 0.0  # if non-zero, cov_M risks not being positive semidefinite...
    kernel = partial(rbf_kernel, gamma=gamma)
    # multiple values per m
    N_ms = []
    cov_ms = []
    cov_Ms = []
    pdfs = []
    for m in range(int(seq_length/2), seq_length-1):
        # first part
        M = seq_length - m
        T_m = T[:m].reshape(m, 1)
        cov_m = A*kernel(T_m.reshape(-1, 1), T_m.reshape(-1, 1))
        cov_ms.append(cov_m)
        # the second part 
        T_M = T[m:].reshape(M, 1)
        cov_mM = kernel(T_M.reshape(-1, 1), T_m.reshape(-1, 1))
        cov_M = sigmasq*(np.eye(M) - lamb*np.dot(np.dot(cov_mM, np.linalg.inv(cov_m)), cov_mM.T))
        cov_Ms.append(cov_M)
    for n in range(num_samples):
        m = m_s[n]
        M = seq_length-m
        # sample the first m
        cov_m = cov_ms[m - int(seq_length/2)]
        Xm = multivariate_normal.rvs(cov=cov_m)
        # generate mean function for second
        Xmin = np.min(Xm)
        initial_val = Xm[-1]
        if Xmin > 1:
            final_val = (1.0 - M/seq_length)*Xmin
        else:
            final_val = (1.0 + M/seq_length)*Xmin
        mu_M = np.linspace(initial_val, final_val, M)
        # sample the rest
        cov_M = cov_Ms[m -int(seq_length/2)]
        XM = multivariate_normal.rvs(mean=mu_M, cov=cov_M)
        # combine the sequence
        # NOTE: just one dimension
        samples[n, :, 0] = np.concatenate([Xm, XM])
    pdf = partial(changepoint_pdf, cov_ms=cov_ms, cov_Ms=cov_Ms)
    return samples, pdf, m_s

def generate_eICU_resampled_patients(resample_factor_in_min=15, 
        upto_in_minutes=None):
    """
    Generates a dataframe with resampled patients. One sample every 15 miniutes

    """
    data_path = 'REDACTED'
    pids = pd.read_csv(data_path)
    pids = np.array(pids).flatten()

    eICU_dir = 'REDACTED'
    variables = ['sao2', 'heartrate', 'respiration', 'systemicmean']

    for num_pat, patient in enumerate(pids):    # have to go patient by patient
        pat_df = pd.read_hdf(eICU_dir + '/vitalPeriodic.h5', 
                             where='patientunitstayid = ' + str(patient), 
                             columns=['observationoffset', 'patientunitstayid'] + variables, 
                             mode='r')
        if not upto_in_minutes is None:
            pat_df = pat_df.loc[0:upto_in_minutes*60]
        # convert the offset to a TimedeltaIndex (necessary for resampling)
        pat_df.observationoffset = pd.TimedeltaIndex(pat_df.observationoffset, unit='m')
        pat_df.set_index('observationoffset', inplace=True)
        pat_df.sort_index(inplace=True)
        # resample by time
        pat_df_resampled = pat_df.resample(str(resample_factor_in_min) + 'T').median()  # pandas ignores NA in median by default
        # rename pid, cast to int
        pat_df_resampled.rename(columns={'patientunitstayid': 'pid'}, inplace=True)
        pat_df_resampled['pid'] = np.int32(pat_df_resampled['pid'])
        # get offsets in minutes from index
        pat_df_resampled['offset'] = np.int32(pat_df_resampled.index.total_seconds()/60)
        # append to file
        if num_pat == 0:
            f = open('REDACTED' + str(resample_factor_in_min) +'min.csv', 'w')
            pat_df_resampled.to_csv(f, header=True, index=False)
        else:
            pat_df_resampled.to_csv(f, header=False, index=False)

        if num_pat % 100 == 0:
            print(num_pat)
            f.flush()
    return True

def get_cohort_of_complete_downsampled_patients(time_in_hours=4, resample_factor_in_min=15):
    """
    Finds the set of patients that have no missing data during the first "time_in_hours".

    """

    resampled_pats = pd.read_csv('REDACTED' + str(resample_factor_in_min) + 'min.csv')

    time_in_minutes = time_in_hours * 60

    # delete patients with any negative offset
    print('Deleting patients with negative offsets...')
    df_posoffset = resampled_pats.groupby('pid').filter(lambda x: np.all(x.offset >= 0))

    # restrict time consideration
    print('Restricting to offsets below', time_in_minutes)
    df = df_posoffset[df_posoffset.offset <= time_in_minutes]

    #variables = ['sao2', 'heartrate', 'respiration', 'systemicmean']
    #variables = ['sao2', 'heartrate', 'respiration']
    
    # patients with no missing values in those variables (this is slow)
    print('Finding patients with no missing values in', ','.join(variables))
    good_patients = df.groupby('pid').filter(lambda x: np.all(x.loc[:, variables].isnull().sum() == 0))
    # extract the pids, save the cohort
    print('Saving...')
    cohort = (set(good_patients.pid))
    f = open('REDACTED' + str(resample_factor_in_min) + 'min.csv', 'w')
    for pid in cohort: f.write(str(int(pid)) + '\n')
    f.close()

    # save the data (only first time_in_hours hours!)
    good_patients.to_csv('REDACTED' + str(resample_factor_in_min) + 'min.csv', index=False)
    
    return True

def get_eICU_with_targets(use_age=False, use_gender=False, save=False):
    """
    Load resampled eICU data and get static prediction targets from demographics
    (patients) file
    """
    if use_age: print('Using age!')
    if use_gender: print('Using gender!')
    if save: print('Save!')
    # load resampled eICU data (the labels are the patientunitstayids)
    samples, pdf, labels = get_data('resampled_eICU', {})

    # load patients static information
    eICU_dir = 'REDACTED'
    pat_dfs = pd.read_hdf(eICU_dir + '/patient.h5', mode='r')

    # keep only static information of patients that are in the resampled table
    pat_dfs = pat_dfs[pat_dfs.patientunitstayid.isin(labels)]
    
    # reordering df to have the same order as samples and labels
    pat_dfs.set_index('patientunitstayid', inplace=True)
    pat_dfs.reindex(labels)
   
    # target variables to keep. For now we don't use hospitaldischargeoffset since it is the only integer variable.
    #target_vars = ['hospitaldischargeoffset', 'hospitaldischargestatus', 'apacheadmissiondx', 'hospitaldischargelocation', 'unittype', 'unitadmitsource']
    real_vars = ['age']
    binary_vars = ['hospitaldischargestatus', 'gender']
    categorical_vars = ['apacheadmissiondx', 'hospitaldischargelocation', 'unittype', 'unitadmitsource']

    target_vars = categorical_vars + ['hospitaldischargestatus']
    if use_age: target_vars += ['age']
    if use_gender: target_vars += ['gender']

    targets_df = pat_dfs.loc[:, target_vars]
    
    # remove patients by criteria
            # missing data in any target
    targets_df.dropna(how='any', inplace=True)
    if use_age:
                # age belonw 18 or above 89
        targets_df = targets_df[targets_df.age != '> 89']       # yes, some ages are strings
        targets_df.age = list(map(int, targets_df.age))
        targets_df = targets_df[targets_df.age >= 18]
    if use_gender:
                # remove non-binary genders (sorry!)
        targets_df['gender'] = targets_df['gender'].replace(['Female', 'Male', 'Other', 'Unknown'], [0, 1, -1, -1])
        targets_df = targets_df[targets_df.gender >= 0]
    # record patients to keep
    keep_indices = [i for (i, pid) in enumerate(labels) if pid in targets_df.index]
    assert len(keep_indices) == targets_df.shape[0]
    new_samples = samples[keep_indices]
    new_labels = np.array(labels)[keep_indices]
    
    # triple check the labels are correct
    assert np.array_equal(targets_df.index, new_labels)
    
    # getn non-one-hot targets (strings)
    targets = targets_df.values

    # one hot encoding of categorical variables
    dummies = pd.get_dummies(targets_df[categorical_vars], dummy_na=True)
    targets_df_oh = pd.DataFrame()
    targets_df_oh[dummies.columns] = dummies
    # convert binary variables to one-hot, too
    targets_df_oh['hospitaldischargestatus']= targets_df['hospitaldischargestatus'].replace(['Alive', 'Expired'],[1, 0])
    if use_gender:
        targets_df_oh['gender'] = targets_df['gender']  # already binarised
    if use_age:
        targets_df_oh['age'] = 2*targets_df['age']/89 - 1     # 89 is max

    # drop dummy columns marking missing data (they should be empty)
    nancols = [col for col in targets_df_oh.columns if col.endswith('nan')]
    assert np.all(targets_df_oh[nancols].sum() == 0)
    targets_df_oh.drop(nancols, axis=1, inplace=True)
    targets_oh = targets_df_oh.values
    
    if save:
        # save!
        # merge with training data, for LR saving
        assert new_samples.shape[0] == targets_df_oh.shape[0]
        flat_samples = new_samples.reshape(new_samples.shape[0], -1)
        features_df = pd.DataFrame(flat_samples)
        features_df.index = targets_df_oh.index
        features_df.columns = ['feature_' + str(i) for i in range(features_df.shape[1])]
        all_data = pd.concat([targets_df_oh, features_df], axis=1)
        all_data.to_csv('./data/eICU_with_targets.csv')

    # do the split
    proportions = [0.6, 0.2, 0.2]
    labels = {'targets': targets, 'targets_oh': targets_oh}
    train_seqs, vali_seqs, test_seqs, labels_split = split(new_samples, proportions, scale=True, labels=labels)
    train_targets, vali_targets, test_targets = labels_split['targets']
    train_targets_oh, vali_targets_oh, test_targets_oh = labels_split['targets_oh']

    return train_seqs, vali_seqs, test_seqs, train_targets, vali_targets, test_targets, train_targets_oh, vali_targets_oh, test_targets_oh


### --- TSTR ---- ####
def generate_synthetic(identifier, epoch, n_train, predict_labels=False):
    """
    - Load a CGAN pretrained model
    - Load its corresponding test data (+ labels)
    - Generate num_examples synthetic training data (+labels)
    - Save to format easy for training classifier on (see Eval)
    """
    settings = json.load(open('./experiments/settings/' + identifier + '.txt', 'r'))
    if not settings['cond_dim'] > 0: 
        assert settings['predict_labels']
        assert predict_labels
    # get the test data
    print('Loading test (real) data for', identifier)
    data_dict = np.load('./experiments/data/' + identifier + '.data.npy').item()
    test_data = data_dict['samples']['test']
    test_labels = data_dict['labels']['test']
    train_data = data_dict['samples']['train']
    train_labels = data_dict['labels']['train']
    print('Loaded', test_data.shape[0], 'test examples')
    print('Sampling', n_train, 'train examples from the model')
    if not predict_labels:
        assert test_data.shape[0] == test_labels.shape[0]
        if 'eICU' in settings['data']:
            synth_labels = train_labels[np.random.choice(train_labels.shape[0], n_train), :]
        else:
            # this doesn't really work for eICU...
            synth_labels = model.sample_C(n_train, settings['cond_dim'], settings['max_val'], settings['one_hot'])
            synth_data = model.sample_trained_model(settings, epoch, n_train, Z_samples=None, cond_dim=settings['cond_dim'], C_samples=synth_labels)
    else:
        assert settings['predict_labels']
        synth_data = model.sample_trained_model(settings, epoch, n_train, Z_samples=None, cond_dim=0)
        # extract the labels
        if 'eICU' in settings['data']:
            n_labels = 7
            synth_labels = synth_data[:, :, -n_labels:]
            train_labels = train_data[:, :, -n_labels:]
            test_labels = test_data[:, :, -n_labels:]
        else:
            n_labels = 6        # mnist
            synth_labels, _ = mode(np.argmax(synth_data[:, :, -n_labels:], axis=2), axis=1)
            train_labels, _ = mode(np.argmax(train_data[:, :, -n_labels:], axis=2), axis=1)
            test_labels, _ = mode(np.argmax(test_data[:, :, -n_labels:], axis=2), axis=1)
        synth_data = synth_data[:, :, :-n_labels]
        train_data = train_data[:, :, :-n_labels]
        test_data = test_data[:, :, :-n_labels]
    # package up, save
    exp_data = dict()
    exp_data['test_data'] = test_data
    exp_data['test_labels'] = test_labels
    exp_data['train_data'] = train_data
    exp_data['train_labels'] = train_labels
    exp_data['synth_data'] = synth_data
    exp_data['synth_labels'] = synth_labels
    # save it all up
    np.save('./experiments/tstr/' + identifier + '_' + str(epoch) + '.data.npy', exp_data)
    return True