quant_iron/components/pauli_string.rs
1use crate::{
2 components::{operator::Operator, operator::Pauli, state::State, gate::Gate},
3 errors::Error,
4};
5use num_complex::Complex;
6use rayon::prelude::*;
7use std::collections::HashMap;
8use std::ops::{Add, Mul};
9
10/// Represents a Pauli string, which is a product of Pauli operators (X, Y, Z) acting on qubits.
11/// Used to represent a term in a Hamiltonian or a quantum operator.
12#[derive(Debug, Clone, PartialEq)]
13pub struct PauliString {
14 /// A mapping from qubit indices to Pauli operators.
15 ops: HashMap<usize, Pauli>,
16 /// The coefficient of the Pauli string, which is a complex number.
17 coefficient: Complex<f64>,
18}
19
20impl PauliString {
21 /// Creates a new Pauli string with the given coefficient and an empty set of operators.
22 ///
23 /// # Arguments
24 /// * `coefficient` - The coefficient of the Pauli string, represented as a complex number.
25 ///
26 /// # Returns
27 /// A new `PauliString` instance with the specified coefficient and no operators.
28 pub fn new(coefficient: Complex<f64>) -> Self {
29 Self {
30 ops: HashMap::new(),
31 coefficient,
32 }
33 }
34
35 /// Returns the length of the Pauli string, defined as the number of operators it contains.
36 ///
37 /// # Returns
38 ///
39 /// * `usize` - The number of operators in the Pauli string.
40 pub fn len(&self) -> usize {
41 self.ops.len()
42 }
43
44 /// Creates a new Pauli string with the given coefficient and a set of operators.
45 ///
46 /// # Arguments
47 /// * `coefficient` - The coefficient of the Pauli string, represented as a complex number.
48 /// * `ops` - A mapping from qubit indices to Pauli operators.
49 ///
50 /// # Returns
51 /// A new `PauliString` instance with the specified coefficient and operators.
52 ///
53 /// Note that the Hashmap ensures uniqueness of operators for each qubit.
54 pub fn with_ops(coefficient: Complex<f64>, ops: HashMap<usize, Pauli>) -> Self {
55 Self { ops, coefficient }
56 }
57
58 /// Adds a Pauli operator to the Pauli string at the specified qubit index.
59 ///
60 /// # Arguments
61 /// * `qubit` - The index of the qubit to which the operator is applied.
62 /// * `op` - The Pauli operator to be added (X, Y, or Z).
63 ///
64 /// # Panics
65 /// This function will panic if an operator for the same qubit index is added more than once.
66 pub fn add_op(&mut self, qubit: usize, op: Pauli) {
67 if self.ops.insert(qubit, op).is_some() {
68 panic!("Duplicate Pauli string operator for qubit: {}", qubit);
69 }
70 }
71
72 /// Adds a Pauli operator to the Pauli string at the specified qubit index and returns the new `PauliString` instance.
73 ///
74 /// # Arguments
75 ///
76 /// * `qubit` - The index of the qubit to which the operator is applied.
77 /// * `op` - The Pauli operator to be added (X, Y, or Z).
78 ///
79 /// # Returns
80 ///
81 /// * `Self` - A new `PauliString` instance with the added operator.
82 ///
83 /// # Panics
84 /// This function will panic if an operator for the same qubit index is added more than once.
85 pub fn with_op(mut self, qubit: usize, op: Pauli) -> Self {
86 self.add_op(qubit, op);
87 self
88 }
89
90 /// Returns the coefficient of the Pauli string.
91 ///
92 /// # Returns
93 ///
94 /// * `Complex<f64>` - The coefficient of the Pauli string, represented as a complex number.
95 pub fn coefficient(&self) -> Complex<f64> {
96 self.coefficient
97 }
98
99 /// Returns a reference to the operators in the Pauli string.
100 ///
101 /// # Returns
102 ///
103 /// * `&HashMap<usize, Pauli>` - A reference to the mapping of qubit indices to Pauli operators.
104 pub fn ops(&self) -> &HashMap<usize, Pauli> {
105 &self.ops
106 }
107
108 /// Returns the list of targets of the Pauli string
109 ///
110 /// # Returns
111 /// * `Vec<usize>` - A vector of qubit indices that the Pauli string acts upon.
112 pub fn get_targets(&self) -> Vec<usize> {
113 let mut keys = self.ops.keys().cloned().collect::<Vec<usize>>();
114 keys.sort();
115 keys
116 }
117
118 /// Converts the Pauli string to a vector of operator gates.
119 ///
120 /// # Returns
121 /// * `Vec<Gate>` - A vector of Gate structs representing the individual Pauli operators.
122 pub fn to_gates(&self) -> Vec<Gate> {
123 self.ops.iter().map(|(qubit, op)| {
124 Gate::Operator(Box::new(*op), vec![*qubit], vec![])
125 }).collect()
126 }
127
128 /// Applies the Pauli string to a given state.
129 ///
130 /// # Arguments
131 /// * `state` - The state to which the Pauli string is applied.
132 ///
133 /// # Returns
134 /// * `Result<State, Error>` - The resulting state after applying the Pauli string, or an error if the operation fails.
135 ///
136 /// # Errors
137 ///
138 /// * Returns an error if the operations in the Pauli string refer to qubits outside the range of the state.
139 pub fn apply(&self, state: &State) -> Result<State, Error> {
140 // If the Pauli string is empty, return the state multiplied by the coefficient
141 if self.ops.is_empty() {
142 return Ok(state.clone() * self.coefficient);
143 }
144
145 // Apply the Pauli string to the state
146 let mut new_state: State = state.clone();
147 for (qubit, op) in &self.ops {
148 new_state = op.apply(&new_state, &[*qubit], &[])?; // Assumes op.apply can take &mut State and modify it or returns a new one
149 }
150 Ok(new_state * self.coefficient)
151 }
152
153 /// Applies the Pauli string to a given state and normalises the new state.
154 ///
155 /// # Arguments
156 /// * `state` - The state to which the Pauli string is applied.
157 ///
158 /// # Returns
159 /// * `Result<State, Error>` - The resulting state after applying the Pauli string, or an error if the operation fails.
160 ///
161 /// # Errors
162 ///
163 /// * Returns an error if the operations in the Pauli string refer to qubits outside the range of the state.
164 /// * Returns an error if the resulting state cannot be normalised (eg., has zero norm).
165 pub fn apply_normalised(&self, state: &State) -> Result<State, Error> {
166 let new_state: State = self.apply_operators(state)?;
167 new_state.normalise()
168 }
169
170 /// Helper function to apply only the operator part of the Pauli string (P_ops) to a state.
171 /// This does not include the PauliString's own coefficient.
172 fn apply_operators(&self, state: &State) -> Result<State, Error> {
173 if self.ops.is_empty() {
174 // If there are no operators, P_ops is effectively Identity.
175 // P_ops |state> = |state>.
176 return Ok(state.clone());
177 }
178
179 let mut current_state = state.clone();
180 for (qubit_idx, pauli_op) in &self.ops {
181 current_state = pauli_op.apply(¤t_state, &[*qubit_idx], &[])?;
182 }
183 Ok(current_state)
184 }
185
186 /// Applies the exponential of the Pauli string to a given state.
187 ///
188 /// # Arguments
189 ///
190 /// * `state` - The state to which the exponential of the Pauli string is applied.
191 ///
192 /// # Returns
193 ///
194 /// * `Result<State, Error>` - The resulting state after applying the exponential of the Pauli string, or an error if the operation fails.
195 ///
196 /// # Errors
197 ///
198 /// * Returns an error if the operations in the Pauli string refer to qubits outside the range of the state.
199 pub fn apply_exp(&self, state: &State) -> Result<State, Error> {
200 // If the Pauli string is empty, return the state multiplied by the exponential of the coefficient
201 let alpha: Complex<f64> = self.coefficient;
202
203 if self.ops.is_empty() {
204 // P_ops is Identity. exp(alpha * I) |state> = exp(alpha) * |state>
205 return Ok(state.clone() * alpha.exp());
206 }
207
208 // 1. Calculate P_ops |state> using the helper
209 let p_ops_psi_state = self.apply_operators(state)?;
210
211 // 2. Calculate scalar coefficients for exp(alpha * P_ops) = cosh(alpha)*I + sinh(alpha)*P_ops
212 let cosh_alpha: Complex<f64> = alpha.cosh();
213 let sinh_alpha: Complex<f64> = alpha.sinh();
214
215 // 3. Calculate cosh(alpha) * |state>
216 let term_identity_part: State = state.clone() * cosh_alpha;
217
218 // 4. Calculate sinh(alpha) * (P_ops |state>)
219 let term_operator_part: State = p_ops_psi_state * sinh_alpha;
220
221 // 5. Add the two terms: cosh(alpha) * |state> + sinh(alpha) * (P_ops |state>)
222 Ok(term_identity_part + term_operator_part)
223 }
224
225 /// Applies the exponential of the Pauli string to a given state with a specified factor.
226 ///
227 /// # Arguments
228 /// * `state` - The state to which the exponential of the Pauli string is applied.
229 /// * `factor` - A complex factor to be multiplied with the coefficient of the Pauli string.
230 ///
231 /// # Returns
232 /// * `Result<State, Error>` - The resulting state after applying the exponential of the Pauli string with the factor, or an error if the operation fails.
233 ///
234 /// # Errors
235 ///
236 /// * Returns an error if the operations in the Pauli string refer to qubits outside the range of the state.
237 pub fn apply_exp_factor(&self, state: &State, factor: Complex<f64>) -> Result<State, Error> {
238 // Calculate the effective coefficient for the exponentiation
239 let alpha: Complex<f64> = self.coefficient * factor;
240
241 if self.ops.is_empty() {
242 // P_ops is Identity. exp(alpha * I) |state> = exp(alpha) * |state>
243 return Ok(state.clone() * alpha.exp());
244 }
245
246 // 1. Calculate P_ops |state> using the helper
247 let p_ops_psi_state = self.apply_operators(state)?;
248 // p_ops_psi_state now holds P_ops |state> (where P_ops is the product of Pauli operators without self.coefficient)
249
250 // 2. Calculate scalar coefficients for exp(alpha * P_ops) = cosh(alpha)*I + sinh(alpha)*P_ops
251 let cosh_alpha: Complex<f64> = alpha.cosh();
252 let sinh_alpha: Complex<f64> = alpha.sinh();
253
254 // 3. Calculate cosh(alpha) * |state>
255 let term_identity_part: State = state.clone() * cosh_alpha;
256
257 // 4. Calculate sinh(alpha) * (P_ops |state>)
258 let term_operator_part: State = p_ops_psi_state * sinh_alpha;
259
260 // 5. Add the two terms: cosh(alpha) * |state> + sinh(alpha) * (P_ops |state>)
261 Ok(term_identity_part + term_operator_part)
262 }
263
264 /// Applies the exponential of the Pauli string to the given state with a negative imaginary factor.
265 /// Represents a time evolution step of the form exp(-i * coefficient * P_ops * dt).
266 /// Guaranteed to return a normalised state if the input state is normalised.
267 ///
268 /// # Arguments
269 ///
270 /// * `state` - The state to which the exponential of the Pauli string is applied.
271 /// * `dt` - The time step for the evolution.
272 ///
273 /// # Returns
274 ///
275 /// * `Result<State, Error>` - The resulting state after applying the exponential of the Pauli string with the negative imaginary factor, or an error if the operation fails.
276 ///
277 /// # Errors
278 ///
279 /// * Returns an error if the operations in the Pauli string refer to qubits outside the range of the state.
280 /// * Returns an error if the coefficient for the Pauli string has an imaginary component.
281 pub fn apply_exp_neg_i_dt(&self, state: &State, dt: f64) -> Result<State, Error> {
282 if self.coefficient.im != 0.0 {
283 return Err(Error::InvalidPauliStringCoefficient(self.coefficient));
284 }
285 let factor: Complex<f64> = Complex::new(0.0, -dt);
286 self.apply_exp_factor(state, factor)
287 }
288
289 /// Returns the Hermitian conjugate of the Pauli string.
290 ///
291 /// # Returns
292 ///
293 /// * `Self` - A new `PauliString` instance representing the Hermitian conjugate of the original Pauli string.
294 pub fn hermitian_conjugate(&self) -> Self {
295 PauliString {
296 ops: self.ops.clone(),
297 coefficient: self.coefficient.conj(),
298 }
299 }
300}
301
302impl Add for PauliString {
303 type Output = SumOp;
304
305 fn add(self, other: Self) -> Self::Output {
306 // Create a new SumOp with the two Pauli strings
307 let terms: Vec<PauliString> = vec![self, other];
308 // Return a new SumOp instance with the terms
309 SumOp::new(terms)
310 }
311}
312
313impl Mul<Complex<f64>> for PauliString {
314 type Output = Self;
315
316 fn mul(self, rhs: Complex<f64>) -> Self::Output {
317 // Create a new Pauli string with the product of the coefficient and the given complex number
318 PauliString {
319 ops: self.ops.clone(),
320 coefficient: self.coefficient * rhs,
321 }
322 }
323}
324
325impl Mul<f64> for PauliString {
326 type Output = Self;
327
328 fn mul(self, rhs: f64) -> Self::Output {
329 // Create a new Pauli string with the product of the coefficient and the given real number
330 PauliString {
331 ops: self.ops.clone(),
332 coefficient: self.coefficient * Complex::new(rhs, 0.0),
333 }
334 }
335}
336
337impl Mul<PauliString> for f64 {
338 type Output = PauliString;
339
340 fn mul(self, rhs: PauliString) -> Self::Output {
341 rhs * self
342 }
343}
344
345impl std::fmt::Display for PauliString {
346 fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result {
347 let coeff_str: String = if self.coefficient.re == 0.0 && self.coefficient.im == 0.0 {
348 "0".to_string()
349 } else if self.coefficient.im == 0.0 {
350 format!("{}", self.coefficient.re)
351 } else if self.coefficient.re == 0.0 {
352 format!("{}i", self.coefficient.im)
353 } else {
354 format!("({} + {}i)", self.coefficient.re, self.coefficient.im)
355 };
356
357 let mut result: String = coeff_str + " * ";
358
359 let mut sorted_ops: Vec<(&usize, &Pauli)> = self.ops.iter().collect();
360 sorted_ops.sort_by(|&(qubit_a, op_a), &(qubit_b, op_b)| {
361 if qubit_a == qubit_b {
362 // If qubits are the same, sort by operator type (X -> Y -> Z)
363 let op_a_str: String = format!("{}", op_a);
364 let op_b_str: String = format!("{}", op_b);
365 match (op_a_str.as_str(), op_b_str.as_str()) {
366 ("X", "Y") | ("X", "Z") | ("Y", "Z") => std::cmp::Ordering::Less,
367 ("Y", "X") | ("Z", "X") | ("Z", "Y") => std::cmp::Ordering::Greater,
368 _ => std::cmp::Ordering::Equal,
369 }
370 } else {
371 // Otherwise, sort by qubit index
372 qubit_a.cmp(qubit_b)
373 }
374 });
375 for (qubit, op) in sorted_ops {
376 result.push_str(&format!("{}[{}] ", op, qubit));
377 }
378
379 write!(f, "{}", result.trim())
380 }
381}
382
383/// A vector of Pauli strings that are summed together.
384/// Useful for representing Hamiltonians or observables.
385///
386/// # Fields
387///
388/// * `terms` - A vector of `PauliString` instances that are summed together.
389#[derive(Debug, Clone, PartialEq)]
390pub struct SumOp {
391 /// A vector of Pauli strings that are summed together.
392 pub terms: Vec<PauliString>,
393}
394
395impl SumOp {
396 /// Creates a new `SumPauliString` instance with the given terms and number of qubits.
397 ///
398 /// # Arguments
399 ///
400 /// * `terms` - A vector of `PauliString` instances that are summed together.
401 ///
402 /// # Returns
403 ///
404 /// * A new `SumPauliString` instance with the specified terms and number of qubits.
405 pub fn new(terms: Vec<PauliString>) -> Self {
406 Self { terms }
407 }
408
409 /// Returns the number of terms in the sum.
410 ///
411 /// # Returns
412 ///
413 /// * `usize` - The number of terms in the sum.
414 pub fn num_terms(&self) -> usize {
415 self.terms.len()
416 }
417
418 /// Adds a new term to the sum.
419 ///
420 /// # Arguments
421 ///
422 /// * `term` - The `PauliString` term to be added to the sum.
423 pub fn add_term(&mut self, term: PauliString) {
424 self.terms.push(term);
425 }
426
427 /// Adds a new term to the sum and returns a new `SumOp` instance.
428 ///
429 /// # Arguments
430 ///
431 /// * `term` - The `PauliString` term to be added to the sum.
432 ///
433 /// # Returns
434 /// * `Self` - A new `SumOp` instance with the added term.
435 pub fn with_term(mut self, term: PauliString) -> Self {
436 self.add_term(term);
437 self
438 }
439
440 /// Applies the sum of Pauli strings to a given state.
441 ///
442 /// # Arguments
443 ///
444 /// * `state` - The state to which the sum of Pauli strings is applied.
445 ///
446 /// # Returns
447 ///
448 /// * `Result<State, Error>` - The resulting state after applying the sum of Pauli strings, or an error if the operation fails.
449 ///
450 /// # Errors
451 ///
452 /// * Returns an error if the operations in the Pauli strings refer to qubits outside the range of the state.
453 pub fn apply(&self, state: &State) -> Result<State, Error> {
454 if self.terms.is_empty() {
455 // An empty sumop is equivalent to the zero operator
456 return Ok(state.clone() * 0.0);
457 }
458
459 Ok(self
460 .terms
461 .par_iter()
462 .map(|term| term.apply(state))
463 .collect::<Result<Vec<State>, Error>>()?
464 .into_iter()
465 .sum())
466 }
467
468 /// Calculates the expectation value <psi|H|psi> = Sum_i <psi|P_i|psi>.
469 ///
470 /// The expectation value is generally real for Hermitian operators and normalised states.
471 /// However, this function returns a `Complex<f64>` as intermediate PauliStrings
472 /// might have complex coefficients or the operator/state might not be strictly physical.
473 ///
474 /// # Arguments
475 /// * `state` - The state |psi> for which to calculate the expectation value.
476 /// For a physically meaningful expectation value, this state should be normalised.
477 ///
478 /// # Returns
479 /// * `Result<Complex<f64>, Error>` - The expectation value, or an error if the operation fails.
480 ///
481 /// # Errors
482 /// * Returns an error if any underlying `PauliString::apply` fails (eg., invalid qubit index).
483 /// * Returns an error if `state.inner_product` fails (eg., mismatched number of qubits,
484 /// though `PauliString::apply` should also catch qubit count issues).
485 pub fn expectation_value(&self, state: &State) -> Result<Complex<f64>, Error> {
486 if self.terms.is_empty() {
487 // The expectation value of a zero operator is zero.
488 return Ok(Complex::new(0.0, 0.0));
489 }
490
491 let expectation_values_per_term: Vec<Complex<f64>> = self
492 .terms
493 .par_iter()
494 .map(|term| {
495 // For each term P_i in the sum H = Sum_i P_i:
496 // 1. Calculate |phi_i> = P_i |psi>
497 let phi_i_state = term.apply(state)?;
498
499 // 2. Calculate <psi | phi_i> = <psi | P_i | psi>
500 state.inner_product(&phi_i_state)
501 })
502 .collect::<Result<Vec<Complex<f64>>, Error>>()?; // Collect into Result<Vec<_>, E>, propagating errors
503
504 // Sum the individual expectation values <psi|P_i|psi>
505 // Complex<f64> from num_complex implements std::iter::Sum.
506 Ok(expectation_values_per_term.into_iter().sum())
507 }
508}
509
510impl std::fmt::Display for SumOp {
511 fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result {
512 let mut result: String = String::new();
513 for (i, term) in self.terms.iter().enumerate() {
514 if i > 0 {
515 result.push_str(" + ");
516 }
517 result.push_str(&format!("{}", term));
518 }
519 write!(f, "{}", result)
520 }
521}
522
523impl Mul<Complex<f64>> for SumOp {
524 type Output = Self;
525
526 fn mul(self, rhs: Complex<f64>) -> Self::Output {
527 // Create a new SumOp with the coefficient multiplied by the given complex number
528 let terms: Vec<PauliString> = self
529 .terms
530 .into_iter()
531 .map(|term| term * rhs)
532 .collect();
533 // Return a new SumOp instance with the modified terms
534 SumOp::new(terms)
535 }
536}
537
538impl Mul<f64> for SumOp {
539 type Output = Self;
540
541 fn mul(self, rhs: f64) -> Self::Output {
542 // Create a new SumOp with the coefficient multiplied by the given real number
543 let terms: Vec<PauliString> = self
544 .terms
545 .into_iter()
546 .map(|term| term * rhs)
547 .collect();
548 // Return a new SumOp instance with the modified terms
549 SumOp::new(terms)
550 }
551}
552
553impl Add for SumOp {
554 type Output = Self;
555
556 fn add(self, other: Self) -> Self::Output {
557 // Create a new SumOp with the two sets of terms
558 let mut terms: Vec<PauliString> = self.terms;
559 terms.extend(other.terms);
560 // Return a new SumOp instance with the combined terms
561 SumOp::new(terms)
562 }
563}
564
565impl Add<PauliString> for SumOp {
566 type Output = Self;
567
568 fn add(self, other: PauliString) -> Self::Output {
569 // Create a new SumOp with the existing terms and the new Pauli string
570 let mut terms: Vec<PauliString> = self.terms;
571 terms.push(other);
572 // Return a new SumOp instance with the combined terms
573 SumOp::new(terms)
574 }
575}