def svgd_kernel(self, h = -1):
        sq_dist = pdist(self.theta)
        pairwise_dists = squareform(sq_dist)**2
        if h < 0: # if h < 0,using median trick
            h = np.median(pairwise_dists)  
            h = np.sqrt(0.5 * h / np.log(self.theta.shape[0]+1))
 
        # compute the rbf kernel
 
        Kxy = np.exp( -pairwise_dists / h**2 / 2)
 
        dxkxy = -np.matmul(Kxy, self.theta)
        sumkxy = np.sum(Kxy, axis=1)
        for i in range(self.theta.shape[1]):
            dxkxy[:, i] = dxkxy[:,i] + np.multiply(self.theta[:,i],sumkxy)
        dxkxy = dxkxy / (h**2)
        return (Kxy, dxkxy)

def KeyGen(**kwargs):
    ''''''
    Appendix B of BLISS paper
    m_bar = m + n
 
    o/p:    
    A: Public Key n x m'' numpy array
    S: Secret Key m''x n numpy array
    ''''''
    q, n, m, alpha = kwargs[''q''], kwargs[''n''], kwargs[''m''], kwargs[''alpha'']
    Aq_bar = util.crypt_secure_matrix(-(q-1)/2, (q-1)/2, m)
    S_bar = util.crypt_secure_matrix(-(2)**alpha, (2)**alpha, n) # alpha is small enough,we need not reduce (modq)
    S = np.vstack((S_bar, np.eye(n, dtype = int))) # dimension is m_bar x n,Elements are in Z mod(2q)
    A = np.hstack((2*Aq_bar, q * np.eye(n, dtype = int) - 2*np.matmul(Aq_bar,S_bar))) # dimension is n x m_bar,Elements are in Z mod(2q)
    #return util.matrix_to_Zq(A,2*q),S,Aq_bar,S_bar
    return util.matrix_to_Zq(A, 2*q), S

def test():
    # Classical SIS parameters
    n, alpha, q = 128, 872, 1, 114356107
    kappa = 20
 
    #discrete Gaussian Parameters
    sd = 300
    eta = 1.2
 
    A, S = KeyGen(q = q,n = n,m = m,alpha = alpha)
    #print np.array(np.matmul(A,S) - q*np.eye(n),dtype=float)/(2*q) #to test AS = q mod(2q)
    z, c = Sign(msg = "Hello Bob",A = A,S = S,sd = sd,q = q,M = 3.0,kappa = kappa)
    print z
    print c
    print Verify(msg = "Hello Bob", A=A, m=m, n=n, sd=sd, q=q, eta=eta, z=z, c=c, kappa = kappa)
    print Verify(msg = "Hello Robert", kappa = kappa)
    print Verify(msg = "Hello Roberto", kappa = kappa)
    print Verify(msg = "Hola Roberto", kappa = kappa)

def Verify(**kwargs):
    ''''''
        Verification for the signature
        i/p:
        msg: the string sent by the sender
        (z,c): vectors in Zq,the signature
        A  : numpy array,Verification Key dimension nxm
        T : the matrix AS mod q,it is used in the Verification of the signature
    ''''''
    msg, z, c, A, T, sd, eta, k, q = kwargs[''msg''], kwargs[''z''], kwargs[''c''], kwargs[''A''], kwargs[''T''], kwargs[''sd''], kwargs[''eta''], kwargs[''k''], kwargs[''q'']
    norm_bound = eta * sd * np.sqrt(m)
    # checks for norm of z being small and that H(Az-Tc mod q,msg) hashes to c
    vec = util.vector_to_Zq(np.array(np.matmul(A,z) - np.matmul(T,c)), q)
    hashedList = util.hash_to_baseb(vec, msg, 3, k)
    print hashedList, c             
    if np.sqrt(z.dot(z)) <= norm_bound and np.array_equal(c, hashedList):
        return True
    else:
        return False

def KeyGen(n, d, q):
    ''''''
        input:
        q : polynomial size prime number
        n,m,k : dimensions specifiers
        d : SIS parameter,hardest instances are where d ~ q^(n/m)
 
        output:
        Signing Key S :  Matrix of dimension mxk with coefficients in [-d.d]
        Verification Key A : Matrix of dimension nxm with coefficients from [-(q-1)/2,(q-1)/2]
        T : the matrix AS,it is used in the Verification of the signature
 
    ''''''
    S = crypt_secure_matrix(d, k)
    A = crypt_secure_matrix((q-1)/2, m)
    T = np.matmul(A, S)
    return S, T

def transfer_color(content, style):
    import scipy.linalg as sl
    # Mean and covariance of content
    content_mean = np.mean(content, axis = (0, 1))
    content_diff = content - content_mean
    content_diff = np.reshape(content_diff, (-1, content_diff.shape[2]))
    content_covariance = np.matmul(content_diff.T, content_diff) / (content_diff.shape[0])
 
    # Mean and covariance of style
    style_mean = np.mean(style, 1))
    style_diff = style - style_mean
    style_diff = np.reshape(style_diff, style_diff.shape[2]))
    style_covariance = np.matmul(style_diff.T, style_diff) / (style_diff.shape[0])
 
    # Calculate A and b
    A = np.matmul(sl.sqrtm(content_covariance), sl.inv(sl.sqrtm(style_covariance)))
    b = content_mean - np.matmul(A, style_mean)
 
    # Construct new style
    new_style = np.reshape(style, style.shape[2])).T
    new_style = np.matmul(A, new_style).T
    new_style = np.reshape(new_style, style.shape)
    new_style = new_style + b
 
    return new_style

def svgd_kernel(self, theta, h = -1):
        sq_dist = pdist(theta)
        pairwise_dists = squareform(sq_dist)**2
        if h < 0: # if h < 0,using median trick
            h = np.median(pairwise_dists)  
            h = np.sqrt(0.5 * h / np.log(theta.shape[0]+1))
 
        # compute the rbf kernel
        Kxy = np.exp( -pairwise_dists / h**2 / 2)
 
        dxkxy = -np.matmul(Kxy, theta)
        sumkxy = np.sum(Kxy, axis=1)
        for i in range(theta.shape[1]):
            dxkxy[:,i] + np.multiply(theta[:, dxkxy)

def run_random_sim(sim, L):
    """
    Run *L* simulations of the state space model specified by *sim*
    (see :func:`setup_random_sim`). Each simulation is added to *sim*
    index by an integer identifier.
    """
    sim[''L''] = L
    for l in range(L):
        sim[l] = defaultdict(list)
        x_i = sim[''mu''] + NP.matmul(sim[''PI_sqrt''], NP.random.randn(sim[''N'']))
        for i in range(sim[''I'']):
            sim[l][''x''].append(x_i)
            # measurement
            v_i = NP.matmul(sim[''R_sqrt''][i], NP.random.randn(sim[''M'']))
            sim[l][''y''].append(NP.matmul(sim[''H''][i], sim[l][''x''][i]) + v_i)
            # time update
            u_i = NP.matmul(sim[''Q_sqrt''][i], NP.random.randn(sim[''N'']))
            x_i = NP.matmul(sim[''F''][i], x_i) + u_i
    return sim

def sqrt_kf_sim(sim):
    """
    Process each simulation trial generated with
    :func:`setup_random_test` with a Kalman filter and return the
    posterior state estimates and error covariances.
    """
    post = defaultdict(dict)
    for l in range(sim[''L'']):
        x_hat_l, P_sqrt_l = sqrt_kalman_filter(sim[l][''y''],
                                               sim[''H''],
                                               sim[''R_sqrt''],
                                               sim[''F''],
                                               sim[''Q_sqrt''],
                                               sim[''mu''],
                                               sim[''PI_sqrt''])
        post[l][''x_hat''] = x_hat_l
        if l == 0:
            post[''P''] = [NP.matmul(x, x.T) for x in P_sqrt_l]
        post[l][''error''] = []
        for x_i, x_hat_i in izip(sim[l][''x''], post[l][''x_hat'']):
            post[l][''error''].append(x_hat_i - x_i)
    return post

def sqrt_kf_tu(x_hat_posterior,
               P_sqrt_posterior,
               F_i,
               Q_sqrt_i,
               z_i=None):
    """
    Square root Kalman filter time update. Given the following:
    - *x_hat_posterior*: posterior state estimate (N)
    - *P_sqrt_posterior*: posterior error covariance square root (NxN)
    - *F_i*: time update operator (NxN)
    - *Q_sqrt_i*: time update noise covariance square root (NxN)
    - *z_i*: optional) systematic time update input (N)
 
    Return the tuple containing the one time step prediction of the
    state and the square root of the error covariance.
    """
    N, _ = F_i.shape
    x_hat_prior = NP.matmul(F_i, x_hat_posterior)
    if z_i is not None:
        x_hat_prior += z_i
    A_T = NP.block([NP.matmul(F_i, P_sqrt_posterior), Q_sqrt_i])
    R_T = NP.linalg.qr(A_T.T, mode=''r'')
    P_sqrt_prior = R_T.T[:, :N]
    return x_hat_prior, P_sqrt_prior

def rotate_ascii_stl(self, rotation_matrix, content, filename):
        """Rotate the mesh array and save as ASCII STL."""
        mesh = np.array(content, dtype=np.float64)
 
        # prefix area vector,if not already done (e.g. in STL format)
        if len(mesh[0]) == 3:
            row_number = int(len(content)/3)
            mesh = mesh.reshape(row_number, 3)
 
        # upgrade numpy with: "pip install numpy --upgrade"
        rotated_content = np.matmul(mesh, rotation_matrix)
 
        v0 = rotated_content[:, 0, :]
        v1 = rotated_content[:, :]
        v2 = rotated_content[:, 2, :]
        normals = np.cross(np.subtract(v1, v0), np.subtract(v2, v0)) \\
            .reshape(int(len(rotated_content)), 3)
        rotated_content = np.hstack((normals, rotated_content))
 
        tweaked = list("solid %s" % filename)
        tweaked += list(map(self.write_facett, list(rotated_content)))
        tweaked.append("\\nendsolid %s\\n" % filename)
        tweaked = "".join(tweaked)
 
        return tweaked

def test_exceptions(self):
        dims = [
            ((1,), (2,)),            # mismatched vector vector
            ((2,         # mismatched matrix vector
            ((2, (1, 2)),          # mismatched vector matrix
            ((1, 2), (3, 1)),        # mismatched matrix matrix
            ((1, ()),              # vector scalar
            ((), (1)),               # scalar vector
            ((1, 1),            # matrix scalar
            ((),            # scalar matrix
            ((2,  # cannot broadcast
            ]
 
        for dt, (dm1, dm2) in itertools.product(self.types, dims):
            a = np.ones(dm1, dtype=dt)
            b = np.ones(dm2, dtype=dt)
            assert_raises(ValueError, self.matmul, a, b)

def test_shapes(self):
        dims = [
            ((1,     # broadcast first argument
            ((2,     # broadcast second argument
            ((2,  # matrix stack sizes match
            ]
 
        for dt, dtype=dt)
            res = self.matmul(a, b)
            assert_(res.shape == (2, 1))
 
        # vector vector returns scalars.
        for dt in self.types:
            a = np.ones((2, dtype=dt)
            b = np.ones((2, dtype=dt)
            c = self.matmul(a, b)
            assert_(np.array(c).shape == ())

def test_numpy_ufunc_override(self):
        # 2016-01-29: NUMPY_UFUNC_disABLED
        return
 
        class A(np.ndarray):
            def __new__(cls, *args, **kwargs):
                return np.array(*args, **kwargs).view(cls)
 
            def __numpy_ufunc__(self, ufunc, method, pos, inputs, **kwargs):
                return "A"
 
        class B(np.ndarray):
            def __new__(cls, **kwargs):
                return NotImplemented
 
        a = A([1, 2])
        b = B([1, 2])
        c = np.ones(2)
        assert_equal(self.matmul(a, b), "A")
        assert_equal(self.matmul(b, a), "A")
        assert_raises(TypeError, b, c)

def get_symmetry_permutation(self):
        """
        This a object function to get the permutation group operators.
        Represented as a table.
        """
        sym_perm = []
        numbers = [i for i in range(self.num_count)]
        sym_mat = spglib.get_symmetry(self._spg_cell, symprec=self.symprec)
        ops = [(r, t) for r, t in zip(sym_mat[''rotations''],
                                      sym_mat[''translations''])]
        for r, t in ops:
            pos_new = np.transpose(np.matmul(r, self._positions.T)) + t
            perm = self._get_new_id_seq(pos_new, numbers)
            sym_perm.append(perm)
 
        return sym_perm

def supercell(self, scale_mat):
        """
        Get the supercell of the origin gcell
        scale_mat is similar as H matrix in superlattice generator
        """
        # return self.__class__(...)
        sarr_lat = np.matmul(scale_mat, self.lattice)
        # coor_conv_pos = np.matmul(self.positions,self.lattice)
        # o_conv_pos = np.matmul(coor_conv_pos,np.linalg.inv(scale_mat))
        o_conv_pos = np.matmul(self.positions, np.linalg.inv(scale_mat))
        o_pos = self.get_frac_from_mat(scale_mat)
 
        l_of_positions = [i for i in map(lambda x: x+o_pos, list(o_conv_pos))]
        pos = np.concatenate(l_of_positions, axis=0)
 
        n = scale_mat.diagonal().prod()
        numbers = np.repeat(self.numbers, n)
 
        return self.__class__(sarr_lat, numbers)

def ams_extractor(x, sr, win_len, shift_len, order):
    from scipy.signal import hilbert
    envelope = np.abs(hilbert(x))
    for i in range(order-1):
        envelope = np.abs(hilbert(envelope))
    envelope = envelope * 1./3.
    frames = (len(envelope) - win_len) // shift_len
    hanning_window = np.hanning(win_len)
    ams_feature = np.zeros(shape=(15, frames))
    wts = cal_triangle_window(0, sr//2, 15, 15.6, 400)
    for i in range(frames):
        one_frame = x[i*shift_len:i*shift_len+win_len]
        one_frame = one_frame * hanning_window
        frame_fft = np.abs(np.fft.fft(one_frame, win_len))
        ams_feature[:,i] = np.matmul(wts, frame_fft)
    return ams_feature

def ams_extractor(x, order=1, decimate_coef=1./4.):
    from scipy.signal import hilbert
    envelope = np.abs(hilbert(x))
    for i in range(order-1):
        envelope = np.abs(hilbert(envelope))
    envelope = envelope * decimate_coef
    frames = (len(envelope) - win_len) // shift_len
    hanning_window = np.hanning(win_len)
    ams_feature = np.zeros(shape=(15, frame_fft)
    return ams_feature

def unkNown_feature_extractor(x, barks, inner_win, inner_shift, win_type, method_version):
    x_spectrum = stft_extractor(x, win_type)
    coef = get_fft_bark_mat(sr, 20, sr//2)
    bark_spect = np.matmul(coef, x_spectrum)
    ams = np.zeros((barks, inner_win//2+1, (bark_spect.shape[1] - inner_win)//inner_shift))
    for i in range(barks):
        channel_stft = stft_extractor(bark_spect[i, :], ''hanning'')
        if method_version == ''v1'':
            ams[i, :, :] = 20 * np.log(np.abs(channel_stft[:inner_win//2+1, :(bark_spect.shape[1] - inner_win)//inner_shift]))
        elif method_version == ''v2'':
            channel_amplitude = np.abs(channel_stft[:inner_win//2+1, :(bark_spect.shape[1] - inner_win)//inner_shift])
            channel_angle = np.angle(channel_stft[:inner_win//2+1, :(bark_spect.shape[1] - inner_win)//inner_shift])
            channel_angle = channel_angle - (np.floor(channel_angle / (2.*np.pi)) * (2.*np.pi))
            ams[i, :] = np.power(channel_amplitude, 1./3.) * channel_angle
        else:
            ams[i, :] = np.abs(channel_stft)
    return ams

def ams_extractor(x, :] = np.abs(channel_stft)
    return ams

def run_test_matmul_aa_ci8_shape(self, shape, transpose=False):
        # **Todo: This currently never triggers the transpose path in the backend
        shape_complex = shape[:-1] + (shape[-1] * 2,)
        # Note: The xGPU-like correlation kernel does not support input values of -128 (only [-127:127])
        a8 = ((np.random.random(size=shape_complex) * 2 - 1) * 127).astype(np.int8)
        a_gold = a8.astype(np.float32).view(np.complex64)
        if transpose:
            a_gold = H(a_gold)
        # Note: np.matmul seems to be slow and inaccurate when there are batch dims
        c_gold = np.matmul(a_gold, H(a_gold))
        triu = np.triu_indices(shape[-2] if not transpose else shape[-1], 1)
        c_gold[..., triu[0], triu[1]] = 0
        a = a8.view(bf.DataType.ci8)
        a = bf.asarray(a, space=''cuda'')
        if transpose:
            a = H(a)
        c = bf.zeros_like(c_gold, space=''cuda'')
        self.linalg.matmul(1, None, c)
        c = c.copy(''system'')
        np.testing.assert_allclose(c, c_gold, RTOL, ATOL)

def run_test_matmul_aa_dtype_shape(self, dtype, axes=None, conj=False):
        a = ((np.random.random(size=shape)) * 127).astype(dtype)
        if axes is None:
            axes = range(len(shape))
        aa = a.transpose(axes)
        if conj:
            aa = aa.conj()
        c_gold = np.matmul(aa, H(aa))
        triu = np.triu_indices(shape[axes[-2]], triu[1]] = 0
        a = bf.asarray(a, space=''cuda'')
        aa = a.transpose(axes)
        if conj:
            aa = aa.conj()
        c = bf.zeros_like(c_gold, aa, ATOL)

def run_test_matmul_ab_ci8_shape(self, transpose=False):
        ashape_complex = shape[:-2] + (shape[-2], k * 2)
        bshape_complex = shape[:-2] + (k, shape[-1] * 2)
        a8 = (np.random.random(size=ashape_complex) * 255).astype(np.int8)
        b8 = (np.random.random(size=bshape_complex) * 255).astype(np.int8)
        a_gold = a8.astype(np.float32).view(np.complex64)
        b_gold = b8.astype(np.float32).view(np.complex64)
        if transpose:
            a_gold, b_gold = H(b_gold), H(a_gold)
        c_gold = np.matmul(a_gold, b_gold)
        a = a8.view(bf.DataType.ci8)
        b = b8.view(bf.DataType.ci8)
        a = bf.asarray(a, space=''cuda'')
        b = bf.asarray(b, space=''cuda'')
        if transpose:
            a, b = H(b), H(a)
        c = bf.zeros_like(c_gold, ATOL)