1use std::cmp::Ordering;
4use std::fmt;
5use std::ops::Mul;
6
7use anyhow::{self, format_err, Context};
8use duplicate::duplicate_item;
9use itertools::Itertools;
10use log;
11use ndarray::{s, Array, Array1, Array2, Axis, Dimension, Ix0, Ix2};
12use ndarray_einsum::*;
13use ndarray_linalg::{solve::Inverse, types::Lapack};
14use num_complex::{Complex, ComplexFloat};
15use num_traits::{ToPrimitive, Zero};
16#[cfg(feature = "python")]
17use pyo3::prelude::*;
18use rayon::prelude::*;
19use serde::{Deserialize, Serialize};
20
21use crate::chartab::chartab_group::CharacterProperties;
22use crate::chartab::{CharacterTable, DecompositionError, SubspaceDecomposable};
23use crate::group::{class::ClassProperties, GroupProperties};
24use crate::io::format::{log_subtitle, qsym2_output};
25use crate::symmetry::symmetry_group::UnitaryRepresentedSymmetryGroup;
26
27#[derive(Clone)]
33pub(crate) struct Metric<'a, T, D: Dimension> {
34 hermitian: &'a Array<T, D>,
39
40 complex_symmetric: Option<&'a Array<T, D>>,
46}
47
48pub trait Overlap<T, D>
61where
62 T: ComplexFloat + fmt::Debug + Lapack,
63 D: Dimension,
64{
65 fn complex_symmetric(&self) -> bool;
68
69 fn overlap(
73 &self,
74 other: &Self,
75 metric: Option<&Array<T, D>>,
76 metric_h: Option<&Array<T, D>>,
77 ) -> Result<T, anyhow::Error>;
78
79 fn overlap_definition(&self) -> String;
82}
83
84pub struct OrbitIterator<'a, G, I>
94where
95 G: GroupProperties,
96{
97 group_iter: <<G as GroupProperties>::ElementCollection as IntoIterator>::IntoIter,
100
101 origin: &'a I,
103
104 action: fn(&G::GroupElement, &I) -> Result<I, anyhow::Error>,
106}
107
108impl<'a, G, I> OrbitIterator<'a, G, I>
109where
110 G: GroupProperties,
111{
112 pub fn new(
124 group: &G,
125 origin: &'a I,
126 action: fn(&G::GroupElement, &I) -> Result<I, anyhow::Error>,
127 ) -> Self {
128 Self {
129 group_iter: group.elements().clone().into_iter(),
130 origin,
131 action,
132 }
133 }
134}
135
136impl<'a, G, I> Iterator for OrbitIterator<'a, G, I>
137where
138 G: GroupProperties,
139{
140 type Item = Result<I, anyhow::Error>;
141
142 fn next(&mut self) -> Option<Self::Item> {
143 self.group_iter
144 .next()
145 .map(|op| (self.action)(&op, self.origin))
146 }
147}
148
149pub trait Orbit<G, I>
155where
156 G: GroupProperties,
157{
158 type OrbitIter: Iterator<Item = Result<I, anyhow::Error>>;
160
161 fn group(&self) -> &G;
163
164 fn origin(&self) -> &I;
166
167 fn iter(&self) -> Self::OrbitIter;
169}
170
171#[derive(Clone, Debug, Serialize, Deserialize, PartialEq)]
182#[cfg_attr(feature = "python", pyclass(eq, eq_int))]
183pub enum EigenvalueComparisonMode {
184 Real,
186
187 Modulus,
189}
190
191impl Default for EigenvalueComparisonMode {
192 fn default() -> Self {
193 Self::Modulus
194 }
195}
196
197impl fmt::Display for EigenvalueComparisonMode {
198 fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
199 match self {
200 EigenvalueComparisonMode::Real => write!(f, "Real part"),
201 EigenvalueComparisonMode::Modulus => write!(f, "Modulus"),
202 }
203 }
204}
205
206pub trait RepAnalysis<G, I, T, D>: Orbit<G, I>
217where
218 T: ComplexFloat + Lapack + fmt::Debug,
219 <T as ComplexFloat>::Real: ToPrimitive,
220 G: GroupProperties + ClassProperties + CharacterProperties,
221 G::GroupElement: fmt::Display,
222 G::CharTab: SubspaceDecomposable<T>,
223 D: Dimension,
224 I: Overlap<T, D> + Clone,
225 Self::OrbitIter: Iterator<Item = Result<I, anyhow::Error>>,
226{
227 fn set_smat(&mut self, smat: Array2<T>);
237
238 #[must_use]
240 fn smat(&self) -> Option<&Array2<T>>;
241
242 #[must_use]
261 fn xmat(&self) -> &Array2<T>;
262
263 #[must_use]
281 fn norm_preserving_scalar_map(&self, i: usize) -> Result<fn(T) -> T, anyhow::Error>;
282
283 #[must_use]
286 fn integrality_threshold(&self) -> <T as ComplexFloat>::Real;
287
288 #[must_use]
291 fn eigenvalue_comparison_mode(&self) -> &EigenvalueComparisonMode;
292
293 fn calc_smat(
308 &mut self,
309 metric: Option<&Array<T, D>>,
310 metric_h: Option<&Array<T, D>>,
311 use_cayley_table: bool,
312 ) -> Result<&mut Self, anyhow::Error> {
313 let order = self.group().order();
314 let mut smat = Array2::<T>::zeros((order, order));
315 let item_0 = self.origin();
316 if let (Some(ctb), true) = (self.group().cayley_table(), use_cayley_table) {
317 log::debug!("Cayley table available. Group closure will be used to speed up overlap matrix computation.");
318 let ovs = self
319 .iter()
320 .map(|item_res| {
321 let item = item_res?;
322 item.overlap(item_0, metric, metric_h)
323 })
324 .collect::<Result<Vec<_>, _>>()?;
325 for (i, j) in (0..order).cartesian_product(0..order) {
326 let jinv = ctb
327 .slice(s![.., j])
328 .iter()
329 .position(|&x| x == 0)
330 .ok_or(format_err!(
331 "Unable to find the inverse of group element `{j}`."
332 ))?;
333 let jinv_i = ctb[(jinv, i)];
334 smat[(i, j)] = self.norm_preserving_scalar_map(jinv)?(ovs[jinv_i]);
335 }
336 } else {
337 log::debug!("Cayley table not available or the use of Cayley table not requested. Overlap matrix will be constructed without group-closure speed-up.");
338 for pair in self
339 .iter()
340 .map(|item_res| item_res.map_err(|err| err.to_string()))
341 .enumerate()
342 .combinations_with_replacement(2)
343 {
344 let (w, item_w_res) = &pair[0];
345 let (x, item_x_res) = &pair[1];
346 let item_w = item_w_res
347 .as_ref()
348 .map_err(|err| format_err!(err.clone()))
349 .with_context(|| "One of the items in the orbit is not available")?;
350 let item_x = item_x_res
351 .as_ref()
352 .map_err(|err| format_err!(err.clone()))
353 .with_context(|| "One of the items in the orbit is not available")?;
354 smat[(*w, *x)] = item_w.overlap(item_x, metric, metric_h).map_err(|err| {
355 log::warn!("{err}");
356 log::warn!(
357 "Unable to calculate the overlap between items `{w}` and `{x}` in the orbit."
358 );
359 err
360 })?;
361 if *w != *x {
362 smat[(*x, *w)] = item_x.overlap(item_w, metric, metric_h).map_err(|err| {
363 log::warn!("{err}");
364 log::warn!(
365 "Unable to calculate the overlap between items `{x}` and `{w}` in the orbit."
366 );
367 err
368 })?;
369 }
370 }
371 }
372 if self.origin().complex_symmetric() {
373 self.set_smat((smat.clone() + smat.t().to_owned()).mapv(|x| x / (T::one() + T::one())))
374 } else {
375 self.set_smat(
376 (smat.clone() + smat.t().to_owned().mapv(|x| x.conj()))
377 .mapv(|x| x / (T::one() + T::one())),
378 )
379 }
380 Ok(self)
381 }
382
383 fn normalise_smat(&mut self) -> Result<&mut Self, anyhow::Error> {
391 let smat = self
392 .smat()
393 .ok_or(format_err!("No orbit overlap matrix to normalise."))?;
394 let norm = smat.diag().mapv(|x| <T as ComplexFloat>::sqrt(x));
395 let nspatial = norm.len();
396 let norm_col = norm
397 .clone()
398 .into_shape_with_order([nspatial, 1])
399 .map_err(|err| format_err!(err))?;
400 let norm_row = norm
401 .into_shape_with_order([1, nspatial])
402 .map_err(|err| format_err!(err))?;
403 let norm_mat = norm_col.dot(&norm_row);
404 let normalised_smat = smat / norm_mat;
405 self.set_smat(normalised_smat);
406 Ok(self)
407 }
408
409 #[must_use]
431 fn calc_tmat(&self, op: &G::GroupElement) -> Result<Array2<T>, anyhow::Error> {
432 let ctb = self
433 .group()
434 .cayley_table()
435 .expect("The Cayley table for the group cannot be found.");
436 let i = self.group().get_index_of(op).unwrap_or_else(|| {
437 panic!("Unable to retrieve the index of element `{op}` in the group.")
438 });
439 let ix = ctb.slice(s![i, ..]).iter().cloned().collect::<Vec<_>>();
440 let twx = self
441 .smat()
442 .ok_or(format_err!("No orbit overlap matrix found."))?
443 .select(Axis(1), &ix);
444 Ok(twx)
445 }
446
447 #[must_use]
467 fn calc_dmat(&self, op: &G::GroupElement) -> Result<Array2<T>, anyhow::Error> {
468 let complex_symmetric = self.origin().complex_symmetric();
469 let xmath = if complex_symmetric {
470 self.xmat().t().to_owned()
471 } else {
472 self.xmat().t().mapv(|x| x.conj())
473 };
474 let smattilde = xmath
475 .dot(
476 self.smat()
477 .ok_or(format_err!("No orbit overlap matrix found."))?,
478 )
479 .dot(self.xmat());
480 let smattilde_inv = smattilde
481 .inv()
482 .expect("The inverse of S~ could not be found.");
483 let dmat = einsum(
484 "ij,jk,kl,lm->im",
485 &[&smattilde_inv, &xmath, &self.calc_tmat(op)?, self.xmat()],
486 )
487 .map_err(|err| format_err!(err))
488 .with_context(|| "Unable to compute the matrix product [(S~)^(-1) X† T X].")?
489 .into_dimensionality::<Ix2>()
490 .map_err(|err| format_err!(err))
491 .with_context(|| {
492 "Unable to convert the matrix product [(S~)^(-1) X† T X] to two dimensions."
493 });
494 dmat
495 }
496
497 #[must_use]
510 fn calc_character(&self, op: &G::GroupElement) -> Result<T, anyhow::Error> {
511 let complex_symmetric = self.origin().complex_symmetric();
512 let xmath = if complex_symmetric {
513 self.xmat().t().to_owned()
514 } else {
515 self.xmat().t().mapv(|x| x.conj())
516 };
517 let smattilde = xmath
518 .dot(
519 self.smat()
520 .ok_or(format_err!("No orbit overlap matrix found."))?,
521 )
522 .dot(self.xmat());
523 let smattilde_inv = smattilde
524 .inv()
525 .expect("The inverse of S~ could not be found.");
526 let chi = einsum(
527 "ij,jk,kl,li",
528 &[&smattilde_inv, &xmath, &self.calc_tmat(op)?, self.xmat()],
529 )
530 .map_err(|err| format_err!(err))
531 .with_context(|| "Unable to compute the trace of the matrix product [(S~)^(-1) X† T X].")?
532 .into_dimensionality::<Ix0>()
533 .map_err(|err| format_err!(err))
534 .with_context(|| "Unable to convert the trace of the matrix product [(S~)^(-1) X† T X] to zero dimensions.")?;
535 chi.into_iter().next().ok_or(format_err!(
536 "Unable to extract the character from the representation matrix."
537 ))
538 }
539
540 #[must_use]
549 fn calc_characters(
550 &self,
551 ) -> Result<Vec<(<G as ClassProperties>::ClassSymbol, T)>, anyhow::Error> {
552 let complex_symmetric = self.origin().complex_symmetric();
553 let xmath = if complex_symmetric {
554 self.xmat().t().to_owned()
555 } else {
556 self.xmat().t().mapv(|x| x.conj())
557 };
558 let smattilde = xmath
559 .dot(
560 self.smat()
561 .ok_or(format_err!("No orbit overlap matrix found."))?,
562 )
563 .dot(self.xmat());
564 let smattilde_inv = smattilde
565 .inv()
566 .expect("The inverse of S~ could not be found.");
567 let chis = (0..self.group().class_number()).map(|cc_i| {
568 let cc = self.group().get_cc_symbol_of_index(cc_i).unwrap();
569 let op = self.group().get_cc_transversal(cc_i).unwrap();
570 let chi = einsum(
571 "ij,jk,kl,li",
572 &[&smattilde_inv, &xmath, &self.calc_tmat(&op)?, self.xmat()],
573 )
574 .map_err(|err| format_err!(err))
575 .with_context(|| "Unable to compute the trace of the matrix product [(S~)^(-1) X† T X].")?
576 .into_dimensionality::<Ix0>()
577 .map_err(|err| format_err!(err))
578 .with_context(|| "Unable to convert the trace of the matrix product [(S~)^(-1) X† T X] to zero dimensions.")?;
579 let chi_val = chi.into_iter().next().ok_or(format_err!(
580 "Unable to extract the character from the representation matrix."
581 ))?;
582 Ok((cc, chi_val))
583 }).collect::<Result<Vec<_>, _>>();
584 chis
585 }
586
587 fn analyse_rep(
599 &self,
600 ) -> Result<
601 <<G as CharacterProperties>::CharTab as SubspaceDecomposable<T>>::Decomposition,
602 DecompositionError,
603 > {
604 let chis = self
605 .calc_characters()
606 .map_err(|err| DecompositionError(err.to_string()))?;
607 let res = self.group().character_table().reduce_characters(
608 &chis.iter().map(|(cc, chi)| (cc, *chi)).collect::<Vec<_>>(),
609 self.integrality_threshold(),
610 );
611 res
612 }
613
614 fn characters_to_string(
625 &self,
626 chis: &[(<G as ClassProperties>::ClassSymbol, T)],
627 integrality_threshold: <T as ComplexFloat>::Real,
628 ) -> String
629 where
630 T: ComplexFloat + Lapack + fmt::Debug,
631 <T as ComplexFloat>::Real: ToPrimitive,
632 {
633 let ndigits = (-integrality_threshold.log10()).to_usize().unwrap_or(6) + 1;
634 let (ccs, characters): (Vec<_>, Vec<_>) = chis
635 .iter()
636 .map(|(cc, chi)| (cc.to_string(), format!("{chi:+.ndigits$}")))
637 .unzip();
638 let cc_width = ccs
639 .iter()
640 .map(|cc| cc.chars().count())
641 .max()
642 .unwrap_or(5)
643 .max(5);
644 let characters_width = characters
645 .iter()
646 .map(|chi| chi.chars().count())
647 .max()
648 .unwrap_or(9)
649 .max(9);
650
651 let divider = "┈".repeat(cc_width + characters_width + 4);
652 let header = format!(" {:<cc_width$} {:<}", "Class", "Character");
653 let body = Itertools::intersperse(
654 chis.iter()
655 .map(|(cc, chi)| format!(" {:<cc_width$} {:<+.ndigits$}", cc.to_string(), chi)),
656 "\n".to_string(),
657 )
658 .collect::<String>();
659
660 Itertools::intersperse(
661 [divider.clone(), header, divider.clone(), body, divider].into_iter(),
662 "\n".to_string(),
663 )
664 .collect::<String>()
665 }
666}
667
668pub trait ProjectionDecomposition<G, I, T, D>: RepAnalysis<G, I, T, D>
675where
676 T: ComplexFloat + Lapack + fmt::Debug,
677 <T as ComplexFloat>::Real: ToPrimitive,
678 G: GroupProperties + ClassProperties + CharacterProperties + Sync + Send,
679 <G as CharacterProperties>::RowSymbol: Sync + Send,
680 G::GroupElement: fmt::Display,
681 G::CharTab: SubspaceDecomposable<T>,
682 D: Dimension,
683 I: Overlap<T, D> + Clone,
684 Self: Sync + Send,
685 Self::OrbitIter: Iterator<Item = Result<I, anyhow::Error>>,
686 for<'a> Complex<f64>: Mul<&'a T, Output = Complex<f64>>,
687{
688 fn calc_projection_compositions(
708 &self,
709 ) -> Result<
710 Vec<(
711 <<G as CharacterProperties>::CharTab as CharacterTable>::RowSymbol,
712 Complex<f64>,
713 )>,
714 anyhow::Error,
715 > {
716 self.group()
717 .character_table()
718 .get_all_rows()
719 .iter()
720 .map(|row_label| {
721 let group = self.group();
722 let chartab = group.character_table();
723 let id_class = group
724 .get_cc_symbol_of_index(0)
725 .ok_or(format_err!("Unable to retrieve the identity class."))?;
726 let dim = chartab
727 .get_character(&row_label, &id_class).complex_value();
728 let group_order = group.order().to_f64().ok_or(
729 DecompositionError("The group order cannot be converted to `f64`.".to_string())
730 )?;
731 let projection: Complex<f64> = (0..group.order())
732 .into_par_iter()
733 .try_fold(|| Complex::<f64>::zero(), |acc, i| {
734 let cc_i = group
735 .get_cc_of_element_index(i)
736 .ok_or(format_err!("Unable to retrieve the conjugacy class index of element index {i}."))?;
737 let cc = group
738 .get_cc_symbol_of_index(cc_i)
739 .ok_or(format_err!("Unable to retrieve the conjugacy class symbol of conjugacy class index {cc_i}."))?;
740 let chi_i_star = group.character_table().get_character(&row_label, &cc).complex_conjugate();
741 let s_0i = self
742 .smat()
743 .ok_or(format_err!("No orbit overlap matrix found."))?
744 .get((0, i))
745 .ok_or(format_err!("Unable to retrieve the overlap matrix element with index `(0, {i})`."))?;
746 Ok::<_, anyhow::Error>(acc + chi_i_star.complex_value() * s_0i)
747 })
748 .try_reduce(|| Complex::<f64>::zero(), |a, s| Ok(a + s))?;
749 Ok::<_, anyhow::Error>((row_label.clone(), &dim * projection / group_order))
750 }).collect::<Result<Vec<_>, anyhow::Error>>()
751 }
752
753 fn projections_to_string(
765 &self,
766 projections: &[(
767 <<G as CharacterProperties>::CharTab as CharacterTable>::RowSymbol,
768 Complex<f64>,
769 )],
770 integrality_threshold: <T as ComplexFloat>::Real,
771 ) -> String {
772 let ndigits = (-integrality_threshold.log10()).to_usize().unwrap_or(6) + 1;
773 let (row_labels, projection_values): (Vec<_>, Vec<_>) = projections
774 .iter()
775 .map(|(row, proj)| (row.to_string(), format!("{proj:+.ndigits$}")))
776 .unzip();
777 let row_width = row_labels
778 .iter()
779 .map(|row| row.chars().count())
780 .max()
781 .unwrap_or(9)
782 .max(9);
783 let proj_width = projection_values
784 .iter()
785 .map(|proj| proj.chars().count())
786 .max()
787 .unwrap_or(10)
788 .max(10);
789
790 let divider = "┈".repeat(row_width + proj_width + 4);
791 let header = format!(" {:<row_width$} {:<}", "Component", "Projection");
792 let body = Itertools::intersperse(
793 projections.iter().map(|(row, proj)| {
794 format!(" {:<row_width$} {:<+.ndigits$}", row.to_string(), proj)
795 }),
796 "\n".to_string(),
797 )
798 .collect::<String>();
799
800 Itertools::intersperse(
801 [divider.clone(), header, divider.clone(), body, divider].into_iter(),
802 "\n".to_string(),
803 )
804 .collect::<String>()
805 }
806}
807
808#[duplicate_item(
811 duplicate!{
812 [ dtype_nested; [f64]; [Complex<f64>] ]
813 [
814 gtype_ [ UnitaryRepresentedSymmetryGroup ]
815 dtype_ [ dtype_nested ]
816 ]
817 }
818)]
819impl<O, I, D> ProjectionDecomposition<gtype_, I, dtype_, D> for O
820where
821 O: RepAnalysis<gtype_, I, dtype_, D>,
822 D: Dimension,
823 I: Overlap<dtype_, D> + Clone,
824 Self: Sync + Send,
825 Self::OrbitIter: Iterator<Item = Result<I, anyhow::Error>>,
826{
827}
828
829macro_rules! fn_calc_xmat_real {
834 ( $(#[$meta:meta])* $vis:vis $func:ident ) => {
835 $(#[$meta])*
836 $vis fn $func(&mut self, preserves_full_rank: bool) -> Result<&mut Self, anyhow::Error> {
837 let thresh = self.linear_independence_threshold;
839 let smat = self
840 .smat
841 .as_ref()
842 .ok_or(format_err!("No overlap matrix found for this orbit."))?;
843 use ndarray_linalg::norm::Norm;
844 if (smat.to_owned() - smat.t()).norm_l2() > thresh {
845 Err(format_err!("Overlap matrix is not symmetric."))
846 } else {
847 let (s_eig, umat) = smat.eigh(UPLO::Lower).map_err(|err| format_err!(err))?;
848 let nonzero_s_indices = match self.eigenvalue_comparison_mode {
849 EigenvalueComparisonMode::Modulus => {
850 s_eig.iter().positions(|x| x.abs() > thresh).collect_vec()
851 }
852 EigenvalueComparisonMode::Real => {
853 s_eig.iter().positions(|x| *x > thresh).collect_vec()
854 }
855 };
856 let nonzero_s_eig = s_eig.select(Axis(0), &nonzero_s_indices);
857 let nonzero_umat = umat.select(Axis(1), &nonzero_s_indices);
858 let nullity = smat.shape()[0] - nonzero_s_indices.len();
859 let xmat = if nullity == 0 && preserves_full_rank {
860 Array2::eye(smat.shape()[0])
861 } else {
862 let s_s = Array2::<f64>::from_diag(&nonzero_s_eig.mapv(|x| 1.0 / x.sqrt()));
863 nonzero_umat.dot(&s_s)
864 };
865 self.smat_eigvals = Some(s_eig);
866 self.xmat = Some(xmat);
867 Ok(self)
868 }
869 }
870 }
871}
872
873macro_rules! fn_calc_xmat_complex {
874 ( $(#[$meta:meta])* $vis:vis $func:ident ) => {
875 $(#[$meta])*
876 $vis fn $func(&mut self, preserves_full_rank: bool) -> Result<&mut Self, anyhow::Error> {
877 let thresh = self.linear_independence_threshold;
879 let smat = self
880 .smat
881 .as_ref()
882 .ok_or(format_err!("No overlap matrix found for this orbit."))?;
883 let (s_eig, umat_nonortho) = smat.eig().map_err(|err| format_err!(err))?;
884
885 let nonzero_s_indices = match self.eigenvalue_comparison_mode {
886 EigenvalueComparisonMode::Modulus => s_eig
887 .iter()
888 .positions(|x| ComplexFloat::abs(*x) > thresh)
889 .collect_vec(),
890 EigenvalueComparisonMode::Real => {
891 if s_eig
892 .iter()
893 .any(|x| Float::abs(ComplexFloat::im(*x)) > thresh)
894 {
895 log::warn!("Comparing eigenvalues using the real parts, but not all eigenvalues are real.");
896 }
897 s_eig
898 .iter()
899 .positions(|x| ComplexFloat::re(*x) > thresh)
900 .collect_vec()
901 }
902 };
903 let nonzero_s_eig = s_eig.select(Axis(0), &nonzero_s_indices);
904 let nonzero_umat_nonortho = umat_nonortho.select(Axis(1), &nonzero_s_indices);
905
906 let nonzero_umat = complex_modified_gram_schmidt(
909 &nonzero_umat_nonortho,
910 self.origin.complex_symmetric(),
911 thresh,
912 )
913 .map_err(
914 |_| format_err!("Unable to orthonormalise the linearly-independent eigenvectors of the overlap matrix.")
915 )?;
916
917 let nullity = smat.shape()[0] - nonzero_s_indices.len();
918 let xmat = if nullity == 0 && preserves_full_rank {
919 Array2::<Complex<T>>::eye(smat.shape()[0])
920 } else {
921 let s_s = Array2::<Complex<T>>::from_diag(
922 &nonzero_s_eig.mapv(|x| Complex::<T>::from(T::one()) / x.sqrt()),
923 );
924 nonzero_umat.dot(&s_s)
925 };
926 self.smat_eigvals = Some(s_eig);
927 self.xmat = Some(xmat);
928 Ok(self)
929 }
930 }
931}
932
933pub(crate) use fn_calc_xmat_complex;
934pub(crate) use fn_calc_xmat_real;
935
936pub(crate) fn log_overlap_eigenvalues<T>(
947 title: &str,
948 eigvals: &Array1<T>,
949 thresh: <T as ComplexFloat>::Real,
950 eigenvalue_comparison_mode: &EigenvalueComparisonMode,
951) where
952 T: std::fmt::LowerExp + ComplexFloat,
953 <T as ComplexFloat>::Real: std::fmt::LowerExp,
954{
955 let mut eigvals_sorted = eigvals.iter().collect::<Vec<_>>();
956 match eigenvalue_comparison_mode {
957 EigenvalueComparisonMode::Modulus => {
958 eigvals_sorted.sort_by(|a, b| a.abs().partial_cmp(&b.abs()).unwrap());
959 }
960 EigenvalueComparisonMode::Real => {
961 eigvals_sorted.sort_by(|a, b| a.re().partial_cmp(&b.re()).unwrap());
962 }
963 }
964 eigvals_sorted.reverse();
965 let eigvals_str = eigvals_sorted
966 .iter()
967 .map(|v| format!("{v:+.3e}"))
968 .collect::<Vec<_>>();
969 log_subtitle(title);
970 qsym2_output!("");
971
972 match eigenvalue_comparison_mode {
973 EigenvalueComparisonMode::Modulus => {
974 qsym2_output!("Eigenvalues are sorted in decreasing order of their moduli.");
975 }
976 EigenvalueComparisonMode::Real => {
977 qsym2_output!("Eigenvalues are sorted in decreasing order of their real parts.");
978 }
979 }
980 let count_length = usize::try_from(eigvals.len().ilog10() + 2).unwrap_or(2);
981 let eigval_length = eigvals_str
982 .iter()
983 .map(|v| v.chars().count())
984 .max()
985 .unwrap_or(20);
986 qsym2_output!("{}", "┈".repeat(count_length + 3 + eigval_length));
987 qsym2_output!("{:>count_length$} Eigenvalue", "#");
988 qsym2_output!("{}", "┈".repeat(count_length + 3 + eigval_length));
989 let mut write_thresh = false;
990 for (i, eigval) in eigvals_str.iter().enumerate() {
991 let cmp = match eigenvalue_comparison_mode {
992 EigenvalueComparisonMode::Modulus => {
993 eigvals_sorted[i].abs().partial_cmp(&thresh).expect(
994 "Unable to compare the modulus of an eigenvalue with the specified threshold.",
995 )
996 }
997 EigenvalueComparisonMode::Real => eigvals_sorted[i].re().partial_cmp(&thresh).expect(
998 "Unable to compare the real part of an eigenvalue with the specified threshold.",
999 ),
1000 };
1001 if cmp == Ordering::Less && !write_thresh {
1002 let comparison_mode_str = match eigenvalue_comparison_mode {
1003 EigenvalueComparisonMode::Modulus => "modulus-based",
1004 EigenvalueComparisonMode::Real => "real-part-based",
1005 };
1006 qsym2_output!(
1007 "{} <-- linear independence threshold ({comparison_mode_str}): {:+.3e}",
1008 "-".repeat(count_length + 3 + eigval_length),
1009 thresh
1010 );
1011 write_thresh = true;
1012 }
1013 qsym2_output!("{i:>count_length$} {eigval}",);
1014 }
1015 qsym2_output!("{}", "┈".repeat(count_length + 3 + eigval_length));
1016}