Python numpy.kron() Examples

def non_redundant_rotation_generators(
        rhf_objective: RestrictedHartreeFockObjective) -> List[FermionOperator]:  # testpragma: no cover
    # coverage: ignore
    """
    Generate the fermionic representation of all non-redundant rotation
    generators for restricted Hartree-fock

    :param rhf_objective: openfermioncirq.experiments.hfvqe.RestrictedHartreeFock object
    :return: list of fermionic generators.
    """
    rotation_generators = []
    for p in range(rhf_objective.nocc * rhf_objective.nvirt):
        grad_params = np.zeros(rhf_objective.nocc * rhf_objective.nvirt)
        grad_params[p] = 1
        kappa_spatial_orbital = rhf_params_to_matrix(grad_params,
                                                     len(rhf_objective.occ) + len(rhf_objective.virt),
                                                     rhf_objective.occ,
                                                     rhf_objective.virt)
        p0 = np.array([[1, 0], [0, 1]])
        kappa_spin_orbital = np.kron(kappa_spatial_orbital, p0)
        fermion_op = get_one_body_fermion_operator(kappa_spin_orbital)
        rotation_generators.append(fermion_op)
    return rotation_generators 

def energy_from_opdm(opdm, constant, one_body_tensor, two_body_tensor):
    """
    Evaluate the energy of an opdm assuming the 2-RDM is opdm ^ opdm

    :param opdm: single spin-component of the full spin-orbital opdm.
    :param constant: constant shift to the Hamiltonian. Commonly this is the
                     nuclear repulsion energy.
    :param one_body_tensor: spatial one-body integrals
    :param two_body_tensor: spatial two-body integrals
    :return:
    """
    spin_opdm = np.kron(opdm, np.eye(2))
    spin_tpdm = 2 * wedge(spin_opdm, spin_opdm, (1, 1), (1, 1))
    molecular_hamiltonian = generate_hamiltonian(constant=constant,
                                                 one_body_integrals=one_body_tensor,
                                                 two_body_integrals=two_body_tensor)
    rdms = InteractionRDM(spin_opdm, spin_tpdm)
    return rdms.expectation(molecular_hamiltonian).real 

def test_simulate_trotter_simulate_controlled(
        hamiltonian, time, initial_state, exact_state, order, n_steps,
        algorithm, result_fidelity):

    n_qubits = openfermion.count_qubits(hamiltonian)
    qubits = cirq.LineQubit.range(n_qubits)

    control = cirq.LineQubit(-1)
    zero = [1, 0]
    one = [0, 1]
    start_state = (numpy.kron(zero, initial_state)
                   + numpy.kron(one, initial_state)) / numpy.sqrt(2)

    circuit = cirq.Circuit(simulate_trotter(
        qubits, hamiltonian, time, n_steps, order, algorithm, control))

    final_state = circuit.final_wavefunction(start_state)
    correct_state = (numpy.kron(zero, initial_state)
                     + numpy.kron(one, exact_state)) / numpy.sqrt(2)
    assert fidelity(final_state, correct_state) > result_fidelity
    # Make sure the time wasn't too small
    assert fidelity(final_state, start_state) < 0.95 * result_fidelity 

def state_to_stdmx(state_vec):
    """
    Convert a state vector into a density matrix.

    Parameters
    ----------
    state_vec : list or tuple
       State vector in the standard (sigma-z) basis.

    Returns
    -------
    numpy.ndarray
        A density matrix of shape (d,d), corresponding to the pure state
        given by the length-`d` array, `state_vec`.
    """
    st_vec = state_vec.view(); st_vec.shape = (len(st_vec), 1)  # column vector
    dm_mx = _np.kron(_np.conjugate(_np.transpose(st_vec)), st_vec)
    return dm_mx  # density matrix in standard (sigma-z) basis 

def unitary_to_process_mx(U):
    """
    Compute the super-operator which acts on (row)-vectorized
    density matrices from a unitary operator (matrix) U which
    acts on state vectors.  This super-operator is given by
    the tensor product of U and conjugate(U), i.e. kron(U,U.conj).

    Parameters
    ----------
    U : numpy array
        The unitary matrix which acts on state vectors.

    Returns
    -------
    numpy array
       The super-operator process matrix.
    """
    # U -> kron(U,Uc) since U rho U_dag -> kron(U,Uc)
    #  since AXB --row-vectorize--> kron(A,B.T)*vec(X)
    return _np.kron(U, _np.conjugate(U)) 

def dictionary_from_transform(transform, n, K, normalized=True, inverse=True):
    """
    Builds a Dictionary matrix from a given transform

    Args:
        transform: A valid transform (e.g. Haar, DCT-II)
        n: number of rows transform dictionary
        K: number of columns transform dictionary
        normalized: If True, the columns will be l2-normalized
        inverse: Uses the inverse transform (as usually needed in applications)

    Returns:
        Dictionary build from the Kronecker-Delta of the transform applied to the identity.
    """
    H = np.zeros((K, n))
    for i in range(n):
        v = np.zeros(n)
        v[i] = 1.0
        H[:, i] = transform(v, sampling_factor=K)
    if inverse:
        H = H.T
    return np.kron(H.T, H.T) 

def random_dictionary(n, K, normalized=True, seed=None):
    """
    Build a random dictionary matrix with K = n
    Args:
        n: square of signal dimension
        K: square of desired dictionary atoms
        normalized: If true, columns will be l2-normalized
        seed: Random seed

    Returns:
        Random dictionary
    """
    if seed:
        np.random.seed(seed)
    H = np.random.rand(n, K) * 255
    if normalized:
        for k in range(K):
            H[:, k] *= 1 / np.linalg.norm(H[:, k])
    return np.kron(H, H) 

def middle_bond_hamiltonian(Jx, Jz, hx, hz, L):
    """" Returns the spin operators sigma_x and sigma_z for L sites."""
    sx = np.array([[0., 1.], [1., 0.]])
    sz = np.array([[1., 0.], [0., -1.]])

    H_bond = Jx * np.kron(sx, sx) + Jz * np.kron(sz, sz)
    H_bond = H_bond + hx / 2 * np.kron(sx, np.eye(2)) + hx / 2 * np.kron(np.eye(2), sx)
    H_bond = H_bond + hz / 2 * np.kron(sz, np.eye(2)) + hz / 2 * np.kron(np.eye(2), sz)
    H_bond = H_bond.reshape(2, 2, 2, 2).transpose(0, 2, 1, 3).reshape(4, 4)  #i1 i2 i1' i2' -->
    U, s, V = np.linalg.svd(H_bond)

    M1 = np.dot(U, np.diag(s)).reshape(2, 2, 1, 4).transpose(2, 3, 0, 1)
    M2 = V.reshape(4, 1, 2, 2)
    M0 = np.tensordot(np.tensordot([1], [1], axes=0), np.eye(2), axes=0)
    W = []

    for i in range(L):
        if i == L / 2 - 1:
            W.append(M1)
        elif i == L / 2:
            W.append(M2)
        else:
            W.append(M0)
    return W 

def init_H_bonds(self):
        """Initialize `H_bonds` hamiltonian.

        Called by __init__().
        """
        sx, sz, id = self.sigmax, self.sigmaz, self.id
        d = self.d
        nbonds = self.L - 1 if self.bc == 'finite' else self.L
        H_list = []
        for i in range(nbonds):
            gL = gR = 0.5 * self.g
            if self.bc == 'finite':
                if i == 0:
                    gL = self.g
                if i + 1 == self.L - 1:
                    gR = self.g
            H_bond = -self.J * np.kron(sx, sx) - gL * np.kron(sz, id) - gR * np.kron(id, sz)
            # H_bond has legs ``i, j, i*, j*``
            H_list.append(np.reshape(H_bond, [d, d, d, d]))
        self.H_bonds = H_list

    # (note: not required for TEBD) 

def test_Kroneker(par):
    """Dot-test, inversion and comparison with np.kron for Kronecker operator
    """
    np.random.seed(10)
    G1 = np.random.normal(0, 10, (par['ny'], par['nx'])).astype(par['dtype'])
    G2 = np.random.normal(0, 10, (par['ny'], par['nx'])).astype(par['dtype'])
    x = np.ones(par['nx']**2) + par['imag']*np.ones(par['nx']**2)

    Kop = Kronecker(MatrixMult(G1, dtype=par['dtype']),
                    MatrixMult(G2, dtype=par['dtype']),
                    dtype=par['dtype'])
    assert dottest(Kop, par['ny']**2, par['nx']**2,
                   complexflag=0 if par['imag'] == 0 else 3)

    xlsqr = lsqr(Kop, Kop * x, damp=1e-20, iter_lim=300, show=0)[0]
    assert_array_almost_equal(x, xlsqr, decimal=2)

    # Comparison with numpy
    assert_array_almost_equal(np.kron(G1, G2), Kop * np.eye(par['nx']**2),
                              decimal=3) 

def graymap(x):

    x = x[..., np.newaxis]
    return np.concatenate([x, x, x], axis=-1)*0.5+0.5

# --------------------------------------
# Visualizing data
# --------------------------------------

# def visualize(x,colormap,name):
#
#     N = len(x)
#     assert(N <= 16)
#
#     x = colormap(x/np.abs(x).max())
#
#     # Create a mosaic and upsample
#     x = x.reshape([1, N, 29, 29, 3])
#     x = np.pad(x, ((0, 0), (0, 0), (2, 2), (2, 2), (0, 0)), 'constant', constant_values=1)
#     x = x.transpose([0, 2, 1, 3, 4]).reshape([1*33, N*33, 3])
#     x = np.kron(x, np.ones([2, 2, 1]))
#
#     PIL.Image.fromarray((x*255).astype('byte'), 'RGB').save(name) 

def plt_vector(x, colormap, num_headers):
    N = len(x)
    assert (N <= 16)
    len_x = 54
    len_y = num_headers
    # size = int(np.ceil(np.sqrt(len(x[0]))))
    length = len_y*len_x
    data = np.zeros((N, length), dtype=np.float64)
    data[:, :x.shape[1]] = x
    data = colormap(data / np.abs(data).max())
    # data = data.reshape([1, N, size, size, 3])
    data = data.reshape([1, N, len_y, len_x, 3])
    # data = np.pad(data, ((0, 0), (0, 0), (2, 2), (2, 2), (0, 0)), 'constant', constant_values=1)
    data = data.transpose([0, 2, 1, 3, 4]).reshape([1 * (len_y), N * (len_x), 3])
    return data
    # data = np.kron(data, np.ones([2, 2, 1])) # scales 

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 cov_params_default(self):  # p.296 (7.2.21)
        # Sigma_co described on p. 287
        beta = self.beta
        if self.det_coef_coint.size > 0:
            beta = vstack((beta, self.det_coef_coint))
        dt = self.deterministic
        num_det = ("co" in dt) + ("lo" in dt)
        num_det += (self.seasons-1) if self.seasons else 0
        if self.exog is not None:
            num_det += self.exog.shape[1]
        b_id = scipy.linalg.block_diag(beta,
                                       np.identity(self.neqs * (self.k_ar-1) +
                                                   num_det))

        y_lag1 = self._y_lag1
        b_y = beta.T.dot(y_lag1)
        omega11 = b_y.dot(b_y.T)
        omega12 = b_y.dot(self._delta_x.T)
        omega21 = omega12.T
        omega22 = self._delta_x.dot(self._delta_x.T)
        omega = np.bmat([[omega11, omega12],
                         [omega21, omega22]]).A

        mat1 = b_id.dot(inv(omega)).dot(b_id.T)
        return np.kron(mat1, self.sigma_u) 

def _orth_cov(self):
        # Lutkepohl 3.7.8

        Ik = np.eye(self.neqs)
        PIk = np.kron(self.P.T, Ik)
        H = self.H

        covs = self._empty_covm(self.periods + 1)
        for i in range(self.periods + 1):
            if i == 0:
                apiece = 0
            else:
                Ci = np.dot(PIk, self.G[i-1])
                apiece = chain_dot(Ci, self.cov_a, Ci.T)

            Cibar = np.dot(np.kron(Ik, self.irfs[i]), H)
            bpiece = chain_dot(Cibar, self.cov_sig, Cibar.T) / self.T

            # Lutkepohl typo, cov_sig correct
            covs[i] = apiece + bpiece

        return covs 

def lr_effect_cov(self, orth=False):
        """
        Returns
        -------

        """
        lre = self.lr_effects
        Finfty = np.kron(np.tile(lre.T, self.lags), lre)
        Ik = np.eye(self.neqs)

        if orth:
            Binf = np.dot(np.kron(self.P.T, np.eye(self.neqs)), Finfty)
            Binfbar = np.dot(np.kron(Ik, lre), self.H)

            return (chain_dot(Binf, self.cov_a, Binf.T) +
                    chain_dot(Binfbar, self.cov_sig, Binfbar.T))
        else:
            return chain_dot(Finfty, self.cov_a, Finfty.T) 

def _var_acf(coefs, sig_u):
    """
    Compute autocovariance function ACF_y(h) for h=1,...,p

    Notes
    -----
    Lütkepohl (2005) p.29
    """
    p, k, k2 = coefs.shape
    assert(k == k2)

    A = util.comp_matrix(coefs)
    # construct VAR(1) noise covariance
    SigU = np.zeros((k*p, k*p))
    SigU[:k, :k] = sig_u

    # vec(ACF) = (I_(kp)^2 - kron(A, A))^-1 vec(Sigma_U)
    vecACF = scipy.linalg.solve(np.eye((k*p)**2) - np.kron(A, A), vec(SigU))

    acf = unvec(vecACF)
    acf = acf[:k].T.reshape((p, k, k))

    return acf 

def test_get_distribution(self):
        # Smoke test
        np.random.seed(34234)
        n = 200
        exog = np.random.normal(size=(n, 2))
        lin_pred = exog.sum(1)
        elin_pred = np.exp(-lin_pred)
        time = -elin_pred * np.log(np.random.uniform(size=n))
        status = np.ones(n)
        status[0:20] = 0
        strata = np.kron(range(5), np.ones(n // 5))

        mod = PHReg(time, exog, status=status, strata=strata)
        rslt = mod.fit()

        dist = rslt.get_distribution()

        fitted_means = dist.mean()
        true_means = elin_pred
        fitted_var = dist.var()
        fitted_sd = dist.std()
        sample = dist.rvs() 

def setup_class(cls):

        vs = Independence()
        family = families.Gaussian()
        np.random.seed(987126)
        Y = np.random.normal(size=100)
        X1 = np.random.normal(size=100)
        X2 = np.random.normal(size=100)
        X3 = np.random.normal(size=100)
        groups = np.kron(np.arange(20), np.ones(5))

        D = pd.DataFrame({"Y": Y, "X1": X1, "X2": X2, "X3": X3})

        md = GEE.from_formula("Y ~ X1 + X2 + X3", groups, D,
                              family=family, cov_struct=vs)
        cls.result1 = md.fit()

        cls.result2 = GLM.from_formula("Y ~ X1 + X2 + X3", data=D).fit() 

def test_margins_gaussian(self):
        # Check marginal effects for a Gaussian GEE fit.  Marginal
        # effects and ordinary effects should be equal.

        n = 40
        np.random.seed(34234)
        exog = np.random.normal(size=(n, 3))
        exog[:, 0] = 1

        groups = np.kron(np.arange(n / 4), np.r_[1, 1, 1, 1])
        endog = exog[:, 1] + np.random.normal(size=n)

        model = sm.GEE(endog, exog, groups)
        result = model.fit(
            start_params=[-4.88085602e-04, 1.18501903, 4.78820100e-02])

        marg = result.get_margeff()

        assert_allclose(marg.margeff, result.params[1:])
        assert_allclose(marg.margeff_se, result.bse[1:])

        # smoke test
        marg.summary() 

def test_constraint_covtype(self):
        # Test constraints with different cov types
        np.random.seed(6432)
        n = 200
        exog = np.random.normal(size=(n, 4))
        endog = exog[:, 0] + exog[:, 1] + exog[:, 2]
        endog += 3 * np.random.normal(size=n)
        group = np.kron(np.arange(n / 4), np.ones(4))
        L = np.array([[1., -1, 0, 0]])
        R = np.array([0., ])
        family = Gaussian()
        va = Independence()
        for cov_type in "robust", "naive", "bias_reduced":
            model = GEE(endog, exog, group, family=family,
                        cov_struct=va, constraint=(L, R))
            result = model.fit(cov_type=cov_type)
            result.standard_errors(cov_type=cov_type)
            assert_allclose(result.cov_robust.shape, np.r_[4, 4])
            assert_allclose(result.cov_naive.shape, np.r_[4, 4])
            if cov_type == "bias_reduced":
                assert_allclose(result.cov_robust_bc.shape, np.r_[4, 4]) 

def test_linear_constrained(self):

        family = Gaussian()

        exog = np.random.normal(size=(300, 4))
        exog[:, 0] = 1
        endog = np.dot(exog, np.r_[1, 1, 0, 0.2]) +\
            np.random.normal(size=300)
        group = np.kron(np.arange(100), np.r_[1, 1, 1])

        vi = Independence()
        ve = Exchangeable()

        L = np.r_[[[0, 0, 0, 1]]]
        R = np.r_[0, ]

        for j, v in enumerate((vi, ve)):
            md = GEE(endog, exog, group, None, family, v,
                     constraint=(L, R))
            mdf = md.fit()
            assert_almost_equal(mdf.params[3], 0, decimal=10) 

def build_indices(self, nelx, nely):
        """ FE: Build the index vectors for the for coo matrix format. """
        self.KE = self.lk()
        self.edofMat = np.zeros((nelx * nely, 8), dtype=int)
        for elx in range(nelx):
            for ely in range(nely):
                el = ely + elx * nely
                n1 = (nely + 1) * elx + ely
                n2 = (nely + 1) * (elx + 1) + ely
                self.edofMat[el, :] = np.array([2 * n1 + 2, 2 * n1 + 3,
                    2 * n2 + 2, 2 * n2 + 3, 2 * n2, 2 * n2 + 1, 2 * n1,
                    2 * n1 + 1])
        # Construct the index pointers for the coo format
        self.iK = np.kron(self.edofMat, np.ones((8, 1))).flatten()
        self.jK = np.kron(self.edofMat, np.ones((1, 8))).flatten() 

def build_indices(self, nelx, nely):
        """ FE: Build the index vectors for the for coo matrix format. """
        self.KE = self.lk()
        self.edofMat = np.zeros((nelx * nely, 8), dtype=int)
        for elx in range(nelx):
            for ely in range(nely):
                el = ely + elx * nely
                n1 = (nely + 1) * elx + ely
                n2 = (nely + 1) * (elx + 1) + ely
                self.edofMat[el, :] = np.array([2 * n1 + 2, 2 * n1 + 3,
                    2 * n2 + 2, 2 * n2 + 3, 2 * n2, 2 * n2 + 1, 2 * n1,
                    2 * n1 + 1])
        # Construct the index pointers for the coo format
        self.iK = np.kron(self.edofMat, np.ones((8, 1))).flatten()
        self.jK = np.kron(self.edofMat, np.ones((1, 8))).flatten() 

def dia(gobj, dm0, gauge_orig=None):
    '''Note the side effects of set_common_origin'''

    if isinstance(dm0, numpy.ndarray) and dm0.ndim == 2: # RHF DM
        return numpy.zeros((3,3))
    mol = gobj.mol
    effspin = mol.spin * .5
    #muB = .5  # Bohr magneton

    dma, dmb = dm0
    #totdm = dma + dmb
    spindm = dma - dmb
    alpha2 = nist.ALPHA ** 2
    #Many choices of qed_fac, see JPC, 101, 3388
    #qed_fac = (nist.G_ELECTRON - 1)
    #qed_fac = nist.G_ELECTRON / 2
    #qed_fac = 1

    assert(not mol.has_ecp())
    if gauge_orig is not None:
        mol.set_common_origin(gauge_orig)

    e11 = numpy.zeros((3,3))
    im, mass_center = rhf_rotg.inertia_tensor(mol)
    for ia in range(mol.natm):
        Z = koseki_charge(mol.atom_charge(ia))
        R = mol.atom_coord(ia) - mass_center
        with mol.with_rinv_origin(R):
            h11 = mol.intor('int1e_drinv', comp=3)  * Z # * mol.atom_charge(ia)
            t1 = numpy.einsum('xij,ij->x', h11, spindm)
            e11 +=  numpy.dot(R, t1)*numpy.eye(3) -  numpy.kron(R, t1).reshape(3,3)

        #GIAO part of dia-magnetic constribution
        #print('kron', numpy.kron(R, t1).reshape(3,3))
            #h22 =  mol.intor('int1e_a01gp', comp=9)
            #e11 -=  Z * numpy.einsum('xij,ij->x', h22, spindm).reshape(3,3)
    gdia = e11 * alpha2 / effspin /  4.
    return gdia 

def basisProductMatrix(sigmaInds, sparse):
    """ Construct the Pauli product matrix from the given `sigmaInds` """
    sigmaVec = (id2x2 / sqrt2, sigmax / sqrt2, sigmay / sqrt2, sigmaz / sqrt2)
    M = _np.identity(1, 'complex')
    for i in sigmaInds:
        M = _np.kron(M, sigmaVec[i])
    return _sps.csr_matrix(M) if sparse else M 

def from_vector(self, v, close=False, nodirty=False):
        """
        Initialize the SPAM vector using a 1D array of parameters.

        Parameters
        ----------
        v : numpy array
            The 1D vector of gate parameters.  Length
            must == num_params()

        Returns
        -------
        None
        """
        if self._prep_or_effect == "prep":
            for sv in self.factors:
                sv.from_vector(v[sv.gpindices], close, nodirty)  # factors hold local indices

        elif all([self.effectLbls[i] == list(povm.keys())[0]
                  for i, povm in enumerate(self.factors)]):
            #then this is the *first* vector in the larger TensorProdPOVM
            # and we should initialize all of the factorPOVMs
            for povm in self.factors:
                local_inds = _modelmember._decompose_gpindices(
                    self.gpindices, povm.gpindices)
                povm.from_vector(v[local_inds], close, nodirty)

        #Update representation, which may be a dense matrix or
        # just fast-kron arrays or a stabilizer state.
        self._update_rep() 

def from_dense_purevec(cls, purevec):
        """
        Create a new StabilizerSPAMVec from a pure-state vector.

        Currently, purevec must be a single computational basis state (it
        cannot be a superpostion of multiple of them).

        Parameters
        ----------
        purevec : numpy.ndarray
            A complex-valued state vector specifying a pure state in the
            standard computational basis.  This vector has length 2^n for
            n qubits.

        Returns
        -------
        StabilizerSPAMVec
        """
        nqubits = int(round(_np.log2(len(purevec))))
        v = (_np.array([1, 0], 'd'), _np.array([0, 1], 'd'))  # (v0,v1)
        for zvals in _itertools.product(*([(0, 1)] * nqubits)):
            testvec = _functools.reduce(_np.kron, [v[i] for i in zvals])
            if _np.allclose(testvec, purevec.flat):
                return cls(nqubits, zvals)
        raise ValueError(("Given `purevec` must be a z-basis product state - "
                          "cannot construct StabilizerSPAMVec")) 

def from_dense_purevec(cls, purevec):
        """
        Create a new StabilizerEffectVec from a pure-state vector.

        Currently, purevec must be a single computational basis state (it
        cannot be a superpostion of multiple of them).

        Parameters
        ----------
        purevec : numpy.ndarray
            A complex-valued state vector specifying a pure state in the
            standard computational basis.  This vector has length 2^n for
            n qubits.

        Returns
        -------
        StabilizerSPAMVec
        """
        nqubits = int(round(_np.log2(len(purevec))))
        v = (_np.array([1, 0], 'd'), _np.array([0, 1], 'd'))  # (v0,v1)
        for zvals in _itertools.product(*([(0, 1)] * nqubits)):
            testvec = _functools.reduce(_np.kron, [v[i] for i in zvals])
            if _np.allclose(testvec, purevec.flat):
                return cls(zvals)
        raise ValueError(("Given `purevec` must be a z-basis product state - "
                          "cannot construct StabilizerEffectVec")) 

def from_dense_vec(cls, vec, evotype):
        """
        Create a new ComputationalSPAMVec from a dense vector.

        Parameters
        ----------
        vec : numpy.ndarray
            A state vector specifying a computational basis state in the
            standard basis.  This vector has length 2^n or 4^n for
            n qubits depending on `evotype`.

        evotype : {"densitymx", "statevec", "stabilizer", "svterm", "cterm"}
            The evolution type of the resulting SPAM vector.  This value
            must be consistent with `len(vec)`, in that `"statevec"` and
            `"stabilizer"` expect 2^n whereas the rest expect 4^n.

        Returns
        -------
        ComputationalSPAMVec
        """
        if evotype in ('stabilizer', 'statevec'):
            nqubits = int(round(_np.log2(len(vec))))
            v0 = _np.array((1, 0), complex)  # '0' qubit state as complex state vec
            v1 = _np.array((0, 1), complex)  # '1' qubit state as complex state vec
        else:
            nqubits = int(round(_np.log2(len(vec)) / 2))
            v0 = 1.0 / _np.sqrt(2) * _np.array((1, 0, 0, 1), 'd')  # '0' qubit state as Pauli dmvec
            v1 = 1.0 / _np.sqrt(2) * _np.array((1, 0, 0, -1), 'd')  # '1' qubit state as Pauli dmvec

        v = (v0, v1)
        for zvals in _itertools.product(*([(0, 1)] * nqubits)):
            testvec = _functools.reduce(_np.kron, [v[i] for i in zvals])
            if _np.allclose(testvec, vec.flat):
                return cls(zvals, evotype)
        raise ValueError(("Given `vec` is not a z-basis product state - "
                          "cannot construct ComputatinoalSPAMVec"))