Python numpy.dot() Examples

The following are code examples for showing how to use numpy.dot(). They are from open source Python projects. You can vote up the examples you like or vote down the ones you don't like.

Example 1
Project: b2ac   Author: hbldh   File: reference.py    MIT License 8 votes vote down vote up
def _calculate_scatter_matrix_py(x, y):
    """Calculates the complete scatter matrix for the input coordinates.

    :param x: The x coordinates.
    :type x: :py:class:`numpy.ndarray`
    :param y: The y coordinates.
    :type y: :py:class:`numpy.ndarray`
    :return: The complete scatter matrix.
    :rtype: :py:class:`numpy.ndarray`

    """
    D = np.ones((len(x), 6), dtype=x.dtype)
    D[:, 0] = x * x
    D[:, 1] = x * y
    D[:, 2] = y * y
    D[:, 3] = x
    D[:, 4] = y

    return D.T.dot(D) 
Example 2
Project: b2ac   Author: hbldh   File: matrix_algorithms.py    MIT License 7 votes vote down vote up
def QR_factorisation_Householder_double(A):
    """Perform QR factorisation in double floating-point precision.

    :param A: The matrix to factorise.
    :type A: :py:class:`numpy.ndarray`
    :returns: The matrix Q and the matrix R.
    :rtype: tuple

    """
    A = np.array(A, 'float')

    n, m = A.shape
    V = np.zeros_like(A, 'float')
    for k in xrange(n):
        V[k:, k] = A[k:, k].copy()
        V[k, k] += np.sign(V[k, k]) * np.linalg.norm(V[k:, k], 2)
        V[k:, k] /= np.linalg.norm(V[k:, k], 2)
        A[k:, k:] -= 2 * np.outer(V[k:, k], np.dot(V[k:, k], A[k:, k:]))
    R = np.triu(A[:n, :n])

    Q = np.eye(m, n)
    for k in xrange((n - 1), -1, -1):
        Q[k:, k:] -= np.dot((2 * (np.outer(V[k:, k], V[k:, k]))), Q[k:, k:])
    return Q, R 
Example 3
Project: b2ac   Author: hbldh   File: double.py    MIT License 7 votes vote down vote up
def _calculate_scatter_matrix_double(x, y):
    """Calculates the complete scatter matrix for the input coordinates.

    :param x: The x coordinates.
    :type x: :py:class:`numpy.ndarray`
    :param y: The y coordinates.
    :type y: :py:class:`numpy.ndarray`
    :return: The complete scatter matrix.
    :rtype: :py:class:`numpy.ndarray`

    """
    D = np.ones((len(x), 6), 'int64')
    D[:, 0] = x * x
    D[:, 1] = x * y
    D[:, 2] = y * y
    D[:, 3] = x
    D[:, 4] = y

    return D.T.dot(D) 
Example 4
Project: tmx2map   Author: joshuaskelly   File: mathhelper.py    MIT License 7 votes vote down vote up
def angle_between(dest, base=(1,0,0)):
    """Returns the angle in positive degrees from base to dest in the
    xy-plane.

    dest: A vector
    base: A vector

    Returns:
        Angle in degrees [0, 360)
    """

    target = dest[0], dest[1]
    p_axis = -base[1], base[0]
    b_axis = base[0], base[1]

    x_proj = numpy.dot(target, b_axis)
    y_proj = numpy.dot(target, p_axis)

    result = math.degrees(math.atan2(y_proj, x_proj))

    return (result + 360) % 360 
Example 5
Project: explirefit   Author: codogogo   File: text_embeddings.py    Apache License 2.0 7 votes vote down vote up
def most_similar(self, embedding, target_lang, num, similarity = True):
		ms = []
		for w in self.lang_vocabularies[target_lang]:
			targ_w_emb = self.get_vector(target_lang, w)
			if len(embedding) != len(targ_w_emb):
				print("Unaligned embedding length: " + w)
			else:
				if similarity:
					nrm = np.linalg.norm(embedding, 2)
					trg_nrm = self.get_norm(target_lang, w)
					sim = np.dot(embedding, targ_w_emb) / (nrm * trg_nrm)
					if (len(ms) < num) or (sim > ms[-1][1]):	
						ms.append((w, sim))
						ms.sort(key = lambda x: x[1], reverse = True)
				else:
					dist = np.linalg.norm(embedding - targ_w_emb)
					if (len(ms) < num) or (dist < ms[-1][1]):	
						ms.append((w, dist))
						ms.sort(key = lambda x: x[1])
				if len(ms) > num: 
						ms.pop() 
		return [ws for ws in ms] 
Example 6
Project: building-boundary   Author: Geodan   File: bounding_box.py    MIT License 7 votes vote down vote up
def rotate_points(points, angle):
    """
    Rotate points in a coordinate system using a rotation matrix based on
    an angle.

    Parameters
    ----------
    points : (Mx2) array
        The coordinates of the points.
    angle : float
        The angle by which the points will be rotated (in radians).

    Returns
    -------
    points_rotated : (Mx2) array
        The coordinates of the rotated points.
    """
    # Compute rotation matrix
    rot_matrix = np.array(((math.cos(angle), -math.sin(angle)),
                           (math.sin(angle), math.cos(angle))))
    # Apply rotation matrix to the points
    points_rotated = np.dot(points, rot_matrix)

    return np.array(points_rotated) 
Example 7
Project: FRIDA   Author: LCAV   File: point_cloud.py    MIT License 7 votes vote down vote up
def classical_mds(self, D):
        ''' 
        Classical multidimensional scaling

        Parameters
        ----------
        D : square 2D ndarray
            Euclidean Distance Matrix (matrix containing squared distances between points
        '''

        # Apply MDS algorithm for denoising
        n = D.shape[0]
        J = np.eye(n) - np.ones((n,n))/float(n)
        G = -0.5*np.dot(J, np.dot(D, J))

        s, U = np.linalg.eig(G)

        # we need to sort the eigenvalues in decreasing order
        s = np.real(s)
        o = np.argsort(s)
        s = s[o[::-1]]
        U = U[:,o[::-1]]

        S = np.diag(s)[0:self.dim,:]
        self.X = np.dot(np.sqrt(S),U.T) 
Example 8
Project: FRIDA   Author: LCAV   File: point_cloud.py    MIT License 7 votes vote down vote up
def flatten(self, ind):
        '''
        Transform the set of points so that the subset of markers given as argument is
        as close as flat (wrt z-axis) as possible.

        Parameters
        ----------
        ind : list of bools
            Lists of marker indices that should be all in the same subspace
        '''

        # Transform references to indices if needed
        ind = [self.key2ind(s) for s in ind]

        # center point cloud around the group of indices
        centroid = self.X[:,ind].mean(axis=1, keepdims=True)
        X_centered = self.X - centroid

        # The rotation is given by left matrix of SVD
        U,S,V = la.svd(X_centered[:,ind], full_matrices=False)

        self.X = np.dot(U.T, X_centered) + centroid 
Example 9
Project: FRIDA   Author: LCAV   File: doa.py    MIT License 7 votes vote down vote up
def compute_mode(self):
        """
        Pre-compute mode vectors from candidate locations (in spherical 
        coordinates).
        """
        if self.num_loc is None:
            raise ValueError('Lookup table appears to be empty. \
                Run build_lookup().')
        self.mode_vec = np.zeros((self.max_bin,self.M,self.num_loc), 
            dtype='complex64')
        if (self.nfft % 2 == 1):
            raise ValueError('Signal length must be even.')
        f = 1.0 / self.nfft * np.linspace(0, self.nfft / 2, self.max_bin) \
            * 1j * 2 * np.pi
        for i in range(self.num_loc):
            p_s = self.loc[:, i]
            for m in range(self.M):
                p_m = self.L[:, m]
                if (self.mode == 'near'):
                    dist = np.linalg.norm(p_m - p_s, axis=1)
                if (self.mode == 'far'):
                    dist = np.dot(p_s, p_m)
                # tau = np.round(self.fs*dist/self.c) # discrete - jagged
                tau = self.fs * dist / self.c  # "continuous" - smoother
                self.mode_vec[:, m, i] = np.exp(f * tau) 
Example 10
Project: FRIDA   Author: LCAV   File: generators.py    MIT License 6 votes vote down vote up
def gen_visibility(alphak, phi_k, pos_mic_x, pos_mic_y):
    """
    generate visibility from the Dirac parameter and microphone array layout
    :param alphak: Diracs' amplitudes
    :param phi_k: azimuths
    :param pos_mic_x: a vector that contains microphones' x coordinates
    :param pos_mic_y: a vector that contains microphones' y coordinates
    :return:
    """
    xk, yk = polar2cart(1, phi_k)
    num_mic = pos_mic_x.size
    visi = np.zeros((num_mic, num_mic), dtype=complex)
    for q in xrange(num_mic):
        p_x_outer = pos_mic_x[q]
        p_y_outer = pos_mic_y[q]
        for qp in xrange(num_mic):
            p_x_qqp = p_x_outer - pos_mic_x[qp]  # a scalar
            p_y_qqp = p_y_outer - pos_mic_y[qp]  # a scalar
            visi[qp, q] = np.dot(np.exp(-1j * (xk * p_x_qqp + yk * p_y_qqp)), alphak)
    return visi 
Example 11
Project: FRIDA   Author: LCAV   File: tools_fri_doa_plane.py    MIT License 6 votes vote down vote up
def cov_mtx_est(y_mic):
    """
    estimate covariance matrix
    :param y_mic: received signal (complex based band representation) at microphones
    :return:
    """
    # Q: total number of microphones
    # num_snapshot: number of snapshots used to estimate the covariance matrix
    Q, num_snapshot = y_mic.shape
    cov_mtx = np.zeros((Q, Q), dtype=complex, order='F')
    for q in range(Q):
        y_mic_outer = y_mic[q, :]
        for qp in range(Q):
            y_mic_inner = y_mic[qp, :]
            cov_mtx[qp, q] = np.dot(y_mic_outer, y_mic_inner.T.conj())
    return cov_mtx / num_snapshot 
Example 12
Project: FRIDA   Author: LCAV   File: tools_fri_doa_plane.py    MIT License 6 votes vote down vote up
def mtx_updated_G(phi_recon, M, mtx_amp2visi_ri, mtx_fri2visi_ri):
    """
    Update the linear transformation matrix that links the FRI sequence to the
    visibilities by using the reconstructed Dirac locations.
    :param phi_recon: the reconstructed Dirac locations (azimuths)
    :param M: the Fourier series expansion is between -M to M
    :param p_mic_x: a vector that contains microphones' x-coordinates
    :param p_mic_y: a vector that contains microphones' y-coordinates
    :param mtx_freq2visi: the linear mapping from Fourier series to visibilities
    :return:
    """
    L = 2 * M + 1
    ms_half = np.reshape(np.arange(-M, 1, step=1), (-1, 1), order='F')
    phi_recon = np.reshape(phi_recon, (1, -1), order='F')
    mtx_amp2freq = np.exp(-1j * ms_half * phi_recon)  # size: (M + 1) x K
    mtx_amp2freq_ri = np.vstack((mtx_amp2freq.real, mtx_amp2freq.imag[:-1, :]))  # size: (2M + 1) x K
    mtx_fri2amp_ri = linalg.lstsq(mtx_amp2freq_ri, np.eye(L))[0]
    # projection mtx_freq2visi to the null space of mtx_fri2amp
    mtx_null_proj = np.eye(L) - np.dot(mtx_fri2amp_ri.T,
                                       linalg.lstsq(mtx_fri2amp_ri.T, np.eye(L))[0])
    G_updated = np.dot(mtx_amp2visi_ri, mtx_fri2amp_ri) + \
                np.dot(mtx_fri2visi_ri, mtx_null_proj)
    return G_updated 
Example 13
Project: DataHack2018   Author: InnovizTech   File: vis_utils.py    BSD 3-Clause "New" or "Revised" License 6 votes vote down vote up
def _zoom(self, value):
        if type(value) is str:
            if value == '-z':
                value = -base.camera.getPos()[2]
            if value == '+z':
                value = base.camera.getPos()[2]
        scale = 0.3
        hpr = base.camera.getHpr()
        rotation_angles = np.array((-hpr[1], 0., -hpr[0])) / 180. * np.pi
        rotation = euler_matrix(rotation_angles[2], rotation_angles[1], rotation_angles[0])
        # inverse
        R = rotation.T
        offset = R.dot(np.array([[0., 1., 0.]]).T)

        new_pos = np.array(base.camera.getPos()) + offset.ravel() * scale * value
        new_pos[2] = max(new_pos[2], 0.)
        base.camera.setPos(new_pos[0], new_pos[1], new_pos[2])
        self._update_location_text() 
Example 14
Project: DataHack2018   Author: InnovizTech   File: math_utils.py    BSD 3-Clause "New" or "Revised" License 6 votes vote down vote up
def quaternion_to_rotation(q):
    """
    Conversion from quaternion vector (w,x,y,z) to a rotation matrix.
    Based on https://en.wikipedia.org/wiki/Rotation_matrix#Quaternion
    """
    w, x, y, z = tuple(q)
    n = np.dot(q, q)
    s = 0. if n == 0. else 2./n
    wx, wy, wz = s*w*x, s*w*y, s*w*z
    xx, xy, xz = s*x*x, s*x*y, s*x*z
    yy, yz, zz = s*y*y, s*y*z, s*z*z

    R = np.array([[1-yy-zz, xy-wz, xz+wy],
                  [xy+wz, 1-xx-zz, yz-wx],
                  [xz-wy, yz+wx, 1-xx-yy]])
    return R 
Example 15
Project: ml_news_popularity   Author: khaledJabr   File: lr.py    MIT License 6 votes vote down vote up
def update_beta(self, data, targets, betas, j):
        n = len(targets)
        x_no_j = np.delete(data, j, 1)
        b_no_j = np.delete(betas, j, 0)
        norm_x_j = np.linalg.norm(data[:, j])
        # predict_no_j = np.dot(b_no_j, x_no_j.T)
        predicted = x_no_j.dot(b_no_j)

        # use normalized targets
        residuals = targets - predicted
        # p_j = data[:,j].dot(residuals)
        z_j = norm_x_j.dot(residuals)/n
        # print("Max Residual = {}".format(np.max(residuals)))
        # print(p_j)
        # update beta_j based on residuals

        # return self.soft_threshold(p_j, self.alpha*n/2)/(norm_x_j**2)
        return self.soft_threshold(z_j, self.alpha) 
Example 16
Project: StructEngPy   Author: zhuoju36   File: dynamic.py    MIT License 6 votes vote down vote up
def Riz_mode(model:Model,n,F):
    """
    Solve the Riz mode of the MDOF system\n
    n: number of modes to extract\n
    F: spacial load pattern
    """
    #            alpha=np.array(mode).T
#            #Grum-Schmidt procedure            
#            beta=[]
#            for i in range(len(mode)):
#                beta_i=alpha[i]
#                for j in range(i):
#                    beta_i-=np.dot(alpha[i],beta[j])/np.dot(alpha[j],alpha[j])*beta[j]
#                beta.append(beta_i)
#            mode_=np.array(beta).T
    pass 
Example 17
Project: StructEngPy   Author: zhuoju36   File: dynamic.py    MIT License 6 votes vote down vote up
def spectrum_analysis(model,n,spec):
    """
    sepctrum analysis
    
    params:
        n: number of modes to use\n
        spec: a list of tuples (period,acceleration response)
    """
    freq,mode=eigen_mode(model,n)
    M_=np.dot(mode.T,model.M)
    M_=np.dot(M_,mode)
    K_=np.dot(mode.T,model.K)
    K_=np.dot(K_,mode)
    C_=np.dot(mode.T,model.C)
    C_=np.dot(C_,mode)
    d_=[]
    for (m_,k_,c_) in zip(M_.diag(),K_.diag(),C_.diag()):
        sdof=SDOFSystem(m_,k_)
        T=sdof.omega_d()
        d_.append(np.interp(T,spec[0],spec[1]*m_))
    d=np.dot(d_,mode)
    #CQC
    return d 
Example 18
Project: StructEngPy   Author: zhuoju36   File: __init__.py    MIT License 6 votes vote down vote up
def resolve_modal_displacement(self,node_id,k): 
        """
        resolve modal node displacement.
        
        params:
            node_id: int.
            k: order of vibration mode.
        return:
            6-array of local nodal displacement.
        """
        if not self.is_solved:
            raise Exception('The model has to be solved first.')
        if node_id in self.__nodes.keys():
            node=self.__nodes[node_id]
            T=node.transform_matrix
            return T.dot(self.mode_[node_id*6:node_id*6+6,k-1]).reshape(6)
        else:
            raise Exception("The node doesn't exists.") 
Example 19
Project: StructEngPy   Author: zhuoju36   File: element.py    MIT License 6 votes vote down vote up
def _BtDB(self,s,r):
        """
        dot product of B^T, D, B
        params:
            s,r:natural position of evalue point.2-array.
        returns:
            3x3 matrix.
        """
        print(self._B(s,r).transpose(2,0,1).shape)
        print(
            np.matmul(
                np.dot(self._B(s,r).T,self._D),
                self._B(s,r).transpose(2,0,1)).shape
            )
        print(self._D.shape)
       
        
        return np.matmul(np.dot(self._B(s,r).T,self._D),self._B(s,r).transpose(2,0,1)).transpose(1,2,0) 
Example 20
Project: StructEngPy   Author: zhuoju36   File: csys.py    MIT License 6 votes vote down vote up
def __init__(self,origin, pt1, pt2, name=None):
        """
        origin: 3x1 vector
        pt1: 3x1 vector
        pt2: 3x1 vector
        """
        self.__origin=origin    
        vec1 = np.array([pt1[0] - origin[0] , pt1[1] - origin[1] , pt1[2] - origin[2]])
        vec2 = np.array([pt2[0] - origin[0] , pt2[1] - origin[1] , pt2[2] - origin[2]])
        cos = np.dot(vec1, vec2)/np.linalg.norm(vec1)/np.linalg.norm(vec2)
        if  cos == 1 or cos == -1:
            raise Exception("Three points should not in a line!!")        
        self.__x = vec1/np.linalg.norm(vec1)
        z = np.cross(vec1, vec2)
        self.__z = z/np.linalg.norm(z)
        self.__y = np.cross(self.z, self.x)
        self.__name=uuid.uuid1() if name==None else name 
Example 21
Project: StructEngPy   Author: zhuoju36   File: csys.py    MIT License 6 votes vote down vote up
def set_by_3pts(self,origin, pt1, pt2):
        """
        origin: tuple 3
        pt1: tuple 3
        pt2: tuple 3
        """
        self.origin=origin    
        vec1 = np.array([pt1[0] - origin[0] , pt1[1] - origin[1] , pt1[2] - origin[2]])
        vec2 = np.array([pt2[0] - origin[0] , pt2[1] - origin[1] , pt2[2] - origin[2]])
        cos = np.dot(vec1, vec2)/np.linalg.norm(vec1)/np.linalg.norm(vec2)
        if  cos == 1 or cos == -1:
            raise Exception("Three points should not in a line!!")        
        self.x = vec1/np.linalg.norm(vec1)
        z = np.cross(vec1, vec2)
        self.z = z/np.linalg.norm(z)
        self.y = np.cross(self.z, self.x) 
Example 22
Project: b2ac   Author: hbldh   File: matrix_algorithms.py    MIT License 5 votes vote down vote up
def QR_factorisation_Givens_double(A):

    n, m = A.shape
    R = np.array(A, dtype='float')
    Q = np.eye(n)
    for i in xrange(m - 1):
        for j in xrange(n - 1, i, -1):
            G = Givens_rotation_matrix_double(R[j - 1, i], R[j, i])
            R[(j - 1):(j + 1), :] = np.dot(G, R[(j - 1):(j + 1), :])
            Q[(j - 1):(j + 1), :] = np.dot(G, Q[(j - 1):(j + 1), :])
    return Q.T, np.triu(R) 
Example 23
Project: b2ac   Author: hbldh   File: reference.py    MIT License 5 votes vote down vote up
def fit_improved_B2AC_int(points):
    """Ellipse fitting in Python with improved B2AC algorithm as described in
    this `paper <http://autotrace.sourceforge.net/WSCG98.pdf>`_.

    This version of the fitting uses int64 storage during calculations and performs the
    eigensolver on an integer array.

    :param points: The [Nx2] array of points to fit ellipse to.
    :type points: :py:class:`numpy.ndarray`
    :return: The conic section array defining the fitted ellipse.
    :rtype: :py:class:`numpy.ndarray`

    """
    S = _calculate_scatter_matrix_c(points[:, 0], points[:, 1])
    S1 = np.array([S[0, 0], S[0, 1], S[0, 2], S[1, 1], S[1, 2], S[2, 2]])
    S3 = np.array([S[3, 3], S[3, 4], S[3, 5], S[4, 4], S[4, 5], S[5, 5]])
    adj_S3, det_S3 = inverse_symmetric_3by3_int(S3)
    S2 = S[:3, 3:]
    T_no_det = - np.dot(np.array(adj_S3.reshape((3, 3)), 'int64'), np.array(S2.T, 'int64'))
    M_term2 = np.dot(np.array(S2, 'int64'), T_no_det) // det_S3
    M = add_symmetric_matrix(M_term2, S1)
    M[[0, 2], :] /= 2
    M[1, :] = -M[1, :]

    e_vals, e_vect = np.linalg.eig(M)

    try:
        elliptical_solution_index = np.where(((4 * e_vect[0, :] * e_vect[2, :]) - ((e_vect[1, :] ** 2))) > 0)[0][0]
    except:
        # No positive eigenvalues. Fit was not ellipse.
        raise ArithmeticError("No elliptical solution found.")
    a = e_vect[:, elliptical_solution_index]
    return np.concatenate((a, np.dot(T_no_det, a) / det_S3)) 
Example 24
Project: b2ac   Author: hbldh   File: int.py    MIT License 5 votes vote down vote up
def fit_improved_B2AC_int(points):
    """Ellipse fitting in Python with improved B2AC algorithm as described in
    this `paper <http://autotrace.sourceforge.net/WSCG98.pdf>`_.

    This version of the fitting uses int64 storage during calculations and performs the
    eigensolver on an integer array.

    :param points: The [Nx2] array of points to fit ellipse to.
    :type points: :py:class:`numpy.ndarray`
    :return: The conic section coefficients array defining the fitted ellipse.
    :rtype: :py:class:`numpy.ndarray`

    """
    e_conds = []
    M, T_no_det, determinant_S3 = _calculate_M_and_T_int64(points)

    e_vals = sorted(QR_algorithm_shift_Givens_int(M)[0])

    a = None
    for ev_ind in [1, 2, 0]:
        # Find the eigenvector that matches this eigenvector.
        eigenvector, e_norm = inverse_iteration_for_eigenvector_int(M, e_vals[ev_ind])
        # See if that eigenvector yields an elliptical solution.
        elliptical_condition = (4 * eigenvector[0] * eigenvector[2]) - (eigenvector[1] ** 2)
        e_conds.append(elliptical_condition)
        if elliptical_condition > 0:
            a = eigenvector
            break

    if a is None:
        raise ArithmeticError("No elliptical solution found.")

    conic_coefficients = np.concatenate((a, np.dot(T_no_det, a) // determinant_S3))
    return conic_coefficients 
Example 25
Project: b2ac   Author: hbldh   File: int.py    MIT License 5 votes vote down vote up
def _calculate_M_and_T_int64(points):
    """Part of the B2AC ellipse fitting algorithm, calculating the M and T
     matrices needed.

     This integer implementation also returns the determinant of the
     scatter matrix, which hasn't been applied to the matrix T yet.

     M is exact in integer values, but is truncated towards zero compared
     to the double version.

    :param points: The [Nx2] array of points to fit ellipse to.
    :type points: :py:class:`numpy.ndarray`
    :return: M, T undivided by determinant and the determinant.
    :rtype: tuple

    """
    S = _calculate_scatter_matrix_c(points[:, 0], points[:, 1])
    # Extract the symmetric parts of the S matrix.
    S1 = np.array([S[0, 0], S[0, 1], S[0, 2], S[1, 1], S[1, 2], S[2, 2]], dtype='int64')
    S3 = np.array([S[3, 3], S[3, 4], S[3, 5], S[4, 4], S[4, 5], S[5, 5]], dtype='int64')

    adj_S3, det_S3 = inverse_symmetric_3by3_int64(S3)
    S2 = S[:3, 3:]

    T_no_det = - np.dot(np.array(adj_S3.reshape((3, 3)), 'int64'), np.array(S2.T, 'int64'))
    T_no_det, det_S3 = scale_T_matrix(T_no_det, det_S3)
    M_term_2 = np.dot(np.array(S2, 'int64'), T_no_det) // det_S3
    M = matrix_add_symmetric(M_term_2, S1)
    M[[0, 2], :] //= 2
    M[1, :] = -M[1, :]

    return M, T_no_det, det_S3 
Example 26
Project: b2ac   Author: hbldh   File: double.py    MIT License 5 votes vote down vote up
def fit_improved_B2AC_double(points):
    """Ellipse fitting in Python with improved B2AC algorithm as described in
    this `paper <http://autotrace.sourceforge.net/WSCG98.pdf>`_.

    This version of the fitting uses float storage during calculations and performs the
    eigensolver on a float array. It only uses `b2ac` package methods for fitting, to
    be as similar to the integer implementation as possible.

    :param points: The [Nx2] array of points to fit ellipse to.
    :type points: :py:class:`numpy.ndarray`
    :return: The conic section array defining the fitted ellipse.
    :rtype: :py:class:`numpy.ndarray`

    """
    e_conds = []
    points = np.array(points, 'float')

    M, T = _calculate_M_and_T_double(points)

    e_vals = sorted(qr.QR_algorithm_shift_Givens_double(M)[0])

    a = None
    for ev_ind in [1, 2, 0]:
        # Find the eigenvector that matches this eigenvector.
        eigenvector = inv_iter.inverse_iteration_for_eigenvector_double(M, e_vals[ev_ind], 5)

        # See if that eigenvector yields an elliptical solution.
        elliptical_condition = (4 * eigenvector[0] * eigenvector[2]) - (eigenvector[1] ** 2)
        e_conds.append(elliptical_condition)
        if elliptical_condition > 0:
            a = eigenvector
            break

    if a is None:
        print("Eigenvalues = {0}".format(e_vals))
        print("Elliptical conditions = {0}".format(e_conds))
        raise ArithmeticError("No elliptical solution found.")

    conic_coefficients = np.concatenate((a, np.dot(T, a)))
    return conic_coefficients 
Example 27
Project: b2ac   Author: hbldh   File: inverse_iteration.py    MIT License 5 votes vote down vote up
def inverse_iteration_for_eigenvector_double(A, eigenvalue, n_iterations=1):
    """Performs a series of inverse iteration steps with a known
    eigenvalue to produce its eigenvector.

    :param A: The 3x3 matrix to which the eigenvalue belongs.
    :type A: :py:class:`numpy.ndarray`
    :param eigenvalue: One eigenvalue of the matrix A.
    :type eigenvalue: float
    :param n_iterations: Number of iterations to perform the multiplication
     with the inverse. For a accurate eigenvalue, one iteration is enough
     for a ~1e-6 correct eigenvector. More than five is usually unnecessary.
    :type n_iterations: int
    :return: The eigenvector of this matrix and eigenvalue combination.
    :rtype: :py:class:`numpy.ndarray`

    """
    A = np.array(A, 'float')
    # Subtract the eigenvalue from the diagonal entries of the matrix.
    # N_POLYPOINTS.B. Also slightly perturb the eigenvalue so the matrix will
    # not be so close to singular!
    for k in xrange(A.shape[0]):
        A[k, k] -= eigenvalue + 0.001
    # Obtain the inverse of the matrix.
    A_inv = mo.inverse_3by3_double(A).reshape((3, 3))
    # Instantiate the eigenvector to iterate with.
    eigenvector = np.ones((A.shape[0], ), 'float')
    eigenvector /= np.linalg.norm(eigenvector)
    # Perform the desired number of iterations.
    for k in xrange(n_iterations):
        eigenvector = np.dot(A_inv, eigenvector)
        eigenvector /= np.linalg.norm(eigenvector)

    if np.any(np.isnan(eigenvector)) or np.any(np.isinf(eigenvector)):
        print("Nan and/or Infs in eigenvector!")

    if (eigenvector[0] < 0) and (eigenvector[2] < 0):
        eigenvector = -eigenvector
    return eigenvector 
Example 28
Project: projection-methods   Author: akshayka   File: polyak.py    GNU General Public License v3.0 5 votes vote down vote up
def solve(self, problem):
        left_set = problem.sets[0]
        right_set = problem.sets[1]

        iterate = (self._initial_iterate if
            self._initial_iterate is not None else np.ones(problem.dimension))
        iterates = [iterate]
        residuals = []

        status = Optimizer.Status.INACCURATE
        for i in xrange(self.max_iters):
            if self.verbose:
                print 'iteration %d' % i
            x_k = iterates[-1]
            x_k_1 = left_set.project(x_k)
            tmp = right_set.project(x_k)

            residuals.append(self._compute_residual(x_k, x_k_1, tmp))
            if self.verbose:
                print '\tresidual: %e' % sum(residuals[-1])
            if self._is_optimal(residuals[-1]):
                status = Optimizer.Status.OPTIMAL
                if not self.do_all_iters:
                    break

            x_k_2 = right_set.project(x_k_1)
            x_k_3 = left_set.project(x_k_2)

            lambda_k = (np.linalg.norm(x_k_1 - x_k_2, ord=2)**2) / (
                np.dot(x_k_1 - x_k_3, x_k_1 - x_k_2))
            x_k_4 = x_k_1 + lambda_k * (x_k_3 - x_k_1)
            
            if self.momentum is not None:
                iterate = heavy_ball_update(
                    iterates=iterates, velocity=x_k_4-x_k,
                    alpha=self.momentum['alpha'],
                    beta=self.momentum['beta'])
            iterates.append(x_k_4)

        return iterates, residuals, status 
Example 29
Project: tmx2map   Author: joshuaskelly   File: tmx2map.py    MIT License 5 votes vote down vote up
def xform_point(point):
                    result = numpy.dot(mat, (*point, 1))
                    return tuple(result.tolist()[:3]) 
Example 30
Project: tmx2map   Author: joshuaskelly   File: mathhelper.py    MIT License 5 votes vote down vote up
def vector_from_angle(degrees):
    """Returns a unit vector in the xy-plane rotated by the given degrees from
    the positive x-axis"""
    radians = math.radians(degrees)
    z_rot_matrix = numpy.identity(4)
    z_rot_matrix[0, 0] = math.cos(radians)
    z_rot_matrix[0, 1] = -math.sin(radians)
    z_rot_matrix[1, 0] = math.sin(radians)
    z_rot_matrix[1, 1] = math.cos(radians)

    return tuple(numpy.dot(z_rot_matrix, (1, 0, 0, 0))[:3]) 
Example 31
Project: explirefit   Author: codogogo   File: text_embeddings.py    Apache License 2.0 5 votes vote down vote up
def word_similarity(self, first_word, second_word, first_language = 'en', second_language = 'en'):	
		if self.do_cache:
			cache_str = min(first_word, second_word) + "-" + max(first_word, second_word)
			if (first_language + "-" + second_language) in self.cache and cache_str in self.cache[first_language + "-" + second_language]:
				return self.cache[first_language + "-" + second_language][cache_str]
		elif (first_word not in self.lang_vocabularies[first_language] and first_word.lower() not in self.lang_vocabularies[first_language]) or (second_word not in self.lang_vocabularies[second_language] and second_word.lower() not in self.lang_vocabularies[second_language]):
			if ((first_word in second_word or second_word in first_word) and first_language == second_language) or first_word == second_word:
					return 1
			else:
					return 0

		index_first = self.lang_vocabularies[first_language][first_word] if first_word in self.lang_vocabularies[first_language] else (self.lang_vocabularies[first_language][first_word.lower()] if first_word.lower() in self.lang_vocabularies[first_language] else -1)
		index_second = self.lang_vocabularies[second_language][second_word] if second_word in self.lang_vocabularies[second_language] else (self.lang_vocabularies[second_language][second_word.lower()] if second_word.lower() in self.lang_vocabularies[second_language] else -1)		

		if index_first >= 0 and index_second >= 0:		
			first_emb = self.lang_embeddings[first_language][index_first]
			second_emb = self.lang_embeddings[second_language][index_second] 

			first_norm = self.lang_emb_norms[first_language][index_first]
			second_norm = self.lang_emb_norms[second_language][index_second]

			score =  np.dot(first_emb, second_emb) / (first_norm * second_norm)
		else:
			score = 0
		
		if self.do_cache:
			if (first_language + "-" + second_language) not in self.cache:
				self.cache[first_language + "-" + second_language] = {}
				if cache_str not in self.cache[first_language + "-" + second_language]:
					self.cache[first_language + "-" + second_language][cache_str] = score		
		return score 
Example 32
Project: explirefit   Author: codogogo   File: simple_stats.py    Apache License 2.0 5 votes vote down vote up
def cosine(vec1, vec2):
	return np.dot(vec1, vec2) / (np.linalg.norm(vec1) * np.linalg.norm(vec2)) 
Example 33
Project: cs207-FinalProject   Author: PYNE-AD   File: AutoDiff.py    MIT License 5 votes vote down vote up
def _calcDerivative(self, der_value, k):
		if k != 0:
			jacobian = (np.array([self.jacobian[0]]))
			der = self._convertNonArray(der_value, k)
			return np.dot(der, jacobian) * self.jacobian**0
		else:
			return self._convertNonArray(der_value, k) 
Example 34
Project: xrft   Author: xgcm   File: test_xrft.py    MIT License 5 votes vote down vote up
def numpy_detrend(da):
    """
    Detrend a 2D field by subtracting out the least-square plane fit.

    Parameters
    ----------
    da : `numpy.array`
        The data to be detrended

    Returns
    -------
    da : `numpy.array`
        The detrended input data
    """
    N = da.shape

    G = np.ones((N[0]*N[1],3))
    for i in range(N[0]):
        G[N[1]*i:N[1]*i+N[1], 1] = i+1
        G[N[1]*i:N[1]*i+N[1], 2] = np.arange(1, N[1]+1)

    d_obs = np.reshape(da.copy(), (N[0]*N[1],1))
    m_est = np.dot(np.dot(spl.inv(np.dot(G.T, G)), G.T), d_obs)
    d_est = np.dot(G, m_est)

    lin_trend = np.reshape(d_est, N)

    return da - lin_trend 
Example 35
Project: FRIDA   Author: LCAV   File: point_cloud.py    MIT License 5 votes vote down vote up
def EDM(self):
        ''' Computes the EDM corresponding to the marker set '''
        if self.X is None:
            raise ValueError('No marker set')

        G = np.dot(self.X.T, self.X)
        return np.outer(np.ones(self.m), np.diag(G)) \
            - 2*G + np.outer(np.diag(G), np.ones(self.m)) 
Example 36
Project: FRIDA   Author: LCAV   File: point_cloud.py    MIT License 5 votes vote down vote up
def align(self, marker, axis):
        '''
        Rotate the marker set around the given axis until it is aligned onto the given marker

        Parameters
        ----------
        marker : int or str
            the index or label of the marker onto which to align the set
        axis : int
            the axis around which the rotation happens
        '''

        index = self.key2ind(marker)
        axis = ['x','y','z'].index(axis) if isinstance(marker, (str, unicode)) else axis

        # swap the axis around which to rotate to last position
        Y = self.X
        if self.dim == 3:
            Y[axis,:], Y[2,:] = Y[2,:], Y[axis,:]

        # Rotate around z to align x-axis to second point
        theta = np.arctan2(Y[1,index],Y[0,index])
        c = np.cos(theta)
        s = np.sin(theta)
        H = np.array([[c, s],[-s, c]])
        Y[:2,:] = np.dot(H,Y[:2,:])

        if self.dim == 3:
            Y[axis,:], Y[2,:] = Y[2,:], Y[axis,:] 
Example 37
Project: FRIDA   Author: LCAV   File: generators.py    MIT License 5 votes vote down vote up
def gen_sig_at_mic(sigmak2_k, phi_k, pos_mic_x,
                   pos_mic_y, omega_band, sound_speed,
                   SNR, Ns=256):
    """
    generate complex base-band signal received at microphones
    :param sigmak2_k: the variance of the circulant complex Gaussian signal
                emitted by the K sources
    :param phi_k: source locations (azimuths)
    :param pos_mic_x: a vector that contains microphones' x coordinates
    :param pos_mic_y: a vector that contains microphones' y coordinates
    :param omega_band: mid-band (ANGULAR) frequency [radian/sec]
    :param sound_speed: speed of sound
    :param SNR: SNR for the received signal at microphones
    :param Ns: number of snapshots used to estimate the covariance matrix
    :return: y_mic: received (complex) signal at microphones
    """
    num_mic = pos_mic_x.size
    xk, yk = polar2cart(1, phi_k)  # source locations in cartesian coordinates
    # reshape to use broadcasting
    xk = np.reshape(xk, (1, -1), order='F')
    yk = np.reshape(yk, (1, -1), order='F')
    pos_mic_x = np.reshape(pos_mic_x, (-1, 1), order='F')
    pos_mic_y = np.reshape(pos_mic_y, (-1, 1), order='F')

    t = np.reshape(np.linspace(0, 10 * np.pi, num=Ns), (1, -1), order='F')
    K = sigmak2_k.size
    sigmak2_k = np.reshape(sigmak2_k, (-1, 1), order='F')

    # x_tilde_k size: K x length_of_t
    # circular complex Gaussian process
    x_tilde_k = np.sqrt(sigmak2_k / 2.) * (np.random.randn(K, Ns) + 1j *
                                           np.random.randn(K, Ns))
    y_mic = np.dot(np.exp(-1j * (xk * pos_mic_x + yk * pos_mic_y) / (sound_speed / omega_band)),
                   x_tilde_k * np.exp(1j * omega_band * t))
    signal_energy = linalg.norm(y_mic, 'fro') ** 2
    noise_energy = signal_energy / 10 ** (SNR * 0.1)
    sigma2_noise = noise_energy / (Ns * num_mic)
    noise = np.sqrt(sigma2_noise / 2.) * (np.random.randn(*y_mic.shape) + 1j *
                                          np.random.randn(*y_mic.shape))
    y_mic_noisy = y_mic + noise
    return y_mic_noisy, y_mic 
Example 38
Project: FRIDA   Author: LCAV   File: tools_fri_doa_plane.py    MIT License 5 votes vote down vote up
def mtx_fri2visi_ri(M, p_mic_x, p_mic_y, D1, D2):
    """
    build the matrix that maps the Fourier series to the visibility in terms of
    REAL-VALUED entries only. (matrix size double)
    :param M: the Fourier series expansion is limited from -M to M
    :param p_mic_x: a vector that contains microphones x coordinates
    :param p_mic_y: a vector that contains microphones y coordinates
    :param D1: expansion matrix for the real-part
    :param D2: expansion matrix for the imaginary-part
    :return:
    """
    return np.dot(cpx_mtx2real(mtx_freq2visi(M, p_mic_x, p_mic_y)),
                  linalg.block_diag(D1, D2)) 
Example 39
Project: FRIDA   Author: LCAV   File: tools_fri_doa_plane.py    MIT License 5 votes vote down vote up
def Tmtx_ri_half_out_half(b_ri, K, D, L, D_coef, mtx_shrink):
    """
    if both b and annihilation filter coefficients are Hermitian symmetric, then
    the output will also be Hermitian symmetric => the effectively output is half the size
    """
    return np.dot(np.dot(mtx_shrink, Tmtx_ri(b_ri, K, D, L)), D_coef) 
Example 40
Project: FRIDA   Author: LCAV   File: tools_fri_doa_plane.py    MIT License 5 votes vote down vote up
def Rmtx_ri(coef_ri, K, D, L):
    coef_ri = np.squeeze(coef_ri)
    coef_r = coef_ri[:K + 1]
    coef_i = coef_ri[K + 1:]
    R_r = linalg.toeplitz(np.concatenate((np.array([coef_r[-1]]),
                                          np.zeros(L - K - 1))),
                          np.concatenate((coef_r[::-1],
                                          np.zeros(L - K - 1)))
                          )
    R_i = linalg.toeplitz(np.concatenate((np.array([coef_i[-1]]),
                                          np.zeros(L - K - 1))),
                          np.concatenate((coef_i[::-1],
                                          np.zeros(L - K - 1)))
                          )
    return np.dot(np.vstack((np.hstack((R_r, -R_i)), np.hstack((R_i, R_r)))), D) 
Example 41
Project: FRIDA   Author: LCAV   File: tools_fri_doa_plane.py    MIT License 5 votes vote down vote up
def Rmtx_ri_half_out_half(coef_half, K, D, L, D_coef, mtx_shrink):
    """
    if both b and annihilation filter coefficients are Hermitian symmetric, then
    the output will also be Hermitian symmetric => the effectively output is half the size
    """
    return np.dot(mtx_shrink, Rmtx_ri(np.dot(D_coef, coef_half), K, D, L)) 
Example 42
Project: FRIDA   Author: LCAV   File: tools_fri_doa_plane.py    MIT License 5 votes vote down vote up
def mtx_updated_G_multiband(phi_recon, M, mtx_amp2visi_ri,
                            mtx_fri2visi_ri, num_bands):
    """
    Update the linear transformation matrix that links the FRI sequence to the
    visibilities by using the reconstructed Dirac locations.
    :param phi_recon: the reconstructed Dirac locations (azimuths)
    :param M: the Fourier series expansion is between -M to M
    :param p_mic_x: a vector that contains microphones' x-coordinates
    :param p_mic_y: a vector that contains microphones' y-coordinates
    :param mtx_freq2visi: the linear mapping from Fourier series to visibilities
    :return:
    """
    L = 2 * M + 1
    ms_half = np.reshape(np.arange(-M, 1, step=1), (-1, 1), order='F')
    phi_recon = np.reshape(phi_recon, (1, -1), order='F')
    mtx_amp2freq = np.exp(-1j * ms_half * phi_recon)  # size: (M + 1) x K
    mtx_amp2freq_ri = np.vstack((mtx_amp2freq.real, mtx_amp2freq.imag[:-1, :]))  # size: (2M + 1) x K
    mtx_fri2amp_ri = linalg.lstsq(mtx_amp2freq_ri, np.eye(L))[0]
    # projection mtx_freq2visi to the null space of mtx_fri2amp
    mtx_null_proj = np.eye(L) - np.dot(mtx_fri2amp_ri.T,
                                       linalg.lstsq(mtx_fri2amp_ri.T, np.eye(L))[0])
    G_updated = np.dot(mtx_amp2visi_ri,
                       linalg.block_diag(*([mtx_fri2amp_ri] * num_bands))
                       ) + \
                np.dot(mtx_fri2visi_ri,
                       linalg.block_diag(*([mtx_null_proj] * num_bands))
                       )
    return G_updated 
Example 43
Project: FRIDA   Author: LCAV   File: tools_fri_doa_plane.py    MIT License 5 votes vote down vote up
def compute_b(G_lst, GtG_lst, beta_lst, Rc0, num_bands, a_ri):
    """
    compute the uniform sinusoidal samples b from the updated annihilating
    filter coeffiients.
    :param GtG_lst: list of G^H G for different subbands
    :param beta_lst: list of beta-s for different subbands
    :param Rc0: right-dual matrix, here it is the convolution matrix associated with c
    :param num_bands: number of bands
    :param L: size of b: L by 1
    :param a_ri: a 2D numpy array. each column corresponds to the measurements within a subband
    :return:
    """
    b_lst = []
    a_Gb_lst = []
    for loop in range(num_bands):
        GtG_loop = GtG_lst[loop]
        beta_loop = beta_lst[loop]
        b_loop = beta_loop - \
                 linalg.solve(GtG_loop,
                              np.dot(Rc0.T,
                                     linalg.solve(np.dot(Rc0, linalg.solve(GtG_loop, Rc0.T)),
                                                  np.dot(Rc0, beta_loop)))
                              )

        b_lst.append(b_loop)
        a_Gb_lst.append(a_ri[:, loop] - np.dot(G_lst[loop], b_loop))

    return np.column_stack(b_lst), linalg.norm(np.concatenate(a_Gb_lst)) 
Example 44
Project: DataHack2018   Author: InnovizTech   File: vis_utils.py    BSD 3-Clause "New" or "Revised" License 5 votes vote down vote up
def euler_matrix(yaw, pitch, roll):
    # angles in radians
    # negating pitch for easier computation
    pi=np.pi+1e-7
    assert -pi <= yaw <= pi and -np.pi <= roll <= np.pi and -np.pi/2 <= pitch <= np.pi/2, \
        "Erroneous yaw, pitch, roll={},{},{}".format(yaw, pitch, roll)
    rotX = np.eye(3)
    rotY = np.eye(3)
    rotZ = np.eye(3)
    rotX[1:, 1:] = rot_mat_2d(roll)
    rotY[::2, ::2] = rot_mat_2d(-pitch)
    rotZ[:2, :2] = rot_mat_2d(yaw)

    return rotX.dot(rotY.dot(rotZ)) 
Example 45
Project: DataHack2018   Author: InnovizTech   File: vis_utils.py    BSD 3-Clause "New" or "Revised" License 5 votes vote down vote up
def _pan_camera(self, movement):
        camera_pos = base.camera.getPos()
        hpr = base.camera.getHpr()
        rot_mat = lambda angle: np.array([[np.cos(angle), -np.sin(angle)],
                                          [np.sin(angle), np.cos(angle)]])
        movement_vec_xy = np.dot(rot_mat(hpr[0] / 180. * np.pi), np.array([movement[0], movement[1]]))
        base.camera.setPos(camera_pos[0] + movement_vec_xy[0],
                           camera_pos[1] + movement_vec_xy[1],
                           camera_pos[2] + movement[2])
        self._update_location_text() 
Example 46
Project: DataHack2018   Author: InnovizTech   File: vis_utils.py    BSD 3-Clause "New" or "Revised" License 5 votes vote down vote up
def _pan_camera_mouse(self, movement):
        camera_transform = np.array(base.camera.getNetTransform().getMat())
        pan_vec = np.array([-movement[0], 0, -movement[1], 1]).dot(camera_transform)
        base.camera.setPos(pan_vec[0],
                           pan_vec[1],
                           pan_vec[2])
        self._update_location_text() 
Example 47
Project: DataHack2018   Author: InnovizTech   File: math_utils.py    BSD 3-Clause "New" or "Revised" License 5 votes vote down vote up
def inverse(self):
        '''
        :return: The inverse transformation
        '''
        _R = self.rotation.T
        _t = -np.dot(self.rotation.T, self.translation)
        return RotationTranslationData(rt=(_R, _t)) 
Example 48
Project: DataHack2018   Author: InnovizTech   File: math_utils.py    BSD 3-Clause "New" or "Revised" License 5 votes vote down vote up
def euler_matrix(yaw, pitch, roll):
    # angles in radians
    # negating pitch for easier computation
    pi=np.pi+1e-7
    assert -pi <= yaw <= pi and -np.pi <= roll <= np.pi and -np.pi/2 <= pitch <= np.pi/2, \
        "Erroneous yaw, pitch, roll={},{},{}".format(yaw, pitch, roll)
    rotX = np.eye(3)
    rotY = np.eye(3)
    rotZ = np.eye(3)
    rotX[1:, 1:] = rot_mat_2d(roll)
    rotY[::2, ::2] = rot_mat_2d(pitch)
    rotZ[:2, :2] = rot_mat_2d(yaw)

    return rotZ.dot(rotY.dot(rotX)) 
Example 49
Project: Automated-Social-Annotation   Author: acadTags   File: data_util.py    MIT License 5 votes vote down vote up
def get_label_sim_matrix(vocabulary_index2word_label,word2vec_model_label_path='../tag-all.bin-300',name_scope='',threshold=0):
    cache_path ='../cache_vocabulary_label_pik/'+ name_scope + "_label_sim_" + str(threshold) + ".pik"
    print("cache_path:",cache_path,"file_exists:",os.path.exists(cache_path))
    if os.path.exists(cache_path):
        with open(cache_path, 'rb') as data_f:
            result=pickle.load(data_f)
            return result
    else:
        model=word2vec.load(word2vec_model_label_path,kind='bin')
        #m = model.vectors.shape[0]-1 #length # the first one is </s>, to be eliminated
        m = len(vocabulary_index2word_label)
        result=np.zeros((m,m))
        count_less_th = 0.0 # count the sim less than the threshold
        for i in range(0,m):
            for j in range(0,m):
                vector_i=model.get_vector(vocabulary_index2word_label[i])
                vector_j=model.get_vector(vocabulary_index2word_label[j])
                #result[i][j] = np.dot(vector_i,vector_j.T) # can be negative here, result in [-1,1]
                result[i][j] = (1+np.dot(vector_i,vector_j.T))/2 # result in [0,1]
                if result[i][j] < threshold:
                    count_less_th = count_less_th + 1
                    result[i][j] = 0
        print("result",result)
        print("result",result.shape)
        print("retained similarities percentage:", str(1-count_less_th/float(m)/float(m)))
        #save to file system if vocabulary of words is not exists.
        if not os.path.exists(cache_path):
            with open(cache_path, 'ab') as data_f:
                pickle.dump(result, data_f)
    return result

# used for other embedding 
Example 50
Project: Automated-Social-Annotation   Author: acadTags   File: data_util.py    MIT License 5 votes vote down vote up
def get_label_sim_matrix(vocabulary_index2word_label,word2vec_model_label_path='../tag-all.bin-300',name_scope='',threshold=0):
    cache_path ='../cache_vocabulary_label_pik/'+ name_scope + "_label_sim_" + str(threshold) + ".pik"
    print("cache_path:",cache_path,"file_exists:",os.path.exists(cache_path))
    if os.path.exists(cache_path):
        with open(cache_path, 'rb') as data_f:
            result=pickle.load(data_f)
            return result
    else:
        model=word2vec.load(word2vec_model_label_path,kind='bin')
        #m = model.vectors.shape[0]-1 #length # the first one is </s>, to be eliminated
        m = len(vocabulary_index2word_label)
        result=np.zeros((m,m))
        count_less_th = 0.0 # count the sim less than the threshold
        for i in range(0,m):
            for j in range(0,m):
                vector_i=model.get_vector(vocabulary_index2word_label[i])
                vector_j=model.get_vector(vocabulary_index2word_label[j])
                #result[i][j] = np.dot(vector_i,vector_j.T) # can be negative here, result in [-1,1]
                result[i][j] = (1+np.dot(vector_i,vector_j.T))/2 # result in [0,1]
                if result[i][j] < threshold:
                    count_less_th = count_less_th + 1
                    result[i][j] = 0
        print("result",result)
        print("result",result.shape)
        print("retained similarities percentage:", str(1-count_less_th/float(m)/float(m)))
        #save to file system if vocabulary of words is not exists.
        if not os.path.exists(cache_path):
            with open(cache_path, 'ab') as data_f:
                pickle.dump(result, data_f)
    return result