diff --git a/ccsrc/lib/simulator/densitymatrix/detail/cpu_common/cpu_densitymatrix_core_policy.cpp b/ccsrc/lib/simulator/densitymatrix/detail/cpu_common/cpu_densitymatrix_core_policy.cpp index bb140a58d..460d5d3a1 100644 --- a/ccsrc/lib/simulator/densitymatrix/detail/cpu_common/cpu_densitymatrix_core_policy.cpp +++ b/ccsrc/lib/simulator/densitymatrix/detail/cpu_common/cpu_densitymatrix_core_policy.cpp @@ -252,12 +252,11 @@ auto CPUDensityMatrixPolicyBase::GetPartialTrace(const qs_ return res; } - template auto CPUDensityMatrixPolicyBase::PureStateVector(const qs_data_p_t& qs, index_t dim) -> py_qs_datas_t { auto p = Purity(qs, dim); - if (std::abs(p - 1) > 1e-8) { + if (1 - p > 1e-6) { throw(std::runtime_error("PureStateVector(): Cannot transform mixed density matrix to vector.")); } if (qs == nullptr) { diff --git a/mindquantum/core/gates/basicgate.py b/mindquantum/core/gates/basicgate.py index 0d98991da..df7baa639 100644 --- a/mindquantum/core/gates/basicgate.py +++ b/mindquantum/core/gates/basicgate.py @@ -70,6 +70,7 @@ class UnivMathGate(NoneParamNonHermMat): def __init__(self, name, matrix_value): """Initialize a UnivMathGate object.""" + _check_input_type('matrix_value', np.ndarray, matrix_value) if len(matrix_value.shape) != 2: raise ValueError(f"matrix_value require shape of 2, but get shape of {matrix_value.shape}") if matrix_value.shape[0] != matrix_value.shape[1]: diff --git a/mindquantum/core/gates/channel.py b/mindquantum/core/gates/channel.py index b661e17e7..cecdf5c36 100644 --- a/mindquantum/core/gates/channel.py +++ b/mindquantum/core/gates/channel.py @@ -621,7 +621,7 @@ def __init__(self, name: str, kraus_op, **kwargs): sum_of_mat = np.zeros((2, 2), 'complex128') for mat in kraus_op: sum_of_mat += np.dot(mat.T.conj(), mat) - if not np.allclose(sum_of_mat, [[1, 0], [0, 1]]): + if not np.allclose(sum_of_mat, [[1, 0], [0, 1]], atol=1e-6): raise ValueError(f"kraus_op need to satisfy the completeness condition, but get {sum_of_mat}") kwargs['name'] = name kwargs['n_qubits'] = 1 diff --git a/mindquantum/simulator/mqsim.py b/mindquantum/simulator/mqsim.py index 1d7bc531b..a8515f6a0 100644 --- a/mindquantum/simulator/mqsim.py +++ b/mindquantum/simulator/mqsim.py @@ -471,7 +471,7 @@ def set_qs(self, quantum_state: np.ndarray): raise ValueError("Wrong quantum state.") self.sim.set_qs(quantum_state / norm_factor) elif len(quantum_state.shape) == 2: - if not np.allclose(quantum_state, quantum_state.T.conj()): + if not np.allclose(quantum_state, quantum_state.T.conj(), atol=1e-6): raise ValueError("density matrix must be hermitian.") if (quantum_state.diagonal() < 0).any(): raise ValueError("the diagonal terms in density matrix cannot be negative.") @@ -528,6 +528,6 @@ def get_pure_state_vector(self) -> np.ndarray: """Get the state vector from a pure density matrix.""" if self.name.startswith('mqvector'): return self.get_qs() - if 1 - self.purity() > 1e-8: + if 1 - self.purity() > 1e-6: raise ValueError("Cannot transform mixed density matrix to vector.") return np.array(self.sim.pure_state_vector()) diff --git a/mindquantum/simulator/simulator.py b/mindquantum/simulator/simulator.py index e0362b9a8..bcae40069 100644 --- a/mindquantum/simulator/simulator.py +++ b/mindquantum/simulator/simulator.py @@ -679,21 +679,21 @@ def fidelity(rho: np.ndarray, sigma: np.ndarray): if np.log2(qs.shape[0]) % 1 != 0: raise ValueError(f"quantum state size {qs.shape[0]} is not power of 2") if qs.ndim == 1: - if not np.allclose(np.sum(np.abs(qs) ** 2), 1): + if not np.allclose(np.sum(np.abs(qs) ** 2), 1, atol=1e-6): raise ValueError("state vector must be normalized.") elif qs.ndim == 2: if qs.shape[0] != qs.shape[1]: raise ValueError("the row of matrix is not equal to column.") - if not np.allclose(qs, qs.T.conj()): + if not np.allclose(qs, qs.T.conj(), atol=1e-6): raise ValueError("density matrix must be hermitian.") if (qs.diagonal() < 0).any(): raise ValueError("the diagonal terms in density matrix cannot be negative.") - if not np.allclose(np.real(np.trace(qs)), 1): + if not np.allclose(np.real(np.trace(qs)), 1, atol=1e-6): raise ValueError("the trace of density matrix must equal to 1.") else: raise ValueError(f"input quantum state requires a vector or matrix, but get shape {rho.shape}.") if rho.ndim == 1 and sigma.ndim == 1: - return np.abs(np.inner(rho, sigma)) ** 2 + return np.abs(np.inner(rho.conj().T, sigma)) ** 2 if rho.ndim == 1 and sigma.ndim == 2: return np.real(rho.conj().T @ sigma @ rho) if rho.ndim == 2 and sigma.ndim == 1: diff --git a/tests/st/test_simulator/test_method_of_mqmatrix.py b/tests/st/test_simulator/test_method_of_mqmatrix.py new file mode 100644 index 000000000..0dde07e8d --- /dev/null +++ b/tests/st/test_simulator/test_method_of_mqmatrix.py @@ -0,0 +1,362 @@ +# Copyright 2021 Huawei Technologies Co., Ltd +# +# Licensed under the Apache License, Version 2.0 (the "License"); +# you may not use this file except in compliance with the License. +# You may obtain a copy of the License at +# +# http://www.apache.org/licenses/LICENSE-2.0 +# +# Unless required by applicable law or agreed to in writing, software +# distributed under the License is distributed on an "AS IS" BASIS, +# wITHOUT wARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +# See the License for the specific language governing permissions and +# limitations under the License. +# ============================================================================ + +# pylint: disable=invalid-name +"""Test simulator.""" + +import numpy as np +import pytest +from scipy.stats import entropy +from scipy.sparse import csr_matrix +from scipy.linalg import logm, sqrtm + +import mindquantum as mq +from mindquantum.core import gates as G +from mindquantum.core.circuit import UN, Circuit +from mindquantum.simulator import Simulator, fidelity +from mindquantum.core.operators import Hamiltonian, QubitOperator +from mindquantum.simulator.available_simulator import SUPPORTED_SIMULATOR +from mindquantum.utils import random_circuit + + +@pytest.mark.level0 +@pytest.mark.platform_x86_gpu_training +@pytest.mark.platform_x86_cpu +@pytest.mark.env_onecard +@pytest.mark.parametrize("virtual_qc", ['mqmatrix']) +@pytest.mark.parametrize("dtype", [mq.complex64, mq.complex128]) +def test_set_qs_and_dm(virtual_qc, dtype): + """ + Description: test setting density matrix + Expectation: success. + """ + qs = np.random.rand(2) + np.random.rand(2) * 1j + sim = Simulator(virtual_qc, 1, dtype=dtype) + sim.set_qs(qs) + qs = qs / np.linalg.norm(qs) + assert np.allclose(sim.get_qs(), np.outer(qs, qs.conj())) + + qs2 = np.random.rand(2) + np.random.rand(2) * 1j + dm = np.outer(qs, qs.conj()) + np.outer(qs2, qs2.conj()) + sim.reset() + sim.set_qs(dm) + dm = dm / np.trace(dm) + assert np.allclose(sim.get_qs(), dm, atol=1e-6) + + +@pytest.mark.level0 +@pytest.mark.platform_x86_gpu_training +@pytest.mark.platform_x86_cpu +@pytest.mark.env_onecard +@pytest.mark.parametrize("virtual_qc", ['mqmatrix']) +@pytest.mark.parametrize("dtype", [mq.complex64, mq.complex128]) +def test_get_partial_trace(virtual_qc, dtype): + """ + Description: test partial trace of density matrix + Expectation: success. + """ + circ = random_circuit(3, 100) + sim = Simulator(virtual_qc, 3, dtype=dtype) + sim.apply_circuit(circ) + mat = sim.get_partial_trace([0, 1]) + qs = sim.get_qs() + ref_mat = np.array( + [[np.trace(qs[0:4, 0:4]), np.trace(qs[0:4, 4:8])], [np.trace(qs[4:8, 0:4]), np.trace(qs[4:8, 4:8])]] + ) + assert np.allclose(mat, ref_mat, atol=1e-6) + + +@pytest.mark.level0 +@pytest.mark.platform_x86_gpu_training +@pytest.mark.platform_x86_cpu +@pytest.mark.env_onecard +@pytest.mark.parametrize("virtual_qc", ['mqmatrix']) +@pytest.mark.parametrize("dtype", [mq.complex64, mq.complex128]) +def test_purity(virtual_qc, dtype): + """ + Description: test purity of density matrix + Expectation: success. + """ + circ = random_circuit(3, 100) + circ = circ.with_noise(G.DepolarizingChannel(0.1)) + sim = Simulator(virtual_qc, 3, dtype=dtype) + sim.apply_circuit(circ) + purity = sim.purity() + ref_purity = np.trace(sim.get_qs() @ sim.get_qs()) + assert np.allclose(purity, ref_purity, atol=1e-6) + + +@pytest.mark.level0 +@pytest.mark.platform_x86_gpu_training +@pytest.mark.platform_x86_cpu +@pytest.mark.env_onecard +@pytest.mark.parametrize("virtual_qc", ['mqmatrix']) +@pytest.mark.parametrize("dtype", [mq.complex64, mq.complex128]) +def test_get_pure_state_vector(virtual_qc, dtype): + """ + Description: test get pure state vector from density matrix + Expectation: success. + """ + init_state = np.random.rand(8) + np.random.rand(8) * 1j + init_state = init_state / np.linalg.norm(init_state) + dm = np.outer(init_state, init_state.conj()) + sim = Simulator(virtual_qc, 3, dtype=dtype) + sim.set_qs(dm) + state_vector = sim.get_pure_state_vector() + assert np.allclose(np.outer(state_vector, state_vector.conj()), sim.get_qs(), atol=1e-6) + + +@pytest.mark.level0 +@pytest.mark.platform_x86_gpu_training +@pytest.mark.platform_x86_cpu +@pytest.mark.env_onecard +@pytest.mark.parametrize("config", list(SUPPORTED_SIMULATOR)) +def test_apply_hamiltonian(config): + """ + Description: test apply hamiltonian + Expectation: success. + """ + virtual_qc, dtype = config + circ = random_circuit(3, 100) + ham0 = Hamiltonian(QubitOperator('X0 Y1'), dtype=dtype) + ham1 = ham0.sparse(3) + ham2 = Hamiltonian(csr_matrix(ham0.hamiltonian.matrix(3)), dtype=dtype) + for ham in (ham0, ham1, ham2): + sim = Simulator(virtual_qc, 3, dtype=dtype) + sim.apply_circuit(circ) + sim.apply_hamiltonian(ham) + qs = sim.get_qs() + sim.reset() + sim.apply_circuit(circ) + sim.apply_gate(G.X.on(0)) + sim.apply_gate(G.Y.on(1)) + ref_qs = sim.get_qs() + assert np.allclose(qs, ref_qs, atol=1e-6) + + +@pytest.mark.level0 +@pytest.mark.platform_x86_gpu_training +@pytest.mark.platform_x86_cpu +@pytest.mark.env_onecard +@pytest.mark.parametrize("config", list(SUPPORTED_SIMULATOR)) +def test_sampling(config): + """ + Description: test sampling + Expectation: success. + """ + virtual_qc, dtype = config + shots = 10000 + qs = np.random.rand(4) + np.random.rand(4) * 1j + sim = Simulator(virtual_qc, 2, dtype=dtype) + sim.set_qs(qs) + res = sim.sampling(Circuit(UN(G.Measure(), [0, 1])), shots=shots) + if virtual_qc.startswith("mqmatrix"): + ref_distribution = sim.get_qs().diagonal().real + else: + ref_distribution = np.abs(qs) ** 2 + nonzero = [] + for key in res.data.keys(): + i = int(key, 2) + nonzero.append(ref_distribution[i]) + difference = entropy(np.array(list(res.data.values())) / shots, nonzero) + assert difference < 1e-3 + + +@pytest.mark.level0 +@pytest.mark.platform_x86_gpu_training +@pytest.mark.platform_x86_cpu +@pytest.mark.env_onecard +@pytest.mark.parametrize("config", list(SUPPORTED_SIMULATOR)) +def test_get_expectation(config): + """ + Description: test get expectation + Expectation: success. + """ + virtual_qc, dtype = config + init_state = np.random.rand(8) + np.random.rand(8) * 1j + init_state = init_state / np.linalg.norm(init_state) + circ = random_circuit(3, 100) + ham0 = Hamiltonian(QubitOperator('X0 Y1'), dtype=dtype) + ham1 = ham0.sparse(3) + ham2 = Hamiltonian(csr_matrix(ham0.hamiltonian.matrix(3)), dtype=dtype) + for ham in (ham0, ham1, ham2): + sim = Simulator(virtual_qc, 3, dtype=dtype) + sim.set_qs(init_state) + f = sim.get_expectation(ham, circ) + ref_f = ( + init_state.T.conj() @ circ.hermitian().matrix() @ ham0.hamiltonian.matrix(3) @ circ.matrix() @ init_state + ) + assert np.allclose(f, ref_f, atol=1e-6) + + +@pytest.mark.level0 +@pytest.mark.platform_x86_gpu_training +@pytest.mark.platform_x86_cpu +@pytest.mark.env_onecard +@pytest.mark.parametrize("config", list(SUPPORTED_SIMULATOR)) +def test_get_expectation_with_grad(config): + """ + Description: test get expectation with gradient + Expectation: success. + """ + virtual_qc, dtype = config + init_state = np.random.rand(8) + np.random.rand(8) * 1j + init_state = init_state / np.linalg.norm(init_state) + circ0 = random_circuit(3, 100) + circ1 = random_circuit(3, 100) + circ = circ0 + G.RX('a').on(0) + circ1 + ham0 = Hamiltonian(QubitOperator('X0 Y1'), dtype=dtype) + ham1 = ham0.sparse(3) + ham2 = Hamiltonian(csr_matrix(ham0.hamiltonian.matrix(3)), dtype=dtype) + for ham in (ham0, ham1, ham2): + sim = Simulator(virtual_qc, 3, dtype=dtype) + sim.set_qs(init_state) + grad_ops = sim.get_expectation_with_grad(ham, circ) + pr = np.random.rand() + f, g = grad_ops([pr]) + ref_f = ( + init_state.T.conj() + @ circ.hermitian().matrix({'a': pr}) + @ ham0.hamiltonian.matrix(3) + @ circ.matrix({'a': pr}) + @ init_state + ) + ref_g = ( + init_state.T.conj() + @ circ.hermitian().matrix({'a': pr}) + @ ham0.hamiltonian.matrix(3) + @ circ1.matrix() + @ np.kron(np.eye(4, 4), G.RX('a').diff_matrix({'a': pr})) + @ circ0.matrix() + @ init_state + ).real * 2 + assert np.allclose(f, ref_f, atol=1e-6) + assert np.allclose(g, ref_g, atol=1e-6) + + +@pytest.mark.level0 +@pytest.mark.platform_x86_gpu_training +@pytest.mark.platform_x86_cpu +@pytest.mark.env_onecard +@pytest.mark.parametrize("virtual_qc", ['mqmatrix']) +@pytest.mark.parametrize("dtype", [mq.complex64, mq.complex128]) +def test_noise_get_expectation_with_grad(virtual_qc, dtype): + """ + Description: test noise circuit get expectation with gradient + Expectation: success. + """ + init_state = np.random.rand(8) + np.random.rand(8) * 1j + init_state = init_state / np.linalg.norm(init_state) + init_dm = np.outer(init_state, init_state.conj()) + circ0 = random_circuit(3, 100, 1.0, 0.0) + circ1 = random_circuit(3, 100, 1.0, 0.0) + circ = circ0 + G.RX('a').on(0) + circ1 + circ = circ.with_noise() + ham0 = Hamiltonian(QubitOperator('X0 Y1'), dtype=dtype) + ham1 = ham0.sparse(3) + ham2 = Hamiltonian(csr_matrix(ham0.hamiltonian.matrix(3)), dtype=dtype) + for ham in (ham0, ham1, ham2): + sim = Simulator(virtual_qc, 3, dtype=dtype) + sim.set_qs(init_dm) + grad_ops = sim.get_expectation_with_grad(ham, circ) + pr = np.random.rand() + f, grad = grad_ops([pr]) + sim.apply_circuit(circ, [pr]) + dm = sim.get_qs() + ref_f = np.trace(ham0.hamiltonian.matrix(3) @ dm) + dm = init_dm + for g in circ: + if g.parameterized: + dm = ( + np.kron(np.eye(4, 4), g.diff_matrix({'a': pr})) + @ dm + @ np.kron(np.eye(4, 4), g.hermitian().matrix({'a': pr})) + ) + elif isinstance(g, G.NoiseGate): + tmp = np.zeros((8, 8), dtype=mq.to_np_type(dtype)) + for m in g.matrix(): + if g.obj_qubits[0] == 0: + big_m = np.kron(np.eye(4, 4), m) + elif g.obj_qubits[0] == 1: + big_m = np.kron(np.kron(np.eye(2, 2), m), np.eye(2, 2)) + else: + big_m = np.kron(m, np.eye(4, 4)) + tmp += big_m @ dm @ big_m.conj().T + dm = tmp + else: + if g.obj_qubits[0] == 0: + big_m = np.kron(np.eye(4, 4), g.matrix()) + elif g.obj_qubits[0] == 1: + big_m = np.kron(np.kron(np.eye(2, 2), g.matrix()), np.eye(2, 2)) + else: + big_m = np.kron(g.matrix(), np.eye(4, 4)) + dm = big_m @ dm @ big_m.conj().T + ref_grad = np.trace(ham0.hamiltonian.matrix(3) @ dm).real * 2 + assert np.allclose(f, ref_f, atol=1e-6) + assert np.allclose(grad, ref_grad, atol=1e-4) + + +@pytest.mark.level0 +@pytest.mark.platform_x86_gpu_training +@pytest.mark.platform_x86_cpu +@pytest.mark.env_onecard +@pytest.mark.parametrize("config", list(SUPPORTED_SIMULATOR)) +def test_entropy(config): + """ + Description: test entropy + Expectation: success. + """ + virtual_qc, dtype = config + circ = random_circuit(3, 100) + circ = circ.with_noise(G.DepolarizingChannel(0.1)) + sim = Simulator(virtual_qc, 3, dtype=dtype) + sim.apply_circuit(circ) + e = sim.entropy() + if virtual_qc.startswith('mqvector'): + ref_entropy = 0 + else: + dm = sim.get_qs() + ref_entropy = -np.trace(dm @ logm(dm)) + assert np.allclose(e, ref_entropy, atol=1e-6) + + +@pytest.mark.level0 +@pytest.mark.platform_x86_gpu_training +@pytest.mark.platform_x86_cpu +@pytest.mark.env_onecard +@pytest.mark.parametrize("config1", list(SUPPORTED_SIMULATOR)) +@pytest.mark.parametrize("config2", list(SUPPORTED_SIMULATOR)) +def test_fidelity(config1, config2): + """ + Description: test fidelity of two quantum states + Expectation: success. + """ + virtual_qc1, dtype1 = config1 + virtual_qc2, dtype2 = config2 + circ1 = random_circuit(3, 100) + circ2 = random_circuit(3, 100) + sim1 = Simulator(virtual_qc1, 3, dtype=dtype1) + sim2 = Simulator(virtual_qc2, 3, dtype=dtype2) + sim1.apply_circuit(circ1) + sim2.apply_circuit(circ2) + qs1 = sim1.get_qs() + qs2 = sim2.get_qs() + f = fidelity(qs1, qs2) + if virtual_qc1.startswith('mqvector'): + qs1 = np.outer(qs1, qs1.conj().T) + if virtual_qc2.startswith('mqvector'): + qs2 = np.outer(qs2, qs2.conj().T) + ref_f = np.trace(sqrtm(sqrtm(qs1) @ qs2 @ sqrtm(qs1))).real ** 2 + assert np.allclose(f, ref_f, atol=1e-6) diff --git a/tests/st/test_simulator/test_simulator.py b/tests/st/test_simulator/test_simulator.py index 94667cab8..0cb56473b 100644 --- a/tests/st/test_simulator/test_simulator.py +++ b/tests/st/test_simulator/test_simulator.py @@ -35,7 +35,7 @@ ) from mindquantum.core.operators import Hamiltonian, QubitOperator from mindquantum.core.parameterresolver import ParameterResolver as PR -from mindquantum.simulator import NoiseBackend, Simulator, inner_product, fidelity +from mindquantum.simulator import NoiseBackend, Simulator, inner_product from mindquantum.simulator.available_simulator import SUPPORTED_SIMULATOR from mindquantum.utils import random_circuit @@ -495,46 +495,6 @@ def test_inner_product(config): assert np.allclose(val_exp, val) -@pytest.mark.level0 -@pytest.mark.platform_x86_gpu_training -@pytest.mark.platform_x86_cpu -@pytest.mark.env_onecard -@pytest.mark.parametrize("config", list(SUPPORTED_SIMULATOR)) -def test_fidelity(config): - """ - Description: test fidelity of two quantum states - Expectation: success. - """ - virtual_qc, dtype = config - sim1 = Simulator(virtual_qc, 1, dtype=dtype) - sim1.apply_gate(G.RX(1.2).on(0)) - sim2 = Simulator(virtual_qc, 1, dtype=dtype) - sim2.apply_gate(G.RY(2.1).on(0)) - val = fidelity(sim1.get_qs(), sim2.get_qs()) - assert np.allclose(val, 0.40853254959045) - - -@pytest.mark.level0 -@pytest.mark.platform_x86_gpu_training -@pytest.mark.platform_x86_cpu -@pytest.mark.env_onecard -@pytest.mark.parametrize("config", list(SUPPORTED_SIMULATOR)) -def test_get_partial_trace(config): - """ - Description: test partial trace of density matrix - Expectation: success. - """ - virtual_qc, dtype = config - if not virtual_qc.startswith('mqmatrix'): - return - circ = Circuit().h(0).x(1, 0) - circ += G.RX(1.2).on(2, 1) - sim = Simulator(virtual_qc, 3, dtype=dtype) - sim.apply_circuit(circ) - mat = sim.get_partial_trace([0, 1]) - assert np.allclose(mat, np.array([[0.84058944, 0.23300977j], [-0.23300977j, 0.15941056]])) - - @pytest.mark.level0 @pytest.mark.platform_x86_gpu_training @pytest.mark.platform_x86_cpu