Python numpy.kron() Examples

The following are 30 code examples for showing how to use numpy.kron(). These examples are extracted from open source projects. You can vote up the ones you like or vote down the ones you don't like, and go to the original project or source file by following the links above each example.

You may check out the related API usage on the sidebar.

You may also want to check out all available functions/classes of the module numpy , or try the search function .

Example 1
Project: OpenFermion-Cirq   Author: quantumlib   File: mfopt.py    License: Apache License 2.0 6 votes vote down vote up
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 
Example 2
Project: OpenFermion-Cirq   Author: quantumlib   File: analysis.py    License: Apache License 2.0 6 votes vote down vote up
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 
Example 3
Project: OpenFermion-Cirq   Author: quantumlib   File: simulate_trotter_test.py    License: Apache License 2.0 6 votes vote down vote up
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 
Example 4
Project: pyGSTi   Author: pyGSTio   File: basistools.py    License: Apache License 2.0 6 votes vote down vote up
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 
Example 5
Project: pyGSTi   Author: pyGSTio   File: optools.py    License: Apache License 2.0 6 votes vote down vote up
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)) 
Example 6
Project: sparselandtools   Author: fubel   File: utils.py    License: MIT License 6 votes vote down vote up
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) 
Example 7
Project: sparselandtools   Author: fubel   File: utils.py    License: MIT License 6 votes vote down vote up
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) 
Example 8
Project: tenpy   Author: tenpy   File: tdvp_numpy.py    License: GNU General Public License v3.0 6 votes vote down vote up
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 
Example 9
Project: tenpy   Author: tenpy   File: b_model.py    License: GNU General Public License v3.0 6 votes vote down vote up
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) 
Example 10
Project: pylops   Author: equinor   File: test_kronecker.py    License: GNU Lesser General Public License v3.0 6 votes vote down vote up
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) 
Example 11
Project: classification-of-encrypted-traffic   Author: SalikLP   File: vis_utils.py    License: MIT License 6 votes vote down vote up
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) 
Example 12
Project: classification-of-encrypted-traffic   Author: SalikLP   File: vis_utils.py    License: MIT License 6 votes vote down vote up
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 
Example 13
Project: D-VAE   Author: muhanzhang   File: test_slinalg.py    License: MIT License 6 votes vote down vote 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) 
Example 14
Project: vnpy_crypto   Author: birforce   File: vecm.py    License: MIT License 6 votes vote down vote up
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) 
Example 15
Project: vnpy_crypto   Author: birforce   File: irf.py    License: MIT License 6 votes vote down vote up
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 
Example 16
Project: vnpy_crypto   Author: birforce   File: irf.py    License: MIT License 6 votes vote down vote up
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) 
Example 17
Project: vnpy_crypto   Author: birforce   File: var_model.py    License: MIT License 6 votes vote down vote up
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 
Example 18
Project: vnpy_crypto   Author: birforce   File: test_phreg.py    License: MIT License 6 votes vote down vote up
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() 
Example 19
Project: vnpy_crypto   Author: birforce   File: test_gee_glm.py    License: MIT License 6 votes vote down vote up
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() 
Example 20
Project: vnpy_crypto   Author: birforce   File: test_gee.py    License: MIT License 6 votes vote down vote up
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() 
Example 21
Project: vnpy_crypto   Author: birforce   File: test_gee.py    License: MIT License 6 votes vote down vote up
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]) 
Example 22
Project: vnpy_crypto   Author: birforce   File: test_gee.py    License: MIT License 6 votes vote down vote up
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) 
Example 23
Project: fenics-topopt   Author: zfergus   File: problem.py    License: MIT License 5 votes vote down vote up
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() 
Example 24
Project: fenics-topopt   Author: zfergus   File: problem.py    License: MIT License 5 votes vote down vote up
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() 
Example 25
Project: pyscf   Author: pyscf   File: uhf.py    License: Apache License 2.0 5 votes vote down vote up
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 
Example 26
Project: pyGSTi   Author: pyGSTio   File: cloudnoisemodel.py    License: Apache License 2.0 5 votes vote down vote up
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 
Example 27
Project: pyGSTi   Author: pyGSTio   File: spamvec.py    License: Apache License 2.0 5 votes vote down vote up
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() 
Example 28
Project: pyGSTi   Author: pyGSTio   File: spamvec.py    License: Apache License 2.0 5 votes vote down vote up
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")) 
Example 29
Project: pyGSTi   Author: pyGSTio   File: spamvec.py    License: Apache License 2.0 5 votes vote down vote up
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")) 
Example 30
Project: pyGSTi   Author: pyGSTio   File: spamvec.py    License: Apache License 2.0 5 votes vote down vote up
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"))