def contrast(spec, sr, a):
	""" 
	contrast() - Compute spectral contrast

		Inputs: 
			spec - spectrogram
			sr - sample rate
			a - percentage of frequencies to consider for peak and valley

		Outputs:
			contrasts - array of centroids for each frame

	"""

	nbins, nframes = spec.shape

	fftSize = (nbins - 1) * 2
	freqRes = sr / 2.0 / fftSize

	contrasts = np.zeros(nframes)
	for i in range(nframes):
		sortedFreqs = np.sort(spec[:,i])
		nbinsToLook = np.round(nbins * a).astype("int")

		valley = np.log(np.sum(sortedFreqs[:nbinsToLook]) + np.spacing(1)) / nbinsToLook
		peak = np.log(np.sum(sortedFreqs[nbins-nbinsToLook:nbins]) + np.spacing(1)) / nbinsToLook

		contrasts[i] = peak - valley

	return contrasts

def add_qualifying_constraint(m,  coefficients, M,  K,  N,  thetaB,  g_flag,  t):
    ''' Add the qualifying constraints to model m.  The qualifying 
        constraint is formed by linearizing the Lagrangian with respect 
        to x about x0. Then linearizing that with respect to theta 
        about thetat.  '''

    qualifying_constraint = []
    x =  [[m.getVarByName("x_%d_%d" % (j,k)) for k in xrange(K)] for j in xrange(M)]

    
    for i in xrange(len(coefficients)):
        #calcuate the qualifying constraint formula
        g_expr = gb.LinExpr()
        g_expr.addConstant(coefficients[i][-1])
        for j in xrange(M*K):
            g_expr.add(x[j/K][j%K]* coefficients[i][j])
        qualifying_constraint.append(g_expr)
        #add the qualifying constraints
        if g_flag[i%K][i/K] != 0:
            ##################### Add constraints: g_expr #########################
            if thetaB[i%K][i/K] == 1:
                m.addConstr(g_expr <= np.spacing(1), name="qc%d_%d_%d" % (t,k,i))
            elif thetaB[i%K][i/K] == 0:
                m.addConstr(g_expr >= -np.spacing(1), name="qc%d_%d_%d" % (t,k,i))
            m.update() 
            

    #return qualifying constraints to calculate lagrange constraint
    return qualifying_constraint

def renyi_information(tfr, timestamps=None, freq=None, alpha=3.0):
    """renyi_information

    :param tfr:
    :param timestamps:
    :param freq:
    :param alpha:
    :type tfr:
    :type timestamps:
    :type freq:
    :type alpha:
    :return:
    :rtype:
    """
    if alpha == 1 and tfr.min().min() < 0:
        raise ValueError("Distribution with negative values not allowed.")
    m, n = tfr.shape
    if timestamps is None:
        timestamps = np.arange(n) + 1
    elif np.allclose(timestamps, np.arange(n)):
        timestamps += 1
    if freq is None:
        freq = np.arange(m)
    freq.sort()
    tfr = tfr / integrate_2d(tfr, timestamps, freq)
    if alpha == 1:
        R = -integrate_2d(tfr * np.log2(tfr + np.spacing(1)), timestamps, freq)
    else:
        R = np.log2(integrate_2d(tfr ** alpha, timestamps, freq) + np.spacing(1))
        R = R / (1 - alpha)
    return R

def colon(r1, inc, r2):
    """
      Matlab's colon operator, althought it doesn't although inc is required

    """

    s = np.sign(inc)

    if s == 0:
        return_value = np.zeros(1)
    elif s == 1:
        n = ((r2 - r1) + 2 * np.spacing(r2 - r1)) // inc
        return_value = np.linspace(r1, r1 + inc * n, n + 1)
    else:  # s == -1:
        # NOTE: I think this is slightly off as we start on the wrong end
        # r1 should be exact, not r2
        n = ((r1 - r2) + 2 * np.spacing(r1 - r2)) // np.abs(inc)
        temp = np.linspace(r2, r2 + np.abs(inc) * n, n + 1)
        return_value = temp[::-1]

    # If the start and steps are whole numbers, we should cast as int
    if(np.equal(np.mod(r1, 1), 0) and
       np.equal(np.mod(s, 1), 0) and
       np.equal(np.mod(r2, 1), 0)):
        return return_value.astype(int)
    else:
        return return_value

 def generateMuonPlusJet(self,nParticles=1,muonEnergy=10):
     particles=["gamma","pi-","pi+","proton","pi0","e+","e-","kaon0L","kaon+","kaon-"]
     config=open("config.txt","w")
     self.theta=self.maxTheta + np.spacing(1.0)
     while(  self.theta>self.maxTheta or self.theta<0 ):
         self.theta=np.random.normal(self.centralTheta,self.sigmaTheta)
         
     eta=-math.log( math.tan(self.theta/2) )
     phi=random.uniform(0,self.maxPhi)
     config.write( "mu-" + " " + str(muonEnergy) + " " + str(eta) + " " + str(phi) )
     for i in range(0,nParticles):
         config.write("\n")
         part=particles[ random.randint(0,len(particles)-1) ]
         self.theta=self.maxTheta + np.spacing(1.0)
         while(  self.theta>self.maxTheta or self.theta<0 ):
             self.theta=np.random.normal(self.centralTheta,self.sigmaTheta)
         
         eta=-math.log( math.tan(self.theta/2) )
         phi=random.uniform(0,self.maxPhi)
         energy=random.uniform(self.minEnergy,self.maxEnergy)
         if( eta<1.5 ):
             print part + " " + str(energy) + " " + str(eta) + " " + str(phi) 
         phi=random.uniform(0,self.maxPhi)
         config.write( part + " " + str(energy) + " " + str(eta) + " " + str(phi) )
     config.close()

    def calculateEntropy(self,Y, mship):
        """
            calculates the split entropy using Y and mship (logical array) telling which 
            child the examples are being split into...

            Input:
            ---------
                Y: a label array
                mship: (logical array) telling which child the examples are being split into, whether
                        each example is assigned to left split or the right one..
            Returns:
            ---------
                entropy: split entropy of the split
        """

        lexam=Y[mship]
        rexam=Y[np.logical_not(mship)]

        pleft= len(lexam) / float(len(Y))
        pright= 1-pleft

        pl= stats.itemfreq(lexam)[:,1] / float(len(lexam)) + np.spacing(1)
        pr= stats.itemfreq(rexam)[:,1] / float(len(rexam)) + np.spacing(1)

        hl= -np.sum(pl*np.log2(pl)) 
        hr= -np.sum(pr*np.log2(pr)) 

        sentropy = pleft * hl + pright * hr

        return sentropy

    def entropy(self, c1, c2):

        total = c1 + c2
        return -1 * (
            (c1 / total) * numpy.log2((c1 / total) + numpy.spacing(1))
            + (c2 / total) * numpy.log2((c2 / total) + numpy.spacing(1))
        )

def klDiv(margDist1, margDist2):
    assert isinstance(margDist1, margDistSCPP) and\
            isinstance(margDist2, margDistSCPP),\
        "margDist1 and margDist2 must be instances of margDistSCPP"
        
    numpy.testing.assert_equal(margDist1.m, 
                               margDist2.m)
    
    kldSum = 0.0
    for idx in xrange(margDist1.m):
        # test that the distributions are over the same bin indices
        # if they are not this calculation is meaninless
        numpy.testing.assert_equal(margDist1.data[idx][1],
                                   margDist2.data[idx][1])
        
        #test that the distributions have the same number of bins
        numpy.testing.assert_equal(len(margDist1.data[idx][0]),
                                   len(margDist2.data[idx][0]))
        
        m1 = margDist1.data[idx][0]
        m1[m1==0] = numpy.spacing(1)
        m2 = margDist2.data[idx][0]
        m2[m2==0] = numpy.spacing(1)
        
             
        kldSum += numpy.sum(margDist1.data[idx][0] * (numpy.log(margDist1.data[idx][0]) 
                                               - numpy.log(margDist2.data[idx][0])))

        
    return kldSum

    def call(self, inputs, mask=None):
        if not isinstance(inputs, list) or len(inputs) <= 1:
            raise TypeError('SpkLifeLongMemory must be called on a list of tensors '
                            '(at least 2). Got: ' + str(inputs))
        # (None(batch), 1), index of speaker
        target_spk_l = inputs[0]
        target_spk_l = K.reshape(target_spk_l, (target_spk_l.shape[0], ))
        if K.dtype(target_spk_l) != 'int32':
            target_spk_l = K.cast(target_spk_l, 'int32')
        # (None(batch), embed_dim)
        spk_vector_l = inputs[1]
        # Start to update life-long memory based on the learned speech vector
        # First do normalization
        spk_vector_eps = K.switch(K.equal(spk_vector_l, 0.), np.spacing(1), spk_vector_l)  # avoid zero
        spk_vector_eps = K.sqrt(K.sum(spk_vector_eps**2, axis=1))
        spk_vector_eps = spk_vector_eps.dimshuffle((0, 'x'))
        spk_vector = T.true_div(spk_vector_l, K.repeat_elements(spk_vector_eps, self.vec_dim, axis=1))
        # Store speech vector into life-long memory according to the speaker identity.
        life_long_mem = T.inc_subtensor(self.life_long_mem[target_spk_l, :], spk_vector)
        # Normalization for memory
        life_long_mem_eps = K.switch(K.equal(life_long_mem, 0.), np.spacing(1), life_long_mem)  # avoid 0
        life_long_mem_eps = K.sqrt(K.sum(life_long_mem_eps**2, axis=1))
        life_long_mem_eps = life_long_mem_eps.dimshuffle((0, 'x'))
        life_long_mem = T.true_div(life_long_mem, K.repeat_elements(life_long_mem_eps, self.vec_dim, axis=1))

        # (None(batch), spk_size, embed_dim)
        return life_long_mem

def solve_master(tree, num_node, Current_node, g_flag, thetaB_list, SUBD,  coefficients, xBar, thetaOpt, lamOpt, muOpt, y,  cons, iteration,  pool=None):
    '''we solve the relaxed master problems based on thetaB_list, then select the infimum of all minimum values.
       Parameters: About the tree : tree, Current_node
                   About the subproblem: SUBD, xBar, thetaOpt, lamOpt, muOpt, y
                   About the boundary: theta_L, theta_U''' 
                   
    (M, N) = np.shape(y)
    K = np.shape(xBar[-1])[0]
    
    x_stor = None
    Q_stor = np.inf
    next_node = -1
    
    #store all the MLBD
    MLBD_stor = [] 

    #store all the MLBD
    if pool == None:
        tree.nodes[Current_node].set_parameters_qualifying_constraint(lamOpt,  thetaOpt, muOpt,  xBar, SUBD,  g_flag,  coefficients)
        #check whether the coefficients are already stored into the parents or not.
        
        print ('\n%d master problems are solving...'  %len(thetaB_list))

        for index in xrange(len(thetaB_list)):
            thetaB = thetaB_list[index].copy()
            status, objVal, xOpt,  thetaB, lagrangian_coefficient= solve_master_s(tree, Current_node, coefficients, thetaOpt, xBar, lamOpt, muOpt, thetaB.copy(), y, g_flag, cons)
            #print objVal, xOpt
            
            if status == 2 and objVal < SUBD - np.spacing(1):
                node = tree.add_node(num_node, 0, 1, Current_node)
                node.set_parameters_thetaB(thetaB,  xOpt, objVal, lagrangian_coefficient)
                MLBD_stor.append(objVal)
                if objVal < Q_stor:
                    Q_stor = objVal
                    next_node = num_node
                    x_stor = xOpt
                num_node = num_node + 1

    else:
        tree.nodes[Current_node].set_parameters_qualifying_constraint(lamOpt,  thetaOpt, muOpt,  xBar, SUBD,  g_flag,  coefficients)
        len_thetaB = len(thetaB_list)
        print ('\n%d master problems are solving...'  %len_thetaB)
        results = [pool.apply_async(solve_master_s,  args = (tree, Current_node, coefficients, thetaOpt, xBar, lamOpt, muOpt, thetaB.copy(), y, g_flag, cons)) for thetaB in thetaB_list]

        #put all the result into the tree.
        for p in results:
            #result = [status, objVal, xOpt, thetaB]
            result = p.get() 
            if result[0] == 2 and result[1] < SUBD - np.spacing(1):
                node = tree.add_node(num_node,  0,  1,  Current_node)
                node.set_parameters_thetaB(result[3],  result[2],  result[1], result[4])
                #node.set_parameter(lamOpt,  thetaOpt, result[3],  muOpt,  xBar,  result[2],  SUBD,  result[1], g_flag,  coefficients)
                MLBD_stor.append(result[1])
                if result[1] < Q_stor:
                    Q_stor = result[1]
                    next_node = num_node
                    x_stor =  result[2]
                num_node += 1

    return x_stor, Q_stor, next_node, num_node, MLBD_stor

def _frank(M, N, alpha):
    if(N<2):
        raise ValueError('Dimensionality Argument [N] must be an integer >= 2')
    elif(N==2):        
        u1 = uniform.rvs(size=M)
        p = uniform.rvs(size=M)
        if abs(alpha) > math.log(sys.float_info.max):
            u2 = (u1 < 0).astype(int) + np.sign(alpha)*u1  # u1 or 1-u1
        elif abs(alpha) > math.sqrt(np.spacing(1)):
            u2 = -1*np.log((np.exp(-alpha*u1)*(1-p)/p + np.exp(-alpha))/(1 + np.exp(-alpha*u1)*(1-p)/p))/alpha
        else:
            u2 = p
        
        U = np.column_stack((u1,u2))
    else:
        # Algorithm 1 described in both the SAS Copula Procedure, as well as the
        # paper: "High Dimensional Archimedean Copula Generation Algorithm"
        if(alpha<=0):
            raise ValueError('For N>=3, alpha >0 in Frank Copula')
            
        U = np.empty((M,N))
        for ii in range(0,M):
            p = -1.0*np.expm1(-1*alpha)
            if(p==1):
                # boundary case protection
                p = 1 - np.spacing(1)
            v = logser.rvs(p, size=1)
            
            # sample N independent uniform random variables
            x_i = uniform.rvs(size=N)
            t = -1*np.log(x_i)/v
            U[ii,:] = -1.0*np.log1p( np.exp(-t)*np.expm1(-1.0*alpha))/alpha
            
    return U

    def test_half_fpe(self):
        oldsettings = np.seterr(all="raise")
        try:
            sx16 = np.array((1e-4,), dtype=float16)
            bx16 = np.array((1e4,), dtype=float16)
            sy16 = float16(1e-4)
            by16 = float16(1e4)

            # Underflow errors
            assert_raises_fpe("underflow", lambda a, b: a * b, sx16, sx16)
            assert_raises_fpe("underflow", lambda a, b: a * b, sx16, sy16)
            assert_raises_fpe("underflow", lambda a, b: a * b, sy16, sx16)
            assert_raises_fpe("underflow", lambda a, b: a * b, sy16, sy16)
            assert_raises_fpe("underflow", lambda a, b: a / b, sx16, bx16)
            assert_raises_fpe("underflow", lambda a, b: a / b, sx16, by16)
            assert_raises_fpe("underflow", lambda a, b: a / b, sy16, bx16)
            assert_raises_fpe("underflow", lambda a, b: a / b, sy16, by16)
            assert_raises_fpe("underflow", lambda a, b: a / b, float16(2.0 ** -14), float16(2 ** 11))
            assert_raises_fpe("underflow", lambda a, b: a / b, float16(-2.0 ** -14), float16(2 ** 11))
            assert_raises_fpe("underflow", lambda a, b: a / b, float16(2.0 ** -14 + 2 ** -24), float16(2))
            assert_raises_fpe("underflow", lambda a, b: a / b, float16(-2.0 ** -14 - 2 ** -24), float16(2))
            assert_raises_fpe("underflow", lambda a, b: a / b, float16(2.0 ** -14 + 2 ** -23), float16(4))

            # Overflow errors
            assert_raises_fpe("overflow", lambda a, b: a * b, bx16, bx16)
            assert_raises_fpe("overflow", lambda a, b: a * b, bx16, by16)
            assert_raises_fpe("overflow", lambda a, b: a * b, by16, bx16)
            assert_raises_fpe("overflow", lambda a, b: a * b, by16, by16)
            assert_raises_fpe("overflow", lambda a, b: a / b, bx16, sx16)
            assert_raises_fpe("overflow", lambda a, b: a / b, bx16, sy16)
            assert_raises_fpe("overflow", lambda a, b: a / b, by16, sx16)
            assert_raises_fpe("overflow", lambda a, b: a / b, by16, sy16)
            assert_raises_fpe("overflow", lambda a, b: a + b, float16(65504), float16(17))
            assert_raises_fpe("overflow", lambda a, b: a - b, float16(-65504), float16(17))
            assert_raises_fpe("overflow", np.nextafter, float16(65504), float16(np.inf))
            assert_raises_fpe("overflow", np.nextafter, float16(-65504), float16(-np.inf))
            assert_raises_fpe("overflow", np.spacing, float16(65504))

            # Invalid value errors
            assert_raises_fpe("invalid", np.divide, float16(np.inf), float16(np.inf))
            assert_raises_fpe("invalid", np.spacing, float16(np.inf))
            assert_raises_fpe("invalid", np.spacing, float16(np.nan))
            assert_raises_fpe("invalid", np.nextafter, float16(np.inf), float16(0))
            assert_raises_fpe("invalid", np.nextafter, float16(-np.inf), float16(0))
            assert_raises_fpe("invalid", np.nextafter, float16(0), float16(np.nan))

            # These should not raise
            float16(65472) + float16(32)
            float16(2 ** -13) / float16(2)
            float16(2 ** -14) / float16(2 ** 10)
            np.spacing(float16(-65504))
            np.nextafter(float16(65504), float16(-np.inf))
            np.nextafter(float16(-65504), float16(np.inf))
            float16(2 ** -14) / float16(2 ** 10)
            float16(-2 ** -14) / float16(2 ** 10)
            float16(2 ** -14 + 2 ** -23) / float16(2)
            float16(-2 ** -14 - 2 ** -23) / float16(2)
        finally:
            np.seterr(**oldsettings)

def interior_point_using_distance(coefficients, nCuts,  dim,  marker, weight, index = -1):
    ''' Get the interior point from the certain region described by marker. The returned point should be far from all the hyperplanes. 
    coefficients: the coefficients of hyperplanes with the same order. For example, (a1, a2, a3, c) in the hyperplane (a1 * x1 + a2 * x2 + a3 * x3 + c = 0)
    nCuts: the number of the hyerplanes.
    Dim: the dimension of the space.
    marker: the sign of the hyperplane.
    weight: calcualte the norm 2 of the coefficients except the constant of all hyperplanes  
    index: the hyperplane we want to trim. Usually index = -1, which means no hyperplane will be trimed.
    
    Here, we add the different weights for different hyperplanes when we calculate the distance.

    maximize z
     subject to  A*x + b + weight*z <= np.spacing(1), for all sign >=0
                       -1*(A*x + b) + weight*z <= np.spacing(1), for all sign <= 0
    '''
    
    #Create a new model
    m = gb.Model('Interior_point')
    
    #Set parameters
    m.setParam('OutputFlag',  False)
    
    x = [0 for i in xrange(dim)]
    for i in xrange(dim):
        x[i] = m.addVar(lb = -1e7,  ub = 1e7,  vtype=gb.GRB.CONTINUOUS,  name = 'x_%d'%(i))
    
    obj = m.addVar(lb = 0.0,  vtype = gb.GRB.CONTINUOUS,  name = 'objective')
    m.update()
    
    for i in xrange(nCuts):
        #hyperplane index is trimed, so there is no constraint for index hyperplane.
        if index != i:
            g_expr = gb.LinExpr()
            g_expr.addConstant(coefficients[i][-1])
            for j in xrange(dim):
                g_expr.add(x[j] * coefficients[i][j])
            if marker[i] == 0:
                m.addConstr(-1 * g_expr + weight[i]*obj <= np.spacing(0),  name = 'qc_%d' %(i))
            elif marker[i] == 1:
                m.addConstr(g_expr + weight[i]*obj <= np.spacing(0),  name = 'qc_%d' %(i))
    m.update()
    
    #Create the objective : maximize obj
    m.setObjective(obj,  gb.GRB.MAXIMIZE)
    m.update()
    
    #Optimize the test problem.
    try:
        m.optimize()
    except gb.GurobiError as e:
        print e.message

    if m.Status == gb.GRB.OPTIMAL:
        xOpt = np.empty(dim)
        for i in xrange(dim):
            xOpt[i] = x[i].x
        return m.Status,  xOpt,  obj.x
    else:
        return m.Status,  np.nan,  np.nan

def ellipse_axis_length( a ):
    b,c,d,f,g,a = a[1]/2, a[2], a[3]/2, a[4]/2, a[5], a[0]
    up = 2*(a*f*f+c*d*d+g*b*b-2*b*d*f-a*c*g)
    down1=(b*b-a*c)*( (c-a)*np.sqrt(1+4*b*b/((a-c)*(a-c)))-(c+a))
    down2=(b*b-a*c)*( (a-c)*np.sqrt(1+4*b*b/((a-c)*(a-c)))-(c+a))
    res1=np.sqrt(up/(np.abs(down1)+np.spacing(1)))
    res2=np.sqrt(up/(np.abs(down2)+np.spacing(1)))
    return np.array([res1, res2])

def t_calc_information(p_x_given_t, PYgivenTs, PXs, PYs):
    """Calculate the MI - I(X;T) and I(Y;T)"""
    Hx = np.nansum(-np.dot(PXs, np.log2(PXs)))
    Hxt = - np.nansum((np.dot(np.multiply(p_x_given_t, np.log2(p_x_given_t)), PXs)))
    Hyt = - np.nansum(np.dot(PYgivenTs*np.log2(PYgivenTs+np.spacing(1)), PTs))
    Hy = np.nansum(-PYs * np.log2(PYs+np.spacing(1)))
    IYT = Hy - Hyt
    ITX = Hx - Hxt
    return ITX, IYT

def lanczos2(x):
    # f = (sin(pi*x) .* sin(pi*x/2) + eps) ./ ((pi^2 * x.^2 / 2) + eps);
    # f = f .* (abs(x) < 2);
    result = (np.sin(np.pi * x) * np.sin(np.pi * x / 2 + np.spacing(1))) / (
        (np.power(np.pi, 2) * np.power(x, 2) / 2) + np.spacing(1)
    )
    print (np.abs(result) < 2)
    result = result * (np.abs(result) < 2)
    return result

    def _evenly_spaced_condition_num_helper(self, p_order):
        import numpy as np

        x_vals = np.linspace(-1, 1, p_order + 1)
        leg_mat = self._call_func_under_test(x_vals)
        kappa2 = np.linalg.cond(leg_mat, p=2)
        # This gives the exponent of kappa2.
        base_exponent = np.log2(np.spacing(1.0))
        return int(np.round(np.log2(np.spacing(kappa2)) - base_exponent))

def estimate_parameters(left_stft, right_stft, stft_window_length):
	#List comprehensions are amazing!
	w_range = [i for j in (range(1, stft_window_length/2+1), range(-stft_window_length/2 + 1, 0)) for i in j]
	w = np.array([2*np.pi*i/(stft_window_length)
		for i in w_range
		if i is not 0])
	delay_estimation = -1/w * np.angle((right_stft + np.spacing(1))/(left_stft + np.spacing(1)))
	attenuation_estimation = np.absolute(right_stft/left_stft) - np.absolute(left_stft/right_stft)
	return attenuation_estimation, delay_estimation

def checkResult(Lbar,eigvec,eigval,k):
	"""
	"input
	"matrix Lbar and k eig values and k eig vectors
	"print norm(Lbar*eigvec[:,i]-lamda[i]*eigvec[:,i])
	"""
	check=[np.dot(Lbar,eigvec[:,i])-eigval[i]*eigvec[:,i] for i in range(0,k)]
	length=[np.linalg.norm(e) for e in check]/np.spacing(1)
	print("Lbar*v-lamda*v are %s*%s" % (length,np.spacing(1)))

def test_segmentation_utils():
    ctx = mx.context.current_context()
    import os
    if not os.path.isdir(os.path.expanduser('~/.mxnet/datasets/voc')):
        return

    transform_fn = transforms.Compose([
        transforms.ToTensor(),
        transforms.Normalize([.485, .456, .406], [.229, .224, .225])
    ])
    # get the dataset
    # TODO FIXME: change it to ADE20K dataset and pretrained model
    dataset = ADE20KSegmentation(split='val')
    # load pretrained net
    net = gluoncv.model_zoo.get_model('fcn_resnet50_ade', pretrained=True, ctx=ctx)
    # count for pixAcc and mIoU
    total_inter, total_union, total_correct, total_label = 0, 0, 0, 0
    np_inter, np_union, np_correct, np_label = 0, 0, 0, 0
    tbar = tqdm(range(10))
    for i in tbar:
        img, mask = dataset[i]
        # prepare data and make prediction
        img = transform_fn(img)
        img = img.expand_dims(0).as_in_context(ctx)
        mask = mask.expand_dims(0)
        pred = net.evaluate(img).as_in_context(mx.cpu(0))
        # gcv prediction
        correct1, labeled1 = batch_pix_accuracy(pred, mask)
        inter1, union1 = batch_intersection_union(pred, mask, dataset.num_class)
        total_correct += correct1
        total_label += labeled1
        total_inter += inter1
        total_union += union1
        pixAcc = 1.0 * total_correct / (np.spacing(1) + total_label)
        IoU = 1.0 * total_inter / (np.spacing(1) + total_union)
        mIoU = IoU.mean()

        # np predicition
        pred2 = np.argmax(pred.asnumpy().astype('int64'), 1) + 1
        mask2 = mask.squeeze().asnumpy().astype('int64') + 1
        _, correct2, labeled2 = pixelAccuracy(pred2, mask2)
        inter2, union2 = intersectionAndUnion(pred2, mask2, dataset.num_class)
        np_correct += correct2
        np_label += labeled2
        np_inter += inter2
        np_union += union2
        np_pixAcc = 1.0 * np_correct / (np.spacing(1) + np_label)
        np_IoU = 1.0 * np_inter / (np.spacing(1) + np_union)
        np_mIoU = np_IoU.mean()
        tbar.set_description('pixAcc: %.3f, np_pixAcc: %.3f, mIoU: %.3f, np_mIoU: %.3f'%\
            (pixAcc, np_pixAcc, mIoU, np_mIoU))

    np.testing.assert_allclose(total_inter, np_inter)
    np.testing.assert_allclose(total_union, np_union)
    np.testing.assert_allclose(total_correct, np_correct)
    np.testing.assert_allclose(total_label, np_label)