我们从Python开源项目中，提取了以下50个代码示例，用于说明如何使用numpy.kron()。
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())
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)))
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])
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
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
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
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
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)
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)
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.")
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
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
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
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
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
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,
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 ")
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)
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
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
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
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)
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)
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
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)
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)
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)
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
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
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 ------------------------------------------------
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])
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)
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:])
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
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
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")
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)
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
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
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
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 ####################################################
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
def normalize2D(x): return x/np.kron(np.ones((1, 2)), utils.mkvc(length2D(x), 2))
def normalize3D(x): return x/np.kron(np.ones((1, 3)), utils.mkvc(length3D(x), 2))
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))
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
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()
