 b/aggmap/_devmap.py
+#!/usr/bin/env python3
+# -*- coding: utf-8 -*-
+"""
+Created on Sun Aug 25 20:29:36 2019
+@author: wanxiang.shen@u.nus.edu
+main aggmap code
+"""
+from aggmap.utils.logtools import print_info, print_warn, print_error
+from aggmap.utils.matrixopt import Scatter2Grid, Scatter2Array, smartpadding
+from aggmap.utils import vismap, summary, calculator
+from sklearn.cluster import AgglomerativeClustering
+from sklearn.manifold import TSNE, MDS
+from sklearn.utils import shuffle
+from joblib import Parallel, delayed, load, dump
+from scipy.spatial.distance import squareform
+from scipy.cluster.hierarchy import fcluster, linkage, dendrogram
+import matplotlib.pylab as plt
+import seaborn as sns
+from umap import UMAP
+from tqdm import tqdm
+from copy import copy
+import pandas as pd
+import numpy as np
+import os
+class Base:
+    def __init__(self):
+        pass
+    def _save(self, filename):
+        return dump(self, filename)
+    def _load(self, filename):
+        return load(filename)
+    def MinMaxScaleClip(self, x, xmin, xmax):
+        scaled = (x - xmin) / ((xmax - xmin) + 1e-8)
+        return scaled
+    def StandardScaler(self, x, xmean, xstd):
+        return (x-xmean) / (xstd + 1e-8)
+def _cluster_model2linkage_matrix(model):
+    counts = np.zeros(model.children_.shape[0])
+    n_samples = len(model.labels_)
+    for i, merge in enumerate(model.children_):
+        current_count = 0
+        for child_idx in merge:
+            if child_idx < n_samples:
+                current_count += 1  # leaf node
+            else:
+                current_count += counts[child_idx - n_samples]
+        counts[i] = current_count
+    linkage_matrix = np.column_stack([model.children_, model.distances_,
+                                      counts]).astype(float)
+    return linkage_matrix
+class AggMap(Base):
+    '''
+    Note: t-SNE initialize method should be changed into 'pca': https://www.nature.com/articles/s41587-020-00809-z
+    >>> mp = AggMap(dfx)
+    >>> mp.fit(emb_method = 'tsne', init = 'pca')
+    '''
+    def __init__(self,
+                 dfx,
+                 metric = 'correlation',
+                ):
+        """
+        paramters
+        -----------------
+        dfx: pandas DataFrame
+        metric: {'cosine', 'correlation', 'euclidean', 'jaccard', 'hamming', 'dice'}, default: 'correlation', measurement of feature distance
+        """
+        assert type(dfx) == pd.core.frame.DataFrame, 'input dfx mush be pandas DataFrame!'
+        super().__init__()
+        self.metric = metric
+        self.isfit = False
+        self.alist = dfx.columns.tolist()
+        self.ftype = 'feature points'
+        ## calculating distance
+        print_info('Calculating distance ...')
+        D = calculator.pairwise_distance(dfx.values, n_cpus=16, method=metric)
+        D = np.nan_to_num(D,copy=False)
+        D_ = squareform(D)
+        self.info_distance = D_.clip(0, np.inf)
+        ## statistic info
+        S = summary.Summary(n_jobs = 10)
+        res= []
+        for i in tqdm(range(dfx.shape[1]), ascii=True):
+            r = S._statistics_one(dfx.values, i)
+            res.append(r)
+        dfs = pd.DataFrame(res, index = self.alist)
+        self.info_scale = dfs
+        print_info('Applying the Agglomerative Clustering ...')
+        cluster_model = AgglomerativeClustering(distance_threshold=0, n_clusters=None)
+        cluster_model.fit(dfx.values.T)
+        self._intrinsic_Z = _cluster_model2linkage_matrix(cluster_model)
+    def _fit_embedding(self,
+                        dist_matrix,
+                        method = 'umap',
+                        n_components = 2,
+                        random_state = 32,
+                        verbose = 2,
+                        n_neighbors = 15,
+                        min_dist = 0.1,
+                        **kwargs):
+        """
+        parameters
+        -----------------
+        dist_matrix: distance matrix to fit
+        method: {'tsne', 'umap', 'mds'}, algorithm to embedd high-D to 2D
+        kwargs: the extra parameters for the conresponding algorithm
+        """
+        if 'metric' in kwargs.keys():
+            metric = kwargs.get('metric')
+            kwargs.pop('metric')
+        else:
+            metric = 'precomputed'
+        if method == 'tsne':
+            embedded = TSNE(n_components=n_components,
+                            random_state=random_state,
+                            metric = metric,
+                            verbose = verbose,
+                            **kwargs)
+        elif method == 'umap':
+            embedded = UMAP(n_components = n_components,
+                            n_neighbors = n_neighbors,
+                            min_dist = min_dist,
+                            verbose = verbose,
+                            random_state=random_state,
+                            metric = metric, **kwargs)
+        elif method =='mds':
+            if 'metric' in kwargs.keys():
+                kwargs.pop('metric')
+            if 'dissimilarity' in kwargs.keys():
+                dissimilarity = kwargs.get('dissimilarity')
+                kwargs.pop('dissimilarity')
+            else:
+                dissimilarity = 'precomputed'
+            embedded = MDS(metric = True,
+                           n_components= n_components,
+                           verbose = verbose,
+                           dissimilarity = dissimilarity,
+                           random_state = random_state, **kwargs)
+        embedded = embedded.fit(dist_matrix)
+        return embedded
+    def fit(self,
+            feature_group_list = [],
+            cluster_channels = 5,
+            var_thr = -1,
+            split_channels = True,
+            fmap_type = 'grid',
+            fmap_shape = None,
+            emb_method = 'umap',
+            min_dist = 0.1,
+            n_neighbors = 15,
+            verbose = 2,
+            random_state = 32,
+            group_color_dict  = {},
+            lnk_method = 'ward',
+            **kwargs):
+        """
+        parameters
+        -----------------
+        feature_group_list: list of the group name for each feature point
+        cluster_channels: int, number of the channels(clusters) if feature_group_list is empty
+        var_thr: float, defalt is -1, meaning that feature will be included only if the conresponding variance larger than this value. Since some of the feature has pretty low variances, we can remove them by increasing this threshold
+        split_channels: bool, if True, outputs will split into various channels using the types of feature
+        fmap_type:{'scatter', 'grid'}, default: 'gird', if 'scatter', will return a scatter mol map without an assignment to a grid
+        fmap_shape: None or tuple, size of molmap, only works when fmap_type is 'scatter', if None, the size of feature map will be calculated automatically
+        emb_method: {'tsne', 'umap', 'mds'}, algorithm to embedd high-D to 2D
+        group_color_dict: dict of the group colors, keys are the group names, values are the colors
+        lnk_method: {'ward','complete', 'average', 'single'}, linkage method
+        kwargs: the extra parameters for the conresponding embedding method
+        """
+        if 'n_components' in kwargs.keys():
+            kwargs.pop('n_components')
+        ## embedding  into a 2d
+        assert emb_method in ['tsne', 'umap', 'mds'], 'No Such Method Supported: %s' % emb_method
+        assert fmap_type in ['scatter', 'grid'], 'No Such Feature Map Type Supported: %s'   % fmap_type
+        self.var_thr = var_thr
+        self.split_channels = split_channels
+        self.fmap_type = fmap_type
+        self.fmap_shape = fmap_shape
+        self.emb_method = emb_method
+        self.lnk_method = lnk_method
+        if fmap_shape != None:
+            assert len(fmap_shape) == 2, "fmap_shape must be a tuple with two elements!"
+        # flist and distance
+        flist = self.info_scale[self.info_scale['var'] > self.var_thr].index.tolist()
+        dfd = pd.DataFrame(squareform(self.info_distance),
+                           index=self.alist,
+                           columns=self.alist)
+        dist_matrix = dfd.loc[flist][flist]
+        self.flist = flist
+        self.x_mean = self.info_scale['mean'].values
+        self.x_std =  self.info_scale['std'].values
+        self.x_min = self.info_scale['min'].values
+        self.x_max = self.info_scale['max'].values
+        #bitsinfo
+        dfb = pd.DataFrame(self.alist, columns = ['IDs'])
+        if feature_group_list != []:
+            self.Z = self._intrinsic_Z
+            assert len(feature_group_list) == len(self.alist), "the length of the input group list is not equal to length of the feature list"
+            self.cluster_channels = len(set(feature_group_list))
+            self.feature_group_list = feature_group_list
+            dfb['Subtypes'] = feature_group_list
+            if set(feature_group_list).issubset(set(group_color_dict.keys())):
+                self.group_color_dict = group_color_dict
+                dfb['colors'] = dfb['Subtypes'].map(group_color_dict)
+            else:
+                unique_types = dfb['Subtypes'].unique()
+                color_list = sns.color_palette("hsv", len(unique_types)).as_hex()
+                group_color_dict = dict(zip(unique_types, color_list))
+                dfb['colors'] = dfb['Subtypes'].map(group_color_dict)
+                self.group_color_dict = group_color_dict
+        else:
+            self.cluster_channels = cluster_channels
+            print_info('applying hierarchical clustering to obtain group information ...')
+            if self.lnk_method != 'ward':
+                Z = linkage(squareform(dfd.values),  lnk_method)
+            else:
+                Z = self._intrinsic_Z
+            labels = fcluster(Z, cluster_channels, criterion='maxclust')
+            feature_group_list = ['cluster_%s' % str(i).zfill(2) for i in labels]
+            dfb['Subtypes'] = feature_group_list
+            dfb = dfb.sort_values('Subtypes')
+            unique_types = dfb['Subtypes'].unique()
+            if not set(unique_types).issubset(set(group_color_dict.keys())):
+                color_list = sns.color_palette("hsv", len(unique_types)).as_hex()
+                group_color_dict = dict(zip(unique_types, color_list))
+            dfb['colors'] = dfb['Subtypes'].map(group_color_dict)
+            self.group_color_dict = group_color_dict
+            self.Z = Z
+            self.feature_group_list = feature_group_list
+        self.bitsinfo = dfb
+        colormaps = dfb.set_index('Subtypes')['colors'].to_dict()
+        colormaps.update({'NaN': '#000000'})
+        self.colormaps = colormaps
+        if fmap_type == 'grid':
+            S = Scatter2Grid()
+        else:
+            if fmap_shape == None:
+                N = len(self.flist)
+                l = np._int(np.sqrt(N))*2
+                fmap_shape = (l, l)
+            S = Scatter2Array(fmap_shape)
+        self._S = S
+        ## 2d embedding first
+        embedded = self._fit_embedding(dist_matrix,
+                                       method = emb_method,
+                                       n_neighbors = n_neighbors,
+                                       random_state = random_state,
+                                       min_dist = min_dist,
+                                       verbose = verbose,
+                                       n_components = 2, **kwargs)
+        self.embedded = embedded
+        df = pd.DataFrame(embedded.embedding_, index = self.flist,columns=['x', 'y'])
+        typemap = self.bitsinfo.set_index('IDs')
+        df = df.join(typemap)
+        df['Channels'] = df['Subtypes']
+        self.df_embedding = df
+        if self.fmap_type == 'scatter':
+            ## naive scatter algorithm
+            print_info('Applying naive scatter feature map...')
+            self._S.fit(self.df_embedding, self.split_channels, channel_col = 'Channels')
+            print_info('Finished')
+        else:
+            ## linear assignment algorithm
+            print_info('Applying grid feature map(assignment), this may take several minutes(1~30 min)')
+            self._S.fit(self.df_embedding, self.split_channels, channel_col = 'Channels')
+            print_info('Finished')
+        ## fit flag
+        self.isfit = True
+        if self.fmap_shape == None:
+            self.fmap_shape = self._S.fmap_shape
+        else:
+            m, n = self.fmap_shape
+            p, q = self._S.fmap_shape
+            assert (m >= p) & (n >=q), "fmap_shape's width must >= %s, height >= %s " % (p, q)
+        return self
+    def transform(self,
+                  arr_1d,
+                  scale = True,
+                  scale_method = 'minmax',
+                  fillnan = 0):
+        """
+        parameters
+        --------------------
+        arr_1d: 1d numpy array feature points
+        scale: bool, if True, we will apply MinMax scaling by the precomputed values
+        scale_method: {'minmax', 'standard'}
+        fillnan: fill nan value, default: 0
+        """
+        if not self.isfit:
+            print_error('please fit first!')
+            return
+        if scale:
+            if scale_method == 'standard':
+                arr_1d = self.StandardScaler(arr_1d, self.x_mean, self.x_std)
+            else:
+                arr_1d = self.MinMaxScaleClip(arr_1d, self.x_min, self.x_max)
+        df = pd.DataFrame(arr_1d).T
+        df.columns = self.alist
+        df = df[self.flist]
+        vector_1d = df.values[0] #shape = (N, )
+        fmap = self._S.transform(vector_1d)
+        p, q, c = fmap.shape
+        if self.fmap_shape != None:
+            m, n = self.fmap_shape
+            if (m > p) | (n > q):
+                fps = []
+                for i in range(c):
+                    fp = smartpadding(fmap[:,:,i], self.fmap_shape)
+                    fps.append(fp)
+                fmap = np.stack(fps, axis=-1)
+        return np.nan_to_num(fmap, nan = fillnan)
+    def batch_transform(self,
+                        array_2d,
+                        scale = True,
+                        scale_method = 'minmax',
+                        n_jobs=4,
+                        fillnan = 0):
+        """
+        parameters
+        --------------------
+        array_2d: 2D numpy array feature points, M(samples) x N(feature ponits)
+        scale: bool, if True, we will apply MinMax scaling by the precomputed values
+        scale_method: {'minmax', 'standard'}
+        n_jobs: number of parallel
+        fillnan: fill nan value, default: 0
+        """
+        if not self.isfit:
+            print_error('please fit first!')
+            return
+        assert type(array_2d) == np.ndarray, 'input must be numpy ndarray!'
+        assert array_2d.ndim == 2, 'input must be 2-D  numpy array!'
+        P = Parallel(n_jobs=n_jobs)
+        res = P(delayed(self.transform)(arr_1d,
+                                        scale,
+                                        scale_method,
+                                        fillnan) for arr_1d in tqdm(array_2d, ascii=True))
+        X = np.stack(res)
+        return X
+    def plot_scatter(self, htmlpath='./', htmlname=None, radius = 2, enabled_data_labels = False):
+        """radius: the size of the scatter, must be int"""
+        df_scatter, H_scatter = vismap.plot_scatter(self,
+                                                    htmlpath=htmlpath,
+                                                    htmlname=htmlname,
+                                                    radius = radius,
+                                                    enabled_data_labels = enabled_data_labels)
+        self.df_scatter = df_scatter
+        return H_scatter
+    def plot_grid(self, htmlpath='./', htmlname=None, enabled_data_labels = False):
+        if self.fmap_type != 'grid':
+            return
+        df_grid, H_grid = vismap.plot_grid(self,
+                                           htmlpath=htmlpath,
+                                           htmlname=htmlname,
+                                           enabled_data_labels = enabled_data_labels)
+        self.df_grid = df_grid
+        return H_grid
+    def plot_tree(self, figsize=(16,8), add_leaf_labels = True, leaf_font_size = 18, leaf_rotation = 90):
+        fig = plt.figure(figsize=figsize)
+        Z = self.Z
+        D_leaf_colors = self.bitsinfo['colors'].to_dict()
+        link_cols = {}
+        for i, i12 in enumerate(Z[:,:2].astype(int)):
+            c1, c2 = (link_cols[x] if x > len(Z) else D_leaf_colors[x] for x in i12)
+            link_cols[i+1+len(Z)] = c1
+        if add_leaf_labels:
+            labels = self.alist
+        else:
+            labels = None
+        P =dendrogram(Z, labels = labels,
+                      leaf_rotation = leaf_rotation,
+                      leaf_font_size = leaf_font_size,
+                      link_color_func=lambda x: link_cols[x])
+        return fig
+    def copy(self):
+        return copy(self)
+    def load(self, filename):
+        return self._load(filename)
+    def save(self, filename):
+        return self._save(filename)