--- a
+++ b/eval.py
@@ -0,0 +1,183 @@
+import argparse
+import metric
+from sklearn.cluster import KMeans
+from sklearn.metrics.cluster import normalized_mutual_info_score, adjusted_rand_score
+from sklearn.metrics.cluster import homogeneity_score, adjusted_mutual_info_score
+import numpy as np
+import random
+import sys,os
+from scipy.io import loadmat
+from sklearn.metrics import confusion_matrix
+import pandas as pd 
+import matplotlib
+matplotlib.use('agg')
+import matplotlib.pyplot as plt
+import seaborn as sns
+sns.set_style("whitegrid", {'axes.grid' : False})
+
+def plot_embedding(X, labels, classes=None, method='tSNE', cmap='tab20', figsize=(8, 8), markersize=15, dpi=300,marker=None,
+                   return_emb=False, save=False, save_emb=False, show_legend=True, show_axis_label=True, **legend_params):
+    if marker is not None:
+        X = np.concatenate([X, marker], axis=0)
+    N = len(labels)
+    matplotlib.rc('xtick', labelsize=20) 
+    matplotlib.rc('ytick', labelsize=20) 
+    matplotlib.rcParams.update({'font.size': 22})
+    if X.shape[1] != 2:
+        if method == 'tSNE':
+            from sklearn.manifold import TSNE
+            X = TSNE(n_components=2, random_state=124).fit_transform(X)
+        if method == 'PCA':
+            from sklearn.decomposition import PCA
+            X = PCA(n_components=2, random_state=124).fit_transform(X)
+        if method == 'UMAP':
+            from umap import UMAP
+            X = UMAP(n_neighbors=15, min_dist=0.1, metric='correlation').fit_transform(X)
+    labels = np.array(labels)
+    plt.figure(figsize=figsize)
+    if classes is None:
+        classes = np.unique(labels)
+    #tab10, tab20, husl, hls
+    if cmap is not None:
+        cmap = cmap
+    elif len(classes) <= 10:
+        cmap = 'tab10'
+    elif len(classes) <= 20:
+        cmap = 'tab20'
+    else:
+        cmap = 'husl'
+    colors = sns.husl_palette(len(classes), s=.8)
+    #markersize = 80
+    for i, c in enumerate(classes):
+        plt.scatter(X[:N][labels==c, 0], X[:N][labels==c, 1], s=markersize, color=colors[i], label=c)
+    if marker is not None:
+        plt.scatter(X[N:, 0], X[N:, 1], s=10*markersize, color='black', marker='*')
+    
+    legend_params_ = {'loc': 'center left',
+                     'bbox_to_anchor':(1.0, 0.45),
+                     'fontsize': 20,
+                     'ncol': 1,
+                     'frameon': False,
+                     'markerscale': 1.5
+                    }
+    legend_params_.update(**legend_params)
+    if show_legend:
+        plt.legend(**legend_params_)
+    sns.despine(offset=10, trim=True)
+    if show_axis_label:
+        plt.xlabel(method+' dim 1', fontsize=12)
+        plt.ylabel(method+' dim 2', fontsize=12)
+
+    if save:
+        plt.savefig(save, format='png', bbox_inches='tight',dpi=dpi)
+
+def cluster_eval(labels_true,labels_infer):
+    purity = metric.compute_purity(labels_infer, labels_true)
+    nmi = normalized_mutual_info_score(labels_true, labels_infer)
+    ari = adjusted_rand_score(labels_true, labels_infer)
+    homogeneity = homogeneity_score(labels_true, labels_infer)
+    ami = adjusted_mutual_info_score(labels_true, labels_infer)
+    print('NMI = {}, ARI = {}, Purity = {},AMI = {}, Homogeneity = {}'.format(nmi,ari,purity,ami,homogeneity))
+    return nmi,ari,homogeneity
+
+def get_best_epoch(exp_dir, dataset, measurement='NMI'):
+    results = []
+    for each in os.listdir('results/%s/%s'%(dataset,exp_dir)):
+        if each.startswith('data'):
+            #print('results/%s/%s/%s'%(dataset,exp_dir,each))
+            data = np.load('results/%s/%s/%s'%(dataset,exp_dir,each))
+            data_x_onehot_,label_y = data['arr_1'],data['arr_2']
+            label_infer = np.argmax(data_x_onehot_, axis=1)
+            nmi,ari,homo = cluster_eval(label_y,label_infer)
+            results.append([each,nmi,ari,homo])
+    if measurement == 'NMI':
+        results.sort(key=lambda a:-a[1])
+    elif measurement == 'ARI':
+        results.sort(key=lambda a:-a[2])
+    elif measurement == 'HOMO':
+        results.sort(key=lambda a:-a[3])
+    else:
+        print('Wrong indicated metric')
+        sys.exit()
+    print('NMI = {}\tARI = {}\tHomogeneity = {}'.format(results[0][1],results[0][2],results[0][3]))
+    return results[0][0]
+
+def save_embedding(emb_feat,save,sep='\t'):
+    index = ['cell%d'%(i+1) for i in range(emb_feat.shape[0])]
+    columns = ['feat%d'%(i+1) for i in range(emb_feat.shape[1])]
+    data_pd = pd.DataFrame(emb_feat,index = index,columns=columns)
+    data_pd.to_csv(save,sep=sep)
+
+def save_clustering(label,save):
+    f = open(save,'w')
+    res_list = ['cell%d\t%s'%(i,str(item)) for i,item in enumerate(label)]
+    f.write('\n'.join(res_list))
+    f.close()
+
+if __name__ == '__main__':
+        parser = argparse.ArgumentParser(description='Simultaneous deep generative modeling and clustering of single cell genomic data')
+        parser.add_argument('--data', '-d', type=str, help='which dataset')
+        parser.add_argument('--timestamp', '-t', type=str, help='timestamp')
+        parser.add_argument('--epoch', '-e', type=int, help='epoch or batch index')
+        parser.add_argument('--train', type=bool, default=False)
+        parser.add_argument('--save', '-s', type=str, help='save latent visualization plot (e.g., t-SNE)')
+        parser.add_argument('--no_label', action='store_true',help='whether the dataset has label')
+        args = parser.parse_args()
+        has_label = not args.no_label
+        if has_label:
+            if args.train:
+                exp_dir = [item for item in os.listdir('results/%s'%args.data) if item.startswith(args.timestamp)][0]
+                if args.epoch is None:
+                    epoch = get_best_epoch(exp_dir,args.data,'ARI')
+                else:
+                    epoch = args.epoch
+                data = np.load('results/%s/%s/%s'%(args.data,exp_dir,epoch))
+                embedding, label_infered_onehot = data['arr_0'],data['arr_1']
+                embedding_before_softmax = embedding[:,-label_infered_onehot.shape[1]:]
+                label_infered = np.argmax(label_infered_onehot, axis=1)
+                label_true = [item.strip() for item  in open('datasets/%s/label.txt'%args.data).readlines()]
+                save_clustering(label_infered,save='results/%s/%s/scDEC_cluster.txt'%(args.data,exp_dir))
+                save_embedding(embedding,save='results/%s/%s/scDEC_embedding.csv'%(args.data,exp_dir),sep='\t')
+                plot_embedding(embedding,label_true,save='results/%s/%s/scDEC_embedding.png'%(args.data,exp_dir))
+            else:
+                if args.data == 'PBMC10k':
+                    data = np.load('results/%s/data_pre.npz'%args.data)
+                    embedding, label_infered_onehot = data['arr_0'],data['arr_1']
+                    embedding_before_softmax = embedding[:,-label_infered_onehot.shape[1]:]
+                    label_infered = np.argmax(label_infered_onehot, axis=1)
+                    barcode2label = {item.split('\t')[0]:item.split('\t')[1].strip() for item in open('datasets/%s/labels_annot.txt'%args.data).readlines()[1:]}
+                    barcodes = [item.strip() for item in open('datasets/%s/barcodes.tsv'%args.data).readlines()]
+                    labels_annot = [barcode2label[item] for i,item in enumerate(barcodes) if item in barcode2label.keys()]
+                    select_idx = [i for i,item in enumerate(barcodes) if item in barcode2label.keys()]
+                    embedding = embedding[select_idx,:] # only evaluated on cells with annotation labels
+                    label_infered = label_infered[select_idx]
+                    uniq_label = list(np.unique(labels_annot))
+                    Y = np.array([uniq_label.index(item) for item in labels_annot])
+                    cluster_eval(Y,label_infered)
+                    save_clustering(label_infered,save='results/%s/scDEC_cluster.txt'%args.data)
+                    save_embedding(embedding,save='results/%s/scDEC_embedding.csv'%args.data,sep='\t')
+                    plot_embedding(embedding,labels_annot,save='results/%s/scDEC_embedding.png'%args.data)
+                else:
+                    data = np.load('results/%s/data_pre.npz'%args.data)
+                    embedding, label_infered_onehot = data['arr_0'],data['arr_1']
+                    embedding_before_softmax = embedding[:,-label_infered_onehot.shape[1]:]
+                    label_infered = np.argmax(label_infered_onehot, axis=1)
+                    label_true = [item.strip() for item  in open('datasets/%s/label.txt'%args.data).readlines()]
+                    save_clustering(label_infered,save='results/%s/scDEC_cluster.txt'%args.data)
+                    save_embedding(embedding,save='results/%s/scDEC_embedding.csv'%args.data,sep='\t')
+                    plot_embedding(embedding,label_true,save='results/%s/scDEC_embedding.png'%args.data)
+        else:
+            if args.epoch is None:
+                print('Provide the epoch or batch index to analyze')
+                sys.exit()
+            else:
+                exp_dir = [item for item in os.listdir('results/%s'%args.data) if item.startswith(args.timestamp)][0]
+                data = np.load('results/%s/%s/data_at_%s.npz'%(args.data,exp_dir,args.epoch))
+                embedding, label_infered_onehot = data['arr_0'],data['arr_1']
+                label_infered = np.argmax(label_infered_onehot, axis=1)
+                save_clustering(label_infered,save='results/%s/%s/scDEC_cluster.txt'%(args.data,exp_dir))
+                
+
+
+    
+