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
python类kron()的实例源码
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 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()
