<jupyter_start><jupyter_text>Mondrian Processes  Various Functions for Mondrian Processes Sampling...<jupyter_code>### SAMPLE MONDRIAN PROCESS ###
def draw_Mondrian(theta_space, budget=5):
    return draw_Mondrian_at_t(theta_space, 0, budget)
    
def draw_Mondrian_at_t(theta_space, t, budget):
    dists = theta_space[:,1] - theta_space[:,0]
    lin_dim = np.sum(dists)
    T = np.random.exponential(scale=1./lin_dim)
    
    if t+T > budget: 
        return (theta_space, None, None)
    
    d = np.argmax(np.random.multinomial(n=1, pvals=dists/lin_dim))
    x = np.random.uniform(low=theta_space[d,0], high=theta_space[d,1])
    
    theta_left = np.copy(theta_space)
    theta_left[d][1] = x 
    M_left = draw_Mondrian_at_t(theta_left, t+T, budget)
    
    theta_right = np.copy(theta_space)
    theta_right[d][0] = x 
    M_right = draw_Mondrian_at_t(theta_right, t+T, budget)
    
    return (theta_space, M_left, M_right)

def comp_log_p_sample(theta_space, data):
    if theta_space[1] == None and theta_space[2] == None:
        if data.shape[0] == 0:
            return 0
        else:
            mu = np.mean(data, axis = 0)
            residual = data - mu
            cov = np.dot(residual.T , residual) / data.shape[0] + np.identity(data.shape[1])*0.001
            return np.log(multivariate_normal.pdf(data, mean=mu, cov=cov)).sum()
    
    # find the dimension and location of first cut
    root_rec = theta_space[0]
    left_rec = theta_space[1][0]
    
    for _ in range(root_rec.shape[0]):
        if root_rec[_,1] != left_rec[_,1]:
            break
    
    dim, pos = _, left_rec[_,1]
    idx_left = data[:,dim] < pos
    idx_right = data[:,dim] >= pos
    log_len_left =  np.log(pos - root_rec[dim,0])
    log_len_right = np.log(root_rec[dim,1] - pos)
    return comp_log_p_sample(theta_space[1], data[idx_left]) + comp_log_p_sample(theta_space[2], data[idx_right])<jupyter_output><empty_output><jupyter_text>Visualization...<jupyter_code>### VISUALIZE 2D MONDRIAN PROCESS ###
def print_partitions(p, trans_level=1., color='k'):
    if not p[1] and not p[2]: 
        plt.plot([p[0][0,0], p[0][0,0]], [p[0][1,0], p[0][1,1]], color+'-', linewidth=3, alpha=trans_level)
        plt.plot([p[0][0,1], p[0][0,1]], [p[0][1,0], p[0][1,1]], color+'-', linewidth=3, alpha=trans_level)
        plt.plot([p[0][0,0], p[0][0,1]], [p[0][1,0], p[0][1,0]], color+'-', linewidth=3, alpha=trans_level)
        plt.plot([p[0][0,0], p[0][0,1]], [p[0][1,1], p[0][1,1]], color+'-', linewidth=3, alpha=trans_level)
    
    else:
        print_partitions(p[1], trans_level, color)
        print_partitions(p[2], trans_level, color)
        
        
### VISUALIZE 2D POSTERIOR WITH DATA###
def print_posterior(data, samples, trans_level=.05, color='k'):

    plt.figure()
    plt.scatter(data[:,0], data[:,1], c='k', edgecolors='k', s=5, alpha=.5)

    #print all samples
    for sample in samples:
        print_partitions(sample, trans_level, color)
        
        
def print_tree_at_leaf(mp_tree, table):

    if mp_tree[1] == None and mp_tree[2] == None: 
        print table.shape
        return 1
    
    
    # find the dimension and location of first cut
    root_rec = mp_tree[0]
    left_rec = mp_tree[1][0]
    
    for _ in range(root_rec.shape[0]):
        if root_rec[_,1] != left_rec[_,1]:
            break
    d, pos = _, left_rec[_,1]
    
    cut_type = ' '.join([str(int(x)) for x in sorted(set(table[table.columns[d]]))]) 
    
    if cut_type in {"-1 0 1", '-1 1'}: 
        idx_table_left = table[table.columns[d]] != 1
        table_left = table.loc[idx_table_left]

        idx_table_right = table[table.columns[d]] != -1
        table_right = table.loc[idx_table_right]
    
    if cut_type == '-1 0':
        idx_table_left = table[table.columns[d]] == -1
        table_left = table.loc[idx_table_left]

        idx_table_right = table[table.columns[d]] == 0
        table_right = table.loc[idx_table_right]
        

    if cut_type == '0 1':
        idx_table_left = table[table.columns[d]] == 0
        table_left = table.loc[idx_table_left]

        idx_table_right = table[table.columns[d]] == 1
        table_right = table.loc[idx_table_right]
    
    return print_tree_at_leaf(mp_tree[1], table_left) + print_tree_at_leaf(mp_tree[2], table_right)<jupyter_output><empty_output><jupyter_text>Mondrian Process Generative ModelWe apply Mondrian Processes (MPs) to flow cytometry data, using the prior information in the table above to guide the axis-aligned cuts.  Instead of uniformly, we draw the cut proportion from $w \sim \text{Beta}(a_{0}, b_{0})$.  Now let's re-implement the MP sampling function, accounting for the prior information...<jupyter_code>### SAMPLE MONDRIAN PROCESS WITH PRIOR INFORMATION ###
def draw_informed_Mondrian(theta_space, table, budget=5):
    
    # INFORMATIVE PRIORS
    upper_cut = (5., 2.)
    lower_cut = (2., 5.)
    middle_cut = (5., 5.)
    neutral_cut = (2., 2.)
    priors_dict = { '-1':lower_cut, '0':neutral_cut, '1':upper_cut, 
                   '-1 0':lower_cut, '-1 1':middle_cut, '0 1':upper_cut,
                   '-1 0 1': middle_cut, '': neutral_cut
                  }
    
    cut_history = [1] * theta_space.shape[0]
    
    return draw_informed_Mondrian_at_t(theta_space, table, priors_dict, cut_history)
    

def draw_informed_Mondrian_at_t(theta_space, table, priors_dict, cut_history):    
    
    if sum(cut_history) == 0 or table.shape[0] == 1:
        return (theta_space, None, None)

    
    types_str = [' '.join([str(int(x)) for x in sorted(set(table[table.columns[d]]))]) 
                 for d in range(table.shape[1])]
    
    if set([types_str[d] for d in range(table.shape[1]) if cut_history[d] == 1]).issubset({'0','1','-1'}):
        return (theta_space, None, None)

    
    low, medium, high, very_high = 0, 1, 100, 1000
    priority_dict = {'-1': low , '0': low, '1': low, 
                   '-1 0': medium, '0 1': medium,
                   '-1 0 1': high, '-1 1':very_high
    }    
        
    types = np.array([priority_dict[_] for _ in types_str])
    

    dists = (theta_space[:,1] - theta_space[:,0])* types    
    lin_dim = np.sum(dists)
    # draw dimension to cut
    dim_probs = ((dists/lin_dim) * np.array(cut_history)) 
    dim_probs /= np.sum(dim_probs)
    d = np.argmax(np.random.multinomial(n=1, pvals=dim_probs))
    cut_history[d] = 0

    prior_type_str = ' '.join([str(int(x)) for x in sorted(set(table[table.columns[d]]))])
    prior_params = priors_dict[prior_type_str]
    
    # make scaled cut
    x = (theta_space[d,1] - theta_space[d,0]) * np.random.beta(prior_params[0], prior_params[1]) + theta_space[d,0]
    
    cut_type = types_str[d]
    
    if cut_type in {"-1 0 1", '-1 1'}: 
        idx_table_left = table[table.columns[d]] != 1
        table_left = table.loc[idx_table_left]

        idx_table_right = table[table.columns[d]] != -1
        table_right = table.loc[idx_table_right]
    
    if cut_type == '-1 0':
        idx_table_left = table[table.columns[d]] == -1
        table_left = table.loc[idx_table_left]

        idx_table_right = table[table.columns[d]] == 0
        table_right = table.loc[idx_table_right]
        

    if cut_type == '0 1':
        idx_table_left = table[table.columns[d]] == 0
        table_left = table.loc[idx_table_left]

        idx_table_right = table[table.columns[d]] == 1
        table_right = table.loc[idx_table_right]

    
    # make lower partition
    theta_left = np.copy(theta_space)
    theta_left[d][1] = x
    M_left = draw_informed_Mondrian_at_t(theta_left, table_left, priors_dict, list(cut_history))
    
    # make upper partition
    theta_right = np.copy(theta_space)
    theta_right[d][0] = x 
    M_right = draw_informed_Mondrian_at_t(theta_right, table_right, priors_dict,list(cut_history))
    
    return (theta_space, M_left, M_right)

def Mondrian_Gaussian_perturbation(theta_space, old_sample, stepsize):
    """
    Input: 
    theta_space: a rectangle
    old_sample: partioned theta_space of a mondrian process
    stepsize: gaussian std
    """
    if old_sample[1] == None and old_sample[2] == None:
        return (theta_space, None, None)
    
    # find the dimension and location of first cut in the old_sample
    for _ in range(old_sample[0].shape[0]):
        if old_sample[0][_,1] > old_sample[1][0][_,1]:
            break    
    dim, pos = _, old_sample[1][0][_,1]
    # propose position of new cut
    good_propose = False
    while good_propose == False:
        new_pos = pos + np.random.normal(0,(old_sample[0][dim,1] - old_sample[0][dim,0])*stepsize,1)[0]
        if new_pos < theta_space[dim,1] and new_pos > theta_space[dim,0]:
            good_propose = True
    
    theta_left = np.copy(theta_space)
    theta_left[dim,1] = new_pos
    theta_right = np.copy(theta_space)
    theta_right[dim,0] = new_pos
    
    new_M_left= Mondrian_Gaussian_perturbation(theta_left, old_sample[1], stepsize)
    new_M_right = Mondrian_Gaussian_perturbation(theta_right, old_sample[2], stepsize)
    
    return (theta_space, new_M_left, new_M_right)

def comp_log_p_prior(theta_space, table, cut_history):
    """
    This function returns prior probability of a Mondrian process theta_space
    """
    if theta_space[1] == None and theta_space[2] == None:
        return 0
    
    log_prior = 0
    

    # INFORMATIVE PRIORS
    upper_cut = (5., 2.)
    lower_cut = (2., 5.)
    middle_cut = (5., 5.)
    neutral_cut = (2., 2.)
    priors_dict = { '-1':lower_cut, '0':neutral_cut, '1':upper_cut, 
                   '-1 0':lower_cut, '-1 1':middle_cut, '0 1':upper_cut,
                   '-1 0 1': middle_cut, '': neutral_cut
                  }
    
    
    # find the dimension and location of first cut
    root_rec = theta_space[0]
    left_rec = theta_space[1][0]
    
    for _ in range(root_rec.shape[0]):
        if root_rec[_,1] != left_rec[_,1]:
            break
    dim  = _
    beta_pos =  (left_rec[_,1] - left_rec[dim,0])/(root_rec[dim,1] - root_rec[dim, 0])
    
    prior_params = priors_dict[' '.join([str(int(x)) \
                                         for x in sorted(set(table[table.columns[dim]]))])]

    # compute the log likelihood of the first cut
    types_str = [' '.join([str(int(x)) for x in sorted(set(table[table.columns[d]]))]) 
                 for d in range(table.shape[1])]
    
    low_priority, medium_priority, high_priority, very_high_priority = 0, 1, 100, 1000
    priority_dict = {'-1': low_priority , '0': low_priority, '1': low_priority, 
                   '-1 0': medium_priority, '0 1': medium_priority,
                   '-1 0 1': high_priority, '-1 1':very_high_priority
    }    
        
    types = np.array([priority_dict[_] for _ in types_str])
    dists = (root_rec[:,1] - root_rec[:,0])* types    
    lin_dim = np.sum(dists)
    
    # probability of dim
    dim_probs = ((dists/lin_dim) * np.array(cut_history)) 
    dim_probs /= np.sum(dim_probs)
    log_prior  += np.log(dim_probs[dim])
    
    
    # probability of pos
    log_prior += np.log(beta.pdf(beta_pos, prior_params[0], prior_params[1]))
    
    # split the table
    cut_history[dim] = 0
    cut_type = types_str[dim]
    
    if cut_type in {"-1 0 1", '-1 1'}: 
        idx_table_left = table[table.columns[dim]] != 1
        table_left = table.loc[idx_table_left]

        idx_table_right = table[table.columns[dim]] != -1
        table_right = table.loc[idx_table_right]
    
    if cut_type == '-1 0':
        idx_table_left = table[table.columns[dim]] == -1
        table_left = table.loc[idx_table_left]

        idx_table_right = table[table.columns[dim]] == 0
        table_right = table.loc[idx_table_right]
        

    if cut_type == '0 1':
        idx_table_left = table[table.columns[dim]] == 0
        table_left = table.loc[idx_table_left]

        idx_table_right = table[table.columns[dim]] == 1
        table_right = table.loc[idx_table_right]
        
        
        
    return log_prior + comp_log_p_prior(theta_space[1], table_left, list(cut_history)) \
        + comp_log_p_prior(theta_space[2], table_right, list(cut_history))<jupyter_output><empty_output><jupyter_text>Classification<jupyter_code>def classify_cells(data, mp_tree, table, cell_type_name2idx):
    Y = np.array([1]*data.shape[0])
    
    if data.shape[0] == 0:
        return Y
    
    if mp_tree[1] == None and mp_tree[2] == None:
        if table.shape[0] > 1:
#            print "more than one clusters, number of data points:", data.shape[0]
            labels = [cell_type_name2idx[table.index[_]] for _ in range(table.shape[0])]
            return np.array(np.random.choice(labels, data.shape[0],replace = True))
        else:
            return Y * cell_type_name2idx[table.index[0]]
    
    
    # find the dimension and location of first cut
    root_rec = mp_tree[0]
    left_rec = mp_tree[1][0]
    
    for _ in range(root_rec.shape[0]):
        if root_rec[_,1] != left_rec[_,1]:
            break
    dim, pos = _, left_rec[_,1]
    
    # find labels that match dim info from table
    idx_table_left = table[table.columns[dim]] != 1
    table_left = table.loc[idx_table_left]

    idx_table_right = table[table.columns[dim]] != -1
    table_right = table.loc[idx_table_right]
    
    # find data INDICIES that go high / low on cut position in dimension dim
    idx_left = data[:,dim]  < pos
    idx_right = data[:,dim]  >= pos
    
    Y[idx_left] = classify_cells(data[idx_left],mp_tree[1],table_left, cell_type_name2idx)
    Y[idx_right] = classify_cells(data[idx_right],mp_tree[2],table_right, cell_type_name2idx)

    return Y<jupyter_output><empty_output><jupyter_text>Running time of Functions<jupyter_code>%%time
# load AML data and table
##### X: np.array, flow cytometry data, arcsin transformed
##### T: table of expert knowledge

np.random.seed(1234)
PATH = '/home/disij/projects/acdc/data/'

### LOAD DATA ###
path = PATH + 'AML_benchmark/'
df = pd.read_csv( path + 'AML_benchmark.csv.gz', sep=',', header = 0, compression = 'gzip', engine='python')
table = pd.read_csv(path + 'AML_table.csv', sep=',', header=0, index_col=0)

### PROCESS: discard ungated events ###
df = df[df.cell_type != 'NotGated']
df = df.drop(['Time', 'Cell_length','file_number', 'event_number', 'DNA1(Ir191)Di',
              'DNA2(Ir193)Di', 'Viability(Pt195)Di', 'subject'], axis = 1)
channels = [item[:item.find('(')] for item in df.columns[:-1]]
df.columns = channels + ['cell_type']
df = df.loc[df['cell_type'] != 'NotDebrisSinglets']

table = table.fillna(0)
X = df[channels].values
table_headers = list(table)
df2 = pd.DataFrame([[0]*table.shape[1]], columns=table.columns)
table.append(df2)

### transform data
data = np.arcsinh((X-1.)/5.)
theta_space = np.array([[data[:,d].min(), data[:,d].max()] for d in range(data.shape[1])])
def comp_depth_tree(sample):
    if sample == None:
        return 0
    else:
        return 1+ max(comp_depth_tree(sample[1]),comp_depth_tree(sample[2]))
%%time
# test efficiency of function draw_informed_Mondrian
sample = draw_informed_Mondrian(theta_space, table)
print comp_depth_tree(sample)
%%time
# test efficiency of function comp_log_p_sample
print comp_log_p_sample(sample, data)
%%time
# test efficiency of function comp_log_p_prior
print comp_log_p_prior(sample, table, [1 for _ in range(theta_space.shape[0])])
%%time
samples = [draw_informed_Mondrian(theta_space, table) for _ in range(100)]
log_p_priors = [comp_log_p_prior(sample, table, [1 for _ in range(theta_space.shape[0])]) for sample in samples]
plt.plot(log_p_priors)
plt.show()<jupyter_output><empty_output><jupyter_text>Flow Cytometry DataLoad AML dataset from [ACDC paper](https://www.ncbi.nlm.nih.gov/pmc/articles/PMC5447237/pdf/btx054.pdf)...<jupyter_code># load AML data and table
##### X: np.array, flow cytometry data, arcsin transformed
##### T: table of expert knowledge

np.random.seed(1234)
PATH = '/home/disij/projects/acdc/data/'

### LOAD DATA ###
path = PATH + 'AML_benchmark/'
df = pd.read_csv( path + 'AML_benchmark.csv.gz', sep=',', header = 0, compression = 'gzip', engine='python')
table = pd.read_csv(path + 'AML_table.csv', sep=',', header=0, index_col=0)
print table.shape

### PROCESS: discard ungated events ###
df = df[df.cell_type != 'NotGated']
df = df.drop(['Time', 'Cell_length','file_number', 'event_number', 'DNA1(Ir191)Di',
              'DNA2(Ir193)Di', 'Viability(Pt195)Di', 'subject'], axis = 1)
channels = [item[:item.find('(')] for item in df.columns[:-1]]
df.columns = channels + ['cell_type']
df = df.loc[df['cell_type'] != 'NotDebrisSinglets']

table = table.fillna(0)
X = df[channels].values
table_headers = list(table)
# df2 = pd.DataFrame([[0]*table.shape[1]], columns=table.columns, index =['unknown'])
# table = table.append(df2)

### transform data
data = np.arcsinh((X-1.)/5.)
theta_space = np.array([[data[:,d].min(), data[:,d].max()] for d in range(data.shape[1])])


cell_type_name2idx = {x:i for i,x in enumerate(table.index)}
Y = np.array([cell_type_name2idx[_] for _ in df.cell_type])
print data.shape<jupyter_output>(104184, 32)<jupyter_text>Experiment -1: Empirical prior.Fit a decision tree to labeled data and extract<jupyter_code>from sklearn import tree
decision_tree = tree.DecisionTreeClassifier()
decision_tree = decision_tree.fit(df[channels], Y)
print channels[21]
import graphviz 
dot_data = tree.export_graphviz(decision_tree, out_file=None) 
graph = graphviz.Source(dot_data)
graph.render("decision_tree_AML")<jupyter_output><empty_output><jupyter_text>Experiment 0: Select pretty 2-D posteriors<jupyter_code>print table.columns[2]
print table.columns[5]
print table.columns[6]
print table.columns[25]
%%time

data_2d = np.hstack([data[:,6][np.newaxis].T, data[:,25][np.newaxis].T])
np.random.shuffle(data_2d)
data_2d = data_2d[:1000,:]
print data_2d.shape

x_min, y_min, x_max, y_max = data_2d[:,0].min(), data_2d[:,1].min(), data_2d[:,0].max(), data_2d[:,1].max()
theta_space_2d = np.array([[x_min, x_max], [y_min, y_max]])
print theta_space_2d

table_2d = table[['CD8','CD3']]
print table_2d

n_mcmc_chain = 20
n_mcmc_sample = 1000
mcmc_gaussin_std = 0.1 # tune step size s.t. acceptance rate ~50%

accepts = [[] for _ in range(n_mcmc_chain)]
rejects = [[] for _ in range(n_mcmc_chain)]


for chain in range(n_mcmc_chain):
    
    if chain % 5 == 0:
        print "Drawing Chain %d ..." % chain
    
    sample = draw_informed_Mondrian(theta_space_2d, table_2d)
    log_p_sample = comp_log_p_sample(sample, data_2d) + \
                     comp_log_p_prior(sample, table_2d, [1 for _ in range(table_2d.shape[1])])
    
    accepts[chain].append(sample)
    
    for idx in xrange(n_mcmc_sample):
        
        new_sample = Mondrian_Gaussian_perturbation(theta_space_2d,sample, mcmc_gaussin_std)

        # perform accept-reject step
        new_log_p_sample = comp_log_p_sample(new_sample, data_2d) + \
                            comp_log_p_prior(new_sample, table_2d, [1 for _ in range(table_2d.shape[1])])
        
        if new_log_p_sample <  log_p_sample and \
            np.log(np.random.uniform(low=0, high=1.)) > new_log_p_sample - log_p_sample:
                rejects[chain].append(new_sample)
        
        else:
            sample = new_sample
            log_p_sample = new_log_p_sample
            accepts[chain].append(sample)
        

        if (idx+1) % 10000 == 0 and (chain+1) % 1 == 0:
            print "Iteration %d, cummulative accepted sample size is %d" %(idx+1, len(accepts[chain]))

print "Total number of accepted samples: %d" %(sum([len(accepts[chain]) for chain in range(n_mcmc_chain)]))
### VISUALIZE 2D MONDRIAN PROCESS ###
def print_partitions(p, trans_level=1., color='k'):
    if not p[1] and not p[2]: 
        plt.plot([p[0][0,0], p[0][0,0]], [p[0][1,0], p[0][1,1]], color+'-', linewidth=3, alpha=trans_level)
        plt.plot([p[0][0,1], p[0][0,1]], [p[0][1,0], p[0][1,1]], color+'-', linewidth=3, alpha=trans_level)
        plt.plot([p[0][0,0], p[0][0,1]], [p[0][1,0], p[0][1,0]], color+'-', linewidth=3, alpha=trans_level)
        plt.plot([p[0][0,0], p[0][0,1]], [p[0][1,1], p[0][1,1]], color+'-', linewidth=3, alpha=trans_level)
    
    else:
        print_partitions(p[1], trans_level, color)
        print_partitions(p[2], trans_level, color)
        
        
### VISUALIZE 2D POSTERIOR WITH DATA###
def print_posterior(data, samples, trans_level=.05, color='k'):

    plt.figure()
    plt.scatter(data[:,0], data[:,1], c='k', edgecolors='k', s=10, alpha=.2)

    #print all samples
    for sample in samples:
        print_partitions(sample, trans_level, color)
        
# print prior
burnt_accepts = [_ for chain in accepts for _ in chain[0:1]]
print_posterior(data_2d, burnt_accepts, trans_level=.2, color='b')

plt.xlabel("CD8")
plt.ylabel("CD3")
plt.title("Mondrian Process Prior Samples")

plt.xlim([x_min, x_max])
plt.ylim([y_min, y_max])
plt.show()


# print posterior
burnt_accepts = [_ for chain in accepts for _ in chain[-2:]]
print_posterior(data_2d, burnt_accepts, trans_level=.1, color='b')

plt.xlabel("CD8")
plt.ylabel("CD3")
#plt.title("Mondrian Process Posterior Samples")

plt.xlim([x_min, x_max])
plt.ylim([y_min, y_max])
plt.show()
%%time

data_2d = np.hstack([data[:,5][np.newaxis].T, data[:,25][np.newaxis].T])
np.random.shuffle(data_2d)
data_2d = data_2d[:1000,:]
print data_2d.shape

x_min, y_min, x_max, y_max = data_2d[:,0].min(), data_2d[:,1].min(), data_2d[:,0].max(), data_2d[:,1].max()
theta_space_2d = np.array([[x_min, x_max], [y_min, y_max]])
print theta_space_2d

table_2d = table[['CD4','CD3']]
print table_2d

n_mcmc_chain = 100
n_mcmc_sample = 100
mcmc_gaussin_std = 0.1 # tune step size s.t. acceptance rate ~50%

accepts = [[] for _ in range(n_mcmc_chain)]
rejects = [[] for _ in range(n_mcmc_chain)]


for chain in range(n_mcmc_chain):
    
    if chain % 5 == 0:
        print "Drawing Chain %d ..." % chain
    
    sample = draw_informed_Mondrian(theta_space_2d, table_2d)
    log_p_sample = comp_log_p_sample(sample, data_2d) + \
                     comp_log_p_prior(sample, table_2d, [1 for _ in range(table_2d.shape[1])])
    
    accepts[chain].append(sample)
    
    for idx in xrange(n_mcmc_sample):
        
        new_sample = Mondrian_Gaussian_perturbation(theta_space_2d,sample, mcmc_gaussin_std)

        # perform accept-reject step
        new_log_p_sample = comp_log_p_sample(new_sample, data_2d) + \
                            comp_log_p_prior(new_sample, table_2d, [1 for _ in range(table_2d.shape[1])])
        
        if new_log_p_sample <  log_p_sample and \
            np.log(np.random.uniform(low=0, high=1.)) > new_log_p_sample - log_p_sample:
                rejects[chain].append(new_sample)
        
        else:
            sample = new_sample
            log_p_sample = new_log_p_sample
            accepts[chain].append(sample)
        

        if (idx+1) % 10000 == 0 and (chain+1) % 1 == 0:
            print "Iteration %d, cummulative accepted sample size is %d" %(idx+1, len(accepts[chain]))

print "Total number of accepted samples: %d" %(sum([len(accepts[chain]) for chain in range(n_mcmc_chain)]))

### VISUALIZE 2D MONDRIAN PROCESS ###
def print_partitions(p, trans_level=1., color='k'):
    if not p[1] and not p[2]: 
        plt.plot([p[0][0,0], p[0][0,0]], [p[0][1,0], p[0][1,1]], color+'-', linewidth=3, alpha=trans_level)
        plt.plot([p[0][0,1], p[0][0,1]], [p[0][1,0], p[0][1,1]], color+'-', linewidth=3, alpha=trans_level)
        plt.plot([p[0][0,0], p[0][0,1]], [p[0][1,0], p[0][1,0]], color+'-', linewidth=3, alpha=trans_level)
        plt.plot([p[0][0,0], p[0][0,1]], [p[0][1,1], p[0][1,1]], color+'-', linewidth=3, alpha=trans_level)
    
    else:
        print_partitions(p[1], trans_level, color)
        print_partitions(p[2], trans_level, color)
        
        
### VISUALIZE 2D POSTERIOR WITH DATA###
def print_posterior(data, samples, trans_level=.05, color='k'):

    plt.figure()
    plt.scatter(data[:,0], data[:,1], c='k', edgecolors='k', s=10, alpha=.3)

    #print all samples
    for sample in samples:
        print_partitions(sample, trans_level, color)

# print prior
burnt_accepts = [_ for chain in accepts for _ in chain[0:1]]
print_posterior(data_2d, burnt_accepts, trans_level=.1, color='b')

plt.xlabel("CD4")
plt.ylabel("CD3")
plt.title("Mondrian Process Prior Samples")

plt.xlim([x_min, x_max])
plt.ylim([y_min, y_max])
plt.show()


# print posterior
burnt_accepts = [_ for chain in accepts for _ in chain[-1:]]
print_posterior(data_2d, burnt_accepts, trans_level=.05, color='b')

# plt.xlabel("CD4")
# plt.ylabel("CD3")
#plt.title("Mondrian Process Posterior Samples")

plt.xlim([x_min, x_max])
plt.ylim([y_min, y_max])
plt.show()

# COMPUTE THE NUMBER OF TREES WHOSE FIRST CUT IS ON CD3

first_cut_counter = dict()
first_cut_counter[0] = 0
first_cut_counter[1] = 0

print len(burnt_accepts)
for tree in burnt_accepts:
    root_rec = tree[0]
    left_rec = tree[1][0]
    for _ in range(2):
        if root_rec[_,1] != left_rec[_,1]:
            break
    d, pos = _, left_rec[_,1]
    first_cut_counter[d] += 1
print first_cut_counter


dim1 = theta_space_2d[0][1] - theta_space_2d[0][0]
dim2 = theta_space_2d[1][1] - theta_space_2d[1][0]
print "dimension of CD4:", dim1
print "dimension of CD3:", dim2
print dim1/dim2<jupyter_output>100
{0: 2, 1: 98}
dimension of CD4: 2.538402261456195
dimension of CD3: 6.281194696224246
0.4041273012890494<jupyter_text>Experiment 1: 2D SimulationLet's now test out the idea on a 2D subset of the data: the *CD4* and *CD8* features.<jupyter_code>data_2d = np.hstack([data[:,5][np.newaxis].T, data[:,6][np.newaxis].T])
np.random.shuffle(data_2d)
data_2d = data_2d[:1000,:]
print data_2d.shape

x_min, y_min, x_max, y_max = data_2d[:,0].min(), data_2d[:,1].min(), data_2d[:,0].max(), data_2d[:,1].max()
theta_space_2d = np.array([[x_min, x_max], [y_min, y_max]])
print theta_space_2d

table_2d = table[['CD4','CD8']]
# test perturbation function
sample = draw_informed_Mondrian(theta_space_2d, table_2d)
plt.scatter(data_2d[:,0], data_2d[:,1], c='k', edgecolors='k', s=10, alpha=.5)
print_partitions(sample)
plt.show()


stepsize = 0.1
propose_sample = Mondrian_Gaussian_perturbation(theta_space_2d, sample, stepsize)

plt.scatter(data_2d[:,0], data_2d[:,1], c='k', edgecolors='k', s=10, alpha=.5)
print_partitions(propose_sample)
plt.show()<jupyter_output><empty_output><jupyter_text>Now run MCMC to collect posterior samples...<jupyter_code>%%time

n_mcmc_chain = 5
n_mcmc_sample = 10000
mcmc_gaussin_std = 0.005 # tune step size s.t. acceptance rate ~50%

accepts = [[] for _ in range(n_mcmc_chain)]
rejects = [[] for _ in range(n_mcmc_chain)]


for chain in range(n_mcmc_chain):
    
    if chain % 5 == 0:
        print "Drawing Chain %d ..." % chain
    
    sample = draw_informed_Mondrian(theta_space_2d, table_2d)
    log_p_sample = comp_log_p_sample(sample, data_2d) + \
                     comp_log_p_prior(sample, table_2d, [1 for _ in range(table_2d.shape[1])])
    
    accepts[chain].append(sample)
    
    for idx in xrange(n_mcmc_sample):
        
        new_sample = Mondrian_Gaussian_perturbation(theta_space_2d,sample, mcmc_gaussin_std)

        # perform accept-reject step
        new_log_p_sample = comp_log_p_sample(new_sample, data_2d) + \
                            comp_log_p_prior(new_sample, table_2d, [1 for _ in range(table_2d.shape[1])])
        
        if new_log_p_sample <  log_p_sample and \
            np.log(np.random.uniform(low=0, high=1.)) > new_log_p_sample - log_p_sample:
                rejects[chain].append(new_sample)
        
        else:
            sample = new_sample
            log_p_sample = new_log_p_sample
            accepts[chain].append(sample)
        

        if (idx+1) % 10000 == 0 and (chain+1) % 1 == 0:
            print "Iteration %d, cummulative accepted sample size is %d" %(idx+1, len(accepts[chain]))

print "Total number of accepted samples: %d" %(sum([len(accepts[chain]) for chain in range(n_mcmc_chain)]))
# get an average model
burnt_accepts = np.array([_ for chain in accepts for _ in chain[len(chain)*10/11:]])
#### plt.figure()

np.random.shuffle(burnt_accepts)
print_posterior(data_2d, burnt_accepts[:100], trans_level=.005, color='b')

plt.xlabel("CD4")
plt.ylabel("CD8")
plt.title("Mondrian Process Posterior Samples")

plt.xlim([x_min, x_max])
plt.ylim([y_min, y_max])
plt.show()
#### plt.figure()


print_posterior(data_2d, accepts[0][:50], trans_level=.05, color='b')

plt.xlabel("CD4")
plt.ylabel("CD8")
plt.title("Mondrian Process Posterior Samples, Chain 0, first 100 acceped samples")

plt.xlim([x_min, x_max])
plt.ylim([y_min, y_max])
plt.show()



print_posterior(data_2d, accepts[0][-50:], trans_level=.05, color='b')

plt.xlabel("CD4")
plt.ylabel("CD8")
plt.title("Mondrian Process Posterior Samples, Chain 0, last 100 acceped samples")

plt.xlim([x_min, x_max])
plt.ylim([y_min, y_max])
plt.show()
#### plt.figure()

print_posterior(data_2d, accepts[1][:50], trans_level=.05, color='b')

plt.xlabel("CD4")
plt.ylabel("CD8")
plt.title("Mondrian Process Posterior Samples, Chain 2, first 150 acceped samples")

plt.xlim([x_min, x_max])
plt.ylim([y_min, y_max])
plt.show()


print_posterior(data_2d, accepts[1][-50:], trans_level=.05, color='b')

plt.xlabel("CD4")
plt.ylabel("CD8")
plt.title("Mondrian Process Posterior Samples, Chain 2, last 100 acceped samples")

plt.xlim([x_min, x_max])
plt.ylim([y_min, y_max])
plt.show()<jupyter_output><empty_output><jupyter_text>Experiment 2: subset AML data, 3 chains<jupyter_code>from sklearn.utils import shuffle

N = 10000
new_df = shuffle(df)[:N]
X_subset = new_df[channels].values
Y_subset = np.array([cell_type_name2idx[_] for _ in new_df.cell_type])
# X_subset = df[channels].values
# Y_subset = np.array([cell_type_name2idx[_] for _ in df.cell_type])
data_subset = np.arcsinh((X_subset-1.)/5.)
N, d = data_subset.shape
print N,d

# rename table header 'HLA-DR' to 'HLADR' to prevent error from '-'
temp_headers = list(table)
temp_headers[29] = "HLADR"
table.columns = temp_headers
print table.columns

emp_bounds = np.array([(data_subset[:,i].min(), data_subset[:,i].max()) for i in range(d)])
%%time

n_mcmc_chain = 10
n_mcmc_sample = 5000

accepts = [[] for _ in range(n_mcmc_chain)]
rejects = [[] for _ in range(n_mcmc_chain)]
logl_accepted_trace = [[] for _ in range(n_mcmc_chain)]
logl_complete_trace = [[] for _ in range(n_mcmc_chain)]
Y_predict_accepted_trace = [[] for _ in range(n_mcmc_chain)]
accuracy_accepted_trace = [[] for _ in range(n_mcmc_chain)]


for chain in range(n_mcmc_chain):
    mcmc_gaussin_std = 0.1 # tune step size s.t. acceptance rate ~50%
    print "Drawing Chain %d ..." % chain
    
    sample = draw_informed_Mondrian(emp_bounds, table)
    log_p_sample = comp_log_p_sample(sample, data_subset)
    
    accepts[chain].append(sample)
    logl_accepted_trace[chain].append(log_p_sample)
    logl_complete_trace[chain].append(log_p_sample)

    
    Y_predict = classify_cells(data_subset, sample, table, cell_type_name2idx)
    accuracy = sum(Y_subset == Y_predict)*1.0/ data_subset.shape[0]
    accuracy_accepted_trace[chain].append(accuracy)
    Y_predict_accepted_trace[chain].append(Y_predict)
    
    for idx in range(n_mcmc_sample):
        
        if idx % (n_mcmc_sample / 10) == 0:
            mcmc_gaussin_std /= 2
        
        new_sample = Mondrian_Gaussian_perturbation(emp_bounds,sample, mcmc_gaussin_std)
        new_log_p_sample = comp_log_p_sample(new_sample, data_subset) 
        logl_complete_trace[chain].append(new_log_p_sample)
        
        if new_log_p_sample <  log_p_sample and \
            np.log(np.random.uniform(low=0, high=1.)) > (new_log_p_sample - log_p_sample):
                rejects[chain].append(new_sample)
        
        else:
            if new_log_p_sample <  log_p_sample:
                print "accepted some bad samples"
                
            sample = new_sample
            log_p_sample = new_log_p_sample
            accepts[chain].append(sample)
            logl_accepted_trace[chain].append(log_p_sample)
        
            Y_predict = classify_cells(data_subset, sample, table, cell_type_name2idx)
            accuracy = sum(Y_subset == Y_predict)*1.0/ data_subset.shape[0]
            accuracy_accepted_trace[chain].append(accuracy)
            Y_predict_accepted_trace[chain].append(Y_predict)
        
        if (idx+1) % 500 == 0:
            print "Iteration %d, cummulative accepted sample size is %d" %(idx+1, len(accepts[chain]))
            
    
    if (chain + 1) % 1 == 0:
        # prediction and visualization
        Y_predict = classify_cells(data_subset, accepts[chain][-1], table, cell_type_name2idx)
        accuracy = sum(Y_subset == Y_predict)*1.0/ data_subset.shape[0]
        print "Chain % d accuracy on subset data: %.3f" % (chain+1,accuracy)

            
print "Total number of accepted samples: %d" %(sum([len(accepts[chain]) for chain in range(n_mcmc_chain)]))
# 3 chain, 5000 samples<jupyter_output><empty_output><jupyter_text>Visualize hierarchical histograms<jupyter_code># generte tree structures for visualization
def find_cuts(list_of_mp_trees):
    if len(list_of_mp_trees) == 0:
        return None
    mp_tree = list_of_mp_trees[0]
    if mp_tree[1] == None and mp_tree[2] == None: 
        return None
    
    # find the dimension and location of first cut
    root_rec = mp_tree[0]
    left_rec = mp_tree[1][0]
    
    for _ in range(root_rec.shape[0]):
        if root_rec[_,1] != left_rec[_,1]:
            break
    d, pos = _, left_rec[_,1]
    
    first_cut = [d]+[mp_tree[1][0][d,1] for mp_tree in list_of_mp_trees]
    
    list_of_left_mp_trees = [mp_tree[1] for mp_tree in list_of_mp_trees]
    list_of_right_mp_trees = [mp_tree[2] for mp_tree in list_of_mp_trees]
    
    return [first_cut, find_cuts(list_of_left_mp_trees),find_cuts(list_of_right_mp_trees)]

def split_data_by_MP(mp_tree, data):
    if mp_tree[1] == None and mp_tree[2] == None:
        return [data,None, None]
    
    # find the dimension and location of first cut
    root_rec = mp_tree[0]
    left_rec = mp_tree[1][0]
    
    for _ in range(root_rec.shape[0]):
        if root_rec[_,1] != left_rec[_,1]:
            break
    d, pos = _, left_rec[_,1]    
    
    data_left = data[data[:,d] < pos]
    data_right = data[data[:,d] >= pos]
    
    return [data,split_data_by_MP(mp_tree[1], data_left),split_data_by_MP(mp_tree[2], data_right)]
    
def print_cuts_on_hist(tree_tructured_data, theta_space_cuts, node_pos, table):
    from scipy.stats import gaussian_kde
    
    if tree_tructured_data == None:
        return
    depth = len(node_pos)
    cuts = theta_space_cuts
    data_leaf = tree_tructured_data
    for _ in node_pos:
        _ = int(_)
        data_leaf = data_leaf[_]

    for _ in node_pos:
        _ = int(_)
        cuts = cuts[_]
        if cuts == None:
            print "leaf node"
            print "size of cell type:,", data_leaf.shape
            return
        
    dim = cuts[0]
    list_pos = cuts[1:]
    
    print data_leaf.shape
    
    # plot histogram    
    n, b, patches = plt.hist(data_leaf[:,dim], bins=30, \
                             alpha = 0.75, label=table.columns[dim],color = '#0343DF')
    plt.legend(loc='upper right',prop={'size': 20})
    
    plt.yscale('log')
    plt.xticks(fontsize=20)
    plt.yticks(fontsize=20)
    
    # plot cuts
    for pos in list_pos:
        plt.plot([pos,pos],[0,n.max()],color='red',linewidth = 5,alpha= 0.1)
    plt.show()
    return
# X_fullset = df[channels].values
# Y_fullset = np.array([cell_type_name2idx[_] for _ in df.cell_type])
# data_fullset = np.arcsinh((X_fullset-1.)/5.)


max_logl_chains =  [sum(_[-1:]) for _ in logl_accepted_trace]
optimal_chain = accepts[max_logl_chains.index(max(max_logl_chains))][-3:]

#optimal_chain = accepts[2][-3:]

optimal_chain = [accepts[7][-1]]
theta_space_cuts = find_cuts(optimal_chain)

#'0','10','20','110','120',''1110','1120','11110',
#'111110','111120', '1111210'
# '11112110',,'11112120'
#'0','10','20','110','120','1110','1120','1210','11110','11120','12120','111120','111210','111220'
node_pos_list = ['0','10','20','110','120','210','1110','1120','11110','111110','111120', '1111210','11112110','11112120']
tree_tructured_data = split_data_by_MP(optimal_chain[-1], data_fullset)
for node_pos in node_pos_list:
    print node_pos
    print_cuts_on_hist(tree_tructured_data,theta_space_cuts,node_pos,table)
n_mcmc_sample = 5000
mcmc_gaussin_std = 0.005

  
sample = accepts_append[-1]
log_p_sample = comp_log_p_sample(sample, data_fullset)

for idx in range(n_mcmc_sample):

    if idx % (n_mcmc_sample / 10) == 0:
        mcmc_gaussin_std /= 2

    new_sample = Mondrian_Gaussian_perturbation(emp_bounds,sample, mcmc_gaussin_std)
    new_log_p_sample = comp_log_p_sample(new_sample, data_fullset) 
    logl_complete_trace[chain].append(new_log_p_sample)

    if new_log_p_sample <  log_p_sample and \
        np.log(np.random.uniform(low=0, high=1.)) > (new_log_p_sample - log_p_sample):
            rejects_append.append(new_sample)

    else:
        if new_log_p_sample <  log_p_sample:
            print "accepted some bad samples"

        sample = new_sample
        log_p_sample = new_log_p_sample
        accepts_append.append(sample)

    if (idx+1) % 500 == 0:
        print "Iteration %d, cummulative accepted sample size is %d" %(idx+1, len(accepts_append))
optimal_chain = accepts_append[-100:]
#optimal_chain = accepts[7][-100:]
theta_space_cuts = find_cuts(optimal_chain)

#'0','10','20','110','120',''1110','1120','11110',
#'111110','111120', '1111210'
# '11112110',,'11112120'
#'0','10','20','110','120','1110','1120','1210','11110','11120','12120','111120','111210','111220'
node_pos_list = ['0','10','20','110','120','210','1110','1120','11110','111110','111120', '1111210','11112110','11112120']
tree_tructured_data = split_data_by_MP(optimal_chain[-1], data_fullset)
for node_pos in node_pos_list:
    print node_pos
    print_cuts_on_hist(tree_tructured_data,theta_space_cuts,node_pos,table)

Y_predict = classify_cells(data_subset, optimal_chain[-1], table, cell_type_name2idx)
accuracy = sum(Y_subset == Y_predict)*1.0/ data_subset.shape[0]
print "Chain % d accuracy on subset data: %.3f" % (chain+1,accuracy)
import matplotlib.style
import matplotlib as mpl
plt.rcParams.update(plt.rcParamsDefault)
optimal_chain = accepts_append[-100:]
#optimal_chain = accepts[7][-100:]
theta_space_cuts = find_cuts(optimal_chain)

#'0','10','20','110','120',''1110','1120','11110',
#'111110','111120', '1111210'
# '11112110',,'11112120'
#'0','10','20','110','120','1110','1120','1210','11110','11120','12120','111120','111210','111220'
node_pos_list = ['0','10','20','110','120','210','1110','1120','11110','111110','111120', '1111210','11112110','11112120']
tree_tructured_data = split_data_by_MP(optimal_chain[-1], data_fullset)
for node_pos in node_pos_list:
    print node_pos
    print_cuts_on_hist(tree_tructured_data,theta_space_cuts,node_pos,table)

Y_predict = classify_cells(data_subset, optimal_chain[-1], table, cell_type_name2idx)
accuracy = sum(Y_subset == Y_predict)*1.0/ data_subset.shape[0]
print "Chain % d accuracy on subset data: %.3f" % (chain+1,accuracy)
PATH_SAMPLES = "/home/disij/projects/sandbox/fc_mondrian/output/AML_optimal_chain_for_visualize"
def write_chains_to_file(accepts, PATH):
    """
    accepts: a list of lists of MP trees
    PATH: directory to write the samples to
    
    each sample is written to a pickle file
    """
    
    ### remove previous files in the path if the directory exists, otherwise create a new empty directory
    if os.path.exists(PATH):
        files = [ f for f in os.listdir(PATH) if f.endswith(".pkl") ]
        for f in files:
            os.remove(PATH + '/' + f)
    else:
        os.makedirs(PATH)
        
    n_mcmc_chain = len(accepts)
    
    for i in range(n_mcmc_chain):
        for j in range(len(accepts[i])):
            filename = '%s/chain_%d_id_%d.pkl' %(PATH,i,j)
            with open(filename, "wb") as fp:
                pickle.dump(accepts[i][j], fp)   
write_chains_to_file([optimal_chain], PATH_SAMPLES)
%%time

n_mcmc_sample = 5000


mcmc_gaussin_std = 0.1 # tune step size s.t. acceptance rate ~50%
print "Drawing Chain %d ..." % chain

sample = draw_informed_Mondrian(emp_bounds, table)
log_p_sample = comp_log_p_sample(sample, data_subset)

accepts[chain].append(sample)
logl_accepted_trace[chain].append(log_p_sample)
logl_complete_trace[chain].append(log_p_sample)


Y_predict = classify_cells(data_subset, sample, table, cell_type_name2idx)
accuracy = sum(Y_subset == Y_predict)*1.0/ data_subset.shape[0]
accuracy_accepted_trace[chain].append(accuracy)
Y_predict_accepted_trace[chain].append(Y_predict)

for idx in range(n_mcmc_sample):

    if idx % (n_mcmc_sample / 10) == 0:
        mcmc_gaussin_std /= 2

    new_sample = Mondrian_Gaussian_perturbation(emp_bounds,sample, mcmc_gaussin_std)
    new_log_p_sample = comp_log_p_sample(new_sample, data_subset) 
    logl_complete_trace[chain].append(new_log_p_sample)

    if new_log_p_sample <  log_p_sample and \
        np.log(np.random.uniform(low=0, high=1.)) > (new_log_p_sample - log_p_sample):
            rejects[chain].append(new_sample)

    else:
        if new_log_p_sample <  log_p_sample:
            print "accepted some bad samples"

        sample = new_sample
        log_p_sample = new_log_p_sample
        accepts[chain].append(sample)
        logl_accepted_trace[chain].append(log_p_sample)

        Y_predict = classify_cells(data_subset, sample, table, cell_type_name2idx)
        accuracy = sum(Y_subset == Y_predict)*1.0/ data_subset.shape[0]
        accuracy_accepted_trace[chain].append(accuracy)
        Y_predict_accepted_trace[chain].append(Y_predict)

    if (idx+1) % 500 == 0:
        print "Iteration %d, cummulative accepted sample size is %d" %(idx+1, len(accepts[chain]))


# prediction and visualization
Y_predict = classify_cells(data_subset, accepts[chain][-1], table, cell_type_name2idx)
accuracy = sum(Y_subset == Y_predict)*1.0/ data_subset.shape[0]
print "Chain % d accuracy on subset data: %.3f" % (chain+1,accuracy)


print "Total number of accepted samples: %d" %(sum([len(accepts[chain]) for chain in range(n_mcmc_chain)]))
# 3 chain, 5000 samples
node_pos_list = ['210','220','1220','11210','11220','12110','111110',\
                 '121210','121220','1111210','1111220','1112110','1112120',\
                 '1112210','1112220']
tree_tructured_data = split_data_by_MP(optimal_chain[-1], data_subset)
for node_pos in node_pos_list:
    print_cuts_on_hist(tree_tructured_data,theta_space_cuts,node_pos,table)
def query_table(table,list_of_condition):
    if len(list_of_condition) == 0:
        return table
    
    dim, label = list_of_condition[0]
    if label == 1:
        table = table[table[dim] >= 0]
    else:
        table = table[table[dim] <= 0]
        
    return query_table(table, list_of_condition[1:])

conditions  = [['CD3',1],['CD4', -1]]
print query_table(table, conditions).index
conditions  = [['CD3',1],['CD4', 1]]
print query_table(table, conditions).index
conditions  = [['CD3',-1],['CD34', 1],['CD123',1]]
print query_table(table, conditions).index
conditions  = [['CD3',-1],['CD34', -1],['CD123',1],['HLADR',-1]]
print query_table(table, conditions).index
conditions  = [['CD3',-1],['CD34', -1],['CD123',1],['HLADR',1]]
print query_table(table, conditions).index
conditions  = [['CD3',-1],['CD34', 1],['CD123',-1],['CD38',-1]]
print query_table(table, conditions).index
conditions = [['CD3',-1],['CD34',-1],['CD123',-1],['CD19',-1],['CD7',-1]]
print query_table(table, conditions).index
conditions  = [['CD3',-1],['CD34', 1],['CD123',-1],['CD38',1],['CD19',-1]]
print query_table(table, conditions).index
conditions  = [['CD3',-1],['CD34', 1],['CD123',-1],['CD38',1],['CD19',1]]
print query_table(table, conditions).index
conditions = [['CD3',-1],['CD34',-1],['CD123',-1],['CD19',-1],['CD7',1],['CD16',-1]]
print query_table(table, conditions).index
conditions = [['CD3',-1],['CD34',-1],['CD123',-1],['CD19',-1],['CD7',1],['CD16',-1]]
print query_table(table, conditions).index
conditions  = [['CD3',-1],['CD34',-1],['CD123',-1],['CD19',1],['HLADR',-1],['CD38',-1]]
print query_table(table, conditions).index
conditions  = [['CD3',-1],['CD34',-1],['CD123',-1],['CD19',1],['HLADR',-1],['CD38',1]]
print query_table(table, conditions).index
conditions  = [['CD3',-1],['CD34',-1],['CD123',-1],['CD19',1],['HLADR',1],['CD38',-1]]
print query_table(table, conditions).index
conditions  = [['CD3',-1],['CD34',-1],['CD123',-1],['CD19',1],['HLADR',1],['CD38',1]]
print query_table(table, conditions).index


print table.loc['CD34+CD38+CD123+ HSPCs']
print table
cell_type_name2idx['Pro B cells']
len(Y[Y == 11]) * 1.0 / len(Y)
Y_predict = classify_cells(data, accepts_append[0][-3], table, cell_type_name2idx)
accuracy = sum(Y == Y_predict)*1.0/ data.shape[0]
print "Accuracy on data: %.3f" % accuracy

bins = table.shape[0]
plt.hist(Y, bins, alpha=0.5, label='Y:cell type')
plt.hist(Y_predict, bins, alpha=0.5, label='Z:prediction')
plt.legend(loc='upper right')
plt.show() 
# plot 5 chains

threshold = -113888.425449/data.shape[0]*N
fig, axs = plt.subplots(5, 3, figsize=(10,10) )
for chain in range(5):
    axs[chain, 0].plot(logl_complete_trace[chain])
    axs[chain, 1].plot(logl_accepted_trace[chain])
    axs[chain, 2].plot(accuracy_accepted_trace[chain])
    axs[chain, 0].set_title('Trace of likelihood Chain %d, all samples' % chain, fontsize=8)
    axs[chain, 1].set_title('Trace of likelihood Chain %d, accepted samples' % chain, fontsize=8)
    axs[chain, 2].set_title('Trace of accuracy Chain %d, accepted samples' % chain, fontsize=8)
    
fig.tight_layout()
print len(accepts)
n_mcmc_chain = 1
n_mcmc_sample = 10000
mcmc_gaussin_std = 0.01 # tune step size s.t. acceptance rate ~50%

accepts_temp = [[] for _ in range(n_mcmc_chain)]
rejects_temp = [[] for _ in range(n_mcmc_chain)]

for chain in range(n_mcmc_chain):
    
    sample = accepts[28][-1]
    log_p_sample = comp_log_p_sample(sample, data_subset)
    accepts_temp[chain].append(sample)
    
    for idx in range(n_mcmc_sample):
        
        new_sample = Mondrian_Gaussian_perturbation(emp_bounds,sample, mcmc_gaussin_std)
        new_log_p_sample = comp_log_p_sample(new_sample, data_subset) 
        logl_complete_trace[chain].append(new_log_p_sample)
        
        if new_log_p_sample <  log_p_sample and \
            np.log(np.random.uniform(low=0, high=1.)) > (new_log_p_sample - log_p_sample):
                rejects_temp[chain].append(new_sample)
        
        else:
            if new_log_p_sample <  log_p_sample:
                print "accepted some bad samples"
                
            sample = new_sample
            log_p_sample = new_log_p_sample
            accepts_temp[chain].append(sample)
        
        if (idx+1) % 500 == 0:
            print "Iteration %d, cummulative accepted sample size is %d" %(idx+1, len(accepts_append[chain]))
                        
print "Total number of accepted samples: %d" %(sum([len(accepts[chain]) for chain in range(n_mcmc_chain)]))
# 3 chain, 5000 samples

optimal_chain = accepts_temp[0][-1:]
theta_space_cuts = find_cuts(optimal_chain)
#'0','10','20','110','120','1110','1120','1210','11110','11120','12120','111120','111210','111220'
node_pos_list = ['0','10','20']
tree_tructured_data = split_data_by_MP(optimal_chain[-1], data_subset)
for node_pos in node_pos_list:
    print_cuts_on_hist(tree_tructured_data,theta_space_cuts,node_pos,table)
subtable = table[table['CD3']<=0]
subtable = subtable[subtable['CD34']>=0]
subtable = subtable[subtable['CD123']<=0]
subtable = subtable[subtable['CD38']<=0]
subtable
print max_logl_chains.index(max(max_logl_chains))
# vote, and compute accuracy

threshold = -113888.425449/data.shape[0]*N  - 10000000000
# keep that samples whose likelihood is above threshold
burnt_samples = []
burnt_predictions = []
for i in range(len(accepts)):
    accepted_chain = accepts[i]
    likelihoods = logl_accepted_trace[i]
    predictions = Y_predict_accepted_trace[i]
    burnt_samples += [accepted_chain[_] for _ in range(len(accepted_chain)) if likelihoods[_] > threshold]
    burnt_predictions += [predictions[_] for _ in range(len(accepted_chain)) if likelihoods[_] > threshold]
    
# vote
votes = np.zeros([data_subset.shape[0], table.shape[0]])
for Y_predict in burnt_predictions:
    for _ in range(len(Y_predict)):
        votes[_,Y_predict[_]] += 1

Y_predict_majority = np.argmax(votes, axis=1)

print votes.sum()

accuracy = sum(Y_subset == Y_predict_majority)*1.0/ data_subset.shape[0]
print "Chain % d accuracy on subset data: %.3f" % (chain+1,accuracy)

bins = table.shape[0]
plt.hist(Y_subset, bins, alpha=0.5, label='Y:cell type')
plt.hist(Y_predict_majority, bins, alpha=0.5, label='Z:prediction')
plt.legend(loc='upper right')
plt.show() 
# emsemble 100 samples from prior

# vote
votes = np.zeros([data_subset.shape[0], table.shape[0]])

for i in range(n_mcmc_chain):
    Y_predict = Y_predict_accepted_trace[i][0]
    for _ in range(len(Y_predict)):
        votes[_,Y_predict[_]] += 1

Y_predict_majority = np.argmax(votes, axis=1)

print votes.sum()

accuracy = sum(Y_subset == Y_predict_majority)*1.0/ data_subset.shape[0]
print "Chain % d accuracy on subset data: %.3f" % (chain+1,accuracy)

bins = table.shape[0]
plt.hist(Y_subset, bins, alpha=0.5, label='Y:cell type')
plt.hist(Y_predict_majority, bins, alpha=0.5, label='Z:prediction')
plt.legend(loc='upper right')
plt.show() 
L = min(len(_) for _ in accuracy_accepted_trace)
res = []
for i in range(L):
    res.append(np.array([_[i] for _ in accuracy_accepted_trace]).mean())
plt.plot(res)
# vote, and compute accuracy

threshold = -113888.425449/data.shape[0]*N  - 10000000000
# keep that samples whose likelihood is above threshold
burnt_samples = []
burnt_predictions = []
for i in range(len(accepts)):
    accepted_chain = accepts[i]
    likelihoods = logl_accepted_trace[i]
    predictions = Y_predict_accepted_trace[i]
    burnt_samples += [accepted_chain[_] for _ in range(22, len(accepted_chain))]
    burnt_predictions +=[classify_cells(data_subset, \
                    sample, table, cell_type_name2idx) for sample in burnt_samples[-1]]
    
# vote
votes = np.zeros([data_subset.shape[0], table.shape[0]])
for Y_predict in burnt_predictions:
    for _ in range(len(Y_predict)):
        votes[_,Y_predict[_]] += 1

Y_predict_majority = np.argmax(votes, axis=1)

print len(burnt_predictions)

accuracy = sum(Y_subset == Y_predict_majority)*1.0/ data_subset.shape[0]
print "Chain % d accuracy on subset data: %.3f" % (chain+1,accuracy)

bins = table.shape[0]
plt.hist(Y_subset, bins, alpha=0.5, label='Y:cell type')
plt.hist(Y_predict_majority, bins, alpha=0.5, label='Z:prediction')
plt.legend(loc='upper right')
plt.show() 
Y_predict = classify_cells(data, optimal_chain[-1], table, cell_type_name2idx)
accuracy = sum(Y == Y_predict)*1.0/ data.shape[0]
print accuracy


correct_table = table.copy()
correct_table.at['Mature B cells','CD38'] = -1.0
print correct_table.loc['Mature B cells']['CD38']
print table.loc['Pre B cells']['CD38']


Y_predict = classify_cells(data_subset, optimal_chain[-1], correct_table, cell_type_name2idx)
accuracy = sum(Y_subset== Y_predict)*1.0/ data_subset.shape[0]
print accuracy
from sklearn.metrics import confusion_matrix

def plot_confusion_matrix(cm, classes,
                          normalize=False,
                          cmap=plt.cm.Blues):
    """
    This function prints and plots the confusion matrix.
    Normalization can be applied by setting `normalize=True`.
    """
    if normalize:
        cm = cm.astype('float') / cm.sum(axis=1)[:, np.newaxis]
        print("Normalized confusion matrix")
    else:
        print('Confusion matrix, without normalization')

    plt.imshow(cm, interpolation='nearest', cmap=cmap)
    plt.colorbar()
    tick_marks = np.arange(len(classes))
    plt.xticks(tick_marks, classes, rotation=90,fontsize ="13")
    plt.yticks(tick_marks, classes,fontsize ="13")

    fmt = '.2f' if normalize else 'd'
    thresh = cm.max() / 2.
#     for i, j in itertools.product(range(cm.shape[0]), range(cm.shape[1])):
#         plt.text(j, i, format(cm[i, j], fmt),
#                  horizontalalignment="center",
#                  color="white" if cm[i, j] > thresh else "black")

    plt.tight_layout()
    plt.ylabel('True label',fontsize ="16")
    plt.xlabel('Predicted label',fontsize ="16")

# Compute confusion matrix
cnf_matrix = confusion_matrix(Y_subset, Y_predict_majority)
class_names = table.index

np.set_printoptions(precision=2)

# Plot non-normalized confusion matrix
# plt.figure()
# plot_confusion_matrix(cnf_matrix, classes=class_names,
#                       title='Confusion matrix, without normalization')

# Plot normalized confusion matrix
plt.figure(figsize=(10,8))
plot_confusion_matrix(cnf_matrix, classes=class_names, normalize=True)
plt.savefig('confusion_matrix.pdf', format='pdf', dpi=1000)
plt.show()
from sklearn.manifold import TSNE
tsne = TSNE(n_components=2, random_state=1234)
Xre = tsne.fit_transform(data_subset)
import matplotlib.pyplot as plt
import matplotlib as mpl   
import matplotlib.patches as mpatches

cmap = mpl.cm.Accent
fig = plt.figure(1, figsize=(11, 11))


idx2ct = [key for idx, key in enumerate(table.index)]
for idx, key in enumerate(idx2ct):
    plt.plot(Xre[Y_subset == idx, 0], Xre[Y_subset == idx, 1], '.',\
                 c=cmap(idx / float(len(idx2ct))), label = key, alpha=0.8)
       
patches = []
for idx, key in enumerate(idx2ct):
    patches.append(mpatches.Patch(color=cmap(idx / float(len(idx2ct))), label=key))
    
plt.legend(handles=patches,loc='upper center', bbox_to_anchor=(0.5, -0.05),\
           ncol=3,prop={'size':15})        
plt.xticks([])
plt.yticks([])
plt.savefig('tsne_true.pdf', format='pdf', dpi=1000)
import matplotlib.pyplot as plt
import matplotlib as mpl   
import matplotlib.patches as mpatches

cmap = mpl.cm.Accent
fig = plt.figure(1, figsize=(11, 11))


idx2ct = [key for idx, key in enumerate(table.index)]
for idx, key in enumerate(idx2ct):
    plt.plot(Xre[Y_predict == idx, 0], Xre[Y_predict == idx, 1], '.',\
                 c=cmap(idx / float(len(idx2ct))), label = key, alpha=0.8)
       
patches = []
for idx, key in enumerate(idx2ct):
    patches.append(mpatches.Patch(color=cmap(idx / float(len(idx2ct))), label=key))
    
plt.legend(handles=patches,loc='upper center', bbox_to_anchor=(0.5, -0.05),\
           ncol=3,prop={'size':15})        
plt.xticks([])
plt.yticks([])
plt.savefig('tsne_predict.pdf', format='pdf', dpi=1000)
print set(Y_predict),Y_predict.shape
print Y_subset.shape
print len(accepts_append)
for i in range(n_mcmc_chain):
    for j in range(len(accepts[i])):
        filename = 'AMML_accepted_samples/chain_%d_id_%d.pkl' %(i,j)
        with open(filename, "wb") as fp:
            pickle.dump(accepts[i][j], fp)<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.15073923875684228, 'fraction_numerical': 0.04240848169633927, 'mean_word_length': 3.819879042860899, 'pattern_counts': {'":': 0, '<': 66, '<?xml version=': 0, '>': 66, 'https://': 1, 'lorem ipsum': 0, 'www.': 1, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 2, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '12', 'dir': 'jupyter-structured-clean-dedup', 'id': '341503', 'license': 'MIT', 'n_tokens_mistral': 21887, 'n_tokens_neox': 20395, 'n_words': 3954, 'path': 'small_run/Flow_Cytometry_Mondrian_Processes-AML-5k.ipynb', 'repo_name': 'disiji/fc_mondrian', 'size': '54989'}