Python numpy.kron() 使用实例

The following are code examples for showing how to use . They are extracted from open source Python projects. You can vote up the examples you like or vote down the exmaples you don’t like. You can also save this page to your account.

Example 1

def apply_gate(self,gate,on_qubit_name):
		on_qubit=self.qubits.get_quantum_register_containing(on_qubit_name)
		if len(on_qubit.get_noop()) > 0:
			print "NOTE this qubit has been measured previously, there should be no more gates allowed but we are reverting that measurement for consistency with IBM's language"
			on_qubit.set_state(on_qubit.get_noop())
			on_qubit.set_noop([])
		if not on_qubit.is_entangled():
			if on_qubit.get_num_qubits()!=1:
				raise Exception("This qubit is not marked as entangled but it has an entangled state")
			on_qubit.set_state(gate*on_qubit.get_state())
		else:
			if not on_qubit.get_num_qubits()>1:
				raise Exception("This qubit is marked as entangled but it does not have an entangled state")
			n_entangled=len(on_qubit.get_entangled())
			apply_gate_to_qubit_idx=[qb.name for qb in on_qubit.get_entangled()].index(on_qubit_name)
			if apply_gate_to_qubit_idx==0:
				entangled_gate=gate
			else:
				entangled_gate=Gate.eye
			for i in range(1,n_entangled):
				if apply_gate_to_qubit_idx==i:
					entangled_gate=np.kron(entangled_gate,gate)
				else:
					entangled_gate=np.kron(entangled_gate,Gate.eye)
			on_qubit.set_state(entangled_gate*on_qubit.get_state()) 

Example 2

def defineSensorLoc(self,XYZ):
        #############################################
        # DEFINE TRANSMITTER AND RECEIVER LOCATIONS
        #
        # XYZ: N X 3 array containing transmitter center locations
        # **NOTE** for this sensor, we know where the receivers are relative to transmitters
        self.TxLoc = XYZ

        dx,dy = np.meshgrid([-0.8,-0.4,0.,0.4,0.8],[-0.8,-0.4,0.,0.4,0.8])
        dx = mkvc(dx)
        dy = mkvc(dy)

        N = np.shape(XYZ)[0]

        X = np.kron(XYZ[:,0],np.ones((25))) + np.kron(np.ones((N)),dx)
        Y = np.kron(XYZ[:,1],np.ones((25))) + np.kron(np.ones((N)),dy)
        Z = np.kron(XYZ[:,2],np.ones((25)))

        self.RxLoc = np.c_[X,Y,Z]

        self.TxID = np.kron(np.arange(1,np.shape(XYZ)[0]+1),np.ones((25)))
        self.RxComp = np.kron(3*np.ones(np.shape(XYZ)[0]),np.ones((25))) 

Example 3

def defineSensorLoc(self,XYZ):
        #############################################
        # DEFINE TRANSMITTER AND RECEIVER LOCATIONS
        #
        # XYZ: N X 3 array containing transmitter center locations
        # **NOTE** for this sensor, we know where the receivers are relative to transmitters
        self.TxLoc = XYZ

        dx = np.kron([-0.18,0.,0.,0.,0.18],np.ones(3))
        dy = np.kron([0.,-0.18,0.,0.18,0.],np.ones(3))

        N = np.shape(XYZ)[0]

        X = np.kron(XYZ[:,0],np.ones((15))) + np.kron(np.ones((N)),dx)
        Y = np.kron(XYZ[:,1],np.ones((15))) + np.kron(np.ones((N)),dy)
        Z = np.kron(XYZ[:,2],np.ones((15)))

        self.RxLoc = np.c_[X,Y,Z]

        self.TxID = np.kron(np.arange(1,np.shape(XYZ)[0]+1),np.ones((15)))
        self.RxComp = np.kron(np.kron(np.ones(np.shape(XYZ)[0]),np.ones((5))),np.r_[1,2,3]) 

Example 4

def toa_calc_d_from_xy(x, y):
    dimx = x.shape
    x_dim = dimx[0]
    m = dimx[1]

    dimy = y.shape
    y_dim = dimy[0]
    n = dimy[1]

    if x_dim > y_dim:
        y[(y_dim + 1):x_dim, ] = np.zeros(x_dim - y_dim, n)
    elif y_dim > x_dim:
        x[(x_dim + 1):y_dim, ] = np.zeros(y_dim - x_dim, n)

    d = sqrt(
        ((np.kron(np.ones(1, n), x) - np.kron(y, np.ones(1, m))) ** 2).sum(
            axis=0))
    d = np.asarray(d)
    d = np.reshape(d, (m, n), order='F')

    return d 

Example 5

def edgeCurl(self):
        """The edgeCurl property."""
        if self.nCy > 1:
            raise NotImplementedError('Edge curl not yet implemented for '
                                      'nCy > 1')
        if getattr(self, '_edgeCurl', None) is None:
            # 1D Difference matricies
            dr = sp.spdiags((np.ones((self.nCx+1, 1))*[-1, 1]).T, [-1, 0],
                            self.nCx, self.nCx, format="csr")
            dz = sp.spdiags((np.ones((self.nCz+1, 1))*[-1, 1]).T, [0, 1],
                            self.nCz, self.nCz+1, format="csr")

            # 2D Difference matricies
            Dr = sp.kron(sp.identity(self.nNz), dr)
            Dz = -sp.kron(dz, sp.identity(self.nCx))

            A = self.area
            E = self.edge
            # Edge curl operator
            self._edgeCurl = (utils.sdiag(1/A)*sp.vstack((Dz, Dr)) *
                              utils.sdiag(E))
        return self._edgeCurl 

Example 6

def aveE2CC(self):
        "Construct the averaging operator on cell edges to cell centers."
        if getattr(self, '_aveE2CC', None) is None:
            # The number of cell centers in each direction
            n = self.vnC
            if self.isSymmetric:
                avR = utils.av(n[0])[:, 1:]
                avR[0, 0] = 1.
                self._aveE2CC = sp.kron(utils.av(n[2]), avR, format="csr")
            else:
                raise NotImplementedError('wrapping in the averaging is not '
                                          'yet implemented')
                # self._aveE2CC = (1./3)*sp.hstack((utils.kron3(utils.av(n[2]),
                #                                               utils.av(n[1]),
                #                                               utils.speye(n[0])),
                #                                   utils.kron3(utils.av(n[2]),
                #                                               utils.speye(n[1]),
                #                                               utils.av(n[0])),
                #                                   utils.kron3(utils.speye(n[2]),
                #                                               utils.av(n[1]),
                #                                               utils.av(n[0]))),
                #                                  format="csr")
        return self._aveE2CC 

Example 7

def f(W,G,y,hparams):
        # f = -1/N*sum_t log(exp(w(yt)'gt)/sum_k exp(wk'gt)) + l*||W||
        # = -1/N*sum_t [w(yt)'*gt - log(sum_k exp(wk'gt))] + l*||W||
        # = -1/N*sum(sum(W(:,y).*G,1),2) + 1/N*sum(log(sumexpWG),2) + l*sum(sum(W.^2));

        #K,l = hparams
        K = hparams['K']
        l = hparams['l']
        d,N = G.shape
        W = W.reshape((d,K))
        
        WG = np.dot(W.T,G) # K x N
        WG -= np.kron(np.ones((K,1)),WG.max(axis=0).reshape(1,N))
        #WG_max = WG.max(axis=0).reshape((1,N))
        #expWG = np.exp(WG-np.kron(np.ones((K,1)),WG_max)) # K x N
        expWG = np.exp(WG) # K x N
        sumexpWG = expWG.sum(axis=0) # N x 1
        WyG = WG[y,range(N)]
        #WyG -= WG_max
        
        fval = -1.0/N*(WyG).sum() \
            + 1.0/N*np.log(sumexpWG).sum() \
            + l*(W**2).sum()#(axis=(0,1))
    
        return fval 

Example 8

def dfdv(W,G,y,hparams):
        # df/dwk = -1/N*sum(x(:,y==k),2) + 1/N*sum_t exp(wk'xt)*xt/(sum_k exp(wk'xt))] + l*2*wk
        K = hparams['K']
        l = hparams['l']
        d,N = G.shape
        shapeW = W.shape
        W = W.reshape((d,K))
        
        WG = np.dot(W.T,G) # K x N
        WG -= np.kron(np.ones((K,1)),WG.max(axis=0).reshape(1,N))
        expWG = np.exp(WG) # K x N
        sumexpWG = expWG.sum(axis=0) # N x 1
        df = np.zeros((d,K))
        for k in range(K):
            indk = np.where(y==k)[0]    
            df[:,k] = -1./N*G[:,indk].sum(axis=1).reshape((d,)) \
                + 1./N*np.dot(G,(expWG[k,:]/sumexpWG).T).reshape((d,)) \
                + 2.*l*W[:,k].reshape((d,))
        
        assert np.isnan(df).any()==False        
        
        return df.reshape(shapeW) 

Example 9

def estimator_cov(self,method):
        """ Creates covariance matrix for the estimators

        Parameters
        ----------
        method : str
            Estimation method

        Returns
        ----------
        A Covariance Matrix
        """         
        
        Y = np.array([reg[self.lags:] for reg in self.data])    
        Z = self._create_Z(Y)
        if method == 'OLS':
            sigma = self.ols_covariance()
        else:           
            sigma = self.custom_covariance(self.latent_variables.get_z_values())
        return np.kron(np.linalg.inv(np.dot(Z,np.transpose(Z))), sigma) 

Example 10

def entangling_mat(gate):
    """
	Helper function to create the entangling gate matrix
	"""
    echoCR = expm(1j * pi / 4 * np.kron(pX, pZ))
    if gate == "CNOT":
        return echoCR
    elif gate == "iSWAP":
        return reduce(lambda x, y: np.dot(y, x),
                      [echoCR, np.kron(C1[6], C1[6]), echoCR])
    elif gate == "SWAP":
        return reduce(lambda x, y: np.dot(y, x),
                      [echoCR, np.kron(C1[6], C1[6]), echoCR, np.kron(
                          np.dot(C1[6], C1[1]), C1[1]), echoCR])
    else:
        raise ValueError("Entangling gate must be one of: CNOT, iSWAP, SWAP.") 

Example 11

def clifford_mat(c, numQubits):
    """
	Return the matrix unitary the implements the qubit clifford C
	"""
    assert numQubits <= 2, "Oops! I only handle one or two qubits"
    if numQubits == 1:
        return C1[c]
    else:
        c = C2Seqs[c]
        mat = np.kron(clifford_mat(c[0][0], 1), clifford_mat(c[0][1], 1))
        if c[1]:
            mat = np.dot(entangling_mat(c[1]), mat)
        if c[2]:
            mat = np.dot(
                np.kron(
                    clifford_mat(c[2][0], 1), clifford_mat(c[2][1], 1)), mat)
        return mat 

Example 12

def proposals(self, X, return_index=False):
        """Retrieve proposals for a video based on its duration

        Parameters
        ----------
        X : ndarray
            m x 1 array with video duration

        Outputs
        -------
        Y : ndarray
            m * n_prop x 2 array with temporal proposals

        """
        if self.priors is None:
            raise ValueError('model has not been trained')

        Y = np.kron(np.expand_dims(X, 1), self.priors)
        idx = np.repeat(np.arange(X.shape[0]), self.n_prop)
        if return_index:
            return Y, idx
        return Y 

Example 13

def proposals(self, X, return_index=False):
        """Retrieve proposals for a video based on its duration

        Parameters
        ----------
        X : ndarray
            m x 1 array with video duration

        Outputs
        -------
        Y : ndarray
            m * n_prop x 2 array with temporal proposals

        """
        if self.priors is None:
            raise ValueError('model has not been trained')

        Y = np.kron(np.expand_dims(X, 1), self.priors)
        idx = np.repeat(np.arange(X.shape[0]), self.n_prop)
        if return_index:
            return Y, idx
        return Y 

Example 14

def apply_gate(self,gate,on_qubit_name):
		on_qubit=self.qubits.get_quantum_register_containing(on_qubit_name)
		if len(on_qubit.get_noop()) > 0:
			print "NOTE this qubit has been measured previously, there should be no more gates allowed but we are reverting that measurement for consistency with IBM's language"
			on_qubit.set_state(on_qubit.get_noop())
			on_qubit.set_noop([])
		if not on_qubit.is_entangled():
			if on_qubit.get_num_qubits()!=1:
				raise Exception("This qubit is not marked as entangled but it has an entangled state")
			on_qubit.set_state(gate*on_qubit.get_state())
		else:
			if not on_qubit.get_num_qubits()>1:
				raise Exception("This qubit is marked as entangled but it does not have an entangled state")
			n_entangled=len(on_qubit.get_entangled())
			apply_gate_to_qubit_idx=[qb.name for qb in on_qubit.get_entangled()].index(on_qubit_name)
			if apply_gate_to_qubit_idx==0:
				entangled_gate=gate
			else:
				entangled_gate=Gate.eye
			for i in range(1,n_entangled):
				if apply_gate_to_qubit_idx==i:
					entangled_gate=np.kron(entangled_gate,gate)
				else:
					entangled_gate=np.kron(entangled_gate,Gate.eye)
			on_qubit.set_state(entangled_gate*on_qubit.get_state()) 

Example 15

def cov(M):
    """
    Compute the sample covariance matrix of a 2D matrix.

    Parameters:
      M: `numpy array`
        2d matrix of HSI data (N x p)

    Returns: `numpy array`
        sample covariance matrix
    """
    N = M.shape[0]
    u = M.mean(axis=0)
    M = M - np.kron(np.ones((N, 1)), u)
    C = np.dot(M.T, M) / (N-1)
    return C 

Example 16

def com(A, d1, d2):
    """Calculation of the center of mass for spatial components
     Inputs:
     A:   np.ndarray
          matrix of spatial components (d x K)
     d1:  int
          number of pixels in x-direction
     d2:  int
          number of pixels in y-direction

     Output:
     cm:  np.ndarray
          center of mass for spatial components (K x 2)
    """
    nr = np.shape(A)[-1]
    Coor = dict()
    Coor['x'] = np.kron(np.ones((d2, 1)), np.expand_dims(range(d1), axis=1))
    Coor['y'] = np.kron(np.expand_dims(range(d2), axis=1), np.ones((d1, 1)))
    cm = np.zeros((nr, 2))        # vector for center of mass
    cm[:, 0] = np.dot(Coor['x'].T, A) / A.sum(axis=0)
    cm[:, 1] = np.dot(Coor['y'].T, A) / A.sum(axis=0)

    return cm 

Example 17

def find_activity_intervals(C,Npeaks = 5, tB=-5, tA = 25, thres = 0.3):
        
    import peakutils
    K,T = np.shape(C)
    L = []
    for i in range(K):
        indexes = peakutils.indexes(C[i,:],thres=thres)        
        srt_ind = indexes[np.argsort(C[i,indexes])][::-1]
        srt_ind = srt_ind[:Npeaks]
        L.append(srt_ind)
        
    LOC = []
    for i in range(K):
        if len(L[i])>0:
            interval = np.kron(L[i],np.ones(tA-tB,dtype=int)) + np.kron(np.ones(len(L[i]),dtype=int),np.arange(tB,tA))                        
            interval[interval<0] = 0
            interval[interval>T-1] = T-1
            LOC.append(np.array(list(set(interval))))        
        else:
            LOC.append(None)                        

        
    return LOC 

Example 18

def paulidouble(i, j, tensor=True):
    pauli_i = pauli(i)
    pauli_j = pauli(j)
    outer = np.zeros((4, 4), dtype=np.complex)
    outer[:2, :2] = pauli_i
    outer[2:, 2:] = pauli_j
    if tensor:
        outer.shape = (2, 2, 2, 2)
    return outer

# def paulitwo_left(i):
    # return np.kron(pauli(i), pauli(0)) 

# def paulitwo_right(i):
    # return np.kron(pauli(0), pauli(i))

# def newrho_DK(Lket, theta_ij):
    # Lbra = np.conjugate(L_before)
    # theta_star = np.conjugate(theta_ij)
    # in_bracket = scon(Lbra, Lket, theta_ij, theta_star,
                # [1], [1], [2, 3, 1, 

Example 19

def to_0xyz_basis(ptm):
    """Transform a Pauli transfer in the 0xy1 basis [1],
    to the the usual 0xyz. The inverse of to_0xy1_basis.

    ptm: The input transfer matrix in 0xy1 basis. Must be 4x4.

    [1] Daniel Greenbaum, Introduction to Quantum Gate Set Tomography, http://arxiv.org/abs/1509.02921v1
    """

    ptm = np.array(ptm)
    if ptm.shape == (4, 4):
        trans_mat = basis_transformation_matrix
        return np.dot(trans_mat, np.dot(ptm, trans_mat))
    elif ptm.shape == (16, 16):
        trans_mat = np.kron(basis_transformation_matrix, basis_transformation_matrix)
        return np.dot(trans_mat, np.dot(ptm, trans_mat))
    else:
        raise ValueError("Dimensions wrong, must be one- or two Pauli transfer matrix ") 

Example 20

def apply_two_ptm(self, bit0, bit1, two_ptm):
        """Apply a two_bit_ptm between bit0 and bit1.
        """
        self.ensure_dense(bit0)
        self.ensure_dense(bit1)

        ptm0 = np.eye(4)
        if bit0 in self.single_ptms_to_do:
            for ptm2 in self.single_ptms_to_do[bit0]:
                ptm0 = ptm2.dot(ptm0)
            del self.single_ptms_to_do[bit0]

        ptm1 = np.eye(4)
        if bit1 in self.single_ptms_to_do:
            for ptm2 in self.single_ptms_to_do[bit1]:
                ptm1 = ptm2.dot(ptm1)
            del self.single_ptms_to_do[bit1]

        full_two_ptm = np.dot(two_ptm, np.kron(ptm1, ptm0))
        self.full_dm.apply_two_ptm(self.idx_in_full_dm[bit0], 
                self.idx_in_full_dm[bit1], full_two_ptm) 

Example 21

def kron_all(op,num,op_2): 
    # returns an addition of sth like xii + ixi + iix for op =x and op_2 =i
    total = np.zeros([len(op)**num,len(op)**num])
    a=op
    for jj in range(num):
        if jj != 0:
            a = op_2
        else:
            a = op
            
        for ii in range(num-1):
            if (jj - ii) == 1:
                
                b = op
            else:
                b = op_2
            a = np.kron(a,b)
        total = total + a
    return a 

Example 22

def outer_product(input_array):
    r'''
    Takes a `NumPy` array and returns the outer (dyadic, Kronecker) product
    with itself. If `input_array` is a vector :math:`\mathbf{x}`, this returns
    :math:`\mathbf{x}\mathbf{x}^T`.
    '''
    la = len(input_array)

    # return outer product as numpy array
    return np.kron(input_array, input_array).reshape(la, la)

##############
# Decorators #
##############


# Main decorator for fit functions 

Example 23

def _upsample_filters(self, filters, rate):
    """Upsamples the filters by a factor of rate along the spatial dimensions.

    Args:
      filters: [h, w, in_depth, out_depth]. Original filters.
      rate: An int, specifying the upsampling rate.

    Returns:
      filters_up: [h_up, w_up, in_depth, out_depth]. Upsampled filters with
        h_up = h + (h - 1) * (rate - 1)
        w_up = w + (w - 1) * (rate - 1)
        containing (rate - 1) zeros between consecutive filter values along
        the filters' spatial dimensions.
    """
    if rate == 1:
      return filters
    # [h, w, in_depth, out_depth] -> [in_depth, out_depth, h, w]
    filters_up = np.transpose(filters, [2, 3, 0, 1])
    ker = np.zeros([rate, rate])
    ker[0, 0] = 1
    filters_up = np.kron(filters_up, ker)[:, :, :-(rate-1), :-(rate-1)]
    # [in_depth, out_depth, h_up, w_up] -> [h_up, w_up, in_depth, out_depth]
    filters_up = np.transpose(filters_up, [2, 3, 0, 1])
    self.assertEqual(np.sum(filters), np.sum(filters_up))
    return filters_up 

Example 24

def test_perform(self):
        if not imported_scipy:
            raise SkipTest('kron tests need the scipy package to be installed')

        for shp0 in [(2,), (2, 3), (2, 3, 4), (2, 3, 4, 5)]:
            x = tensor.tensor(dtype='floatX',
                              broadcastable=(False,) * len(shp0))
            a = numpy.asarray(self.rng.rand(*shp0)).astype(config.floatX)
            for shp1 in [(6,), (6, 7), (6, 7, 8), (6, 7, 8, 9)]:
                if len(shp0) + len(shp1) == 2:
                    continue
                y = tensor.tensor(dtype='floatX',
                                  broadcastable=(False,) * len(shp1))
                f = function([x, y], kron(x, y))
                b = self.rng.rand(*shp1).astype(config.floatX)
                out = f(a, b)
                # Newer versions of scipy want 4 dimensions at least,
                # so we have to add a dimension to a and flatten the result.
                if len(shp0) + len(shp1) == 3:
                    scipy_val = scipy.linalg.kron(
                        a[numpy.newaxis, :], b).flatten()
                else:
                    scipy_val = scipy.linalg.kron(a, b)
                utt.assert_allclose(out, scipy_val) 

Example 25

def testCholesky(self):
    # Tests the cholesky function
    np.random.seed(8)

    # generating two symmetric positive-definite tt-cores
    L_1 = np.tril(np.random.normal(scale=2., size=(2, 2)))
    L_2 = np.tril(np.random.normal(scale=2., size=(3, 3)))
    K_1 = L_1.dot(L_1.T)
    K_2 = L_2.dot(L_2.T)
    K = np.kron(K_1, K_2)
    initializer = tensor_train.TensorTrain([K_1[None, :, :, None], 
                                            K_2[None, :, :, None]], 
                                            tt_ranks=7*[1])
    kron_mat = variables.get_variable('kron_mat', initializer=initializer)
    init_op = tf.global_variables_initializer()
    with self.test_session() as sess:
      sess.run(init_op)
      desired = np.linalg.cholesky(K)
      actual = ops.full(kr.cholesky(kron_mat)).eval()
      self.assertAllClose(desired, actual) 

Example 26

def _iter_contrasts(n_subjects, factor_levels, effect_picks):
    """ Aux Function: Setup contrasts """
    from scipy.signal import detrend
    sc = []
    n_factors = len(factor_levels)
    # prepare computation of Kronecker products
    for n_levels in factor_levels:
        # for each factor append
        # 1) column vector of length == number of levels,
        # 2) square matrix with diagonal == number of levels

        # main + interaction effects for contrasts
        sc.append([np.ones([n_levels, 1]),
                   detrend(np.eye(n_levels), type='constant')])

    for this_effect in effect_picks:
        contrast_idx = _get_contrast_indices(this_effect + 1, n_factors)
        c_ = sc[0][contrast_idx[n_factors - 1]]
        for i_contrast in range(1, n_factors):
            this_contrast = contrast_idx[(n_factors - 1) - i_contrast]
            c_ = np.kron(c_, sc[i_contrast][this_contrast])
        df1 = matrix_rank(c_)
        df2 = df1 * (n_subjects - 1)
        yield c_, df1, df2 

Example 27

def test_homoskedastic_direct(cov_data, debias):
    x, z, eps, sigma = cov_data
    cov = HomoskedasticCovariance(x, eps, sigma, sigma, debiased=debias)
    k = len(x)
    nobs = x[0].shape[0]
    big_x = []
    for i in range(k):
        row = []
        for j in range(k):
            if i == j:
                row.append(x[i])
            else:
                row.append(np.zeros((nobs, x[j].shape[1])))
        big_x.append(np.concatenate(row, 1))
    big_x = np.concatenate(big_x, 0)
    xeex = big_x.T @ np.kron(sigma, np.eye(nobs)) @ big_x / nobs
    xpxi = _xpxi(x)
    direct = xpxi @ xeex @ xpxi / nobs
    direct = (direct + direct.T) / 2
    assert_allclose(np.diag(direct), np.diag(cov.cov))
    s = np.sqrt(np.diag(direct))[:, None]
    r_direct = direct / (s @ s.T)
    s = np.sqrt(np.diag(cov.cov))[:, None]
    c_direct = direct / (s @ s.T)
    assert_allclose(r_direct, c_direct, atol=1e-5) 

Example 28

def test_inner_product_short_circuit(data):
    y, x, sigma = data
    sigma = np.eye(len(x))
    efficient = blocked_inner_prod(x, sigma)
    nobs = x[0].shape[0]
    omega = np.kron(sigma, np.eye(nobs))
    k = len(x)
    bigx = []
    for i in range(k):
        row = []
        for j in range(k):
            if i == j:
                row.append(x[i])
            else:
                row.append(np.zeros((nobs, x[j].shape[1])))
        bigx.append(np.hstack(row))
    bigx = np.vstack(bigx)
    expected = bigx.T @ omega @ bigx
    assert_allclose(efficient, expected) 

Example 29

def test_diag_product(data):
    y, x, sigma = data
    efficient = blocked_diag_product(x, sigma)
    nobs = x[0].shape[0]
    omega = np.kron(sigma, np.eye(nobs))
    k = len(x)
    bigx = []
    for i in range(k):
        row = []
        for j in range(k):
            if i == j:
                row.append(x[i])
            else:
                row.append(np.zeros((nobs, x[j].shape[1])))
        bigx.append(np.hstack(row))
    bigx = np.vstack(bigx)
    expected = omega @ bigx
    assert_allclose(efficient, expected) 

Example 30

def test_quantumnumbers_ordinary():
    print 'test_quantumnumbers_ordinary'
    a=QuantumNumbers('C',([SQN(-1.0),SQN(0.0),SQN(1.0)],[1,1,1]),protocol=QuantumNumbers.COUNTS)
    print 'a: %s'%a
    print 'copy(a): %s'%copy(a)
    print 'deepcopy(a): %s'%deepcopy(a)

    b,permutation=QuantumNumbers.kron([a]*2).sort(history=True)
    print 'b: ',b
    print 'permutation of b:%s'%permutation
    print 'b.reorder(permutation,protocol="EXPANSION"): \n%s'%b.reorder(permutation,protocol="EXPANSION")
    print 'b.reorder([4,3,2,1,0],protocol="CONTENTS"): \n%s'%b.reorder([4,3,2,1,0],protocol="CONTENTS")

    c=b.to_ordereddict(protocol=QuantumNumbers.COUNTS)
    print 'c(b.to_ordereddict(protocol=QuantumNumbers.COUNTS)):\n%s'%('\n'.join('%s: %s'%(key,value) for key,value in c.iteritems()))
    print 'QuantumNumbers.from_ordereddict(SQN,c,protocol=QuantumNumbers.COUNTS):\n%s'%QuantumNumbers.from_ordereddict(SQN,c,protocol=QuantumNumbers.COUNTS)

    d=b.to_ordereddict(protocol=QuantumNumbers.INDPTR)
    print 'd(b.to_ordereddict(protocol=QuantumNumbers.INDPTR)):\n%s'%('\n'.join('%s: %s'%(key,value)for key,value in d.iteritems()))
    print 'QuantumNumbers.from_ordereddict(SQN,d,protocol=QuantumNumbers.INDPTR):\n%s'%QuantumNumbers.from_ordereddict(SQN,d,protocol=QuantumNumbers.INDPTR)
    print 

Example 31

def hmm_trans_matrix(self):
        # NOTE: more general version, allows different delays, o/w we could
        # construct with np.kron
        if self._hmm_trans_matrix is None:
            ps, delays = map(np.array,zip(*[(d.p,d.delay) for d in self.dur_distns]))
            starts, ends = cumsum(delays,strict=True), cumsum(delays,strict=False)
            trans_matrix = self._hmm_trans_matrix = np.zeros((ends[-1],ends[-1]))

            for (i,j), Aij in np.ndenumerate(self.trans_matrix):
                block = trans_matrix[starts[i]:ends[i],starts[j]:ends[j]]
                if i == j:
                    block[:-1,1:] = np.eye(block.shape[0]-1)
                    block[-1,-1] = 1-ps[i]
                else:
                    block[-1,0] = ps[j]*Aij

        return self._hmm_trans_matrix 

Example 32

def _initR(X, A, lmbdaR):
	_log.debug('Initializing R (SVD) lambda R: %s' % str(lmbdaR))
	rank = A.shape[1]
	U, S, Vt = svd(A, full_matrices=False)
	Shat = kron(S, S)
	Shat = (Shat / (Shat ** 2 + lmbdaR)).reshape(rank, rank)
	R = []
	ep = 1e-9
	for i in range(len(X)): # parallelize
		Rn = Shat * dot(U.T, X[i].dot(U))
		Rn = dot(Vt.T, dot(Rn, Vt))
 
		negativeVal = Rn.min()
		Rn.__iadd__(-negativeVal+ep)
		# if Rn.min() < 0 :
		# 	print("Negative Rn!")
		# 	raw_input("Press Enter: ")
		# Rn = np.eye(rank)
		# Rn = dot(A.T,A)

		R.append(Rn)
	print('Initialized R')
	return R

# ------------------ Update V ------------------------------------------------ 

Example 33

def concurrence(self, rho):
		"""Compute the concurrence of the density matrix.

		:param numpy_array rho: Density matrix

		:return: The concurrence, see https://en.wikipedia.org/wiki/Concurrence_(quantum_computing).
		:rtype: complex
		"""

		rhoTilde = np.dot( np.dot(np.kron(self.PAULI[1], self.PAULI[1]) , rho.conj()) , np.kron(self.PAULI[1], self.PAULI[1]))

		rhoRoot = scipy.linalg.fractional_matrix_power(rho, 1/2.0)

		R = scipy.linalg.fractional_matrix_power(np.dot(np.dot(rhoRoot,rhoTilde),rhoRoot),1/2.0)

		eigValues, eigVecors = np.linalg.eig(R)
		sortedEigValues = np.sort(eigValues)

		con = sortedEigValues[3]-sortedEigValues[2]-sortedEigValues[1]-sortedEigValues[0]
		return np.max([con,0]) 

Example 34

def pauli_mpp(nr_sites, local_dim):
    r"""Pauli POVM tensor product as MP-POVM

    The resulting MP-POVM will contain all tensor products of the
    elements of the local Pauli POVM from :func:`mpp.pauli_povm()`.

    :param int nr_sites: Number of sites of the returned MP-POVM
    :param int local_dim: Local dimension
    :rtype: MPPovm

    For example, for two qubits the `(1, 3)` measurement outcome is
    `minus X` on the first and `minus Y` on the second qubit:

    >>> nr_sites = 2
    >>> local_dim = 2
    >>> pauli = pauli_mpp(nr_sites, local_dim)
    >>> xy = np.kron([1, -1], [1, -1j]) / 2
    >>> xyproj = np.outer(xy, xy.conj())
    >>> proj = pauli.get([1, 3], astype=mp.MPArray) \
    ...             .to_array_global().reshape((4, 4))
    >>> abs(proj - xyproj / 3**nr_sites).max() <= 1e-10
    True

    The prefactor `1 / 3**nr_sites` arises because X, Y and Z are in a
    single POVM.

    """
    from mpnum.povm import pauli_povm
    return MPPovm.from_local_povm(pauli_povm(local_dim), nr_sites) 

Example 35

def mkron(*args):
    """np.kron() with an arbitrary number of n >= 1 arguments"""
    if len(args) == 1:
        return args[0]
    return mkron(np.kron(args[0], args[1]), *args[2:]) 

Example 36

def test_partialdot(nr_sites, local_dim, rank, rgen, dtype):
    # Only for at least two sites, we can apply an operator to a part
    # of a chain.
    if nr_sites < 2:
        return
    part_sites = nr_sites // 2
    start_at = min(2, nr_sites // 2)

    mpo = factory.random_mpa(nr_sites, (local_dim, local_dim), rank,
                             randstate=rgen, dtype=dtype)
    op = mpo.to_array_global().reshape((local_dim**nr_sites,) * 2)
    mpo_part = factory.random_mpa(part_sites, (local_dim, local_dim), rank,
                                  randstate=rgen, dtype=dtype)
    op_part = mpo_part.to_array_global().reshape((local_dim**part_sites,) * 2)
    op_part_embedded = np.kron(
        np.kron(np.eye(local_dim**start_at), op_part),
        np.eye(local_dim**(nr_sites - part_sites - start_at)))

    prod1 = np.dot(op, op_part_embedded)
    prod2 = np.dot(op_part_embedded, op)
    prod1_mpo = mp.partialdot(mpo, mpo_part, start_at=start_at)
    prod2_mpo = mp.partialdot(mpo_part, mpo, start_at=start_at)
    prod1_mpo = prod1_mpo.to_array_global().reshape((local_dim**nr_sites,) * 2)
    prod2_mpo = prod2_mpo.to_array_global().reshape((local_dim**nr_sites,) * 2)

    assert_array_almost_equal(prod1, prod1_mpo)
    assert_array_almost_equal(prod2, prod2_mpo)
    assert prod1_mpo.dtype == dtype
    assert prod2_mpo.dtype == dtype 

Example 37

def state_from_string(qubit_state_string):
		if not all(x in '01' for x in qubit_state_string):
			raise Exception("Description must be a string in binary")
		state=None
		for qubit in qubit_state_string:
			if qubit=='0':
				new_contrib=State.zero_state
			elif qubit=='1':
				new_contrib=State.one_state
			if state==None:
				state=new_contrib
			else:
				state=np.kron(state,new_contrib)
		return state 

Example 38

def separate_state(qubit_state):
		"""This only works if the state is fully separable at present

		Throws exception if not a separable state"""
		n_entangled=QuantumRegister.num_qubits(qubit_state)
		if list(qubit_state.flat).count(1)==1:
			separated_state=[]
			idx_state=list(qubit_state.flat).index(1)
			add_factor=0
			size=qubit_state.shape[0]
			while(len(separated_state)<n_entangled):
				size=size/2
				if idx_state<(add_factor+size):
					separated_state+=[State.zero_state]
					add_factor+=0
				else:
					separated_state+=[State.one_state]
					add_factor+=size
			return separated_state
		else:
			# Try a few naive separations before giving up
			cardinal_states=[State.zero_state,State.one_state,State.plus_state,State.minus_state,State.plusi_state,State.minusi_state]
			for separated_state in itertools.product(cardinal_states, repeat=n_entangled):
				candidate_state=reduce(lambda x,y:np.kron(x,y),separated_state)
				if np.allclose(candidate_state,qubit_state):
					return separated_state
			# TODO: more general separation methods
			raise StateNotSeparableException("TODO: Entangled qubits not represented yet in quantum computer implementation. Can currently do manual calculations; see TestBellState for Examples") 

Example 39

def apply_two_qubit_gate_CNOT(self,first_qubit_name,second_qubit_name):
		""" Should work for all combination of qubits"""
		first_qubit=self.qubits.get_quantum_register_containing(first_qubit_name)
		second_qubit=self.qubits.get_quantum_register_containing(second_qubit_name)
		if len(first_qubit.get_noop())>0 or len(second_qubit.get_noop())>0:
			raise Exception("Control or target qubit has been measured previously, no more gates allowed")
		if not first_qubit.is_entangled() and not second_qubit.is_entangled():
			combined_state=np.kron(first_qubit.get_state(),second_qubit.get_state())
			if first_qubit.get_num_qubits()!=1 or second_qubit.get_num_qubits()!=1:
				raise Exception("Both qubits are marked as not entangled but one or the other has an entangled state")
			new_state=Gate.CNOT2_01*combined_state
			if State.is_fully_separable(new_state):
				second_qubit.set_state(State.get_second_qubit(new_state))
			else:
				self.qubits.entangle_quantum_registers(first_qubit,second_qubit)
				first_qubit.set_state(new_state)
		else:
			if not first_qubit.is_entangled_with(second_qubit):
				# Entangle the state
				combined_state=np.kron(first_qubit.get_state(),second_qubit.get_state())
				self.qubits.entangle_quantum_registers(first_qubit,second_qubit)
			else:
				# We are ready to do the operation
				combined_state=first_qubit.get_state()
			# Time for more meta programming!
			# Select gate based on indices
			control_qubit_idx,target_qubit_idx=first_qubit.get_indices(second_qubit)
			gate_size=QuantumRegister.num_qubits(combined_state)
			try:
				exec 'gate=Gate.CNOT%d_%d%d' %(gate_size,control_qubit_idx,target_qubit_idx) 
			except:
				print 'gate=Gate.CNOT%d_%d%d' %(gate_size,control_qubit_idx,target_qubit_idx)
				raise Exception("Unrecognized combination of number of qubits")
			first_qubit.set_state(gate*combined_state) 

Example 40

def computeMisfit2(self,q):

        assert self.r0 is not None, "Must have current estimate of polarizations"
        assert self.dunc is not None, "Must have set uncertainties"
        assert self.dobs is not None, "Must have observed data"

        dunc = mkvc(self.dunc)
        dobs = mkvc(self.dobs)
        r0 = self.r0

        Hp = self.computeHp(r0=r0,update=False)
        Brx = self.computeBrx(r0=r0,update=False)
        P = self.computeP(Hp,Brx)

        N = np.size(dobs)
        M = len(self.times)

        Px = np.kron(np.diag(np.ones(M)),P)

        dpre = np.dot(Px,q)

        v = (dpre-dobs)/dunc

        Phi = np.dot(v.T,v)

        return Phi/N 

Example 41

def updatePolarizationsQP(self,r0,UB):

        # Set operator and solution array
        Hp = self.computeHp(r0=r0)
        Brx = self.computeBrx(r0=r0)
        P = self.computeP(Hp,Brx)
        dunc = self.dunc
        dobs = self.dobs
        
        K = np.shape(dobs)[1]
        q = np.zeros((6,K))
        
        # Bounds
        ub = UB*np.ones(6)
        e = matrix(np.r_[np.zeros(9),ub])

        # Constraints
        C1 = np.r_[np.c_[-1,0,0,0,0,0],np.c_[0,0,0,-1,0,0],np.c_[0,0,0,0,0,-1]]
        C2 = np.r_[np.c_[-0.5,1,0,-0.5,0,0],np.c_[-0.5,0,1,0,0,-0.5],np.c_[0,0,0,-0.5,1,-0.5]]
        C3 = np.r_[np.c_[-0.5,-1,0,-0.5,0,0],np.c_[-0.5,0,-1,0,0,-0.5],np.c_[0,0,0,-0.5,-1,-0.5]]

        # G = sp.kron(sp.diags([-np.ones(K),np.ones(K)],[0,1],shape=(K-1,K)),sp.diags(np.ones(6)))
        I = np.diag(np.ones(6))
        G = np.r_[C1,C2,C3,I]
        G = matrix(G)

        solvers.options['show_progress'] = False
        
        for kk in range(0,K):
            
            LHS = P/np.c_[dunc[:,kk],dunc[:,kk],dunc[:,kk],dunc[:,kk],dunc[:,kk],dunc[:,kk]]
            RHS = dobs[:,kk]/dunc[:,kk]

            A = matrix(np.dot(LHS.T,LHS))
            b = matrix(np.dot(LHS.T,-RHS))
            
            sol = solvers.qp(A, b, G=G, h=e)

            q[:,kk] = np.reshape(sol['x'],(6))

        self.q = q 

Example 42

def updatePolarizationsQP(self,r0,UB):

        # Set operator and solution array
        Hp = self.computeHp(r0=r0)
        Brx = self.computeBrx(r0=r0)
        P = self.computeP(Hp,Brx)
        dunc = self.dunc
        dobs = self.dobs
        
        K = np.shape(dobs)[1]
        q = np.zeros((6,K))
        
        # Bounds
        ub = UB*np.ones(6)
        e = matrix(np.r_[np.zeros(9),ub])

        # Constraints
        C1 = np.r_[np.c_[-1,0,0,0,0,0],np.c_[0,0,0,-1,0,0],np.c_[0,0,0,0,0,-1]]
        C2 = np.r_[np.c_[-0.5,1,0,-0.5,0,0],np.c_[-0.5,0,1,0,0,-0.5],np.c_[0,0,0,-0.5,1,-0.5]]
        C3 = np.r_[np.c_[-0.5,-1,0,-0.5,0,0],np.c_[-0.5,0,-1,0,0,-0.5],np.c_[0,0,0,-0.5,-1,-0.5]]

        # G = sp.kron(sp.diags([-np.ones(K),np.ones(K)],[0,1],shape=(K-1,K)),sp.diags(np.ones(6)))
        I = np.diag(np.ones(6))
        G = np.r_[C1,C2,C3,I]
        G = matrix(G)

        solvers.options['show_progress'] = False
        
        for kk in range(0,K):
            
            LHS = P/np.c_[dunc[:,kk],dunc[:,kk],dunc[:,kk],dunc[:,kk],dunc[:,kk],dunc[:,kk]]
            RHS = dobs[:,kk]/dunc[:,kk]

            A = matrix(np.dot(LHS.T,LHS))
            b = matrix(np.dot(LHS.T,-RHS))
            
            sol = solvers.qp(A, b, G=G, h=e)

            q[:,kk] = np.reshape(sol['x'],(6))

        self.q = q 

Example 43

def area(self):
        """Face areas"""
        if getattr(self, '_area', None) is None:
            if self.nCy > 1:
                raise NotImplementedError('area not yet implemented for 3D '
                                          'cyl mesh')
            areaR = np.kron(self.hz, 2*pi*self.vectorNx)
            areaZ = np.kron(np.ones_like(self.vectorNz), pi*(self.vectorNx**2 -
                            np.r_[0, self.vectorNx[:-1]]**2))
            self._area = np.r_[areaR, areaZ]
        return self._area 

Example 44

def vol(self):
        """Volume of each cell"""
        if getattr(self, '_vol', None) is None:
            if self.nCy > 1:
                raise NotImplementedError('vol not yet implemented for 3D '
                                          'cyl mesh')
            az = pi*(self.vectorNx**2 - np.r_[0, self.vectorNx[:-1]]**2)
            self._vol = np.kron(self.hz, az)
        return self._vol

    ####################################################
    # Operators
    #################################################### 

Example 45

def aveF2CC(self):
        "Construct the averaging operator on cell faces to cell centers."
        if getattr(self, '_aveF2CC', None) is None:
            n = self.vnC
            if self.isSymmetric:
                avR = utils.av(n[0])[:, 1:]
                avR[0, 0] = 1.
                self._aveF2CC = ((0.5)*sp.hstack((sp.kron(utils.speye(n[2]),
                                                          avR),
                                                  sp.kron(utils.av(n[2]),
                                                          utils.speye(n[0]))),
                                                 format="csr"))
            else:
                raise NotImplementedError('wrapping in the averaging is not '
                                          'yet implemented')
                # self._aveF2CC = (1./3.)*sp.hstack((utils.kron3(utils.speye(n[2]),
                #                                                utils.speye(n[1]),
                #                                                utils.av(n[0])),
                #                                    utils.kron3(utils.speye(n[2]),
                #                                                utils.av(n[1]),
                #                                                utils.speye(n[0])),
                #                                    utils.kron3(utils.av(n[2]),
                #                                                utils.speye(n[1]),
                #                                                utils.speye(n[0]))),
                #                                   format="csr")
        return self._aveF2CC 

Example 46

def normalize2D(x):
    return x/np.kron(np.ones((1, 2)), utils.mkvc(length2D(x), 2)) 

Example 47

def normalize3D(x):
    return x/np.kron(np.ones((1, 3)), utils.mkvc(length3D(x), 2)) 

Example 48

def test_kron_matrix(self, level=rlevel):
        # Ticket #71
        x = np.matrix('[1 0; 1 0]')
        assert_equal(type(np.kron(x, x)), type(x)) 

Example 49

def kron(a, b):
    """Returns the kronecker product of two arrays.

    Args:
        a (~cupy.ndarray): The first argument.
        b (~cupy.ndarray): The second argument.

    Returns:
        ~cupy.ndarray: Output array.

    .. seealso:: :func:`numpy.kron`

    """
    a_ndim = a.ndim
    b_ndim = b.ndim
    if a_ndim == 0 or b_ndim == 0:
        return cupy.multiply(a, b)

    ndim = b_ndim
    a_shape = a.shape
    b_shape = b.shape
    if a_ndim != b_ndim:
        if b_ndim > a_ndim:
            a_shape = (1,) * (b_ndim - a_ndim) + a_shape
        else:
            b_shape = (1,) * (a_ndim - b_ndim) + b_shape
            ndim = a_ndim

    axis = ndim - 1
    out = core.tensordot_core(a, b, None, a.size, b.size, 1, a_shape + b_shape)
    for _ in six.moves.range(ndim):
        out = core.concatenate_method(out, axis=axis)

    return out 

Example 50

def selftest1():
    
    # Generate data
    D0 = 5
    K1 = 2
    K2 = 2
    NperClass = 500
    N = NperClass*K1*K2
    #l = 1.0e-3
    X = np.zeros((D0,NperClass,K1,K2))
    y1 = np.zeros((NperClass,K1,K2),dtype=int)
    y2 = np.zeros((NperClass,K1,K2),dtype=int)
    bias1 = np.random.normal(scale=1.0,size=(D0,K1))
    bias2 = np.random.normal(scale=1.0,size=(D0,K2))    
    for k1 in range(K1):
        for k2 in range(K2):
            X[:,:,k1,k2] = \
                np.random.normal(scale=0.25, size=(D0,NperClass)) \
                + np.kron(np.ones((1,NperClass)),bias1[:,k1].reshape((D0,1))) \
                + np.kron(np.ones((1,NperClass)),bias2[:,k2].reshape((D0,1)))
            y1[:,k1,k2] = k1*np.ones((NperClass,))
            y2[:,k1,k2] = k2*np.ones((NperClass,))

    X = X.reshape((D0,N))
    y1 = y1.reshape((N,))
    y2 = y2.reshape((N,))
    
    E,dd = run(X,y1,y2)
    print dd

    plt.figure(1)
    plt.clf()
    plt.subplot(1,2,1)
    plt.plot(dd)
    plt.subplot(1,2,2)
    plt.imshow(E, aspect='auto', interpolation='none')
    plt.colorbar()
    plt.show() 
点赞

发表评论

电子邮件地址不会被公开。 必填项已用*标注