diff --git a/ITISCFunction.m b/ITISCFunction.m
new file mode 100644
index 0000000..4977ead
--- /dev/null
+++ b/ITISCFunction.m
@@ -0,0 +1,27 @@
+function history= ITISCFunction(filename)
+load(filename,'initCenters','data','T1','T2','outputFile')
+x0=initCenters;
+history.x = [];
+history.fval = [];
+disp('Init center is');
+disp(x0) 
+
+options = optimoptions(@fminunc,'OutputFcn',@outfun,'Display','iter','Algorithm','quasi-newton','HessianApproximation','bfgs');
+[xsol,fval,exitflag,output] = fminunc(@objfun,x0,options) 
+ function stop = outfun(x,optimValues,state)
+     stop = false;
+     switch state
+         case 'iter'
+           history.fval = [history.fval; optimValues.fval];
+           history.x = cat(3,history.x, x); 
+         otherwise
+     end
+ end
+ 
+ function f = objfun(x)
+          f =T2.*log(sum(sum((pdist2(data, x).^2).^(-1/T1),2).^(-T1/T2)));
+ end
+
+historyCenters = permute(history.x,[3 1 2]);
+save(outputFile,"historyCenters","exitflag") 
+end
diff --git a/ITISC_AO.py b/ITISC_AO.py
new file mode 100644
index 0000000..8194aaa
--- /dev/null
+++ b/ITISC_AO.py
@@ -0,0 +1,278 @@
+# ITISC Alternative Optimization algorithm
+import math
+import random
+import numpy as np
+from scipy.linalg import norm
+from scipy.spatial.distance import cdist # l2: default
+from math import sqrt
+import matplotlib.pyplot as plt
+import imageio
+
+class ITISC_AO:
+    def __init__(self, args, initCenters=None):
+        self.T1 = args.T1
+        self.T2 = args.T2
+        self.max_iter = args.max_iter
+        assert self.T1 > 0
+        assert self.T2 > 0
+        self.T2_list = [] 
+        self.u, self.centers = None, None
+        self.n_clusters = args.C
+        self.max_iter = args.max_iter
+        self.error = args.error
+        self.itisc_random_state = args.itisc_random_state
+        self.args = args
+
+        # intermediate step
+        self.all_W = []
+        self.all_centers = []
+        self.all_pred = []
+        self.convergeOrNot = None 
+        
+        # save gif
+        self.gif = args.gif
+        self.gif_gap = args.gif_gap
+        self.prob = args.prob
+        self.image_list = []
+        self.initCenters = initCenters
+        self.init_method = args.init_method
+        self.plot_boundary = args.plot_boundary # isoprobability curve
+
+    def fit(self, X):
+        self.n_samples = X.shape[0]
+        self.r = np.random.RandomState(self.itisc_random_state)
+        self.U = self.r.rand(self.n_samples,self.n_clusters)
+        self.U = self.U / np.tile(self.U.sum(axis=1)[np.newaxis].T, self.n_clusters)
+        self.W = self.r.rand(self.n_samples,1)
+        self.W = self.W / np.sum(self.W) # normalize
+        self.all_W = [self.W]
+        
+        for iteration in range(self.max_iter):
+            if iteration == 0:
+               self.centers = self.init_centers(X)
+               self.all_centers = [self.centers]
+
+            U_old = self.U.copy()
+            centers_old = self.centers.copy()
+            self.U = self.next_u(X)
+            
+            self.W = self.next_w(X)
+            self.all_W = self.all_W + [self.W]
+
+            self.centers = self.next_centers(X, iteration)
+            self.all_centers = self.all_centers + [self.centers]
+            self.all_pred = self.all_pred + [self.predict(X)]
+            # Stopping rule
+            if not np.isnan(self.centers).any():
+               # if norm(self.U - U_old) < self.error:
+               if norm(self.centers - centers_old) < self.error: # same as ITISC_R
+                   print('Centers Y not changed, Break iteration={}'.format(iteration))
+                   break
+            else:
+                print('Centers Y contains nan value, break at iteration={}'.format(iteration))
+                break
+
+            ## save intermediate step
+            if self.gif:
+                if iteration % self.gif_gap == 0: 
+                    image, fig= self.get_gif(X, iteration, plot_boundary=self.plot_boundary)
+                    fig.clf()
+                    self.image_list = self.image_list + [image]
+        if self.gif:
+           imageio.mimsave('./gif/ITISC_AO_T1={}_T2={}.gif'.format(self.T1, self.T2), self.image_list, fps=1)  
+ 
+        print('ITISC_AO break at iteration:{}, T1={:.2f}, T2={:.2f}'.format(iteration, self.T1, self.T2))
+        print('iteration={}, centers_diff={}'.format(iteration, norm(self.centers - centers_old)))
+        print('center norm={}'.format(norm(self.centers - centers_old)))
+        self.convergeOrNot = norm(self.centers - centers_old) < self.error
+        print('The ITISC_AO algorithm converge={}'.format(self.convergeOrNot))
+        return self.T1, self.T2, self.convergeOrNot
+      
+    def init_centers(self, X):
+        # provided init Centers: may come from another clustering algorithm
+        if not self.initCenters is None: 
+           initCenters = self.initCenters
+        elif self.init_method == 'random': # random initialization from data range
+            random.seed(self.args.itisc_random_state) 
+            xmax = np.max(X, axis=0)
+            xmin = np.min(X, axis=0)
+            initCenters = []
+            for c in range(self.n_clusters):
+                tmp = [random.uniform(tmp_min, tmp_max) for tmp_min, tmp_max in zip(xmin, xmax)]
+                initCenters.append(tmp)
+            initCenters = np.array(initCenters)
+        elif self.init_method == 'kmeans':
+            from sklearn.cluster import KMeans
+            kmeans = KMeans(n_clusters=self.n_clusters, random_state=self.args.random_state).fit(X)
+            initCenters = kmeans.cluster_centers_ 
+        elif self.init_method == 'fcm':
+            from fcm_github import FCM
+            fcm = FCM(self.args)
+            fcm.fit(X)
+            initCenters = fcm.centers
+        else:
+            exit('Unrecognized init center method')
+        return initCenters
+
+    def next_centers(self, X, iteration):
+        # based on updated u_{ij} and w_{i}
+        # xi: observation
+        # yj: current center
+        # i : index for observation，total N
+        # j : index for center，total C
+        # U : [N,C] # numerator:[N,C] denominator:[N,1] --> update u_{ij} (公式1)
+        # W : [N,1] # numerator:[N,1] denominator:[1,1] --> update w_{i}  (公式2)
+        # X:  [N,R] # R: feature dimension
+        W = self.W ** (1-self.T2) 
+        U = self.U ** (1+self.T1)
+        D = cdist(X, self.centers)**2 # [N,C] 
+        
+        ### xi overlap with cluster center
+        if 0 in D:
+            zero_ind = np.argwhere(D==0) 
+            W = np.delete(W, zero_ind[:,0], axis=0)
+            U = np.delete(U, zero_ind[:,0], axis=0)
+            D = np.delete(D, zero_ind[:,0], axis=0)
+            X = np.delete(X, zero_ind[:,0], axis=0)
+            print('Updating centers, center and data overlap, len={}'.format(len(zero_ind)))
+                
+        denom = (U.T) @ W # [C,N] * [N,1] --> [C,1]
+        W = W.repeat(X.shape[-1], axis=1) # [N,1] --> [N,R]
+        numer = (U.T) @ (W * X) # [C,N] * [N,R] --> [C,R]
+
+        new_center = numer / denom
+        return new_center
+    
+    def next_u(self, X):
+        alpha = 1/self.T1
+        temp = cdist(X, self.centers)**2  # [400,3] --> [400,1,3] -->[400,3,3]
+        
+        if 0 in temp:
+           zero_ind_exists = True
+           zero_ind = np.argwhere(temp == 0)
+           print('Updating U, center and data overlap, len={}'.format(len(zero_ind)))
+           for ind in zero_ind:
+               temp[ind[0], ind[1]] = 1
+        else:
+            zero_ind_exists = False
+           
+        denominator = temp.reshape((X.shape[0], 1, -1)).repeat(temp.shape[-1], axis=1)
+        denominator = temp[:, :, np.newaxis] / denominator # [400,3,1]/[400,3,3] --> [400,3,3]
+        denominator = denominator ** alpha
+        res = 1/denominator.sum(2)
+        if zero_ind_exists:
+           for ind in zero_ind:
+               res[ind[0], ...] = 0
+               res[ind[0],ind[1]] = 1
+        return res
+    
+    ## Trick from FCM for stable computation: put numer to denom
+    def next_w(self, X):
+        # based on updated u_{ij} and w_{i}
+        # xi: observation
+        # yj: current center
+        # i : index for observation，total N
+        # j : index for center，total C
+        # U : [N,C] # numerator:[N,C] denominator:[N,1] --> update u_{ij} (公式1)
+        # W : [N,1] # numerator:[N,1] denominator:[1,1] --> update w_{i}  (公式2)
+        # X:  [N,R] # R: feature dimension
+        alpha = -1 / self.T1
+        beta = self.T1 / self.T2
+        temp = cdist(X, self.centers)**2 
+        
+        if 0 in temp:
+           zero_ind_exists = True
+           zero_ind = np.argwhere(temp == 0)
+           print('Updating W, center and data overlap, len={}'.format(len(zero_ind)))
+           for ind in zero_ind:
+               temp[ind[0], ind[1]] = 1
+        else:
+            zero_ind_exists = False
+
+        temp = temp ** alpha
+        temp = np.sum(temp, axis=1) # [N,C]
+        denom = temp.reshape((X.shape[0], -1)).repeat(temp.shape[0], axis=1)
+        res = (temp / denom) ** beta
+        res = 1 / np.sum(res, axis=0)
+
+        if zero_ind_exists:
+           for ind in zero_ind:
+               res[ind[0], ...] = 0
+        return res[:,np.newaxis]
+
+    # predict: original data
+    def predict(self, X):
+        if len(X.shape) == 1:
+            X = np.expand_dims(X, axis=0)
+
+        U = self.next_u(X) 
+        return np.argmax(U, axis=-1)
+    
+    # predict: new data
+    def predict_new(self, X):
+        if len(X.shape) == 1:
+            X = np.expand_dims(X, axis=0)
+        return self.next_u(X) 
+
+    def get_U(self):
+        return self.U
+
+    def get_W(self):
+        return self.W
+
+    def get_gif(self, X, iter, plot_boundary=False):
+        fig, ax = plt.subplots(1,1, figsize=(8,8))
+        colors_tmp = ['g', 'r', 'b', 'c', 'm', 'k','y','cyan','pink','gray']
+        point_size = 3
+        alpha = 0.8
+        color_tmp = colors_tmp[:self.n_clusters]
+        pred = self.predict(X)
+
+        if  plot_boundary:
+            rtol = 1e-2
+            mesh_grid_pred = get_meshgrid(X)
+            mesh_grid_U_pred = self.predict_new(mesh_grid_pred)
+            cls_ind_list = []
+            for ind in range(self.n_clusters):
+                cls_ind_tmp = get_boundary_index(mesh_grid_U_pred, ind, self.prob, rtol=rtol)
+                cls_ind_list = cls_ind_list + [cls_ind_tmp]
+        
+        ax.plot(self.centers[:,0],  self.centers[:,1], "*" , markersize = 10, c = 'blue') 
+
+        for ind in range(self.n_clusters):
+            ax.scatter(X[pred==ind,0], X[pred==ind,1], s=5, alpha=0.5, c=colors_tmp[ind])
+            if plot_boundary:
+               ax.scatter(mesh_grid_pred[cls_ind_list[ind],0],mesh_grid_pred[cls_ind_list[ind],1], c = colors_tmp[ind], s=point_size, alpha=alpha)
+
+        ax.set_title('iter={}, T1={:.2f}, T2={:.2f}'.format(iter, self.T1, self.T2))
+        fig.canvas.draw()       
+        image = np.frombuffer(fig.canvas.tostring_rgb(), dtype='uint8')
+        image  = image.reshape(fig.canvas.get_width_height()[::-1] + (3,))
+        return image, fig
+
+
+if __name__ == '__main__':    
+    from hyperparam import args
+    from utility_code import get_data, generate_data_ind, get_meshgrid, get_boundary_index
+
+    # data
+    args.T2 = 0.7
+    args.C = 3
+    args.R = 1
+    args = get_data(args)
+    data, labels, true_centers, data_mean, dataDistDesInd = generate_data_ind(args)
+    
+    # model
+    model = ITISC_AO(args)
+    model.fit(data)
+    centers = model.centers
+    pred = model.predict(data)
+    print('centers\n', centers)
+    
+    if not args.gif: 
+        plt.scatter(data[:,0],data[:,1], c=pred, s=10)
+        plt.scatter(centers[:,0],centers[:,1], c='k', s=30)
+        plt.show()      
+
+
+
diff --git a/ITISC_R.py b/ITISC_R.py
new file mode 100644
index 0000000..7a37749
--- /dev/null
+++ b/ITISC_R.py
@@ -0,0 +1,208 @@
+### ITISC-R: MATLAB
+import os 
+import random
+import imageio
+import numpy as np
+import scipy.io
+from scipy.spatial.distance import cdist
+from scipy.linalg import norm
+import matlab.engine # import after scipy
+import matplotlib.pyplot as plt
+
+class ITISC_R:
+    def __init__(self, args, initCenters=None):
+        self.T1 = args.T1
+        self.T2 = args.T2
+        assert self.T1 > 0
+        assert self.T2 > 0
+        self.C = args.C
+        self.centers = None
+        self.historyCenters = None
+        self.itisc_random_state = args.itisc_random_state
+        self.kmeans_random_state = args.kmeans_random_state
+        self.fcm_random_state = args.fcm_random_state
+        self.init_method = args.init_method
+        self.initCenters = initCenters
+        
+        # 初始化
+        self.itisc_random_state = args.itisc_random_state
+        self.kmeans_random_state = args.kmeans_random_state
+        self.args = args # FCM
+
+        # 画图
+        self.gif = args.gif
+        self.plot_boundary = args.plot_boundary # isoprobability curve
+        self.prob = args.prob
+        
+
+    def init_centers(self, X):
+        if not self.initCenters is None:
+           initCenters = self.initCenters
+        elif self.init_method == 'random':
+            random.seed(self.itisc_random_state)
+            xmax = np.max(X, axis=0)
+            xmin = np.min(X, axis=0)
+            initCenters = []
+            for c in range(self.C):
+                tmp = [random.uniform(tmp_min, tmp_max) for tmp_min, tmp_max in zip(xmin, xmax)]
+                initCenters.append(tmp)
+            self.initCenters = np.array(initCenters)
+        elif self.init_method == 'kmeans':
+            from sklearn.cluster import KMeans
+            kmeans = KMeans(n_clusters=self.C, random_state=self.kmeans_random_state).fit(X)
+            self.initCenters = kmeans.cluster_centers_ 
+        elif self.init_method == 'fcm':
+            from fcm_github import FCM
+            fcm = FCM(self.args)
+            fcm.fit(X)
+            self.initCenters = fcm.centers
+        return self.initCenters
+            
+    def fit(self, X):
+        initCenters = self.init_centers(X)  
+        matlabFileName = 'matlabCenters/MultiGaussian.mat'
+        pythonFileName = 'matlabCenters/matlabResult.mat'
+        matlabDir = os.path.join(os.getcwd(),'models')
+        outputFile = 'matlabCenters/matlabResult'
+        scipy.io.savemat(matlabFileName, dict(initCenters = initCenters, data = X, T1 = self.T1, T2 = self.T2, outputFile=outputFile)) # 将数据转化为matlab的格式
+        eng = matlab.engine.start_matlab()
+        eng.addpath(matlabDir, nargout=0) # add dir to the matlab search directory 
+        eng.ITISCFunction(matlabFileName)  # Matlab quasi-newton
+        mat = scipy.io.loadmat(pythonFileName) # load matlab data
+        self.historyCenters = mat['historyCenters']
+        # exitflag=0:other status，exitflag=1:converge
+        self.exitflag = mat['exitflag'][0][0] # matlab exit flag
+        self.centers = self.historyCenters[-1]
+        self.U = self.next_u(X) # according to current centers
+        self.W = self.next_w(X) # according to current centers 
+        if self.gif:
+           self.get_gif(X)
+        return self.centers
+      
+    def next_u(self, X):
+        # based on updated u_{ij} and w_{i}
+        # xi: observation
+        # yj: current center
+        # i : index for observation，total N
+        # j : index for center，total C
+        # U : [N,C] # numerator:[N,C] denominator:[N,1] --> update u_{ij} (公式1)
+        # W : [N,1] # numerator:[N,1] denominator:[1,1] --> update w_{i}  (公式2)
+        # X:  [N,R] # R: feature dimension
+        alpha = 1/self.T1
+        temp = cdist(X, self.centers)**2  # [400,3] --> [400,1,3] -->[400,3,3]
+        denominator = temp.reshape((X.shape[0], 1, -1)).repeat(temp.shape[-1], axis=1)
+        denominator = temp[:, :, np.newaxis] / denominator # [400,3,1]/[400,3,3] --> [400,3,3]
+        denominator = denominator ** alpha
+        res = 1/denominator.sum(2)
+        return res
+    
+    def next_w(self, X):
+        alpha = -1 / self.T1
+        beta = self.T1 / self.T2
+        temp = cdist(X, self.centers)**2 
+        temp = temp ** alpha
+        temp = np.sum(temp, axis=1) 
+        denom = temp.reshape((X.shape[0], -1)).repeat(temp.shape[0], axis=1)
+        res = (temp / denom) ** beta 
+        res = 1 / np.sum(res, axis=0)
+        return res[:,np.newaxis]
+
+    # predict original dataset
+    def predict(self, X):
+        if len(X.shape) == 1:
+            X = np.expand_dims(X, axis=0)
+        U = self.next_u(X) 
+        return np.argmax(U, axis=-1)
+
+    # predict new dataset
+    def predict_new(self, X):
+        if len(X.shape) == 1:
+            X = np.expand_dims(X, axis=0)
+        return self.next_u(X)
+
+    def get_U(self):
+        return self.U
+
+    def get_W(self):
+        return self.W
+
+    def get_gif_iter(self, X, iter):
+        def u_given_center(X, centers):
+            if len(X.shape) == 1:
+                X = np.expand_dims(X, axis=0)
+            # same as next_u with given center
+            alpha = 1/self.T1
+            temp = cdist(X, centers)**2 
+            denominator = temp.reshape((X.shape[0], 1, -1)).repeat(temp.shape[-1], axis=1)
+            denominator = temp[:, :, np.newaxis] / denominator # [N,C,1]/[N,C,C] --> [N,C,C]
+            denominator = denominator ** alpha
+            res = 1/denominator.sum(2)
+            return res
+        
+        fig, ax = plt.subplots(1,1, figsize=(8,8))
+        point_size = 3
+        centersize = 50
+        alpha = 0.8
+        colors_tmp = ['g', 'r', 'b', 'c', 'm', 'k','y','cyan','pink','gray']
+        color_tmp = colors_tmp[:self.C]
+        
+        x_U_pred = u_given_center(X, self.historyCenters[iter,...]) # data pred at current iteration
+        pred = np.argmax(x_U_pred, axis=-1)
+        
+        if  self.plot_boundary:
+            rtol = 1e-2
+            mesh_grid_pred = get_meshgrid(X)
+            mesh_grid_U_pred = u_given_center(mesh_grid_pred, self.historyCenters[iter,...])
+            cls_ind_list = []
+            for ind in range(self.C):
+                cls_ind_tmp = get_boundary_index(mesh_grid_U_pred, ind, self.prob, rtol=rtol)
+                cls_ind_list = cls_ind_list + [cls_ind_tmp]
+                
+        ax.scatter(self.historyCenters[iter,:,0],  self.historyCenters[iter,:,1], marker="*" , s = centersize, c = 'blue', label='ITISC-R')
+        for ind in range(self.C):
+            ax.scatter(X[pred==ind,0], X[pred==ind,1], s=5, alpha=0.5, c=colors_tmp[ind])
+            if self.plot_boundary:
+               ax.scatter(mesh_grid_pred[cls_ind_list[ind],0],mesh_grid_pred[cls_ind_list[ind],1], c = colors_tmp[ind], s=point_size, alpha=alpha)
+               
+        ax.set_title('iter={},T1={},T2={}'.format(iter,self.T1, self.T2))
+        ax.legend()
+        fig.canvas.draw()
+        image = np.frombuffer(fig.canvas.tostring_rgb(), dtype='uint8')
+        image  = image.reshape(fig.canvas.get_width_height()[::-1] + (3,))
+        return image, fig
+    
+    def get_gif(self, X):
+        image_list = []
+        for iteration in range(self.historyCenters.shape[0]):
+            image, fig= self.get_gif_iter(X, iteration)
+            image_list = image_list + [image]
+            fig.clf()
+        imageio.mimsave('./gif/ITISC_R_T1={}_T2={}.gif'.format(self.T1, self.T2), image_list, fps=1)
+
+
+if __name__ == '__main__':
+    from hyperparam import args
+    from utility_code import get_data, generate_data_ind, get_meshgrid, get_boundary_index
+    
+    # data
+    args.T2 = 0.1
+    args.C = 3
+    args.R = 1
+    args = get_data(args)
+    data, labels, true_centers, data_mean, dataDistDesInd = generate_data_ind(args)
+    
+    # model
+    model = ITISC_R(args)
+    model.fit(data)
+    centers = model.centers
+    pred = model.predict(data)
+    print('centers\n', centers)
+    
+    if not args.gif:
+        plt.scatter(data[:,0],data[:,1], c=pred, s=10)
+        plt.scatter(centers[:,0],centers[:,1], c='k', s=30)
+        plt.show()  
+
+
+
+
diff --git a/fcm_github.py b/fcm_github.py
new file mode 100644
index 0000000..09f743b
--- /dev/null
+++ b/fcm_github.py
@@ -0,0 +1,234 @@
+# code-from：https://github.com/omadson/fuzzy-c-means/blob/master/fcmeans/fcm.py
+import argparse
+import imageio
+import numpy as np
+from scipy.linalg import norm
+from scipy.spatial.distance import cdist # l2: default
+import matplotlib.pyplot as plt
+import pandas as pd
+
+class FCM:
+    """Fuzzy C-means
+    
+    Parameters
+    ----------
+    n_clusters: int, optional (default=10)
+        The number of clusters to form as well as the number of
+        centroids to generate
+
+    max_iter: int, optional (default=150)
+        Hard limit on iterations within solver.
+
+    m: float, optional (default=2.0)
+        Exponent for the fuzzy partition matrix, specified as a
+        scalar greater than 1.0. This option controls the amount of
+        fuzzy overlap between clusters, with larger values indicating
+        a greater degree of overlap.
+
+
+    error: float, optional (default=1e-5)
+        Tolerance for stopping criterion.
+
+    random_state: int, optional (default=42)
+        Determines random number generation for centroid initialization. Use
+        an int to make the randomness deterministic.
+
+    Attributes
+    ----------
+    n_samples: int
+        Number of examples in the data set
+
+    n_features: int
+        Number of features in samples of the data set
+
+    u: array, shape = [n_samples, n_clusters]
+        Fuzzy partition array, returned as an array with n_samples rows
+        and n_clusters columns. Element u[i,j] indicates the degree of
+        membership of the jth data point in the ith cluster. For a given
+        data point, the sum of the membership values for all clusters is one.
+
+    centers: array, shape = [n_class-1, n_SV]
+        Final cluster centers, returned as an array with n_clusters rows
+        containing the coordinates of each cluster center. The number of
+        columns in centers is equal to the dimensionality of the data being
+        clustered.
+
+    r: 
+    Container for the Mersenne Twister pseudo-random number generator.
+    
+    Methods
+    -------
+    fit(X)
+        fit the data
+
+    _predict(X)
+        use fitted model and output cluster memberships
+
+    predict(X)
+        use fitted model and output 1 cluster for each sample
+
+    References
+    ----------
+    .. [1] `Pattern Recognition with Fuzzy Objective Function Algorithms
+        <https://doi.org/10.1007/978-1-4757-0450-1>`_
+    .. [2] `FCM: The fuzzy c-means clustering algorithm
+        <https://doi.org/10.1016/0098-3004(84)90020-7>`_
+
+    """
+    def __init__(self, args):
+        self.u, self.centers = None, None
+        self.n_clusters = args.C
+        self.max_iter = args.max_iter
+        self.m = args.m
+        self.error = args.error
+        self.random_state = args.fcm_random_state
+        assert self.m > 1
+
+        # save gif
+        self.gif = args.gif
+        self.gif_gap = args.gif_gap
+        self.plot_boundary = args.plot_boundary
+        self.prob = args.prob
+        self.image_list = []
+
+    def fit(self, X):
+        """Compute fuzzy C-means clustering.
+
+        Parameters
+        ----------
+        X : array-like, shape = [n_samples, n_features]
+            Training instances to cluster.
+        """
+        self.n_samples = X.shape[0]
+        r = np.random.RandomState(self.random_state)
+        self.u = r.rand(self.n_samples,self.n_clusters)
+        # u matrix initialization: random 
+        self.u = self.u / np.tile(self.u.sum(axis=1)[np.newaxis].T, self.n_clusters)
+
+        for iteration in range(self.max_iter):
+            u_old = self.u.copy()
+            self.centers = self.next_centers(X)
+            self.u = self._predict(X)
+            if norm(self.u - u_old) < self.error: # Stopping rule
+               break
+        
+            # save intermediate result
+            if self.gif:
+               if iteration % self.gif_gap == 0: 
+                  image, fig= self.get_gif(X, iteration, plot_boundary=self.plot_boundary)
+                  fig.clf()
+                  self.image_list = self.image_list + [image]
+
+        if self.gif:
+            imageio.mimsave('./gif/FCM_m={}.gif'.format(self.m), self.image_list, fps=1)
+        print('FCM break at iteration:{}'.format(iteration))
+
+    def next_centers(self, X):
+        """Update cluster centers"""
+        um = self.u ** self.m 
+        return (X.T @ um / np.sum(um, axis=0)).T 
+    
+    ### U_update
+    def _predict(self, X):
+        """ 
+        Parameters
+        ----------
+        X : array, shape = [n_samples, n_features]
+            New data to predict.
+
+        Returns
+        -------
+        u: array, shape = [n_samples, n_clusters]
+            Fuzzy partition array, returned as an array with n_samples rows
+            and n_clusters columns.
+
+        """
+        power = float(2 / (self.m - 1))
+        temp = cdist(X, self.centers) ** power 
+        denominator_ = temp.reshape((X.shape[0], 1, -1)).repeat(temp.shape[-1], axis=1)
+        denominator_ = temp[:, :, np.newaxis] / denominator_
+        res = 1 / denominator_.sum(2)
+        return res
+
+    def predict(self, X):
+        """Predict the closest cluster each sample in X belongs to.
+        
+        Parameters
+        ----------
+        X : array, shape = [n_samples, n_features]
+            New data to predict.
+
+        Returns
+        -------
+        labels : array, shape = [n_samples,]
+            Index of the cluster each sample belongs to.
+
+        """
+        if len(X.shape) == 1:
+            X = np.expand_dims(X, axis=0)
+        u = self._predict(X)
+        return np.argmax(u, axis=-1)
+    
+
+    def predict_new(self, X):
+        if len(X.shape) == 1:
+            X = np.expand_dims(X, axis=0)
+        u = self._predict(X)
+        return u
+
+    def get_gif(self, X, iter, plot_boundary=False):
+        fig, ax = plt.subplots()
+        colors_tmp = ['g', 'r', 'b', 'c', 'm', 'k','y','cyan','pink','gray']
+        point_size = 3
+        alpha = 0.8
+        color_tmp = colors_tmp[:self.n_clusters]
+        pred = self.predict(X)
+
+        if  plot_boundary:
+            rtol = 1e-2
+            mesh_grid_pred = get_meshgrid(X)
+            mesh_grid_U_pred = self.predict_new(mesh_grid_pred)
+            cls_ind_list = []
+            for ind in range(self.n_clusters):
+                cls_ind_tmp = get_boundary_index(mesh_grid_U_pred, ind, self.prob, rtol=rtol)
+                cls_ind_list = cls_ind_list + [cls_ind_tmp]
+        
+        ax.plot(self.centers[:,0],  self.centers[:,1], "*" , markersize = 10, c = 'blue') # DA centers
+
+        for ind in range(self.n_clusters):
+            ax.scatter(X[pred==ind,0], X[pred==ind,1], s=5, alpha=0.5, c=colors_tmp[ind]) # 数据点
+            if plot_boundary:
+               ax.scatter(mesh_grid_pred[cls_ind_list[ind],0],mesh_grid_pred[cls_ind_list[ind],1], c = colors_tmp[ind], s=point_size, alpha=alpha)
+
+        ax.set_title('iter={}, m={}'.format(iter, self.m))
+        fig.canvas.draw()      
+        image = np.frombuffer(fig.canvas.tostring_rgb(), dtype='uint8')
+        image  = image.reshape(fig.canvas.get_width_height()[::-1] + (3,))
+        return image, fig
+
+
+if __name__ == '__main__':    
+    from hyperparam import args
+    from utility_code import get_data, generate_data_ind, get_meshgrid, get_boundary_index
+    
+    # data
+    args.T2 = 0.7
+    args.C = 3
+    args.R = 1
+    args = get_data(args)
+    data, labels, true_centers, data_mean, dataDistDesInd = generate_data_ind(args)
+    
+    # model
+    model = FCM(args)
+    model.fit(data)
+    centers = model.centers
+    pred = model.predict(data)
+    print('centers\n', centers)
+    
+    if not args.gif:
+        plt.scatter(data[:,0],data[:,1], c=pred, s=10)
+        plt.scatter(centers[:,0],centers[:,1], c='k', s=30)
+        plt.show()  
+
+
+
diff --git a/hierarchical_clustering.py b/hierarchical_clustering.py
new file mode 100644
index 0000000..780a1d1
--- /dev/null
+++ b/hierarchical_clustering.py
@@ -0,0 +1,48 @@
+from sklearn import cluster
+from sklearn.neighbors import NearestCentroid
+from scipy.spatial.distance import cdist # l2: default
+import matplotlib.pyplot as plt
+
+class HC:
+    def __init__(self, C, linkage):
+        self.n_clusters = C
+        self.linkage = linkage
+        self.model = None 
+        self.centers = None
+        self.pred = None
+        self.clf = None
+
+    def fit(self,X):
+        self.model = cluster.AgglomerativeClustering(n_clusters = self.n_clusters, linkage=self.linkage)
+        self.model = self.model.fit(X)
+        self.pred = self.model.labels_ 
+        self.clf = NearestCentroid() # hierarchical clustering centers 
+        self.clf.fit(X, self.pred)
+        self.centers = self.clf.centroids_ 
+    
+    # hierarchical clustering predict on new dataset
+    def predict(self,X):
+        pred_new = self.clf.predict(X)
+        return pred_new
+
+if __name__ == '__main__':
+    from hyperparam import args
+    from utility_code import get_data, generate_data_ind
+    
+    # data
+    args.T2 = 0.1
+    args.C = 3
+    args.R = 1
+    args = get_data(args)
+    data, labels, true_centers, data_mean, dataDistDesInd = generate_data_ind(args)
+    
+    # model
+    model = HC(C=args.C, linkage=args.linkage)
+    model.fit(data)
+    centers = model.centers
+    pred = model.predict(data)
+    print('centers\n', centers)
+    
+    plt.scatter(data[:,0],data[:,1], c=pred, s=10)
+    plt.scatter(centers[:,0],centers[:,1], c='k', s=30)
+    plt.show()  
diff --git a/hyperparam.py b/hyperparam.py
new file mode 100644
index 0000000..53a857b
--- /dev/null
+++ b/hyperparam.py
@@ -0,0 +1,74 @@
+### hyperparameter for ITISC
+import argparse
+parser = argparse.ArgumentParser(description='ITISC algorithms')
+parser.add_argument('--T1',default = 1.0, type=float, help='parameter T1')
+parser.add_argument('--T2',default = 0.5, type=float, help='parameter T2')
+parser.add_argument('--m',default = 2.0, type=float, help='parameter m for FCM algorithm')
+
+parser.add_argument('--N1',default = 200, type=int, help='# of points in cluster 1')
+parser.add_argument('--N2',default = 200, type=int, help='# of points in cluster 2')
+parser.add_argument('--N3',default = 200, type=int, help='# of points in cluster 3')
+parser.add_argument('--N4',default = 200, type=int, help='# of points in cluster 4')
+parser.add_argument('--N5',default = 200, type=int, help='# of points in cluster 5')
+parser.add_argument('--N6',default = 200, type=int, help='# of points in cluster 6')
+parser.add_argument('--N',default = 50, type=int, help='# uniform random noise')
+
+parser.add_argument('--max_iter',default = 150, type=int, help='maximum number of iterations')
+parser.add_argument('--random_state',default = 0, type=int, help='random seed') # data generation
+parser.add_argument('--itisc_random_state',default = 0, type=int, help='random seed for itisc initialization')   # itisc random seed
+parser.add_argument('--fcm_random_state',default = 0, type=int, help='random seed for fcm initialization')       # fcm random state
+parser.add_argument('--kmeans_random_state',default = 0, type=int, help='random seed for kmeans initialization') # kmeans random state
+parser.add_argument('--error',default = 1e-5, type=float, help='convergence error')
+
+### ITISC initialization
+parser.add_argument('--init_method', default='random', type=str, 
+                    help='method to initialize clustering centers(default:random)',
+                    choices=['kmeans','fcm', 'random','given','ITISC-R'])
+
+# argment for generating synthetic 2D Gaussian generation
+parser.add_argument('--C', default = 3, type=int, help='modeled cluster numbers') # required=True
+parser.add_argument('--mean1', nargs='+', default=[1,0], type=float, help='true mean value for the 1-th cluster')
+parser.add_argument('--mean2', nargs='+', default=[2,0], type=float, help='true mean value for the 2-th cluster')
+parser.add_argument('--mean3', nargs='+', default=[3,0], type=float, help='true mean value for the 3-th cluster')
+parser.add_argument('--mean4', nargs='+', default=[4,0], type=float, help='true mean value for the 4-th cluster')
+parser.add_argument('--mean5', nargs='+', default=[5,0], type=float, help='true mean value for the 5-th cluster')
+parser.add_argument('--mean6', nargs='+', default=[6,0], type=float, help='true mean value for the 6-th cluster')
+parser.add_argument('--cov1', nargs='+', default=[0.2,0,0,0.2], type=float, help='1-th covariance matrix')
+parser.add_argument('--cov2', nargs='+', default=[0.2,0,0,0.2], type=float, help='2-th covariance matrix')
+parser.add_argument('--cov3', nargs='+', default=[0.2,0,0,0.2], type=float, help='3-th covariance matrix')
+parser.add_argument('--cov4', nargs='+', default=[0.2,0,0,0.2], type=float, help='4-th covariance matrix')
+parser.add_argument('--cov5', nargs='+', default=[0.2,0,0,0.2], type=float, help='5-th covariance matrix')
+parser.add_argument('--cov6', nargs='+', default=[0.2,0,0,0.2], type=float, help='6-th covariance matrix')
+parser.add_argument('--R', default = 1.0, type=float, help='Gaussian synthetic dataset generation hyperparameter')
+# M-boundaryDist
+parser.add_argument('--max_num', default = 10, type=int, help='Number of marginal points')
+# analysis 
+parser.add_argument('--model', default='itisc', type=str, help='clustering method', choices=['kmeans','fcm', 'itisc','dafcm','isfcm'])
+parser.add_argument('--save_plot',action='store_true',help='whether to save plots')
+# plot
+parser.add_argument('--prob',default = 1/3, type=float, help='probability for isoprobability curve')
+parser.add_argument('--plot',action='store_true',help='test code, plot')
+parser.add_argument('--plot_boundary',action='store_true',help='whether to isoprobability curve')
+parser.add_argument('--gif',action='store_true',help='whether to draw gif in ITISC')
+parser.add_argument('--gif_gap', default = 1, type=int, help='record interval between gifs')
+### hierarchical clustering 
+parser.add_argument('--linkage', default='ward', type=str, help='linkage type for hierarchical clustering')
+
+# python arg.py -l 1234 2345 3456 4567
+# check python or jupyter
+def is_notebook() -> bool:
+    try:
+        shell = get_ipython().__class__.__name__
+        if shell == 'ZMQInteractiveShell':
+            return True   # Jupyter notebook or qtconsole
+        elif shell == 'TerminalInteractiveShell':
+            return False  # Terminal running IPython
+        else:
+            return False  # Other type (?)
+    except NameError:
+        return False
+
+if is_notebook():
+   args = parser.parse_args(args=[]) # jupyter
+else:
+   args = parser.parse_args() # python
diff --git a/utility_code.py b/utility_code.py
new file mode 100644
index 0000000..93ea413
--- /dev/null
+++ b/utility_code.py
@@ -0,0 +1,294 @@
+### utility function for ITISC
+import matplotlib
+import matplotlib.pyplot as plt
+import pickle
+import os
+import csv
+import numpy as np
+import random
+from math import pi, cos, sin
+from scipy.spatial.distance import cdist # l2: default
+from math import pi, cos, sin
+from hyperparam import args
+
+def get_color_list():
+    c1 = np.array([1,86,153])/255
+    c2 = np.array([250,192,15])/255
+    c3 = np.array([243,118,74])/255
+    c4 = np.array([95,198,201])/255
+    c5 = np.array([79,89,100])/255
+    c6 = np.array([176,85,42])/255
+    color_list = [c1,c2,c3,c4,c5,c6]
+    return color_list
+
+def get_color_tmp():
+    colors_tmp = ['g', 'r', 'b', 'c', 'm', 'k','y','cyan','pink','gray']
+    return colors_tmp
+
+def save_obj(obj, name, dir):
+    save_path = os.path.join(dir, name + '.pkl')
+    with open(save_path, 'wb') as f:
+        pickle.dump(obj, f, pickle.HIGHEST_PROTOCOL)
+
+def load_obj(name, dir):
+    save_path = os.path.join(dir, name + '.pkl')
+    with open(save_path, 'rb') as f:
+        return pickle.load(f)
+
+def get_kmeans(args, data):
+    from sklearn.cluster import KMeans
+    kmeans = KMeans(n_clusters=args.C, random_state=args.kmeans_random_state).fit(data) # kmeans++ initialization
+    kmeans_centers = kmeans.cluster_centers_
+    kmeans_pred = kmeans.predict(data) # data:[N,R]
+    print('Kmeans-centers\n', kmeans_centers)
+    return kmeans, kmeans_centers, kmeans_pred
+
+def get_fcm(args,data):
+    from fcm_github import FCM
+    fcm = FCM(args)
+    fcm.fit(data)
+    fcm_centers = fcm.centers
+    fcm_pred = fcm.predict(data)
+    print('fcm_centers\n', fcm_centers)
+    return fcm, fcm_centers, fcm_pred
+
+def get_hc(args, data):
+    from hierarchical_clustering import HC
+    hc = HC(C=args.C, linkage=args.linkage)
+    hc.fit(data)
+    hc_centers = hc.centers
+    hc_pred = hc.pred
+    print('hc_centers\n', hc_centers)
+    return hc, hc_centers, hc_pred
+
+### ITISC_AO
+def get_itisc_ao(args, data):
+    from ITISC_AO import ITISC_AO
+    itisc = ITISC_AO(args)
+    _ = itisc.fit(data)
+    itisc_centers = itisc.centers
+    itisc_pred = itisc.predict(data)
+    print('itisc_ao_centers', itisc_centers)
+    return itisc, itisc_centers, itisc_pred
+
+### ITISC-R: matlab: log d(x,y)
+def get_itisc_r(args, data):
+    from ITISC_R import ITISC_R
+    itisc = ITISC_R(args)
+    _ = itisc.fit(data)
+    itisc_centers = itisc.centers
+    itisc_pred = itisc.predict(data)
+    print('itisc_r_centers', itisc_centers)
+    return itisc, itisc_centers, itisc_pred
+
+
+### KL divergence between two Gaussian
+def get_KL_two_multivariate_gaussian(mu_1, mu_2, sigma_1, sigma_2):
+    n = sigma_1.shape[0]
+    det_sigma_1 = np.linalg.det(sigma_1)
+    det_sigma_2 = np.linalg.det(sigma_2)
+    sigma_2_inv = np.linalg.inv(sigma_2) # inverse matrix of sigma_2
+    result = np.log(det_sigma_2/det_sigma_1) - n + np.matrix.trace(sigma_2_inv @ sigma_1) + \
+    np.transpose(mu_2 - mu_1) @ sigma_2_inv @ (mu_2 - mu_1)
+    return 0.5 * result
+
+### mesh grid for isoprobability curve
+def get_meshgrid(data, h=0.025):
+    # h: meshgrid gap 
+    x_min, x_max = data[:, 0].min() - 1, data[:, 0].max() + 1
+    y_min, y_max = data[:, 1].min() - 1, data[:, 1].max() + 1
+    xx, yy = np.meshgrid(np.arange(x_min, x_max, h), np.arange(y_min, y_max, h))
+    data_pred = np.c_[xx.ravel(), yy.ravel()] # Obtain labels for each point in mesh. Use last trained model.
+    return data_pred
+
+def get_boundary_index(U_pred, cls_ind, prob, rtol=1e-2):
+    # U_pred:  [N_new, C]
+    # cls_ind: cluster center index
+    # prob:    probability for isoprobability curve
+    # rtol:    allowed error 1e-2
+    cls_num = np.isclose(U_pred[:,cls_ind], prob, rtol)
+    cls_num_ind  = np.argwhere(cls_num == 1)
+    return cls_num_ind
+
+
+def get_rotate_matrix(cov,rotate_degree):  
+    cov = np.array(cov).reshape(2,2)
+    theta = pi * rotate_degree
+    rotate_matrix = [cos(theta), -sin(theta), sin(theta), cos(theta)]   # clockwise rotate
+    # rotate_matrix = [cos(theta), sin(theta), -sin(theta), cos(theta)] # counter clockwise rotate
+    rotate_matrix = np.array(rotate_matrix).reshape(2,2)
+    tmp = rotate_matrix @ cov @ rotate_matrix.T
+    tmp = tmp.flatten()
+    return tmp
+
+# extreme dataset: C=3
+def get_extreme_3_clusters(args):
+    args.mean1 = [1,0]
+    args.mean2 = [8,0]
+    args.mean3 = [4,8]
+    args.cov1 = [0.8,0.4,0.4,0.8]
+    args.cov2 = [0.8,0.4,0.4,0.8]
+    args.cov3 = [0.8,-0.4,-0.4,0.8]
+    args.N1 = 2
+    args.N2 = 100
+    args.N3 = 2
+    return args 
+
+# extreme dataset: C=2
+def get_extreme_2_clusters(args):
+    args.C = 2
+    args.mean1 = [1,0]
+    args.mean2 = [5,0]
+    args.cov1 = [1,0.5,0.5,1]
+    args.cov2 = [1,-0.5,-0.5,1]
+    return args
+
+    
+def get_data(args):
+    # args.R = 1.5
+    N = 200
+    args.N1 = N
+    args.N2 = N
+    args.N3 = N
+    args.N4 = N
+    args.N5 = N
+    args.N6 = N
+    ### C = 2
+    if args.C == 2:
+        args.mean1 = [args.R,0]
+        args.mean2 = [-args.R,0]
+        cov = [1,0,0,0.3]
+        args.cov1 = get_rotate_matrix(cov,0.25)
+        args.cov2 = get_rotate_matrix(cov,0.75)   
+    ### C=3
+    elif args.C == 3:
+        args.mean1 = [args.R,0]
+        args.mean2 = [-1 * args.R / np.sqrt(3),-1*args.R]
+        args.mean3 = [-1 * args.R / np.sqrt(3),args.R]
+        args.cov1 = [1,0,0,0.3]
+        args.cov2 = get_rotate_matrix(args.cov1,1/3) # cov2 
+        args.cov3 = get_rotate_matrix(args.cov1,2/3) # cov3
+    ### C=4
+    elif args.C == 4:
+        args.mean1 = [args.R,args.R]
+        args.mean2 = [args.R,-args.R]
+        args.mean3 = [-args.R,-args.R]
+        args.mean4 = [-args.R,args.R]
+        cov = [1,0,0,0.1]
+        rotate_degree = 0.5 
+        args.cov1 = get_rotate_matrix(cov,0.25) # cov2
+        args.cov2 = get_rotate_matrix(cov,0.75) # cov2
+        args.cov3 = get_rotate_matrix(cov,1.25) # cov3
+        args.cov4 = get_rotate_matrix(cov,1.75) # cov4
+    ### C=6
+    elif args.C == 6:
+        tmp = np.sqrt(3)
+        args.mean1 = [args.R/2,args.R/2*tmp] 
+        args.mean2 = [-args.R/2,args.R/2*tmp] 
+        args.mean3 = [-args.R,0] 
+        args.mean4 = [-args.R/2,-args.R/2*tmp] 
+        args.mean5 = [args.R/2,-args.R/2*tmp] 
+        args.mean6 = [args.R,0]
+        cov = [1,0,0,0.3]
+        rotate_degree = 1/3 
+        args.cov1 = get_rotate_matrix(cov,rotate_degree) # cov2
+        args.cov2 = get_rotate_matrix(cov,rotate_degree*2) 
+        args.cov3 = get_rotate_matrix(cov,rotate_degree*3) 
+        args.cov4 = get_rotate_matrix(cov,rotate_degree*4) 
+        args.cov5 = get_rotate_matrix(cov,rotate_degree*5) 
+        args.cov6 = get_rotate_matrix(cov,rotate_degree*6)
+    ### 
+    print('mean1', args.mean1)
+    print('mean2', args.mean2)
+    print('mean3', args.mean3)
+    print('mean4', args.mean4)
+    print('mean5', args.mean5)
+    print('mean6', args.mean6)
+
+    print('cov1', args.cov1)
+    print('cov2', args.cov2)
+    print('cov3', args.cov3)
+    print('cov4', args.cov4)
+    print('cov5', args.cov5)
+    print('cov6', args.cov6)
+    return args
+
+
+### 2D data generation
+def generate_data(args):
+    rs = np.random.RandomState(args.random_state)
+    mean_list = []
+    cov_list = []
+    N_list = []
+
+    for i in range(args.true_C):
+        mean_list = mean_list + [eval('args.mean{}'.format(i+1))]
+        cov_list = cov_list + [eval('args.cov{}'.format(i+1))]
+        N_list = N_list + [eval('args.N{}'.format(i+1))]
+
+    labels = []
+    true_centers = []
+
+    for i in range(args.true_C):
+        mean_tmp = mean_list[i]
+        cov_tmp = cov_list[i]
+        N_tmp = N_list[i]
+
+        cov_tmp = np.array(cov_tmp).reshape(2,2)
+        x_tmp, y_tmp = rs.multivariate_normal(mean_tmp, cov_tmp, N_tmp).T
+        if i == 0:   
+           x = x_tmp 
+           y = y_tmp
+        else:
+           x = np.concatenate((x, x_tmp), axis=0)
+           y = np.concatenate((y, y_tmp), axis=0)
+
+        labels = labels + [i] * N_tmp
+        true_centers = true_centers + [mean_tmp]
+
+    # random points
+    labels = np.array(labels)
+    x = x[:,np.newaxis]
+    y = y[:,np.newaxis]
+    data = np.concatenate((x,y),1)
+    return data, labels, np.array(true_centers)
+
+
+def generate_data_ind(args):
+    rs = np.random.RandomState(args.random_state)
+    mean_list = []
+    cov_list = []
+    N_list = []
+    for i in range(args.C):
+        mean_list = mean_list + [eval('args.mean{}'.format(i+1))]
+        cov_list = cov_list + [eval('args.cov{}'.format(i+1))]
+        N_list = N_list + [eval('args.N{}'.format(i+1))]
+    labels = []
+    true_centers = []
+    for i in range(args.C):
+        mean_tmp = mean_list[i]
+        cov_tmp = cov_list[i]
+        N_tmp = N_list[i]
+        cov_tmp = np.array(cov_tmp).reshape(2,2)
+        x_tmp, y_tmp = rs.multivariate_normal(mean_tmp, cov_tmp, N_tmp).T
+        if i == 0: 
+           x = x_tmp 
+           y = y_tmp
+        else:
+           x = np.concatenate((x, x_tmp), axis=0)
+           y = np.concatenate((y, y_tmp), axis=0)
+        labels = labels + [i] * N_tmp
+        true_centers = true_centers + [mean_tmp]
+    # random points
+    labels = np.array(labels)
+    x = x[:,np.newaxis]
+    y = y[:,np.newaxis]
+    data = np.concatenate((x,y),1)
+    ### 
+    data_mean = np.mean(data, axis=0, keepdims=True)
+    print('data_mean', data_mean) # data centroid
+    dataMeanDist = cdist(data, data_mean)**2 # [N,1] 
+    dataDistDesInd = dataMeanDist.squeeze().argsort()[::-1] # Data Distance Descending Index
+    return data, labels, np.array(true_centers), data_mean, dataDistDesInd
+
+