diff --git a/Cargo.lock b/Cargo.lock index 5b16fb34..9a684dcc 100644 --- a/Cargo.lock +++ b/Cargo.lock @@ -1056,7 +1056,7 @@ dependencies = [ [[package]] name = "rosy" -version = "1.3.1" +version = "1.4.0" dependencies = [ "anyhow", "clap", diff --git a/rosy/Cargo.toml b/rosy/Cargo.toml index 2adde353..87c2d009 100644 --- a/rosy/Cargo.toml +++ b/rosy/Cargo.toml @@ -1,6 +1,6 @@ [package] name = "rosy" -version = "1.3.1" +version = "1.4.0" edition = "2024" [lib] diff --git a/rosy/src/program/statements/da/daprv/rosy_output.txt b/rosy/src/program/statements/da/daprv/rosy_output.txt index 3e182af2..fde4d03c 100644 --- a/rosy/src/program/statements/da/daprv/rosy_output.txt +++ b/rosy/src/program/statements/da/daprv/rosy_output.txt @@ -1,4 +1,3 @@ - I COEFFICIENT 1 ORDER EXPONENTS - 1 2.000000000000000 0 0 0 - 2 1.000000000000000 1 1 0 ------------------------------------------------------- + 2.000000000000 00 + 1.000000000000 10 + ------------------------------------------------------------------------------ diff --git a/rosy/src/rosy_lib/core/da_ops.rs b/rosy/src/rosy_lib/core/da_ops.rs index 59863da4..f0bd75c2 100644 --- a/rosy/src/rosy_lib/core/da_ops.rs +++ b/rosy/src/rosy_lib/core/da_ops.rs @@ -203,6 +203,10 @@ pub fn rosy_dafset(template_da: Vec) -> Result<()> { /// is used as a universal mask). pub fn rosy_dafilt(input: &Vec, result: &mut Vec) -> Result<()> { let filter = get_filter_da()?; + let epsilon = { + let rt = get_runtime().context("DAFILT requires DA initialized")?; + rt.config.epsilon + }; match filter { None => { @@ -217,7 +221,7 @@ pub fn rosy_dafilt(input: &Vec, result: &mut Vec) -> Result<()> { // Keep only monomials that are nonzero in the mask let mut new_da = DA::zero(); for &k in &src.nonzero { - if mask.coeffs[k as usize].abs() > 0.0 { + if mask.coeffs[k as usize].abs() > epsilon { new_da.coeffs[k as usize] = src.coeffs[k as usize]; new_da.nonzero.push(k); } @@ -234,9 +238,9 @@ pub fn rosy_dafilt(input: &Vec, result: &mut Vec) -> Result<()> { /// `sparsity` is the fraction of monomials that will be set nonzero /// (0.0 = all zero, 1.0 = all filled). Uses the global Rosy RNG. pub fn rosy_daran(da: &mut Vec, sparsity: f64) -> Result<()> { - let num_monomials = { + let (num_monomials, epsilon) = { let rt = get_runtime().context("DARAN requires DA to be initialized (call DAINI first)")?; - rt.num_monomials + (rt.num_monomials, rt.config.epsilon) }; let sparsity = sparsity.clamp(0.0, 1.0); @@ -246,8 +250,10 @@ pub fn rosy_daran(da: &mut Vec, sparsity: f64) -> Result<()> { for k in 0..num_monomials { if crate::rosy_lib::core::rng::rng_f64() < sparsity { let val = crate::rosy_lib::core::rng::rng_f64_symmetric(); - da_el.coeffs[k] = val; - da_el.nonzero.push(k as u32); + if val.abs() > epsilon { + da_el.coeffs[k] = val; + da_el.nonzero.push(k as u32); + } } } } diff --git a/rosy/src/rosy_lib/core/dapew.rs b/rosy/src/rosy_lib/core/dapew.rs index 99ff3577..e372c3cb 100644 --- a/rosy/src/rosy_lib/core/dapew.rs +++ b/rosy/src/rosy_lib/core/dapew.rs @@ -10,7 +10,7 @@ use anyhow::{Context, Result}; use crate::rosy_lib::core::display::RosyDisplay; -use crate::rosy_lib::taylor::{DA, MAX_VARS}; +use crate::rosy_lib::taylor::{DA, MAX_VARS, get_runtime}; use crate::rosy_lib::taylor::Monomial; // ============================================================================ @@ -50,49 +50,42 @@ fn decode_transport_id(id: u64) -> [u8; MAX_VARS] { pub fn rosy_darea(unit: u64, da: &mut Vec, num_vars: usize) -> Result<()> { use crate::rosy_lib::core::file_io::rosy_read_from_unit; - // Skip the header line - let _header = rosy_read_from_unit(unit).context("Failed to read header line in DAREA")?; - - // Collect data lines until a separator (all dashes) - let mut lines: Vec = Vec::new(); - loop { - let line = rosy_read_from_unit(unit).context("Failed to read data line in DAREA")?; - if line.trim().chars().all(|c| c == '-') && !line.trim().is_empty() { - break; - } - lines.push(line); - } - - // Ensure the output array has at least one element + // Ensure the output array has at least one element and is zeroed while da.is_empty() { da.push(DA::zero()); } da[0] = DA::zero(); - for line in &lines { + let nv = num_vars.min(MAX_VARS); + + // Read COSY DAPRV format: coefficient lines until separator (all dashes) + loop { + let line = rosy_read_from_unit(unit).context("Failed to read data line in DAREA")?; let trimmed = line.trim(); + + // Separator line ends the block + if trimmed.chars().all(|c| c == '-') && !trimmed.is_empty() { + break; + } if trimmed.is_empty() { continue; } - // DAPRV row layout (1-component): idx coeff order exp1 exp2 ... + // COSY DAPRV format: `{coeff:17.12} {concatenated_exps}` let tokens: Vec<&str> = trimmed.split_whitespace().collect(); - if tokens.len() < 3 { + if tokens.len() < 2 { continue; } - let coeff: f64 = tokens[1].parse().unwrap_or(0.0); + let coeff: f64 = tokens[0].parse().unwrap_or(0.0); if coeff.abs() <= 1e-15 { continue; } - // Exponents start at token index 3 (after idx, coeff, order) - let exp_start = 3; + // Exponents are concatenated single digits, e.g. "10" = x1=1, x2=0 let mut exponents = [0u8; MAX_VARS]; - for i in 0..num_vars.min(MAX_VARS) { - if exp_start + i < tokens.len() { - exponents[i] = tokens[exp_start + i].parse().unwrap_or(0); - } + for (i, ch) in tokens[1].chars().enumerate().take(nv) { + exponents[i] = ch.to_digit(10).unwrap_or(0) as u8; } let monomial = Monomial::new(exponents); @@ -122,13 +115,19 @@ pub fn rosy_dapew(unit: u64, da: &Vec, var_i: usize, order_n: u32) -> Result return Ok(()); } + let num_vars = get_runtime() + .context("DAPEW requires DA initialized")? + .config.num_vars; + let da0 = &da[0]; let mut output = String::new(); + // COSY DAPEW header format output.push_str(&format!( - " I COEFFICIENT ORDER EXPONENTS (DAPEW var={} order={})\n", - var_i, order_n + " ORDER{:>4} IN COLUMN{:>5}\n", + order_n, var_i )); + output.push_str(" I COEFFICIENT ORDER EXPONENTS\n"); let mut terms: Vec<(Monomial, f64)> = da0 .coeffs_iter() @@ -142,7 +141,7 @@ pub fn rosy_dapew(unit: u64, da: &Vec, var_i: usize, order_n: u32) -> Result }) .collect(); - // Sort by total order for deterministic output + // Sort by total order ascending then lexicographic terms.sort_by(|(a, _), (b, _)| { a.total_order .cmp(&b.total_order) @@ -151,12 +150,12 @@ pub fn rosy_dapew(unit: u64, da: &Vec, var_i: usize, order_n: u32) -> Result for (idx, (monomial, coeff)) in terms.iter().enumerate() { let order = monomial.total_order; - let exp_parts: Vec = (0..MAX_VARS) + let exp_parts: Vec = (0..num_vars.min(MAX_VARS)) .map(|i| format!("{:>2}", monomial.exponents[i])) .collect(); let exp_str = exp_parts.join(" "); output.push_str(&format!( - "{:>3} {} {:>5} {}\n", + " {:>4} {} {:>5} {}\n", idx + 1, coeff.rosy_display(), order, @@ -164,11 +163,9 @@ pub fn rosy_dapew(unit: u64, da: &Vec, var_i: usize, order_n: u32) -> Result )); } - if terms.is_empty() { - output.push_str(" (no matching terms)\n"); - } - - output.push_str(&"-".repeat(50)); + // COSY separator: 5 spaces + 35 dashes + output.push_str(" "); + output.push_str(&"-".repeat(35)); output.push('\n'); if unit == 6 { diff --git a/rosy/src/rosy_lib/core/daprv.rs b/rosy/src/rosy_lib/core/daprv.rs index a528492b..e030000a 100644 --- a/rosy/src/rosy_lib/core/daprv.rs +++ b/rosy/src/rosy_lib/core/daprv.rs @@ -5,10 +5,8 @@ use anyhow::{Result, Context, bail}; -use crate::rosy_lib::taylor::{DA, get_config}; -use crate::rosy_lib::taylor::da::DACoefficient; +use crate::rosy_lib::taylor::{DA, get_config, get_runtime}; use crate::rosy_lib::taylor::Monomial; -use crate::rosy_lib::core::display::RosyDisplay; /// Write an array of DA vectors in COSY INFINITY DAPRV format. /// @@ -40,16 +38,25 @@ pub fn rosy_daprv( Ok(()) } -/// Format DAPRV output in COSY-compatible format. +/// Format DAPRV output in COSY INFINITY-compatible format. +/// +/// COSY format (per component block): +/// - No header line +/// - Each non-zero term: `{coeff:17.12} {exponents_concatenated}\n` +/// - Separator: ` ` + 78 dashes + `\n` fn format_daprv( array: &Vec, num_components: usize, _max_vars: usize, current_vars: usize, ) -> Result { + let epsilon = get_runtime() + .context("DAPRV requires DA to be initialized (call OV first)")? + .config.epsilon; + let mut output = String::new(); - // Collect all unique monomials from all components + // Collect all unique monomials across all components let mut all_monomials: Vec = Vec::new(); for i in 0..num_components.min(array.len()) { for (m, _) in array[i].coeffs_iter() { @@ -59,7 +66,7 @@ fn format_daprv( } } - // Sort monomials: first by total order, then by exponents + // Sort by total order ascending, then reverse-lexicographic on exponents all_monomials.sort_by(|m1, m2| { m1.total_order.cmp(&m2.total_order) .then_with(|| { @@ -73,40 +80,23 @@ fn format_daprv( }) }); - // If no monomials, print a zero entry - if all_monomials.is_empty() { - all_monomials.push(Monomial::constant()); - } - - // Print header - output.push_str(&format!( - " I COEFFICIENT ")); - for comp in 1..=num_components.min(array.len()) { - output.push_str(&format!(" {:>2} ", comp)); - } - output.push_str("ORDER EXPONENTS\n"); - - // Print each monomial row - for (idx, monomial) in all_monomials.iter().enumerate() { - let order = monomial.total_order; - let exp_str = build_exp_str(&monomial.exponents, current_vars); - - output.push_str(&format!("{:>3} ", idx + 1)); - - // Print coefficient for each component - for i in 0..num_components.min(array.len()) { - let coeff = array[i].get_coeff(monomial); - output.push_str(&format!("{} ", coeff.rosy_display())); + // One block per component (single-column COSY format) + let nv = current_vars.min(6); + for comp_idx in 0..num_components.min(array.len()) { + for monomial in &all_monomials { + let coeff = array[comp_idx].get_coeff(monomial); + if coeff.abs() <= epsilon { + continue; + } + let exp_str: String = monomial.exponents[..nv] + .iter() + .map(|&e| char::from_digit(e as u32, 10).unwrap_or('?')) + .collect(); + output.push_str(&format!("{:17.12} {}\n", coeff, exp_str)); } - - output.push_str(&format!("{:>5} {}\n", order, exp_str)); + output.push_str(&format!(" {}\n", "-".repeat(78))); } - // Print separator - let sep_len = 30 + num_components.min(array.len()) * 24; - output.push_str(&"-".repeat(sep_len.min(132))); - output.push('\n'); - Ok(output) } @@ -125,85 +115,46 @@ pub fn rosy_darev( current_vars: usize, unit: u64, ) -> Result<()> { - // Read lines from the file - let mut lines = Vec::new(); - - // Read the header line - let _header = crate::rosy_lib::core::file_io::rosy_read_from_unit(unit) - .context("Failed to read header line in DAREV")?; - - // Read coefficient lines until we hit the separator - loop { - let line = crate::rosy_lib::core::file_io::rosy_read_from_unit(unit) - .context("Failed to read line in DAREV")?; - - // Check if this is a separator line (all dashes) - if line.trim().chars().all(|c| c == '-') && !line.trim().is_empty() { - break; - } - - lines.push(line); - } - - // Ensure array is big enough + // Ensure array is big enough and zeroed while array.len() < num_components { array.push(DA::zero()); } - - // Zero out the components we're reading into for i in 0..num_components.min(array.len()) { array[i] = DA::zero(); } - // Parse each line - for line in &lines { - let trimmed = line.trim(); - if trimmed.is_empty() { - continue; - } + let nv = current_vars.min(6); - // Parse the line: index, coefficients, order, exponents - let tokens: Vec<&str> = trimmed.split_whitespace().collect(); - if tokens.len() < 2 + num_components { - continue; // Skip malformed lines - } + // Read one block per component; each block ends with a separator line (all dashes) + for comp_idx in 0..num_components.min(array.len()) { + loop { + let line = crate::rosy_lib::core::file_io::rosy_read_from_unit(unit) + .context("Failed to read line in DAREV")?; + let trimmed = line.trim(); - // First token is the index (1-based), skip it - // Next num_components tokens are coefficients - // Then order - // Then exponents - let mut coeffs = Vec::new(); - for i in 0..num_components { - if let Ok(coeff) = tokens[1 + i].parse::() { - coeffs.push(coeff); - } else { - coeffs.push(0.0); + // Separator line (all dashes) ends this component's block + if trimmed.chars().all(|c| c == '-') && !trimmed.is_empty() { + break; + } + if trimmed.is_empty() { + continue; } - } - // Order is after the coefficients - let order_idx = 1 + num_components; - if order_idx >= tokens.len() { - continue; - } - - // Exponents start after order - let exp_start = order_idx + 1; - let mut exponents = [0u8; 6]; - for i in 0..current_vars.min(6) { - if exp_start + i < tokens.len() { - if let Ok(exp) = tokens[exp_start + i].parse::() { - exponents[i] = exp; - } + let tokens: Vec<&str> = trimmed.split_whitespace().collect(); + if tokens.len() < 2 { + continue; } - } - let monomial = Monomial::new(exponents); + let coeff: f64 = tokens[0].parse().unwrap_or(0.0); + // Exponents are concatenated single digits per variable, e.g. "10" = x1=1, x2=0 + let mut exponents = [0u8; 6]; + for (i, ch) in tokens[1].chars().enumerate().take(nv) { + exponents[i] = ch.to_digit(10).unwrap_or(0) as u8; + } + let monomial = Monomial::new(exponents); - // Set coefficients for each component - for (i, &coeff) in coeffs.iter().enumerate() { - if i < array.len() && coeff.abs() > 1e-15 { - array[i].set_coeff(monomial, coeff); + if coeff.abs() > 1e-15 { + array[comp_idx].set_coeff(monomial, coeff); } } } @@ -302,19 +253,6 @@ pub fn rosy_datrn( Ok(()) } -/// Build exponent string for DAPRV display. -fn build_exp_str(exponents: &[u8], num_vars: usize) -> String { - let mut result = String::new(); - for i in 0..num_vars.min(exponents.len()) { - if i % 2 == 0 { - result.push_str(&format!("{:>2}", exponents[i])); - } else { - result.push_str(&format!("{:>2} ", exponents[i])); - } - } - result.trim_end().to_string() -} - /// DAPLU: Replace independent variable xi by constant C in a DA vector. /// /// For each term c·x₁^e₁·…·xᵢ^eᵢ·…·xₙ^eₙ, the result accumulates diff --git a/rosy/src/rosy_lib/core/display.rs b/rosy/src/rosy_lib/core/display.rs index 57b8f302..1ea3ef58 100644 --- a/rosy/src/rosy_lib/core/display.rs +++ b/rosy/src/rosy_lib/core/display.rs @@ -1,4 +1,5 @@ use crate::rosy_lib::{RE, ST, LO, CM, VE, DA, CD}; +use crate::rosy_lib::taylor::get_runtime; fn sci(x: f64) -> (f64, i32) { if x == 0.0 { @@ -83,10 +84,9 @@ fn display_re ( ) } } -fn build_exp_str ( - exps: &[u8], -) -> String { - exps.iter() +fn build_exp_str(exps: &[u8], num_vars: usize) -> String { + exps[..num_vars.min(exps.len())] + .iter() .enumerate() .fold(String::new(), |mut acc, (i, exp)| { if i % 2 == 0 { @@ -170,13 +170,12 @@ impl RosyDisplay for &DA { for (idx, (monomial, coeff)) in sorted.iter().enumerate() { let order = monomial.total_order; let exp_str = { - // For 6 exponents, should match: '1 0 1 0 0 0' let exps = &monomial.exponents; - - build_exp_str(exps) + let nv = get_runtime().map(|rt| rt.config.num_vars).unwrap_or(exps.len()); + build_exp_str(exps, nv) }; output.push_str(&format!( - "{} {} {} {}\n", + "{} {} {} {}\n", idx + 1, coeff.rosy_display(), format!("{:>3}", order), @@ -238,15 +237,14 @@ impl RosyDisplay for &CD { let imag_coeff = imag_part.get_coeff(monomial); let order = monomial.total_order; let exp_str = { - // For 6 exponents, should match: '1 0 1 0 0 0' let exps = &monomial.exponents; - - build_exp_str(exps) + let nv = get_runtime().map(|rt| rt.config.num_vars).unwrap_or(exps.len()); + build_exp_str(exps, nv) }; output.push_str(&format!( " {} {} {} {:>3} {}\n", - idx + 1, - real_coeff.rosy_display(), + idx + 1, + real_coeff.rosy_display(), imag_coeff.rosy_display(), order, exp_str.trim_end() diff --git a/rosy/src/rosy_lib/core/rkco.rs b/rosy/src/rosy_lib/core/rkco.rs index 9deb4ca7..f6162aa3 100644 --- a/rosy/src/rosy_lib/core/rkco.rs +++ b/rosy/src/rosy_lib/core/rkco.rs @@ -38,20 +38,21 @@ pub fn rosy_rkco() -> Result<(Vec, Vec, Vec, Vec, Vec)> 1.0, ]; - // --- b: 8th-order weights --- + // --- b: 8th-order weights (Hairer, Nørsett, Wanner "Solving ODEs I", dop853.f) --- + // Sum = 1.0 exactly; b[1..4] and b[12] = 0 (unused stages) let b: Vec = vec![ 5.42937341165687296e-2, 0.0, 0.0, 0.0, 0.0, - 4.45031289275240888e-1, - 1.89237478148923991e-1, - -2.72937341165687296e-2, - 3.05326994405566566e-2, - 1.79592280957798019e-2, - 2.49919795974755027e-3, - 0.0, // b[11] = 0 in DOP853 + 4.45031289275240888e+0, + 1.89151789931450038e+0, + -5.80120396001058478e+0, + 3.11164366957819894e-1, + -1.52160949662516079e-1, + 2.01365400804030348e-1, + 4.47106157277725905e-2, 0.0, ]; @@ -63,8 +64,8 @@ pub fn rosy_rkco() -> Result<(Vec, Vec, Vec, Vec, Vec)> 0.0, 0.0, 0.0, - -0.1225156446376204440e+1, - -0.4957589496572501915e+0, + -0.1225156446376204440e-1, + -0.4957589496572501915e-1, 0.1664377182454986536e+0, -0.3558496486701148929e+0, 0.9340847839611065608e+0, @@ -113,11 +114,11 @@ pub fn rosy_rkco() -> Result<(Vec, Vec, Vec, Vec, Vec)> -2.00087205822486249e-2, 1.57982909820588250e1, 2.57112430717927171e1, -4.05313840176771403e1, -1.37316482655824625e1, 2.13374040065074902e1, 2.93930402093266800e0, - // row 13 (12 elements) — same as b (FSAL property, last row = b weights) + // row 13 (12 elements) — FSAL: mirrors b weights exactly 5.42937341165687296e-2, 0.0, 0.0, 0.0, 0.0, - 4.45031289275240888e-1, 1.89237478148923991e-1, -2.72937341165687296e-2, - 3.05326994405566566e-2, 1.79592280957798019e-2, 2.49919795974755027e-3, - 0.0, + 4.45031289275240888e+0, 1.89151789931450038e+0, -5.80120396001058478e+0, + 3.11164366957819894e-1, -1.52160949662516079e-1, 2.01365400804030348e-1, + 4.47106157277725905e-2, ]; Ok((c, b, e, a1, a2)) diff --git a/rosy/src/rosy_lib/taylor/da.rs b/rosy/src/rosy_lib/taylor/da.rs index a94a3dc9..a0181712 100644 --- a/rosy/src/rosy_lib/taylor/da.rs +++ b/rosy/src/rosy_lib/taylor/da.rs @@ -337,9 +337,10 @@ impl DA { /// O(1) amortized. Used by Horner's method. pub fn add_constant_in_place(&mut self, value: T) { - let was_nz = self.coeffs[0].abs() > 1e-15; + let epsilon = get_runtime().map(|rt| rt.config.epsilon).unwrap_or(1e-15); + let was_nz = self.coeffs[0].abs() > epsilon; self.coeffs[0] = self.coeffs[0] + value; - let is_nz = self.coeffs[0].abs() > 1e-15; + let is_nz = self.coeffs[0].abs() > epsilon; if is_nz && !was_nz { self.nonzero.push(0); @@ -358,7 +359,8 @@ impl DA { /// this avoids the RwLock acquisition and HashMap lookup in `set_coeff`. pub fn make_prime(&self) -> Self { let mut prime = self.clone(); - if prime.coeffs[0].abs() > 1e-15 { + let epsilon = get_runtime().map(|rt| rt.config.epsilon).unwrap_or(1e-15); + if prime.coeffs[0].abs() > epsilon { prime.coeffs[0] = T::zero(); // Remove 0 from nonzero list if let Some(pos) = prime.nonzero.iter().position(|&i| i == 0) {