sparse_ir_minipole/minipole/
mini_pole_dlr.rs

1//! MPM-DLR: minimal poles from DLR coefficients (port of
2//! `mini_pole/mini_pole_dlr.py`).
3//!
4//! Ported from Green-Phys/MiniPole (commit 15e4a54, MIT License,
5//! Copyright (c) 2024 lzphy); see `LICENSE-THIRD-PARTY`.
6
7use super::con_map::{ConMap, ConMapGapless, ConMapGeneric};
8use super::{MiniPoleResult, assemble};
9use crate::basis_trait::Basis;
10use crate::dlr::DiscreteLehmannRepresentation;
11use crate::error::{ArrayRole, Error, Result};
12use crate::esprit::{ErrType, Esprit, EspritParams};
13use crate::traits::StatisticsType;
14use num_complex::Complex;
15use tenferro_tensor::{TensorScalar, TypedTensor};
16
17type C64 = Complex<f64>;
18
19/// Parameters of [`mini_pole_dlr`], with the defaults of the reference.
20#[derive(Debug, Clone, PartialEq)]
21pub struct MiniPoleDlrParams {
22    /// Index of the lowest Matsubara frequency `ω_{n0} = (2 n0 + 1)π/β` of
23    /// the contour, typically in `0..10`.
24    pub n0: usize,
25    /// Cutoff `nmax` of the upper end `(2 nmax + 1)π/β` of the contour when
26    /// `symmetry` is false; `None` uses `β`.
27    pub nmax: Option<f64>,
28    /// Error tolerance of ESPRIT. One of `err` and `m` is required.
29    pub err: Option<f64>,
30    /// Interpretation of `err` (default absolute).
31    pub err_type: ErrType,
32    /// Number of poles.
33    pub m: Option<usize>,
34    /// Impose up-down symmetry (the gapless map).
35    pub symmetry: bool,
36    /// Maximum number of moments (default 200).
37    pub k_max: usize,
38    /// Ratio `L/N` of ESPRIT (default 0.4).
39    pub lfactor: f64,
40}
41
42impl MiniPoleDlrParams {
43    /// Parameters with contour start `n0`, tolerance `err` and defaults
44    /// otherwise.
45    pub fn new(n0: usize, err: f64) -> Self {
46        Self {
47            n0,
48            nmax: None,
49            err: Some(err),
50            err_type: ErrType::Abs,
51            m: None,
52            symmetry: false,
53            k_max: 200,
54            lfactor: 0.4,
55        }
56    }
57}
58
59/// Minimal pole representation of `G(z) = Σ_l A_l / (z - x_l)`.
60///
61/// `al` has shape `[r, ...]` (residues along axis 0, trailing channel axes),
62/// `xl` the `r` real pole locations.
63///
64/// # Errors
65/// [`Error::ShapeMismatch`] if `al` does not have `len(xl)` rows,
66/// [`Error::InvalidParameter`] for invalid parameters, and ESPRIT errors.
67pub fn mini_pole_dlr<T>(
68    al: &TypedTensor<T>,
69    xl: &[f64],
70    beta: f64,
71    params: &MiniPoleDlrParams,
72) -> Result<MiniPoleResult>
73where
74    T: TensorScalar + Copy + Into<C64>,
75{
76    let shape = al.shape().to_vec();
77    let r = xl.len();
78    if shape.first() != Some(&r) {
79        let mut expected = shape.clone();
80        match expected.first_mut() {
81            Some(n) => *n = r,
82            None => expected.push(r),
83        }
84        return Err(Error::ShapeMismatch {
85            which: ArrayRole::Input,
86            expected,
87            actual: shape,
88        });
89    }
90    if !(beta.is_finite() && beta > 0.0) {
91        return Err(Error::InvalidParameter {
92            name: "beta",
93            value: format!("{beta:?}"),
94            reason: "must be finite and positive".to_string(),
95        });
96    }
97    let d: usize = shape[1..].iter().product();
98    let al: Vec<C64> = al.host_data()?.iter().map(|&v| v.into()).collect();
99
100    let pi = std::f64::consts::PI;
101    let w_n0 = (2.0 * params.n0 as f64 + 1.0) * pi / beta;
102    let generic;
103    let gapless;
104    let map: &dyn ConMap = if !params.symmetry {
105        let nmax = params.nmax.unwrap_or(beta);
106        let w_nmax = (2.0 * nmax + 1.0) * pi / beta;
107        generic = ConMapGeneric {
108            w_m: 0.5 * (w_n0 + w_nmax),
109            dw_h: 0.5 * (w_nmax - w_n0),
110        };
111        if generic.dw_h.is_nan() || generic.dw_h <= 0.0 {
112            return Err(Error::InvalidParameter {
113                name: "nmax",
114                value: format!("{nmax:?}"),
115                reason: format!("must exceed n0 = {}", params.n0),
116            });
117        }
118        &generic
119    } else {
120        gapless = ConMapGapless { w_min: w_n0 };
121        &gapless
122    };
123
124    // h_k = Σ_l xl_p^k A_l / z'(xl_p), k < min(int((r + 1) / Lfactor), k_max).
125    let xl_p: Vec<C64> = xl.iter().map(|&x| map.w(C64::new(x, 0.0))).collect();
126    let n = (((r + 1) as f64 / params.lfactor) as usize).min(params.k_max);
127    let mut h = vec![C64::new(0.0, 0.0); n * d];
128    for (l, &q) in xl_p.iter().enumerate() {
129        let f = map.dz(q).inv();
130        let mut p = C64::new(1.0, 0.0);
131        for k in 0..n {
132            for c in 0..d {
133                h[k + n * c] += p * (al[l + r * c] * f);
134            }
135            p *= q;
136        }
137    }
138
139    let esprit = Esprit::new(
140        &h,
141        n,
142        d,
143        &EspritParams {
144            err: params.err,
145            err_type: params.err_type,
146            m: params.m,
147            lfactor: params.lfactor,
148            ..EspritParams::default()
149        },
150    )?;
151    let keep: Vec<usize> = (0..esprit.gamma.len())
152        .filter(|&j| esprit.gamma[j].norm() < 1.0)
153        .collect();
154    let location: Vec<C64> = keep.iter().map(|&j| map.z(esprit.gamma[j])).collect();
155    let mut weight = vec![C64::new(0.0, 0.0); keep.len() * d];
156    for (jn, &j) in keep.iter().enumerate() {
157        let f = map.dz(esprit.gamma[j]);
158        for c in 0..d {
159            weight[jn + keep.len() * c] = esprit.omega[j + esprit.m * c] * f;
160        }
161    }
162    let mut hshape = shape;
163    hshape[0] = n;
164    assemble(
165        location,
166        weight,
167        vec![C64::new(0.0, 0.0); d],
168        h,
169        hshape,
170        esprit,
171        params.n0,
172        None,
173    )
174}
175
176/// [`mini_pole_dlr`] of DLR coefficients: `A_l = g_l w_l` at the DLR poles
177/// `x_l`, with `w_l` the pole weights of the DLR.
178///
179/// # Errors
180/// See [`mini_pole_dlr`].
181pub fn mini_pole_dlr_from<S, T>(
182    dlr: &DiscreteLehmannRepresentation<S>,
183    g_dlr: &TypedTensor<T>,
184    params: &MiniPoleDlrParams,
185) -> Result<MiniPoleResult>
186where
187    S: StatisticsType + 'static,
188    T: TensorScalar + Copy + Into<C64>,
189{
190    let shape = g_dlr.shape().to_vec();
191    let r = dlr.poles().len();
192    if shape.first() != Some(&r) {
193        let mut expected = shape.clone();
194        match expected.first_mut() {
195            Some(n) => *n = r,
196            None => expected.push(r),
197        }
198        return Err(Error::ShapeMismatch {
199            which: ArrayRole::Input,
200            expected,
201            actual: shape,
202        });
203    }
204    let weights = dlr.pole_weights();
205    let g = g_dlr.host_data()?;
206    let al: Vec<C64> = g
207        .iter()
208        .enumerate()
209        .map(|(i, &v)| v.into() * weights[i % r])
210        .collect();
211    let al = TypedTensor::from_vec_col_major(shape, al)?;
212    mini_pole_dlr(&al, dlr.poles(), dlr.beta(), params)
213}
sparse_ir_minipole/minipole/mini_pole_dlr.rs

sparse_ir_minipole/minipole/
mini_pole_dlr.rs
