def __init__(self,m,n,l1_weight,eval_func, eval_func_for_test_set = None, output_func = None, tol = None):
        self.M = m
        self.N = n
        
        d = np.ones((n))
       
        self.diagH = sparse.diags(d,0)
        # to reference the element, use diag.data[0,0]

        self.iter = 0
        
        self.l1_weight = l1_weight
        self.eval_func = eval_func
        
        self.eval_func_for_test_set  = eval_func_for_test_set
        self.output_func = output_func

        self.x = scipy.mat([0]*n).T
        self.grad = scipy.mat([0]*n).T
        self.dir = scipy.mat([0]*n).T
        self.loss = 0
        
        if tol == None:
            self.tol = 1e-4
        else:
            self.tol = tol

def calcInvFisher(sigma, invSigma=None, factorSigma=None):
    """ Efficiently compute the exact inverse of the FIM of a Gaussian.
    Returns a list of the diagonal blocks. """
    if invSigma == None:
        invSigma = inv(sigma)
    if factorSigma == None:
        factorSigma = cholesky(sigma)
    dim = sigma.shape[0]

    invF = [mat(1 / (invSigma[-1, -1] + factorSigma[-1, -1] ** -2))]
    invD = 1 / invSigma[-1, -1]
    for k in reversed(list(range(dim - 1))):
        v = invSigma[k + 1:, k]
        w = invSigma[k, k]
        wr = w + factorSigma[k, k] ** -2
        u = dot(invD, v)
        s = dot(v, u)
        q = 1 / (w - s)
        qr = 1 / (wr - s)
        t = -(1 + q * s) / w
        tr = -(1 + qr * s) / wr
        invF.append(blockCombine([[qr, tr * u], [mat(tr * u).T, invD + qr * outer(u, u)]]))
        invD = blockCombine([[q , t * u], [mat(t * u).T, invD + q * outer(u, u)]])

    invF.append(sigma)
    invF.reverse()
    return invF

def steepest_descent(A, b, x0, tol=1e-8):
    """
    Uses the steepest descent method to find the x that satisfies Ax = b.

    Inputs:
        A: An m x n NumPy array
        b: An m x 1 NumPy array
        x0: An n x 1 NumPy array that represents the initial guess at a
            solution.
        tol (optional): The tolerance level for convergence. This is compared
                        against the norm(x_n+1 - x_n) each iteration.

    Outputs:
        x: The x that satisfies the equation.
    """
    A = sp.mat(A)
    b = sp.reshape(sp.mat(b),(b.size,1))


    def grad(A, b, x):
        """
        Find the gradient of ||Ax - b||
        Inputs:
            A: An m x n NumPy matrix.
            b: An m x 1 NumPy matrix.
            x: An n x a NumPy matrix.

        Outputs:
            grad: A NumPy matrix representing the gradient of ||Ax - b||
        """
        return np.mat(2  * A.T*(A*x - b))

    def solve_alpha_k(A, b, x):
        """
        Solves for alpha in the steepest descent algorithm
        x_n+1 = x_n - alpha * grad(x_n)

        Inputs:
            A: An m x n NumPy array
            b: An m x 1 NumPy array
            x: The x value where you want alpha to be defined for.

        Outputs:
            alpha: The alpha satisfying the algorithm above.
        """

        gradient = grad(A, b, x)
        return np.array(
            (gradient.T * gradient)/(2 * gradient.T * A.T * A * gradient))[0]



    xold = sp.reshape(sp.mat(x0),(x0.size,1))
    xnew = xold - grad(A, b, xold) * solve_alpha_k(A,b,xold)

    while la.norm(xold - xnew) > tol:
        xold = xnew
        xnew = xold - grad(A, b, xold) * solve_alpha_k(A,b,xold)

    return xnew

def dsimul(sys,u):
    """Simulate the discrete system sys
    Only for discrete systems!!!

    Call:
    y=dsimul(sys,u)

    Parameters
    ----------
    sys : Discrete System in State Space form
    u   : input vector
    Returns
    -------
    y: ndarray
    Simulation results

    """
    a=mat(sys.A)
    b=mat(sys.B)
    c=mat(sys.C)
    d=mat(sys.D)
    nx=shape(a)[0]
    ns=shape(u)[1]
    xk=zeros((nx,1))
    for i in arange(0,ns):
        uk=u[:,i]
        xk_1=a*xk+b*uk
        yk=c*xk+d*uk
        xk=xk_1
        if i==0:
            y=yk
        else:
            y=hstack((y,yk))
    y=array(y).T
    return y

def qhull(V, qstring):
    """
    Use qhull to determine convex hull / volume / normals.
     V - [matrix] vertices
     qstring - [string] arguments to pass to qhull
    """
    try:
        qhullp = subprocess.Popen(["qhull", qstring],
                              stdin=subprocess.PIPE, stdout=subprocess.PIPE)
        Vc = qhullp.communicate(qhullstr(V))[0] #qhull output to Vc
        
        if qstring == "FS": #calc area and volume
            ks = Vc.split('\n')[-2]
            Vol = float(ks.split(' ')[-2]) #get volume of D-hull
            return Vol
        elif qstring == "Ft": #calc vertices and facets
            ks = Vc.split('\n')
            fms = int(ks[1].split(' ')[1]) #get size of facet matrix
            fmat = ks[-fms-1:-1]
            fmat = mat(';'.join(fmat)) #generate matrix
            fmatv = fmat[:, 1:] #vertices on facets
            return array(fmatv)
        elif qstring == "n": #calc convex hull and get normals
            ks = ';'.join(Vc.split('\n')[2:]) #remove leading dimension output
            k = mat(ks[:-1]) #convert to martrix with vertices
            return array(k)
        else:
            exit(1)
    except:
        raise NameError('QhullError')

def con2vert(A, b):
    """
    Convert sets of constraints to a list of vertices (of the feasible region).
    If the shape is open, con2vert returns False for the closed property.
    """
    # Python implementation of con2vert.m by Michael Kleder (July 2005),
    #  available: http://www.mathworks.com/matlabcentral/fileexchange/7894
    #  -con2vert-constraints-to-vertices
    # Author: Michael Kelder (Original)
    #         Andre Campher (Python implementation)
    c = linalg.lstsq(mat(A), mat(b))[0]
    btmp = mat(b)-mat(A)*c
    D = mat(A)/matlib.repmat(btmp, 1, A.shape[1])

    fmatv = qhull(D, "Ft") #vertices on facets

    G  = zeros((fmatv.shape[0], D.shape[1]))
    for ix in range(0, fmatv.shape[0]):
        F = D[fmatv[ix, :], :].squeeze()
        G[ix, :] = linalg.lstsq(F, ones((F.shape[0], 1)))[0].transpose()

    V = G + matlib.repmat(c.transpose(), G.shape[0], 1)
    ux = uniqm(V)

    eps = 1e-13
    Av = dot(A, ux.T)
    bv = tile(b, (1, ux.shape[0]))
    closed = sciall(Av - bv <= eps)

    return ux, closed

def qnwcheb1(n, a, b):
    """ Univariate Gauss-Chebyshev quadrature nodes and weights

    Parameters
    -----------
    n : int
        number of nodes
    a : float
        left endpoint
    b : float
        right endpoint

    Returns
    ---------
    x : array, shape (n,)
        nodes
    x : array, shape (n,)
        weights

    Notes
    ---------
    
    Port of the qnwcheb1 function in the compecon matlab toolbox.
    """
    x = ((b + a) / 2 - (b - a) / 2
         * sp.cos(sp.pi / n * sp.arange(0.5, n + 0.5, 1)))
    w2 =  sp.r_[1, -2. / (sp.r_[1:(n - 1):2] * sp.r_[3:(n + 1):2])]
    w1 = (sp.cos(sp.pi / n * sp.mat((sp.r_[0:n] + 0.5)).T *
                 sp.mat((sp.r_[0:n:2]))).A)
    w0 = (b - a) / n
    w = w0 * sp.dot(w1, w2)
    return x, w

def get_projection_matrix(u,z,v,k=TOP_NUM_SINGULAR_VALUES):
    """
        generate the projection matrix which contains a 
        score vector for the patterns for each word pair 
    """ 
    #determine the best patterns for comparison   
    column_indexes=get_top_k_column_indexes(z)
    
    #using the column indexes of the best patterns recreate the 
    #u & z matrices containing only those correspoding columns
    
    #creating the uk matrix
    uk=[]
    for r in range(len(u)):
        uk.append([])
        for index in column_indexes:
            uk[r].append(u[r][index])
        r+=1

    #creating the zk matrix
    zk=[]
    for index in column_indexes:
        zk.append([])
        for col in range(len(v)):
            if (col==index):
                zk[len(zk)-1].append(z[index])
            else: 
                zk[len(zk)-1].append(0)
    
    #calcualte the projecttion matrix by u.z
    return mat(uk)*mat(zk)

def arnoldi(A, v0, k):
    """
    Arnoldi algorithm (Krylov approximation of a matrix)
        input: 
            A: matrix to approximate
            v0: initial vector (should be in matrix form) 
            k: number of Krylov steps 
        output: 
            V: matrix (large, N*k) containing the orthogonal vectors
            H: matrix (small, k*k) containing the Krylov approximation of A

    Author: Vasile Gradinaru, 14.12.2007 (Rennes)
    """
    #print 'ARNOLDI METHOD'
    inputtype = A.dtype.type
    V = mat( v0.copy() / norm(v0), dtype=inputtype)
    H = mat( zeros((k+1,k), dtype=inputtype) )
    for m in xrange(k):
        vt = A*V[ :, m]
        for j in xrange( m+1):
            H[ j, m] = (V[ :, j].H * vt )[0,0]
            vt -= H[ j, m] * V[:, j]
        H[ m+1, m] = norm(vt);
        if m is not k-1:
            V =  hstack( (V, vt.copy() / H[ m+1, m] ) ) 
    return V,  H

 def genInitSigmaFactor(self):
     """ depending on the algorithm settings, we start out with in identity matrix, or perturb it """
     if self.perturbedInitSigma:
         res = mat(eye(self.xdim)*self.initSigmaCoeff+randn(self.xdim, self.xdim)*self.initSigmaRandCoeff)            
     else:
         res = mat(eye(self.xdim)*self.initSigmaCoeff)
     return res   

def minreal(sys):
    """Minimal representation for state space systems

    Usage
    =====
    [sysmin]=minreal[sys]

    Inputs
    ------

    sys: system in ss or tf form

    Outputs
    -------
    sysfin: system in state space form
    """
    a=mat(sys.A)
    b=mat(sys.B)
    c=mat(sys.C)
    d=mat(sys.D)
    nx=shape(a)[0]
    ni=shape(b)[1]
    no=shape(c)[0]

    out=tb03ad(nx,no,ni,a,b,c,d,'R')

    nr=out[3]
    A=out[0][:nr,:nr]
    B=out[1][:nr,:ni]
    C=out[2][:no,:nr]
    sysf=ss(A,B,C,sys.D,sys.Tsamp)
    return sysf

def accuracy(p, colname='ClassLabel'):
    # find out how many classes are in the experiment
    numSamples = {}      # this is a list containing list of labels
    labels = p.column(col=2)[1]   # get the column of the actual labels
    for l in labels:
        if not(l in numSamples.keys()):
            numSamples.update({l:0})
    numLabels = len(numSamples.keys()) 
    # count number of samples per class
    labels = p.column(col=2)[1]
    for l in labels:
        numSamples[l] = numSamples[l] + 1  
    numLabels = len(numSamples.keys())
    confusionMatrix = scipy.mat( numpy.zeros( (numLabels,numLabels) ) ) 
    hdr_actual = '%s-actual' % colname
    hdr_predic = '%s-prediction' % colname
    mapLabels = {}
    cnt = 0
    for l in numSamples.keys():
       mapLabels.update({l:cnt})
       cnt = cnt + 1
    confusionMatrix = scipy.mat( numpy.zeros( (numLabels,numLabels) ) )  
    hdr_actual = '%s-actual' % colname
    hdr_predic = '%s-prediction' % colname
    actualLabels = p.column(hdr_actual)
    predictLabels = p.column(hdr_predic)
    for cnt in range(0,len(actualLabels[0])):
        pLabel = mapLabels[predictLabels[1][cnt]]     # predicted label
        aLabel = mapLabels[actualLabels[1][cnt]]      # actual Label
        if ( not(  predictLabels[0][cnt] ==   actualLabels[0][cnt]   ) ):     # it is just a sanity check, this event should happen ever. Just to make sure that it is comparing the same subject
              assert False, "This event should NOT happend ever !!! Are you sure you are using the correct Pyxel version??? I am comparing labels of two different subjects !! " 
        confusionMatrix[aLabel,pLabel] = confusionMatrix[aLabel,pLabel] + 1.0
    return confusionMatrix

def dcgain(sys):
    """Return the steady state value of the step response os sys

    Usage
    =====
    dcgain=dcgain(sys)

    Inputs
    ------

    sys: system

    Outputs
    -------
    dcgain : steady state value
    """

    a=mat(sys.A)
    b=mat(sys.B)
    c=mat(sys.C)
    d=mat(sys.D)
    nx=shape(a)[0]
    if sys.Tsamp!=0.0:
        a=a-eye(nx,nx)
    r=rank(a)
    if r<nx:
        gm=[]
    else:
        gm=-c*inv(a)*b+d
    return array(gm)

def full_obs(sys,poles):
    """Full order observer of the system sys

    Call:
    obs=full_obs(sys,poles)

    Parameters
    ----------
    sys : System in State Space form
    poles: desired observer poles

    Returns
    -------
    obs: ss
    Observer

    """
    if isinstance(sys, TransferFunction):
        "System must be in state space form"
        return
    a=mat(sys.A)
    b=mat(sys.B)
    c=mat(sys.C)
    d=mat(sys.D)
    L=place(a.T,c.T,poles)
    L=mat(L).T
    Ao=a-L*c
    Bo=hstack((b-L*d,L))
    n=shape(Ao)
    m=shape(Bo)
    Co=eye(n[0],n[1])
    Do=zeros((n[0],m[1]))
    obs=StateSpace(Ao,Bo,Co,Do,sys.dt)
    return obs

def ned2ecef(lat, lon, alt, n, e, d):
    X0, Y0, Z0 = coord.geodetic2ecef(lat, lon, alt)
    lat, lon = radians(lat), radians(lon)
    
    pitch = math.pi/2 + lat
    yaw = -lon 
    
    my = mat('[%f %f %f; %f %f %f; %f %f %f]' %
        (cos(pitch), 0, -sin(pitch),
         0,1,0,
         sin(pitch), 0, cos(pitch)))
    
    mz = mat('[%f %f %f; %f %f %f; %f %f %f]' %
        (cos(yaw), sin(yaw),0,
         -sin(yaw),cos(yaw),0,
         0,0,1))
    
    mr = mat('[%f %f %f; %f %f %f; %f %f %f]' %
        (-cos(lon)*sin(lat), -sin(lon), -cos(lat) * cos(lon), 
         -sin(lat)*sin(lon), cos(lon), -sin(lon)*cos(lat),
         cos(lat), 0, -sin(lat)))
    
    geo = mat('[%f; %f; %f]' % (X0, Y0, Z0))
    ned = mat('[%f; %f; %f]' % (n, e, d))
    res = mr*ned + geo
    return res[0], res[1], res[2]  

def acker(A,B,poles):
    """Pole placemenmt using Ackermann method

    Call:
    k=acker(A,B,poles)

    Parameters
    ----------
    A, B : State and input matrix of the system
    poles: desired poles

    Returns
    -------
    k: matrix
    State feedback gains

    """
    a=mat(A)
    b=mat(B)
    p=real(poly(poles))
    ct=ctrb(A,B)
    if det(ct)==0:
        k=0
        print "Pole placement invalid"
    else:
        n=size(p)
        pmat=p[n-1]*a**0
        for i in arange(1,n):
            pmat=pmat+p[n-i-1]*a**i
        k=inv(ct)*pmat
        k=k[-1][:]
    return k

def nonna_lsq(target, aux, idx=(), names=(), order=2):
	"""
	This function returns the coefficients of the least square prediction of the target
	signal, using the auxiliary signals and their powers, as specified by the order argument.
	
	Input arguments:
	target = target signal
	aux    = matrix of auxiliary signals
	idx    = boolean vector to select a subset of the data for the LSQ fit
	order  = order of the polynomial of aux signals to be used in the fit, default is 2
	names  = list of the auxiliary signal names
	
	Output:
	p      = list of coefficients
	X      = matrix of the signals used in the reconstruction
	cnames = list of the corresponding signals
	
	Note that the mean will be removed from the auxiliary signals. 
	"""
	# number of auxiliary channels
	naux = scipy.shape(aux[1])
	
	if len(names) == 0:
		# since the user didn't provide signal names, let's build some
		names = map(lambda x: 'S'+str(x), scipy.arange(naux)+1)
		
	if len(idx) == 0:
		# no index means use all
		idx = numpy.array(target, dtype=bool)
		idx[:] = True
	
	##### PREPARE CHANNELS FOR LSQ PREDICTION 

	# prepare channels and their squared values
	X = scipy.zeros((scipy.shape(aux)[0], order*scipy.shape(aux)[1]+1))
	cnames = []
	for i in range(scipy.shape(aux)[1]):
		for j in range(order):
			# add the (j+1)th power of the signal after removing the mean
			X[:,order*i+j] = numpy.power((aux[:,i] - scipy.mean(aux[idx,i])), j+1)
			# then remove the mean of the result
			X[:,order*i+j] = X[:,order*i+j] - scipy.mean(X[idx,order*i+j])
			# save the name, including the power
			if j==0:
				cnames.append(names[i])
			else:
				cnames.append(names[i]+'^'+str(j+1))
				
	# add a constant at the end of the list
	X[:,-1] = 1
	cnames.append('1')
	# convert to matrix object for simpler manipulation
	X = scipy.mat(X)
	
	##### best estimate of coefficients to minimize the squared error
	p = scipy.linalg.inv(X[idx,:].T * X[idx,:]) * X[idx,:].T * scipy.mat(target[idx]).T

	# return all the results
	return p, X, cnames

def kNN(segment_list, k, emb_dict):
  
  pred_matrix = []
  gold_word_list = []
  
  '''
  # TEST
  segment_list = []
  counter = 0
  for key, value in emb_dict.items():
    segment_list.append(key)
    counter += 1
    if counter == 5:
      break
  '''
  
  for seg in segment_list:
    print(len(emb_dict[seg]))
    print(seg)
    pred_matrix.append(emb_dict[seg])
    gold_word_list.append(seg)
    
  # Obtain now the closest neighbors.
  word_list = []
  emb_matrix = []
  for key, value in emb_dict.items():
    if containsPunctuation(key + " "):
      continue
    word_list.append(key)
    emb_matrix.append(value)
 
  pred_matrix=mat(pred_matrix)
  emb_matrix = mat(emb_matrix)
  rows=len(segment_list) + 1

  print('Calculating simi_matrix...')
  print(pred_matrix.shape)
  print(emb_matrix.shape)
  simi_matrix=1-cdist(pred_matrix, emb_matrix, 'cosine')
  print('...simi_matrix done!')
    
  max_index_matrix=simi_matrix.argsort()[:,-k-1:]
  #max_index_matrix=simi_matrix.argsort()[:,k+1:]

  pred_word_matrix=[]
  for row in range(max_index_matrix.shape[0]):
    pred_list=[word_list[i] for i in max_index_matrix[row]] 
    pred_word_matrix.append(pred_list)

  for i in range(len(segment_list)):
    print('The original segment: ' + gold_word_list[i])
    #for j in [4,3,2,1,0]: 
    for ij in range(k):
            j = k - ij
            #if pred_word_matrix[i][j] == gold_word_list[i]:
            #    continue
            print('Predicted neighbor segment: ' + pred_word_matrix[i][j])
            #print('Distance: ' + str(cdist([word_list[i]], [emb_dict[pred_word_matrix[i][j]]], 'cosine')))
            print('Similarity: ' + str(1-cdist([emb_dict[gold_word_list[i]]], [emb_dict[pred_word_matrix[i][j]]], 'cosine')))

def scipy_cosine_similarity(tf_idf1, tf_idf2):
  # input; vector 1, 2 (vector의 길이가 다른 경우 없는 key에 대해서는 0을 채워서 같은 길이로 만든 vector임)
  # output; cosine similarity
  # scipy library를 사용해 계산

  # http://stackoverflow.com/questions/21980644/calculate-cosine-similarity-of-two-matrices-python
  a, b = mat(tf_idf1), mat(tf_idf2)
  c = dot(a, b.T) / linalg.norm(a) / linalg.norm(b)
  return c.A1[0]    # http://stackoverflow.com/questions/3337301/numpy-matrix-to-array

 def distanceState(self, belief1, belief2):
     '''
     Calculate distance between two beliefs by 1 - (cosine similarity)
     '''
     if len(belief1) != len(belief2):
         return 2.0
     b1 = mat(belief1)
     b2 = mat(belief2)
     cosSim = dot(b1,b2.T)/linalg.norm(b1)/linalg.norm(b2)
     return 1.0 - cosSim[0,0]