Python 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 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 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 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 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