GitVita
New repositoryImport repositoryNew issue
Your repositories Your dashboard Sign in

microsoft/qdk

Public

mirrored fromhttps://github.com/microsoft/qdkAvailable

Watch0Fork0Star0
CodeCommitsIssuesPull requestsActionsInsightsSecurity
billti/num2-sim

Branches

Tags

  • No tags available.
0Branches0Tags
Go to file
Add file
Code

Clone

HTTPS

Download ZIP

qdk/source/noisy_simulator/src

source/noisy_simulator/src/density_matrix_simulator.rs

370lines · modecode

RawDownload
Latest commit unavailable.
1// Copyright (c) Microsoft Corporation.
2// Licensed under the MIT License.
3
4//! This module contains two structs: the `DensityMatrixSimulator` and its
5//! internal `DensityMatrix` state.
6
7#[cfg(test)]
8mod tests;
9
10use crate::{
11 ComplexVector, Error, NoisySimulator, SquareMatrix, TOLERANCE, handle_error,
12 instrument::Instrument, kernel::apply_kernel, operation::Operation,
13};
14use num_complex::Complex;
15use rand::{Rng, SeedableRng, rngs::StdRng};
16
17/// A vectorized density matrix.
18#[derive(Debug, Clone)]
19pub struct DensityMatrix {
20 /// Dimension of the matrix. E.g.: If the matrix is 5 x 5, then dimension is 5.
21 dimension: usize,
22 /// Number of qubits in the system.
23 number_of_qubits: usize,
24 /// Theoretical change in trace due to operations that have been applied so far.
25 trace_change: f64,
26 /// Vector storing the entries of the density matrix.
27 data: ComplexVector,
28}
29
30impl DensityMatrix {
31 fn new(number_of_qubits: usize) -> Self {
32 let dimension = 1 << number_of_qubits;
33 let mut data = ComplexVector::zeros(dimension * dimension);
34 data[0].re = 1.0;
35 Self {
36 dimension,
37 number_of_qubits,
38 trace_change: 1.0,
39 data,
40 }
41 }
42
43 /// Builds a `DensityMatrix` from its raw fields. Returns `None` if
44 /// the provided args don't represent a valid `DensityMatrix`.
45 ///
46 /// This method is to be used from the `PyO3` wrapper.
47 pub fn try_from(
48 dimension: usize,
49 number_of_qubits: usize,
50 trace_change: f64,
51 data: ComplexVector,
52 ) -> Result<Self, Error> {
53 if 1 << number_of_qubits != dimension {
54 return Err(Error::DensityMatrixTryFromError(format!(
55 "the system has {number_of_qubits} qubits and density matrix has dimension {dimension} but 2 ^ {number_of_qubits} != {dimension}"
56 )));
57 }
58
59 if data.len() != dimension * dimension {
60 return Err(Error::DensityMatrixTryFromError(format!(
61 "matrix dimension is {} but data has {} entries, {} x {} != {}",
62 dimension,
63 data.len(),
64 dimension,
65 dimension,
66 data.len()
67 )));
68 }
69
70 Ok(Self {
71 dimension,
72 number_of_qubits,
73 trace_change,
74 data,
75 })
76 }
77
78 /// Returns a reference to the vector containing the density matrix's data.
79 #[must_use]
80 pub fn data(&self) -> &ComplexVector {
81 &self.data
82 }
83
84 /// Returns dimension of the matrix. E.g.: If the matrix is 5 x 5, then dimension is 5.
85 #[must_use]
86 pub fn dimension(&self) -> usize {
87 self.dimension
88 }
89
90 /// Returns the number of qubits in the system.
91 #[must_use]
92 pub fn number_of_qubits(&self) -> usize {
93 self.number_of_qubits
94 }
95
96 /// Returns `true` if the matrix is Hermitian.
97 fn is_hermitian(&self) -> bool {
98 for row in 0..self.dimension {
99 for col in 0..self.dimension {
100 let elt = self.data[self.dimension * row + col];
101 let mirror_elt = self.data[self.dimension * col + row];
102 if (elt.re - mirror_elt.re).abs() > TOLERANCE
103 || (elt.im + mirror_elt.im).abs() > TOLERANCE
104 {
105 return false;
106 }
107 }
108 }
109 true
110 }
111
112 /// Returns `true` if the trace of the matrix is 1.
113 fn is_normalized(&self) -> Result<bool, Error> {
114 Ok((self.trace()? - 1.0).abs() <= TOLERANCE)
115 }
116
117 /// Returns the trace of the matrix. The trace is the sum of the diagonal entries of a matrix.
118 fn trace(&self) -> Result<f64, Error> {
119 let mut trace: Complex<f64> = Complex::ZERO;
120 for idx in 0..self.dimension {
121 trace += self.data[(self.dimension + 1) * idx];
122 }
123
124 if trace.im > TOLERANCE {
125 Err(Error::TraceIsNotReal(trace.im))
126 } else {
127 Ok(trace.re)
128 }
129 }
130
131 /// Return theoretical change in trace due to operations that have been applied so far.
132 /// In reality, the density matrix is always renormalized after instruments / operations
133 /// have been applied.
134 #[must_use]
135 pub fn trace_change(&self) -> f64 {
136 self.trace_change
137 }
138
139 /// Renormalizes the matrix such that the trace is 1.
140 fn renormalize(&mut self) -> Result<(), Error> {
141 self.renormalize_with_trace(self.trace()?)
142 }
143
144 /// Renormalizes the matrix such that the trace is 1. Uses a precomputed `trace`.
145 fn renormalize_with_trace(&mut self, trace: f64) -> Result<(), Error> {
146 if trace < TOLERANCE {
147 return Err(Error::ProbabilityZeroEvent);
148 }
149 self.trace_change *= trace;
150 let renormalization_factor = 1.0 / trace;
151 self.data.scale_mut(renormalization_factor);
152 Ok(())
153 }
154
155 /// Applies the operation matrix to the target qubits.
156 fn apply_operation_matrix(
157 &mut self,
158 operation_matrix: &SquareMatrix,
159 qubits: &[usize],
160 ) -> Result<(), Error> {
161 // TODO [Research]: Figure out why they do this qubits_expanded thing.
162 let mut qubits_expanded = Vec::with_capacity(2 * qubits.len());
163 for id in qubits {
164 qubits_expanded.push(*id);
165 }
166 for id in qubits {
167 qubits_expanded.push(*id + self.number_of_qubits());
168 }
169 apply_kernel(&mut self.data, operation_matrix, &qubits_expanded)
170 }
171}
172
173/// A quantum circuit simulator using a density matrix.
174///
175/// All the simulator methods return a `Result<_, Error>`. If the simulator reaches an
176/// invalid state due to a numerical error, it will return that last error from there on.
177pub struct DensityMatrixSimulator {
178 /// A `DensityMatrix` representing the current state of the quantum system.
179 state: Result<DensityMatrix, Error>,
180 /// Dimension of the density matrix. We need this field to verify the size of the
181 /// quantum system in the `set_state` method in the case that `self.state == Err(...)`.
182 dimension: usize,
183 /// Random number generator used for probabilistic operations.
184 rng: StdRng,
185}
186
187impl DensityMatrixSimulator {
188 fn check_out_of_bounds_qubits(&self, qubits: &[usize]) -> Result<(), Error> {
189 let number_of_qubits = self.state.as_ref()?.number_of_qubits;
190 if let Some(id) = qubits.iter().find(|id| **id >= number_of_qubits) {
191 Err(Error::QubitIdOutOfBounds(*id))
192 } else {
193 Ok(())
194 }
195 }
196}
197
198impl NoisySimulator for DensityMatrixSimulator {
199 type State = DensityMatrix;
200
201 /// Creates a new `DensityMatrixSimulator`.
202 fn new(number_of_qubits: usize) -> Self {
203 let density_matrix = DensityMatrix::new(number_of_qubits);
204 let dimension = density_matrix.dimension();
205 Self {
206 state: Ok(density_matrix),
207 dimension,
208 rng: StdRng::from_entropy(),
209 }
210 }
211
212 /// Creates a new `DensityMatrixSimulator` with a given seed for its random number generator.
213 fn new_with_seed(number_of_qubits: usize, seed: u64) -> Self {
214 let density_matrix = DensityMatrix::new(number_of_qubits);
215 let dimension = density_matrix.dimension();
216 Self {
217 state: Ok(density_matrix),
218 dimension,
219 rng: StdRng::seed_from_u64(seed),
220 }
221 }
222
223 /// Apply an operation to the given qubit ids.
224 fn apply_operation(&mut self, operation: &Operation, qubits: &[usize]) -> Result<(), Error> {
225 self.check_out_of_bounds_qubits(qubits)?;
226
227 self.state
228 .as_mut()?
229 .apply_operation_matrix(operation.matrix(), qubits)?;
230 if let Err(err) = self.state.as_mut()?.renormalize() {
231 handle_error!(self, err);
232 }
233 Ok(())
234 }
235
236 /// Apply non selective evolution to the given qubit ids.
237 fn apply_instrument(&mut self, instrument: &Instrument, qubits: &[usize]) -> Result<(), Error> {
238 self.check_out_of_bounds_qubits(qubits)?;
239
240 self.state
241 .as_mut()?
242 .apply_operation_matrix(instrument.non_selective_operation_matrix(), qubits)?;
243 if let Err(err) = self.state.as_mut()?.renormalize() {
244 handle_error!(self, err);
245 }
246 Ok(())
247 }
248
249 /// Performs selective evolution under the given instrument.
250 /// Returns the index of the observed outcome.
251 ///
252 /// Use this method to perform measurements on the quantum system.
253 fn sample_instrument(
254 &mut self,
255 instrument: &Instrument,
256 qubits: &[usize],
257 ) -> Result<usize, Error> {
258 let sample = self.rng.r#gen();
259 self.sample_instrument_with_distribution(instrument, qubits, sample)
260 }
261
262 /// Performs selective evolution under the given instrument using a custom random distribution.
263 /// Returns the index of the observed outcome.
264 ///
265 /// This method is used for testing and debugging purposes.
266 fn sample_instrument_with_distribution(
267 &mut self,
268 instrument: &Instrument,
269 qubits: &[usize],
270 random_sample: f64,
271 ) -> Result<usize, Error> {
272 self.check_out_of_bounds_qubits(qubits)?;
273
274 let mut tmp_state = self.state.clone()?;
275 apply_kernel(
276 &mut tmp_state.data,
277 instrument.total_effect_transposed(),
278 qubits,
279 )?;
280 let total_effect_trace = tmp_state.trace()?;
281 if total_effect_trace < TOLERANCE {
282 let err = Error::ProbabilityZeroEvent;
283 handle_error!(self, err);
284 }
285 let mut last_non_zero_trace_outcome: usize = 0;
286 let mut last_non_zero_trace: f64 = 0.0;
287 let mut summed_probability: f64 = 0.0;
288
289 for outcome in 0..instrument.num_operations() {
290 if summed_probability > random_sample {
291 break;
292 }
293 tmp_state = self.state.clone()?;
294 apply_kernel(
295 &mut tmp_state.data,
296 instrument.operation(outcome).effect_matrix_transpose(),
297 qubits,
298 )?;
299 let outcome_trace = tmp_state.trace()?;
300 summed_probability += outcome_trace / total_effect_trace;
301 if outcome_trace >= TOLERANCE {
302 last_non_zero_trace_outcome = outcome;
303 last_non_zero_trace = outcome_trace;
304 }
305 }
306
307 if summed_probability + TOLERANCE <= random_sample || last_non_zero_trace < TOLERANCE {
308 let err = Error::FailedToSampleInstrumentOutcome;
309 handle_error!(self, err);
310 }
311
312 if let Err(err) = self.state.as_mut()?.apply_operation_matrix(
313 instrument.operation(last_non_zero_trace_outcome).matrix(),
314 qubits,
315 ) {
316 handle_error!(self, err);
317 }
318 if let Err(err) = self
319 .state
320 .as_mut()?
321 .renormalize_with_trace(last_non_zero_trace)
322 {
323 handle_error!(self, err);
324 }
325 Ok(last_non_zero_trace_outcome)
326 }
327
328 /// Returns the `DensityMatrix` if the simulator is in a valid state.
329 fn state(&self) -> Result<&DensityMatrix, &Error> {
330 self.state.as_ref()
331 }
332
333 /// Set state of the quantum system.
334 fn set_state(&mut self, new_state: DensityMatrix) -> Result<(), Error> {
335 if self.dimension != new_state.dimension() {
336 return Err(Error::InvalidState(format!(
337 "the provided state should have the same dimensions as the quantum system's state, {} != {}",
338 self.dimension,
339 new_state.dimension(),
340 )));
341 }
342 if !new_state.is_normalized()? {
343 return Err(Error::InvalidState(format!(
344 "`state` is not normalized, trace is {}",
345 new_state.trace()?
346 )));
347 }
348 if !new_state.is_hermitian() {
349 return Err(Error::InvalidState("`state` is not Hermitian".to_string()));
350 }
351 self.state = Ok(new_state);
352 Ok(())
353 }
354
355 /// Return theoretical change in trace due to operations that have been applied so far
356 /// In reality, the density matrix is always renormalized after instruments/operations
357 /// have been applied.
358 fn trace_change(&self) -> Result<f64, Error> {
359 Ok(self.state.as_ref()?.trace_change())
360 }
361
362 /// Set the trace of the quantum system.
363 fn set_trace(&mut self, trace: f64) -> Result<(), Error> {
364 if trace < TOLERANCE || (trace - 1.) > TOLERANCE {
365 return Err(Error::NotNormalized(trace));
366 }
367 self.state.as_mut()?.trace_change = trace;
368 Ok(())
369 }
370}
371