microsoft/qdk
Publicmirrored fromhttps://github.com/microsoft/qdkAvailable
source/noisy_simulator/src/density_matrix_simulator.rs
373lines · modecode
| 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)] |
| 8 | mod tests; |
| 9 | |
| 10 | use crate::{ |
| 11 | ComplexVector, Error, NoisySimulator, SquareMatrix, TOLERANCE, handle_error, |
| 12 | instrument::Instrument, kernel::apply_kernel, operation::Operation, |
| 13 | }; |
| 14 | use num_complex::Complex; |
| 15 | use rand::{Rng, SeedableRng, rngs::StdRng}; |
| 16 | |
| 17 | /// A vectorized density matrix. |
| 18 | #[derive(Debug, Clone)] |
| 19 | pub 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 | |
| 30 | impl 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. |
| 177 | pub 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 | |
| 187 | impl 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 | /// Apply non selective evolution to the given qubit ids. |
| 198 | pub fn apply_instrument( |
| 199 | &mut self, |
| 200 | instrument: &Instrument, |
| 201 | qubits: &[usize], |
| 202 | ) -> Result<(), Error> { |
| 203 | self.check_out_of_bounds_qubits(qubits)?; |
| 204 | self.state |
| 205 | .as_mut()? |
| 206 | .apply_operation_matrix(instrument.non_selective_operation_matrix(), qubits)?; |
| 207 | if let Err(err) = self.state.as_mut()?.renormalize() { |
| 208 | handle_error!(self, err); |
| 209 | } |
| 210 | Ok(()) |
| 211 | } |
| 212 | } |
| 213 | |
| 214 | impl NoisySimulator for DensityMatrixSimulator { |
| 215 | type State = DensityMatrix; |
| 216 | |
| 217 | /// Creates a new `DensityMatrixSimulator`. |
| 218 | fn new(number_of_qubits: usize) -> Self { |
| 219 | let density_matrix = DensityMatrix::new(number_of_qubits); |
| 220 | let dimension = density_matrix.dimension(); |
| 221 | Self { |
| 222 | state: Ok(density_matrix), |
| 223 | dimension, |
| 224 | rng: StdRng::from_entropy(), |
| 225 | } |
| 226 | } |
| 227 | |
| 228 | /// Creates a new `DensityMatrixSimulator` with a given seed for its random number generator. |
| 229 | fn new_with_seed(number_of_qubits: usize, seed: u64) -> Self { |
| 230 | let density_matrix = DensityMatrix::new(number_of_qubits); |
| 231 | let dimension = density_matrix.dimension(); |
| 232 | Self { |
| 233 | state: Ok(density_matrix), |
| 234 | dimension, |
| 235 | rng: StdRng::seed_from_u64(seed), |
| 236 | } |
| 237 | } |
| 238 | |
| 239 | /// Apply an operation to the given qubit ids. |
| 240 | fn apply_operation(&mut self, operation: &Operation, qubits: &[usize]) -> Result<(), Error> { |
| 241 | self.check_out_of_bounds_qubits(qubits)?; |
| 242 | |
| 243 | self.state |
| 244 | .as_mut()? |
| 245 | .apply_operation_matrix(operation.matrix(), qubits)?; |
| 246 | if let Err(err) = self.state.as_mut()?.renormalize() { |
| 247 | handle_error!(self, err); |
| 248 | } |
| 249 | Ok(()) |
| 250 | } |
| 251 | |
| 252 | /// Performs selective evolution under the given instrument. |
| 253 | /// Returns the index of the observed outcome. |
| 254 | /// |
| 255 | /// Use this method to perform measurements on the quantum system. |
| 256 | fn sample_instrument( |
| 257 | &mut self, |
| 258 | instrument: &Instrument, |
| 259 | qubits: &[usize], |
| 260 | ) -> Result<usize, Error> { |
| 261 | let sample = self.rng.r#gen(); |
| 262 | self.sample_instrument_with_distribution(instrument, qubits, sample) |
| 263 | } |
| 264 | |
| 265 | /// Performs selective evolution under the given instrument using a custom random distribution. |
| 266 | /// Returns the index of the observed outcome. |
| 267 | /// |
| 268 | /// This method is used for testing and debugging purposes. |
| 269 | fn sample_instrument_with_distribution( |
| 270 | &mut self, |
| 271 | instrument: &Instrument, |
| 272 | qubits: &[usize], |
| 273 | random_sample: f64, |
| 274 | ) -> Result<usize, Error> { |
| 275 | self.check_out_of_bounds_qubits(qubits)?; |
| 276 | |
| 277 | let mut tmp_state = self.state.clone()?; |
| 278 | apply_kernel( |
| 279 | &mut tmp_state.data, |
| 280 | instrument.total_effect_transposed(), |
| 281 | qubits, |
| 282 | )?; |
| 283 | let total_effect_trace = tmp_state.trace()?; |
| 284 | if total_effect_trace < TOLERANCE { |
| 285 | let err = Error::ProbabilityZeroEvent; |
| 286 | handle_error!(self, err); |
| 287 | } |
| 288 | let mut last_non_zero_trace_outcome: usize = 0; |
| 289 | let mut last_non_zero_trace: f64 = 0.0; |
| 290 | let mut summed_probability: f64 = 0.0; |
| 291 | |
| 292 | for outcome in 0..instrument.num_operations() { |
| 293 | if summed_probability > random_sample { |
| 294 | break; |
| 295 | } |
| 296 | tmp_state = self.state.clone()?; |
| 297 | apply_kernel( |
| 298 | &mut tmp_state.data, |
| 299 | instrument.operation(outcome).effect_matrix_transpose(), |
| 300 | qubits, |
| 301 | )?; |
| 302 | let outcome_trace = tmp_state.trace()?; |
| 303 | summed_probability += outcome_trace / total_effect_trace; |
| 304 | if outcome_trace >= TOLERANCE { |
| 305 | last_non_zero_trace_outcome = outcome; |
| 306 | last_non_zero_trace = outcome_trace; |
| 307 | } |
| 308 | } |
| 309 | |
| 310 | if summed_probability + TOLERANCE <= random_sample || last_non_zero_trace < TOLERANCE { |
| 311 | let err = Error::FailedToSampleInstrumentOutcome; |
| 312 | handle_error!(self, err); |
| 313 | } |
| 314 | |
| 315 | if let Err(err) = self.state.as_mut()?.apply_operation_matrix( |
| 316 | instrument.operation(last_non_zero_trace_outcome).matrix(), |
| 317 | qubits, |
| 318 | ) { |
| 319 | handle_error!(self, err); |
| 320 | } |
| 321 | if let Err(err) = self |
| 322 | .state |
| 323 | .as_mut()? |
| 324 | .renormalize_with_trace(last_non_zero_trace) |
| 325 | { |
| 326 | handle_error!(self, err); |
| 327 | } |
| 328 | Ok(last_non_zero_trace_outcome) |
| 329 | } |
| 330 | |
| 331 | /// Returns the `DensityMatrix` if the simulator is in a valid state. |
| 332 | fn state(&self) -> Result<&DensityMatrix, &Error> { |
| 333 | self.state.as_ref() |
| 334 | } |
| 335 | |
| 336 | /// Set state of the quantum system. |
| 337 | fn set_state(&mut self, new_state: DensityMatrix) -> Result<(), Error> { |
| 338 | if self.dimension != new_state.dimension() { |
| 339 | return Err(Error::InvalidState(format!( |
| 340 | "the provided state should have the same dimensions as the quantum system's state, {} != {}", |
| 341 | self.dimension, |
| 342 | new_state.dimension(), |
| 343 | ))); |
| 344 | } |
| 345 | if !new_state.is_normalized()? { |
| 346 | return Err(Error::InvalidState(format!( |
| 347 | "`state` is not normalized, trace is {}", |
| 348 | new_state.trace()? |
| 349 | ))); |
| 350 | } |
| 351 | if !new_state.is_hermitian() { |
| 352 | return Err(Error::InvalidState("`state` is not Hermitian".to_string())); |
| 353 | } |
| 354 | self.state = Ok(new_state); |
| 355 | Ok(()) |
| 356 | } |
| 357 | |
| 358 | /// Return theoretical change in trace due to operations that have been applied so far |
| 359 | /// In reality, the density matrix is always renormalized after instruments/operations |
| 360 | /// have been applied. |
| 361 | fn trace_change(&self) -> Result<f64, Error> { |
| 362 | Ok(self.state.as_ref()?.trace_change()) |
| 363 | } |
| 364 | |
| 365 | /// Set the trace of the quantum system. |
| 366 | fn set_trace(&mut self, trace: f64) -> Result<(), Error> { |
| 367 | if trace < TOLERANCE || (trace - 1.) > TOLERANCE { |
| 368 | return Err(Error::NotNormalized(trace)); |
| 369 | } |
| 370 | self.state.as_mut()?.trace_change = trace; |
| 371 | Ok(()) |
| 372 | } |
| 373 | } |
| 374 | |