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',
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',
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()
# 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)
size=qubit_state.shape[0]
while(len(separated_state)<n_entangled):
size=size/2
separated_state+=[State.zero_state]
else:
separated_state+=[State.one_state]
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()