dace_rs/
io.rs

1//! Text and binary I/O of [`Da`] values: the `daceWrite` display format and
2//! its parser (`daceRead`), the binary blob format (`daceExportBlob` /
3//! `daceImportBlob`), and the code-generation formatter
4//! (`DASimpleFormatter`).
5//!
6//! The blob format preserves the C layout: magic `0x1E304144`, then `no`,
7//! `nv1`, `nv2`, `len` as little-endian `u32`s, then `len` packed
8//! `{i1, i2, cc}` records (`u32, u32, f64`, little-endian). All supported
9//! targets are little-endian, matching the C library on the reference
10//! platforms.
11
12use std::fmt;
13
14use crate::da::Da;
15
16const BEGSTR: &str = "     I  COEFFICIENT              ORDER EXPONENTS";
17const ENDSTR: &str = "------------------------------------------------";
18const ZEROSTR: &str = "        ALL COEFFICIENTS ZERO";
19
20/// Format a coefficient exactly like C's `%24.16e`: 16 fractional digits and
21/// a signed two-digit exponent (e.g. `1.0000000000000000e+00`).
22fn fmt_c_e(c: f64) -> String {
23    let s = format!("{c:.16e}");
24    let (mantissa, exponent) = s.split_once('e').expect("scientific notation");
25    let exp: i32 = exponent.parse().unwrap_or(0);
26    format!(
27        "{mantissa}e{}{:02}",
28        if exp < 0 { '-' } else { '+' },
29        exp.abs()
30    )
31}
32
33impl fmt::Display for Da {
34    fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
35        if self.terms.is_empty() {
36            writeln!(f, "{ZEROSTR}")?;
37            writeln!(f, "{ENDSTR}")?;
38            return Ok(());
39        }
40
41        writeln!(f, "{BEGSTR}")?;
42        let ctx = &self.ctx;
43        let mut jj = vec![0u32; ctx.nvmax as usize];
44        let mut iout: usize = 1;
45        for ioa in 0..=ctx.nomax {
46            for t in &self.terms {
47                if ctx.order_of(t.idx) != ioa {
48                    continue;
49                }
50                ctx.decode_into(t.idx, &mut jj);
51                write!(f, "{iout:6}  {:>24}", fmt_c_e(t.c))?;
52                write!(f, "{ioa:4} ")?;
53                for &e in &jj {
54                    write!(f, " {e:2}")?;
55                }
56                writeln!(f)?;
57                iout += 1;
58            }
59        }
60        writeln!(f, "{ENDSTR}")?;
61        Ok(())
62    }
63}
64
65impl std::str::FromStr for Da {
66    type Err = crate::error::DaceError;
67
68    /// Parse the text format emitted by [`Display`][fmt::Display] (the
69    /// `daceRead` parser, including its tolerance: line-number and order
70    /// mismatches warn, order/truncation-excess rows are skipped, duplicate
71    /// monomials accumulate).
72    ///
73    /// # Panics
74    ///
75    /// Panics with [`crate::DaceError`] if DACE has not been initialized.
76    fn from_str(s: &str) -> Result<Da, crate::error::DaceError> {
77        let lines: Vec<&str> = s.lines().collect();
78        if lines.is_empty() {
79            return Err(crate::error::DaceError::new(
80                634,
81                "Not enough lines to read",
82            ));
83        }
84        let first = lines[0];
85        if first.starts_with(ZEROSTR) || first.starts_with("        ALL COMPONENTS ZERO") {
86            return Ok(Da::new());
87        }
88        // DACE and COSY header differ only in the coefficient field width;
89        // dace-rs emits only the DACE one but accepts both.
90        let cosy = first.starts_with("     I  COEFFICIENT            ORDER EXPONENTS");
91        let dace = first.starts_with(BEGSTR);
92        if !cosy && !dace {
93            return Err(crate::error::DaceError::new(632, "Unknown format"));
94        }
95        let coefflen = if cosy { 22 } else { 24 };
96
97        let ctx = crate::context::Context::current();
98        let (_eps, nocut) = crate::context::eps_nocut();
99        let mut cc = vec![0.0; ctx.nmmax as usize];
100        let mut jj = vec![0u32; ctx.nvmax as usize];
101
102        for (iin, line) in lines.iter().enumerate().skip(1) {
103            if line.len() < 4 {
104                return Err(crate::error::DaceError::new(632, "Unknown format"));
105            }
106            if line[4..].starts_with(ENDSTR) {
107                break;
108            }
109            let b = line.as_bytes();
110            if line.len() < 37 {
111                return Err(crate::error::DaceError::new(632, "Unknown format"));
112            }
113            // line number (columns 0-5)
114            let ii: u32 = std::str::from_utf8(&b[..6])
115                .ok()
116                .and_then(|t| t.trim().parse().ok())
117                .unwrap_or(0);
118            // coefficient (columns 8..8+coefflen)
119            let c: f64 = std::str::from_utf8(&b[8..8 + coefflen])
120                .ok()
121                .and_then(|t| t.trim().parse().ok())
122                .unwrap_or(0.0);
123            // order (columns 32-35)
124            let io1: u32 = std::str::from_utf8(&b[32..36])
125                .ok()
126                .and_then(|t| t.trim().parse().ok())
127                .unwrap_or(0);
128            // exponents: walk sequentially (1 space, then 2 digits per
129            // variable; in COSY format only every other variable has the
130            // leading space)
131            let mut pos = 37usize;
132            for (i, slot) in jj.iter_mut().enumerate() {
133                *slot = 0;
134                if line.len() > pos && (!cosy || i % 2 == 0) {
135                    pos += 1;
136                }
137                if line.len() >= pos + 2 {
138                    *slot = std::str::from_utf8(&b[pos..pos + 2])
139                        .ok()
140                        .and_then(|t| t.trim().parse().ok())
141                        .unwrap_or(0);
142                    pos += 2;
143                }
144            }
145
146            // check line numbers (informational in C)
147            if ii != iin as u32 {
148                log::warn!("DACE info 164: numbering out of order while reading");
149            }
150            // check order and hence number of variables
151            let io: u32 = jj.iter().sum();
152            if io != io1 {
153                log::warn!("DACE info 165: inaccurate estimate while reading; line skipped");
154                continue;
155            }
156            // check cutoff order
157            if io > nocut {
158                continue;
159            }
160            let icc = ctx.encode(&jj).expect("validated order");
161            if cc[icc as usize] != 0.0 {
162                log::warn!("DACE info 166: duplicate monomial while reading");
163            }
164            cc[icc as usize] += c;
165        }
166
167        Ok(crate::eval::pack(&ctx, &mut cc))
168    }
169}
170
171/// Binary magic `0x1E304144` (little-endian "DA0" + record separator).
172const DACE_BINARY_MAGIC: u32 = 0x1E30_4144;
173
174impl Da {
175    /// Export in the C binary blob format (`daceExportBlob`): magic, `no`,
176    /// `nv1`, `nv2`, `len` as little-endian `u32`s, then `len` packed
177    /// `{i1, i2, cc}` records. As in C, the buffer always includes one
178    /// record slot (zeroed when the DA is empty).
179    ///
180    /// All supported targets are little-endian, so blobs interoperate with
181    /// the C library on the reference platforms.
182    pub fn to_blob(&self) -> Vec<u8> {
183        let ctx = &self.ctx;
184        let len = self.terms.len();
185        let mut out = Vec::with_capacity(20 + 16 * len.max(1));
186        out.extend_from_slice(&DACE_BINARY_MAGIC.to_le_bytes());
187        out.extend_from_slice(&ctx.nomax.to_le_bytes());
188        out.extend_from_slice(&ctx.nv1.to_le_bytes());
189        out.extend_from_slice(&ctx.nv2.to_le_bytes());
190        out.extend_from_slice(&(len as u32).to_le_bytes());
191        for t in &self.terms {
192            out.extend_from_slice(&ctx.ie1[t.idx as usize].to_le_bytes());
193            out.extend_from_slice(&ctx.ie2[t.idx as usize].to_le_bytes());
194            out.extend_from_slice(&t.c.to_le_bytes());
195        }
196        // C always reserves one record slot in the header struct.
197        if len == 0 {
198            out.extend_from_slice(&[0u8; 16]);
199        }
200        out
201    }
202
203    /// Import from the C binary blob format (`daceImportBlob`): silently
204    /// truncates orders above the current maximum computation order and any
205    /// extra variables present, as in C.
206    ///
207    /// # Errors
208    ///
209    /// Returns [`crate::DaceError`] code 631 ("Invalid data") when the magic is
210    /// wrong or the buffer is too short.
211    pub fn from_blob(blob: &[u8]) -> Result<Da, crate::error::DaceError> {
212        let invalid = || crate::error::DaceError::new(631, "Invalid data");
213        if blob.len() < 20 {
214            return Err(invalid());
215        }
216        let magic = u32::from_le_bytes(blob[0..4].try_into().unwrap());
217        if magic != DACE_BINARY_MAGIC {
218            return Err(invalid());
219        }
220        let no = u32::from_le_bytes(blob[4..8].try_into().unwrap());
221        let nv1 = u32::from_le_bytes(blob[8..12].try_into().unwrap());
222        let nv2 = u32::from_le_bytes(blob[12..16].try_into().unwrap());
223        let len = u32::from_le_bytes(blob[16..20].try_into().unwrap()) as usize;
224        if blob.len() < 20 + 16 * len {
225            return Err(invalid());
226        }
227
228        let ctx = crate::context::Context::current();
229        let nv = nv1 + nv2;
230        let mut p = vec![0u32; nv.max(ctx.nvmax) as usize];
231        let mut cc = vec![0.0; ctx.nmmax as usize];
232
233        for i in 0..len {
234            let off = 20 + 16 * i;
235            let i1 = u32::from_le_bytes(blob[off..off + 4].try_into().unwrap());
236            let i2 = u32::from_le_bytes(blob[off + 4..off + 8].try_into().unwrap());
237            let c = f64::from_le_bytes(blob[off + 8..off + 16].try_into().unwrap());
238
239            // decode with the blob's parameters
240            let base = no + 1;
241            let mut order = 0u32;
242            let mut ic = i1;
243            for slot in p[..nv1 as usize].iter_mut() {
244                *slot = ic % base;
245                ic /= base;
246                order += *slot;
247            }
248            let mut ic = i2;
249            for slot in p[nv1 as usize..nv as usize].iter_mut() {
250                *slot = ic % base;
251                ic /= base;
252                order += *slot;
253            }
254
255            // order of variables outside the current setup
256            let extravar: u32 = p[ctx.nvmax as usize..nv as usize].iter().sum();
257
258            if order <= ctx.nomax && extravar == 0 {
259                let idx = ctx.encode(&p[..ctx.nvmax as usize]).expect("order checked");
260                cc[idx as usize] = c;
261            }
262        }
263
264        Ok(crate::eval::pack(&ctx, &mut cc))
265    }
266}
267
268// ---------------------------------------------------------------------------
269// Code-generation formatter (C++ DASimpleFormatter)
270// ---------------------------------------------------------------------------
271
272/// The elements of a simple code-generation format (C++
273/// `DASimpleFormat`): strings used around signs, products, variables, and
274/// powers, plus line wrapping.
275#[derive(Debug, Clone)]
276pub struct SimpleFormat {
277    /// Prefix of a positive term (e.g. `"+"`).
278    pub pos: String,
279    /// Prefix of a negative term with the coefficient negated (e.g. `"-"`).
280    pub neg: String,
281    /// Multiplication separator (e.g. `"*"`).
282    pub mul: String,
283    /// Power prefix (e.g. `"pow("`).
284    pub pre_pow: String,
285    /// Variable name prefix (e.g. `"p"`).
286    pub var: String,
287    /// Before the variable index (e.g. `"["`).
288    pub pre_var: String,
289    /// After the variable index (e.g. `"]"`).
290    pub post_var: String,
291    /// Power separator (e.g. `","`).
292    pub pow: String,
293    /// Power suffix (e.g. `")"`).
294    pub post_pow: String,
295    /// Line break inserted between wrapped lines.
296    pub linebreak: String,
297    /// Offset added to the 0-based variable position.
298    pub first_var: i64,
299    /// Offset added to each exponent.
300    pub first_pow: i64,
301    /// Monomials per line before wrapping.
302    pub monperline: u32,
303    /// Skip the power syntax for exponent 1.
304    pub shorten: bool,
305}
306
307impl SimpleFormat {
308    // A 14-field constructor mirroring the C++ aggregate initializer.
309    #[allow(clippy::too_many_arguments)]
310    fn new(
311        pos: &str,
312        neg: &str,
313        mul: &str,
314        pre_pow: &str,
315        var: &str,
316        pre_var: &str,
317        post_var: &str,
318        pow: &str,
319        post_pow: &str,
320        linebreak: &str,
321        first_var: i64,
322        first_pow: i64,
323        monperline: u32,
324        shorten: bool,
325    ) -> SimpleFormat {
326        SimpleFormat {
327            pos: pos.to_string(),
328            neg: neg.to_string(),
329            mul: mul.to_string(),
330            pre_pow: pre_pow.to_string(),
331            var: var.to_string(),
332            pre_var: pre_var.to_string(),
333            post_var: post_var.to_string(),
334            pow: pow.to_string(),
335            post_pow: post_pow.to_string(),
336            linebreak: linebreak.to_string(),
337            first_var,
338            first_pow,
339            monperline,
340            shorten,
341        }
342    }
343
344    /// C array-output preset (`p[0]*p[1]`).
345    pub fn c() -> SimpleFormat {
346        Self::new(
347            "+", "-", "*", "", "p", "[", "", "][", "]", " \\\n\t", 0, -1, 20, false,
348        )
349    }
350
351    /// C preset using `pow(...)`.
352    pub fn c_pow() -> SimpleFormat {
353        Self::new(
354            "+", "-", "*", "pow(", "x", "[", "]", ",", ")", " \\\n\t", 0, 0, 20, true,
355        )
356    }
357
358    /// Fortran preset (`p(1)*p(2)`).
359    pub fn fortran() -> SimpleFormat {
360        Self::new(
361            "+",
362            "-",
363            "*",
364            "",
365            "p",
366            "(",
367            "",
368            ",",
369            ")",
370            " &\n     &",
371            1,
372            0,
373            20,
374            false,
375        )
376    }
377
378    /// Fortran preset using `**( )`.
379    pub fn fortran_pow() -> SimpleFormat {
380        Self::new(
381            "+",
382            "-",
383            "*",
384            "",
385            "x",
386            "(",
387            ")",
388            "**(",
389            ")",
390            " &\n     &",
391            1,
392            0,
393            20,
394            true,
395        )
396    }
397
398    /// MATLAB preset (`p(1).*p(2)`).
399    pub fn matlab() -> SimpleFormat {
400        Self::new(
401            "+", "-", ".*", "", "p", "(", "", ",", ")", " ...\n\t", 1, 0, 20, false,
402        )
403    }
404
405    /// MATLAB preset using `.( )`.
406    pub fn matlab_pow() -> SimpleFormat {
407        Self::new(
408            "+", "-", ".*", "", "x", "(", ")", ".^(", ")", " ...\n\t", 1, 0, 20, true,
409        )
410    }
411
412    /// LaTeX preset (`x_{1} \cdot x_{2}^{3}`).
413    pub fn latex() -> SimpleFormat {
414        Self::new(
415            " +", " -", " \\cdot ", "", "x", "_{", "}", "^{", "}", " \n\t", 1, 0, 20, true,
416        )
417    }
418}
419
420/// Format a coefficient like a C++ ostream with `precision(16)` (16
421/// significant digits, trailing zeros trimmed, scientific notation outside
422/// `1e-4..1e16`).
423fn fmt_g16(c: f64) -> String {
424    if c == 0.0 {
425        return "0".to_string();
426    }
427    let sci = format!("{:.15e}", c);
428    let (mantissa, exponent) = sci.split_once('e').expect("scientific notation");
429    let exp: i32 = exponent.parse().unwrap_or(0);
430    if !(-4..16).contains(&exp) {
431        let mut m = mantissa
432            .trim_end_matches('0')
433            .trim_end_matches('.')
434            .to_string();
435        if m.is_empty() {
436            m = "0".to_string();
437        }
438        format!("{m}e{}{:02}", if exp < 0 { '-' } else { '+' }, exp.abs())
439    } else {
440        let decimals = (15 - exp).max(0) as usize;
441        let mut s = format!("{c:.decimals$}");
442        if s.contains('.') {
443            s = s.trim_end_matches('0').trim_end_matches('.').to_string();
444        }
445        s
446    }
447}
448
449/// Format a DA as source code in the given [`SimpleFormat`]
450/// (`DASimpleFormatter::format`).
451pub fn format_da(da: &Da, sf: &SimpleFormat) -> String {
452    let monomials = da.iter_monomials();
453    let all: Vec<_> = monomials.collect();
454    let mut res = String::new();
455    for (i, m) in all.iter().enumerate() {
456        if m.c < 0.0 {
457            res.push_str(&sf.neg);
458            res.push_str(&fmt_g16(-m.c));
459        } else {
460            res.push_str(&sf.pos);
461            res.push_str(&fmt_g16(m.c));
462        }
463        for (j, &e) in m.jj.iter().enumerate() {
464            if e == 0 {
465                continue;
466            } else if sf.shorten && e == 1 {
467                res.push_str(&sf.mul);
468                res.push_str(&sf.var);
469                res.push_str(&sf.pre_var);
470                res.push_str(&(j as i64 + sf.first_var).to_string());
471                res.push_str(&sf.post_var);
472            } else {
473                res.push_str(&sf.mul);
474                res.push_str(&sf.pre_pow);
475                res.push_str(&sf.var);
476                res.push_str(&sf.pre_var);
477                res.push_str(&(j as i64 + sf.first_var).to_string());
478                res.push_str(&sf.post_var);
479                res.push_str(&sf.pow);
480                res.push_str(&(e as i64 + sf.first_pow).to_string());
481                res.push_str(&sf.post_pow);
482            }
483        }
484        if (i + 1) % sf.monperline as usize == 0 && i + 1 < all.len() {
485            res.push_str(&sf.linebreak);
486        }
487    }
488    res
489}
490
491/// Format a vector of DAs, one per line (`DASimpleFormatter::format`).
492pub fn format_das(das: &[Da], sf: &SimpleFormat) -> String {
493    das.iter().map(|da| format_da(da, sf) + "\n").collect()
494}
495
496#[cfg(test)]
497mod tests {
498    use super::*;
499    use crate::test_support::CONTEXT_LOCK;
500
501    #[test]
502    fn display_matches_c_format() {
503        let _g = CONTEXT_LOCK.lock();
504        crate::context::init(3, 2).unwrap();
505        let x = Da::variable(1);
506        let y = Da::variable(2);
507        let f = 1.0 + 2.0 * x.clone() - 0.5 * y.clone() * y.clone();
508        let s = f.to_string();
509        let lines: Vec<&str> = s.lines().collect();
510        assert_eq!(lines[0], "     I  COEFFICIENT              ORDER EXPONENTS");
511        assert_eq!(lines[1], "     1    1.0000000000000000e+00   0   0  0");
512        assert_eq!(lines[2], "     2    2.0000000000000000e+00   1   1  0");
513        assert_eq!(lines[3], "     3   -5.0000000000000000e-01   2   0  2");
514        assert_eq!(lines[4], "------------------------------------------------");
515        assert_eq!(lines.len(), 5);
516
517        // Zero DA prints the special line plus footer.
518        let z = Da::new().to_string();
519        assert_eq!(
520            z,
521            "        ALL COEFFICIENTS ZERO\n------------------------------------------------\n"
522        );
523    }
524    #[test]
525    fn from_str_roundtrip() {
526        let _g = CONTEXT_LOCK.lock();
527        crate::context::init(4, 2).unwrap();
528        let x = Da::variable(1);
529        let y = Da::variable(2);
530        let f = 0.75 - 1.25 * x.clone() + 3.5 * (x.clone() * y.clone()) - 0.125 * y.clone();
531        let parsed: Da = f.to_string().parse().expect("parses");
532        assert_eq!(parsed.size(), f.size());
533        for m in f.iter_monomials() {
534            assert!(
535                (parsed.get_coefficient(&m.jj) - m.c).abs() == 0.0,
536                "{:?}: {} vs {}",
537                m.jj,
538                parsed.get_coefficient(&m.jj),
539                m.c
540            );
541        }
542
543        let z: Da = Da::new().to_string().parse().unwrap();
544        assert_eq!(z.size(), 0);
545
546        assert!("garbage".parse::<Da>().is_err());
547        assert!("".parse::<Da>().is_err());
548    }
549
550    #[test]
551    fn blob_roundtrip_and_layout() {
552        let _g = CONTEXT_LOCK.lock();
553        crate::context::init(5, 3).unwrap();
554        let x = Da::variable(1);
555        let y = Da::variable(2);
556        let z = Da::variable(3);
557        let f = 1.5 + 0.5 * (x.clone() * z.clone()) - 2.0 * y.clone();
558
559        let blob = f.to_blob();
560        // header (20) + one slot per term; empty DAs still reserve a slot
561        assert_eq!(blob.len(), 20 + 16 * f.size().max(1));
562        assert_eq!(
563            u32::from_le_bytes(blob[0..4].try_into().unwrap()),
564            0x1E304144
565        );
566        assert_eq!(u32::from_le_bytes(blob[4..8].try_into().unwrap()), 5);
567        assert_eq!(u32::from_le_bytes(blob[8..12].try_into().unwrap()), 2);
568        assert_eq!(u32::from_le_bytes(blob[12..16].try_into().unwrap()), 1);
569        assert_eq!(
570            u32::from_le_bytes(blob[16..20].try_into().unwrap()) as usize,
571            f.size()
572        );
573
574        let back = Da::from_blob(&blob).expect("imports");
575        for m in f.iter_monomials() {
576            assert!((back.get_coefficient(&m.jj) - m.c).abs() == 0.0);
577        }
578
579        // empty DA roundtrip
580        let zb = Da::new().to_blob();
581        assert_eq!(zb.len(), 36);
582        assert_eq!(Da::from_blob(&zb).unwrap().size(), 0);
583
584        // error paths
585        assert!(Da::from_blob(&[0u8; 36]).is_err());
586        assert!(Da::from_blob(&[0u8; 10]).is_err());
587
588        // truncation on import: order-4 term dropped at nomax=3 context
589        crate::context::init(3, 3).unwrap();
590        let g = Da::monomial(&[2, 1, 0], 2.0) + Da::monomial(&[1, 1, 1], 1.0);
591        let gblob = g.to_blob();
592        let gback = Da::from_blob(&gblob).unwrap();
593        assert_eq!(gback.get_coefficient(&[1, 1, 1]), 1.0);
594        // extra variable: blob made with 3 vars, read with fewer
595    }
596
597    #[test]
598    fn formatter_presets() {
599        let _g = CONTEXT_LOCK.lock();
600        crate::context::init(3, 2).unwrap();
601        let x = Da::variable(1);
602        let y = Da::variable(2);
603        let f =
604            1.5 + 2.0 * x.clone() - 0.5 * (x.clone() * y.clone()) + 0.25 * (x.clone() * x.clone());
605
606        let c = format_da(&f, &SimpleFormat::c());
607        assert_eq!(c, "+1.5+2*p[0][0]-0.5*p[0][0]*p[1][0]+0.25*p[0][1]");
608
609        let cp = format_da(&f, &SimpleFormat::c_pow());
610        assert_eq!(cp, "+1.5+2*x[0]-0.5*x[0]*x[1]+0.25*pow(x[0],2)");
611
612        let fo = format_da(&f, &SimpleFormat::fortran());
613        assert_eq!(fo, "+1.5+2*p(1,1)-0.5*p(1,1)*p(2,1)+0.25*p(1,2)");
614
615        let fp = format_da(&f, &SimpleFormat::fortran_pow());
616        assert_eq!(fp, "+1.5+2*x(1)-0.5*x(1)*x(2)+0.25*x(1)**(2)");
617
618        let ml = format_da(&f, &SimpleFormat::matlab());
619        assert_eq!(ml, "+1.5+2.*p(1,1)-0.5.*p(1,1).*p(2,1)+0.25.*p(1,2)");
620
621        let lx = format_da(&f, &SimpleFormat::latex());
622        assert_eq!(
623            lx,
624            " +1.5 +2 \\cdot x_{1} -0.5 \\cdot x_{1} \\cdot x_{2} +0.25 \\cdot x_{1}^{2}"
625        );
626        let multi = format_das(&[f.clone(), f.clone()], &SimpleFormat::c());
627        assert_eq!(multi.matches('\n').count(), 2);
628    }
629}
dace_rs/io.rs

dace_rs/
io.rs
