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 
Example 51
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 
Example 52
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 53
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 
Example 54
Project: ml_news_popularity   Author: khaledJabr   File: lr.py    MIT License 5 votes vote down vote up
def predict(self, data):
        data = self.std_scaler.transform(data)
        predicted = data.dot(self.betas)
        return predicted

        # norm_data = np.zeros((len(data), len(data[0])))
        # # normalize features
        # for j in range(len(data[0])):
        #     norm_data[:,j] = self.normalize(data[:,j], self.beta_mins[j], self.beta_maxes[j])

        # norm_predicted = np.dot(self.betas, norm_data.T)
        # return self.denormalize_zero_max(norm_predicted, self.max_target)

    # def mse(self):
    #     return mean_squared_error(self.targets, self.predicted) 
Example 55
Project: Adversarial-Face-Attack   Author: ppwwyyxx   File: face_attack.py    GNU General Public License v3.0 5 votes vote down vote up
def distance_to_victim(self, img):
        emb = self.eval_embeddings([img])
        dist = np.dot(emb, self.victim_embeddings.T).flatten()
        stats = np.percentile(dist, [10, 30, 50, 70, 90])
        return stats 
Example 56
Project: RandomFourierFeatures   Author: tiskw   File: PyRFF.py    MIT License 5 votes vote down vote up
def conv(self, X):
        ts = np.dot(X, self.W)
        cs = np.cos(ts)
        ss = np.sin(ts)
        return np.bmat([cs, ss]) 
Example 57
Project: RandomFourierFeatures   Author: tiskw   File: PyRFF.py    MIT License 5 votes vote down vote up
def conv(self, X):
        ts = np.dot(X, self.W)
        cs = np.cos(ts)
        ss = np.sin(ts)
        return np.bmat([cs, ss]) 
Example 58
Project: RandomFourierFeatures   Author: tiskw   File: PyRFF.py    MIT License 5 votes vote down vote up
def predict(self, X, **args):
        svc = RFFSVC(self.dim, self.std, self.W, **args)
        return np.argmax(np.dot(svc.conv(X), self.coef.T) + self.icpt.flatten(), axis = 1) 
Example 59
Project: StructEngPy   Author: zhuoju36   File: dynamic.py    MIT License 5 votes vote down vote up
def Newmark_beta(model:Model,T,F,u0,v0,a0,beta=0.25,gamma=0.5):
    """
    beta,gamma: parameters.\n
    u0,v0,a0: initial state.\n
    T: time list with uniform interval.\n
    F: list of time-dependent force vectors.
    """
    dt=T[1]-T[0]
    b1=1/(beta*dt*dt)
    b2=-1/(beta*dt)
    b3=1/2/beta-1
    b4=gamma*dt*b1
    b5=1+gamma*dt*b2
    b6=dt*(1+gamma*b3-gamma)
    K_=self.K+b4*self.C+b1*self.M 
    u=[u0]
    v=[v0]
    a=[a0]
    tt=[0]
    for (t,ft) in zip(T,F):
        ft_=ft+np.dot(self.M,b1*u[-1]-b2*v[-1]-b3*a[-1])+np.dot(self.C,b4*u[-1]-b5*v[-1]-b6*a[-1])
        ut=np.linalg.solve(K_,ft_)
        vt=b4*(ut-u[-1])+b5*v[-1]+b6*a[-1]
        at=b1*(ut-u[-1])+b2*v[-1]+b3*a[-1]
        u.append(ut)
        v.append(vt)
        a.append(at)
        tt.append(t)
    df=pd.DataFrame({'t':tt,'u':u,'v':v,'a':a})
    return df 
Example 60
Project: StructEngPy   Author: zhuoju36   File: dynamic.py    MIT License 5 votes vote down vote up
def Wilson_theta(model:Model,T,F,u0=0,v0=0,a0=0,beta=0.25,gamma=0.5,theta=1.4):
    """
    beta,gamma,theta: parameters.\n
    u0,v0,a0: initial state.\n
    T: time list with uniform interval.\n
    F: list of time-dependent force vectors.
    """
    dt=T[1]-T[0]
    dt_=theta*dt
    b1=1/(beta*dt_**2)
    b2=-1/(beta*dt_)
    b3=(1/2-beta)/beta
    b4=gamma*dt_*b1
    b5=1+gamma*dt_*b2
    b6=dt_*(1+gamma*b3-gamma)
    K_=self.K+b4*self.C+b1*self.M 
    u=[u0]
    v=[v0]
    a=[a0]
    tt=[0]

    R_=[F[0]]
    for i in range(len(F)-1):
        R_.append(F[i]+theta*(F[i+1]-F[i]))
        
    for (t,ft) in zip(T,R_):
        ft_=ft+np.dot(self.M,b1*u[-1]-b2*v[-1]-b3*a[-1])+np.dot(self.C,b4*u[-1]-b5*v[-1]-b6*a[-1])
        ut_=np.linalg.solve(K_,ft_)
        vt_=b4*(ut_-u[-1])+b5*v[-1]+b6*a[-1]
        at_=b1*(ut_-u[-1])+b2*v[-1]+b3*a[-1]          
        
        at=a[-1]+1/theta*(at_-a[-1])
        vt=v[-1]+((1-gamma)*a[-1]+gamma*at)*dt
        ut=u[-1]+v[-1]*dt+(1/2-beta)*a[-1]*dt**2+beta*at*dt**2

        u.append(ut)
        v.append(vt)
        a.append(at)
        tt.append(t)
    df=pd.DataFrame({'t':tt,'u':u,'v':v,'a':a})
    return df 
Example 61
Project: StructEngPy   Author: zhuoju36   File: __init__.py    MIT License 5 votes vote down vote up
def resolve_node_disp(self,node_id):
        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.d_[node_id*6:node_id*6+6]).reshape(6)
        else:
            raise Exception("The node doesn't exists.") 
Example 62
Project: StructEngPy   Author: zhuoju36   File: __init__.py    MIT License 5 votes vote down vote up
def resolve_node_reaction(self,node_id):
        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.r_[node_id*6:node_id*6+6,0]).reshape(6)
        else:
            raise Exception("The node doesn't exists.") 
Example 63
Project: StructEngPy   Author: zhuoju36   File: __init__.py    MIT License 5 votes vote down vote up
def resolve_beam_force(self,beam_id):
        if not self.is_solved:
            raise Exception('The model has to be solved first.')
        if beam_id in self.__beams.keys():
            beam=self.__beams[beam_id]
            i=beam.nodes[0].hid
            j=beam.nodes[1].hid
            T=beam.transform_matrix
            ue=np.vstack([
                        self.d_[i*6:i*6+6],
                        self.d_[j*6:j*6+6]
                        ])   
            return (beam.Ke_.dot(T.dot(ue))+beam.re_).reshape(12)
        else:
            raise Exception("The element doesn't exists.") 
Example 64
Project: StructEngPy   Author: zhuoju36   File: element.py    MIT License 5 votes vote down vote up
def __init__(self,node_i,node_j,node_k,t,E,mu,rho,dof,name=None,tol=1e-6):
        super(Tri,self).__init__(2,dof,name)
        self._nodes=[node_i,node_j,node_k]
        #Initialize local CSys
        o=[(node_i.x+node_j.x+node_k.x)/3,
            (node_i.y+node_j.y+node_k.y)/3,
            (node_i.z+node_j.z+node_k.z)/3]
        pt1 = [ node_j.x, node_j.y, node_j.z ]
        pt2 = [ node_i.x, node_i.y, node_i.z ]
        self._local_csys = Cartisian(o, pt1, pt2) 

        self._area=0.5*np.linalg.det(np.array([[1,1,1],
                                    [node_j.x-node_i.x,node_j.y-node_i.y,node_j.z-node_i.z],
                                    [node_k.x-node_i.x,node_k.y-node_i.y,node_k.z-node_i.z]]))
        self._t=t
        self._E=E
        self._mu=mu

        E=self._E
        mu=self._mu
        D0=E/(1-mu**2)
        self._D=np.array([[1,mu,0],
                    [mu,1,0],
                    [0,0,(1-mu)/2]])*D0
        #3D to local 2D
        V=self._local_csys.transform_matrix
        x3D=np.array([[node_i.x,node_i.y,node_i.z],
                    [node_j.x,node_j.y,node_j.z],
                    [node_k.x,node_k.y,node_k.z]])
        x2D=np.dot(x3D,V.T)
        self._x2D=x2D[:,:2] 
Example 65
Project: StructEngPy   Author: zhuoju36   File: element.py    MIT License 5 votes vote down vote up
def _N(self,x):
        """
        interpolate function.
        return: 3x1 array represent x,y
        """
        x,y=x[0],x[1]
        L=np.array((3,1))
        L[0]=self._abc(1,2).dot(np.array([1,x,y]))/2/self._area
        L[1]=self._abc(2,0).dot(np.array([1,x,y]))/2/self._area
        L[2]=self._abc(0,1).dot(np.array([1,x,y]))/2/self._area
        return L.reshape(3,1) 
Example 66
Project: StructEngPy   Author: zhuoju36   File: element.py    MIT License 5 votes vote down vote up
def _x(self,L):
        """
        convert csys from L to x
        return: 2x1 array represent x,y
        """
        return np.dot(np.array(L).reshape(1,3),self._x0).reshape(2,1) 
Example 67
Project: StructEngPy   Author: zhuoju36   File: element.py    MIT License 5 votes vote down vote up
def _S(self,s,r):
        """
        stress matrix
        """
        return np.dot(self._D,self._B(s,r)) 
Example 68
Project: b2ac   Author: hbldh   File: matrix_algorithms.py    MIT License 4 votes vote down vote up
def convert_to_Hessenberg_double(A):
    """ Tridiagonalize a square matrix to upper Hessenberg form using Householder reflections.

    Costs 10/3 * n^3 + O(n^2) operations, where n is the size of the matrix.

    .. code-block:: matlab

        function A = mytridiag(A)

        [m,n] = size(A);
        if (m ~= n)
            error('Input matrix is not square.')
        end

        for k = 1:(m - 2)
            vk = A((k+1):m,k);
            vk(1) = vk(1) + sign(vk(1)) * norm(vk);
            vk = vk / norm(vk);
            A((k+1):m,k:m) = A((k+1):m,k:m) - ...
                                2 * vk * (vk' * A((k+1):m,k:m));
            A(1:m,(k+1):m) = A(1:m,(k+1):m) - ...
                               2 * (A(1:m,(k+1):m) * vk) * vk';
        end

        end

    :param A: The matrix to convert.
    :type A: :py:class:`numpy.ndarray`
    :return: The Hessenberg matrix.
    :rtype: :py:class:`numpy.ndarray`

    """
    m, n = A.shape
    A = np.array(A, 'float')

    for k in xrange(m - 1):
        vk = A[(k + 1):, k].copy()
        vk[0] += np.sign(vk[0]) * np.linalg.norm(vk, 2)
        vk /= np.linalg.norm(vk, 2)
        A[(k + 1):, k:] -= 2 * np.outer(vk, np.dot(vk, A[(k + 1):, k:]))
        A[:, (k + 1):] -= 2 * np.outer(np.dot(A[:, (k + 1):], vk), vk)

    return np.triu(A, -1) 
Example 69
Project: b2ac   Author: hbldh   File: matrix_algorithms.py    MIT License 4 votes vote down vote up
def convert_to_Hessenberg_symmetric_double(A):
    """Tridiagonalize a symmetric square matric to Hessenberg form using Householder reflections.

    Costs 4/3 * n^3 + O(n^2) operations.

    .. code-block:: matlab

        function A = mytridiagturbo(A)

        [m,n] = size(A);
        if (m ~= n)
            error('Input matrix is not square.')
        end

        for k = 1:(m - 2)
            vk = A((k+1):m,k);
            beta = sign(vk(1))*norm(vk);
            vk(1) = vk(1) + beta;
            vk = vk / norm(vk);

            % Calculate A_hat for current iteration.
            u = -2*(A((k+1:m),(k+1:m)) * vk);
            d = -vk' * u;
            w = u + d * vk;
            wv = w * vk';
            A((k+1:m),(k+1:m)) = A((k+1:m),(k+1:m)) + wv + wv';

            % Adjust tridiagonality.
            A((k+1):m,k) = [-beta; zeros(m-k-1,1)];
            A(k,(k+1):m) = A((k+1):m,k)';
         end

    :param A: The matrix to convert.
    :type A: :py:class:`numpy.ndarray`
    :return: The Hessenberg matrix.
    :rtype: :py:class:`numpy.ndarray`

    """

    m, n = A.shape
    A = np.array(A, 'float')

    for k in xrange(m - 1):
        vk = A[(k + 1):, k].copy()
        beta = np.sign(vk[0]) * np.linalg.norm(vk)
        vk[0] += beta
        vk = vk / np.linalg.norm(vk)

        # Calculate A_hat for current iteration.
        u = -2 * np.dot(A[(k + 1):, (k + 1):], vk)
        d = np.dot(-vk, u)
        w = u + d * vk
        wv = np.outer(w, vk)
        A[(k + 1):, (k + 1):] += wv + wv.T

        # Adjust tridiagonality.
        t = np.zeros((m - (k + 1), ), dtype='float')
        t[0] = -beta
        A[(k + 1):, k] = t
        A[k, (k + 1):] = t
    return A 
Example 70
Project: b2ac   Author: hbldh   File: reference.py    MIT License 4 votes vote down vote up
def fit_improved_B2AC_numpy(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 simply applies NumPy:s methods for calculating
    the conic section, modelled after the Matlab code in the paper:

    .. code-block::

        function a = fit_ellipse(x, y)

        D1 = [x .ˆ 2, x .* y, y .ˆ 2]; % quadratic part of the design matrix
        D2 = [x, y, ones(size(x))]; % linear part of the design matrix
        S1 = D1’ * D1; % quadratic part of the scatter matrix
        S2 = D1’ * D2; % combined part of the scatter matrix
        S3 = D2’ * D2; % linear part of the scatter matrix
        T = - inv(S3) * S2’; % for getting a2 from a1
        M = S1 + S2 * T; % reduced scatter matrix
        M = [M(3, :) ./ 2; - M(2, :); M(1, :) ./ 2]; % premultiply by inv(C1)
        [evec, eval] = eig(M); % solve eigensystem
        cond = 4 * evec(1, :) .* evec(3, :) - evec(2, :) .ˆ 2; % evaluate a’Ca
        a1 = evec(:, find(cond > 0)); % eigenvector for min. pos. eigenvalue
        a = [a1; T * a1]; % ellipse coefficients

    :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`

    """
    x = points[:, 0]
    y = points[:, 1]

    D1 = np.vstack([x ** 2, x * y, y ** 2]).T
    D2 = np.vstack([x, y, np.ones((len(x), ), dtype=x.dtype)]).T
    S1 = D1.T.dot(D1)
    S2 = D1.T.dot(D2)
    S3 = D2.T.dot(D2)
    T = -np.linalg.inv(S3).dot(S2.T)
    M = S1 + S2.dot(T)
    M = np.array([M[2, :] / 2, -M[1, :], M[0, :] / 2])
    eval, evec = np.linalg.eig(M)
    cond = (4 * evec[:, 0] * evec[:, 2]) - (evec[:, 1] ** 2)
    I = np.where(cond > 0)[0]
    a1 = evec[:, I[np.argmin(cond[I])]]
    return np.concatenate([a1, T.dot(a1)]) 
Example 71
Project: b2ac   Author: hbldh   File: overlap_functions.py    MIT License 4 votes vote down vote up
def quickhull(sample):
    """Calculates the convex hull of an arbitrary 2D point cloud.

    This is a pure Python version of the Quick Hull algorithm.
    It's based on the version of ``literateprograms``, but fixes some
    old-style Numeric function calls.

    This version works with numpy version > 1.2.1

    References:

    * `Literateprograms <http://en.literateprograms.org/Quickhull_(Python,_arrays)>`_
    * `Wikipedia <http://en.wikipedia.org/wiki/QuickHull>`_

    Code adapted from:

    `<http://members.home.nl/wim.h.bakker/python/quickhull2d.py>`_

    :param sample: Points to which the convex hull is desired to be found.
    :type sample: :py:class:`numpy.ndarray`
    :return: The convex hull of the points.
    :rtype: :py:class:`numpy.ndarray`

    """

    def calculate_convex_hull(sample):
        link = lambda a, b: np.concatenate((a, b[1:]))
        edge = lambda a, b: np.concatenate(([a], [b]))

        def dome(sample, base):
            h, t = base
            dists = np.dot(sample-h, np.dot(((0, -1), (1, 0)), (t - h)))
            outer = np.repeat(sample, dists > 0, axis=0)

            if len(outer):
                pivot = sample[np.argmax(dists)]
                return link(dome(outer, edge(h, pivot)),
                            dome(outer, edge(pivot, t)))
            else:
                return base

        if len(sample) > 2:
            axis = sample[:, 0]
            base = np.take(sample, [np.argmin(axis), np.argmax(axis)], axis=0)
            return link(dome(sample, base),
                        dome(sample, base[::-1]))
        else:
            return sample

    # Perform a reversal of points here to get points ordered clockwise instead of
    # counter clockwise that the QuickHull above returns.
    return calculate_convex_hull(sample)[::-1, :] 
Example 72
Project: xrft   Author: xgcm   File: xrft.py    MIT License 4 votes vote down vote up
def detrendn(da, axes=None):
    """
    Detrend by subtracting out the least-square plane or least-square cubic fit
    depending on the number of axis.

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

    Returns
    -------
    da : `numpy.array`
        The detrended input data
    """
    N = [da.shape[n] for n in axes]
    M = []
    for n in range(da.ndim):
        if n not in axes:
            M.append(da.shape[n])

    if len(N) == 2:
        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)
        if type(da) == xr.DataArray:
            d_obs = np.reshape(da.copy().values, (N[0]*N[1],1))
        else:
            d_obs = np.reshape(da.copy(), (N[0]*N[1],1))
    elif len(N) == 3:
        if type(da) == xr.DataArray:
            if da.ndim > 3:
                raise NotImplementedError("Cubic detrend is not implemented "
                                         "for 4-dimensional `xarray.DataArray`."
                                         " We suggest converting it to "
                                         "`dask.array`.")
            else:
                d_obs = np.reshape(da.copy().values, (N[0]*N[1]*N[2],1))
        else:
            d_obs = np.reshape(da.copy(), (N[0]*N[1]*N[2],1))

        G = np.ones((N[0]*N[1]*N[2],4))
        G[:,3] = np.tile(np.arange(1,N[2]+1), N[0]*N[1])
        ys = np.zeros(N[1]*N[2])
        for i in range(N[1]):
            ys[N[2]*i:N[2]*i+N[2]] = i+1
        G[:,2] = np.tile(ys, N[0])
        for i in range(N[0]):
            G[len(ys)*i:len(ys)*i+len(ys),1] = i+1
    else:
        raise NotImplementedError("Detrending over more than 4 axes is "
                                 "not implemented.")

    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, da.shape)

    return da - lin_trend 
Example 73
Project: comet-commonsense   Author: atcbosselut   File: demo_bilinear.py    Apache License 2.0 4 votes vote down vote up
def score(term1, term2, words, We, rel, Rel, Weight, Offset, evaType):
    v1 = getVec(We, words, term1)
    v2 = getVec(We, words, term2)
    result = {}

    for k, v in rel.items():
        v_r = Rel[rel[k], :]
        gv1 = np.tanh(np.dot(v1, Weight)+Offset)
        gv2 = np.tanh(np.dot(v2, Weight)+Offset)

        temp1 = np.dot(gv1,  v_r)
        score = np.inner(temp1, gv2)
        result[k] = (sigmoid(score))

    if(evaType.lower() == 'max'):
        result = sorted(result.items(),  key=lambda x: x[1],  reverse = True)
        for k, v in result[:1]:
            print k,  'score:',  v
        return result[:1]
    if(evaType.lower() == 'topfive'):
        result = sorted(result.items(),  key=lambda x: x[1],  reverse = True)
        for k, v in result[:5]:
            print k,  'score:',  v
        return result[:5]
    if(evaType.lower() == 'sum'):
        result = sorted(result.items(),  key=lambda x: x[1],  reverse = True)
        total = 0
        for i in result:
            total = total + i[1]
        print 'total score is:', total
        return total
    if(evaType.lower() == 'all'):
        result = sorted(result.items(),  key=lambda x: x[1],  reverse = True)
        for k, v in result[:]:
            print k,  'score:',  v
        return result
    else:
        tar_rel = evaType.lower()
        if result.get(tar_rel) == None:
            print 'illegal relation,  please re-enter a valid relation'
            return 'None'
        else:
            return result.get(tar_rel) 
Example 74
Project: FRIDA   Author: LCAV   File: point_cloud.py    MIT License 4 votes vote down vote up
def normalize(self, refs=None):
        '''
        Reposition points such that x0 is at origin, x1 lies on c-axis
        and x2 lies above x-axis, keeping the relative position to each other.
        The z-axis is defined according to right hand rule by default.

        Parameters:
        -----------
        refs : list of 3 ints or str
            The index or label of three markers used to define (origin, x-axis, y-axis)
        left_hand : bool, optional (default False)
            Normally the z-axis is defined using right-hand rule, this flag allows to override this behavior
        '''

        if refs is None:
            refs = [0,1,2,3]

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

        if self.dim == 2 and len(refs) < 3:
            raise ValueError('In 2D three reference points are needed to define a reference frame')
        elif self.dim == 3 and len(refs) < 4:
            raise ValueError('In 3D four reference points are needed to define a reference frame')

        # set first point to origin
        X0 = self.X[:,refs[0],None]
        Y = self.X - X0

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

        if self.dim == 2:
            # set third point to lie above x-axis
            if Y[1,refs[2]] < 0:
                Y[1,:] *= -1

        elif self.dim == 3:
            # In 3D we also want to make sur z-axis points up
            theta = np.arctan2(Y[2,refs[2]],Y[1,refs[2]])
            c = np.cos(theta)
            s = np.sin(theta)
            H = np.array([[1, 0, 0], [0, c, s],[0, -s, c]])
            Y = np.dot(H,Y)

        # Flip the z-axis if requested
        if self.dim == 3 and Y[2,refs[3]] < 0:
            Y[2,:] *= -1

        self.X = Y 
Example 75
Project: ieml   Author: IEMLdev   File: tree_graph.py    GNU General Public License v3.0 4 votes vote down vote up
def __init__(self, list_transitions):
        """
        Transitions list must be the (start, end, data) the data will be stored as the transition tag
        :param list_transitions:
        """
        # transitions : dict
        #
        self.transitions = defaultdict(list)
        for t in list_transitions:
            self.transitions[t[0]].append((t[1], t))

        self.nodes = sorted(set(self.transitions) | {e[0] for l in self.transitions.values() for e in l})

        # sort the transitions
        for s in self.transitions:
            self.transitions[s].sort(key=lambda t: self.nodes.index(t[0]))

        self.nodes_index = {n: i for i, n in enumerate(self.nodes)}
        _count = len(self.nodes)
        self.array = numpy.zeros((len(self.nodes), len(self.nodes)), dtype=bool)

        for t in self.transitions:
            for end in self.transitions[t]:
                self.array[self.nodes_index[t]][self.nodes_index[end[0]]] = True

        # checking
        # root checking, no_parent hold True for each index where the node has no parent
        parents_count = numpy.dot(self.array.transpose().astype(dtype=int), numpy.ones((_count,), dtype=int))
        no_parents = parents_count == 0
        roots_count = numpy.count_nonzero(no_parents)

        if roots_count == 0:
            raise InvalidTreeStructure('No root node found, the graph has at least a cycle.')
        elif roots_count > 1:
            raise InvalidTreeStructure('Several root nodes found.')

        self.root = self.nodes[no_parents.nonzero()[0][0]]

        if (parents_count > 1).any():
            raise InvalidTreeStructure('A node has several parents.')

        def __stage():
            current = [self.root]
            while current:
                yield current
                current = [child[0] for parent in current for child in self.transitions[parent]]

        self.stages = list(__stage()) 
Example 76
Project: ml_news_popularity   Author: khaledJabr   File: lr.py    MIT License 4 votes vote down vote up
def fit(self, data, targets):
        self.training = data
        self.targets = targets

        self.std_scaler = preprocessing.StandardScaler().fit(data)
        data = self.std_scaler.transform(data)
        targets = preprocessing.scale(targets, with_mean=True, with_std=False)

        # self.max_target = np.max(targets)
        # self.beta_mins = [np.min(data[:,j]) for j in range(len(data[0]))]
        # # print(self.beta_mins)
        # self.beta_maxes = [np.max(data[:,j]) for j in range(len(data[0]))]
        # # print(self.beta_maxes)
        # self.norm_targets = self.normalize_zero_max(self.targets, self.max_target)

        # norm = np.zeros((len(data), len(data[0])))
        # # normalize features
        # for j in range(len(data[0])):
        #     norm[:,j] = self.normalize(data[:,j], self.beta_mins[j], self.beta_maxes[j])
        
        # initialize weights
        self.betas = np.zeros(len(data[0]))
        # norm = normalize(self.training)
        # norm = self.training
        # iterate till not converged
        for iter in range(self.max_iter):
            # print("=== Iter {}".format(iter))
            # iterate over all features (j=0,1...M)            
            for j in range(len(self.betas)):
                # remove feature j to determine effect on residuals
                # b_no_j = self.betas.copy()
                # b_no_j[j] = 0
                # x_no_j = self.training
                self.betas[j] = self.update_beta(data, targets, self.betas, j)
                # x_no_j = np.delete(self.training, j, 1)
                # b_no_j = np.delete(self.betas, j, 0)
                # # predict_no_j = np.dot(b_no_j, x_no_j.T)
                # predict_no_j = x_no_j.dot(b_no_j)

                # # use normalized targets
                # residuals = self.norm_targets - predict_no_j
                # p_j = np.dot(norm[:,j], residuals)
                # print("Max Residual = {}".format(np.max(residuals)))
                # print(p_j)
                # # update beta_j based on residuals
                # self.betas[j] = self.soft_threshold(p_j, self.alpha / 2);
            # if iter % 100 == 0:
                # print(self.betas)
                
            # print(self.betas) 
Example 77
Project: Tacotron   Author: ElwynWang   File: signal_process.py    GNU General Public License v3.0 4 votes vote down vote up
def get_spectrograms(fpath):
    '''Returns normalized log(melspectrogram) and log(magnitude) from `sound_file`.
    Args:
      sound_file: A string. The full path of a sound file.

    Returns:
      mel: A 2d array of shape (T, n_mels) <- Transposed
      mag: A 2d array of shape (T, 1+n_fft/2) <- Transposed
    '''
    # Loading sound file
    y, sr = librosa.load(fpath, sr=Hp.sample_rate)

    # Trimming
    y, _ = librosa.effects.trim(y)

    # Preemphasis
    y = np.append(y[0], y[1:] - Hp.preemphasis * y[:-1])

    # stft
    linear = librosa.stft(y=y,
                          n_fft=Hp.num_fft,
                          hop_length=Hp.hop_length,
                          win_length=Hp.win_length)

    # magnitude spectrogram
    mag = np.abs(linear)  # (1+n_fft//2, T)

    # mel spectrogram
    mel_basis = librosa.filters.mel(Hp.sample_rate, Hp.num_fft, Hp.num_mels)  # (n_mels, 1+n_fft//2)
    mel = np.dot(mel_basis, mag)  # (n_mels, t)

    # to decibel
    mel = 20 * np.log10(np.maximum(1e-5, mel))
    mag = 20 * np.log10(np.maximum(1e-5, mag))

    # normalize
    mel = np.clip((mel - Hp.ref_db + Hp.max_db) / Hp.max_db, 1e-8, 1)
    mag = np.clip((mag - Hp.ref_db + Hp.max_db) / Hp.max_db, 1e-8, 1)

    # Transpose
    mel = mel.T.astype(np.float32)  # (T, n_mels)
    mag = mag.T.astype(np.float32)  # (T, 1+n_fft//2)

    return mel, mag 
Example 78
Project: dc_tts   Author: Kyubyong   File: utils.py    Apache License 2.0 4 votes vote down vote up
def get_spectrograms(fpath):
    '''Parse the wave file in `fpath` and
    Returns normalized melspectrogram and linear spectrogram.

    Args:
      fpath: A string. The full path of a sound file.

    Returns:
      mel: A 2d array of shape (T, n_mels) and dtype of float32.
      mag: A 2d array of shape (T, 1+n_fft/2) and dtype of float32.
    '''
    # Loading sound file
    y, sr = librosa.load(fpath, sr=hp.sr)

    # Trimming
    y, _ = librosa.effects.trim(y)

    # Preemphasis
    y = np.append(y[0], y[1:] - hp.preemphasis * y[:-1])

    # stft
    linear = librosa.stft(y=y,
                          n_fft=hp.n_fft,
                          hop_length=hp.hop_length,
                          win_length=hp.win_length)

    # magnitude spectrogram
    mag = np.abs(linear)  # (1+n_fft//2, T)

    # mel spectrogram
    mel_basis = librosa.filters.mel(hp.sr, hp.n_fft, hp.n_mels)  # (n_mels, 1+n_fft//2)
    mel = np.dot(mel_basis, mag)  # (n_mels, t)

    # to decibel
    mel = 20 * np.log10(np.maximum(1e-5, mel))
    mag = 20 * np.log10(np.maximum(1e-5, mag))

    # normalize
    mel = np.clip((mel - hp.ref_db + hp.max_db) / hp.max_db, 1e-8, 1)
    mag = np.clip((mag - hp.ref_db + hp.max_db) / hp.max_db, 1e-8, 1)

    # Transpose
    mel = mel.T.astype(np.float32)  # (T, n_mels)
    mag = mag.T.astype(np.float32)  # (T, 1+n_fft//2)

    return mel, mag 
Example 79
Project: StructEngPy   Author: zhuoju36   File: dynamic.py    MIT License 4 votes vote down vote up
def response_spectrum(model:Model,spec,mdd,n=60,comb='CQC'):
    """
    spec: a {'T':period,'a':acceleration} dictionary of spectrum\n
    mdd: a list of modal damping ratio\n
    comb: combination method, 'CQC' or 'SRSS'
    """
    K=model.K_
    M=model.M_
    DOF=model.DOF
    w,f,T,mode=eigen_mode(model,DOF)
    mode[n:,:]=np.zeros((DOF-n,DOF))#use n modes only.
    mode[:,n:]=np.zeros((DOF,DOF-n))
    M_=np.dot(np.dot(mode.T,M),mode)#generalized mass
    K_=np.dot(np.dot(mode.T,K),mode)#generalized stiffness
    L_=np.dot(np.diag(M),mode)
    px=[]
    Vx=[]
    Xm=[]
    gamma=[]
    mx=np.diag(M)
    for i in range(len(mode)):
        #mass participate factor
        px.append(-np.dot(mode[:,i].T,mx))
        Vx.append(px[-1]**2)
        Xm.append(Vx[-1]/3/m)
        #modal participate factor
        gamma.append(L_[i]/M_[i,i])    
    S=np.zeros((DOF,mode.shape[0]))
    

    for i in range(mode.shape[1]):        
        xi=mdd[i]
        y=np.interp(T[i],spec['T'],spec['a'])
        y/=w[i]**2
        S[:,i]=gamma[i]*y*mode[:,i]

    if comb=='CQC':
        cqc=0    
        rho=np.diag(np.ones(mode.shape[1]))
        for i in range(mode.shape[1]):
            for j in range(mode.shape[1]):
                if i!=j:
                    r=T[i]/T[j]
                    rho[i,j]=8*xi**2*(1+r)*r**1.5/((1-r**2)**2+4*xi**2*r*(1+r)**2)
                cqc+=rho[i,j]*S[:,i]*S[:,j]
        cqc=np.sqrt(cqc)
        print(cqc)
    elif comb=='SRSS':
        srss=0
        for i in range(mode.shape[1]):
            srss+=S[:,i]**2
        srss=np.sqrt(srss)
        print(srss) 
Example 80
Project: StructEngPy   Author: zhuoju36   File: element.py    MIT License 4 votes vote down vote up
def __init__(self,node_i, node_j, node_k, t, E, mu, rho, name=None):
        """
        params:
            node_i,node_j,node_k: Node, corners of triangle.
            t: float, thickness
            E: float, elastic modulus
            mu: float, Poisson ratio
            rho: float, mass density
        """
        super(Membrane3,self).__init__(node_i,node_j,node_k,t,E,mu,rho,6,name)

        x0=np.array([(node.x,node.y,node.z) for node in self._nodes])
        V=self._local_csys.transform_matrix
        o=self._local_csys.origin
        self._x0=(x0-np.array(o)).dot(V.T)[:,:2]
        
        D=self._D

        #calculate strain matrix
        abc0=self._abc(1,2)
        abc1=self._abc(2,0)
        abc2=self._abc(0,1)
        B0= np.array([[abc0[1],      0],
                      [      0,abc0[2]],
                      [abc0[2],abc0[1]]])
        B1= np.array([[abc1[1],     0],
                      [      0,abc1[2]],
                      [abc1[2],abc1[1]]])
        B2= np.array([[abc2[1],      0],
                      [      0,abc2[2]],
                      [abc2[2],abc2[1]]])
        self._B=np.hstack([B0,B1,B2])/2/self.area

        _Ke_=np.dot(np.dot(self._B(0).T,D),self._B(0))*self.area*self._t

        row=[a for a in range(0*2,0*2+2)]+\
            [a for a in range(1*2,1*2+2)]+\
            [a for a in range(2*2,2*2+2)]
        col=[a for a in range(0*6,0*6+2)]+\
            [a for a in range(1*6,1*6+2)]+\
            [a for a in range(2*6,2*6+2)]
        elm_node_count=3
        elm_dof=2
        data=[1]*(elm_node_count*elm_dof)
        G=sp.sparse.csr_matrix((data,(row,col)),shape=(elm_node_count*elm_dof,elm_node_count*6))
        self._Ke=G.transpose()*_Ke_*G

        #Concentrated mass matrix, may be wrong
        self._Me=np.eye(18)*rho*self.area*t/3
        
        self._re =np.zeros((18,1))