File size: 2,529 Bytes
224e773 | 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 | //! Density Matrix for Jordan algebra evolution
use hyperkitty_core::{GOLDEN_RATIO_INV, GOLDEN_RATIO_INV_SQ};
#[derive(Debug, Clone)]
pub struct DensityMatrix {
pub data: Vec<Vec<f64>>,
}
impl DensityMatrix {
pub fn new(n: usize) -> Self {
Self { data: vec![vec![0.0; n]; n] }
}
pub fn identity(n: usize) -> Self {
let mut d = Self::new(n);
for i in 0..n {
d.data[i][i] = 1.0;
}
d
}
pub fn size(&self) -> usize {
self.data.len()
}
pub fn multiply(&self, other: &Self) -> Self {
let n = self.size();
let mut result = Self::new(n);
for i in 0..n {
for j in 0..n {
let mut sum = 0.0;
for k in 0..n {
sum += self.data[i][k] * other.data[k][j];
}
result.data[i][j] = sum;
}
}
result
}
pub fn scale(&self, s: f64) -> Self {
let mut result = self.clone();
for row in &mut result.data {
for val in row {
*val *= s;
}
}
result
}
pub fn add(&self, other: &Self) -> Self {
let n = self.size();
let mut result = Self::new(n);
for i in 0..n {
for j in 0..n {
result.data[i][j] = self.data[i][j] + other.data[i][j];
}
}
result
}
pub fn transpose(&self) -> Self {
let n = self.size();
let mut result = Self::new(n);
for i in 0..n {
for j in 0..n {
result.data[i][j] = self.data[j][i];
}
}
result
}
pub fn evolve(&self, u: &Self) -> Self {
let ut = u.transpose();
let urho = u.multiply(self);
let urhou = urho.multiply(&ut);
let term1 = urhou.scale(GOLDEN_RATIO_INV);
let term2 = self.scale(GOLDEN_RATIO_INV_SQ);
term1.add(&term2)
}
pub fn frobenius_norm(&self) -> f64 {
self.data.iter().flat_map(|row| row).map(|x| x * x).sum::<f64>().sqrt()
}
pub fn distance_to(&self, other: &Self) -> f64 {
self.add(&other.scale(-1.0)).frobenius_norm()
}
pub fn trace(&self) -> f64 {
let n = self.size();
(0..n).map(|i| self.data[i][i]).sum()
}
}
#[cfg(test)]
mod tests {
use super::*;
#[test]
fn identity_trace() {
let id = DensityMatrix::identity(3);
assert!((id.trace() - 3.0).abs() < 1e-10);
}
}
|