Skip to content
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension


Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
2 changes: 1 addition & 1 deletion Cargo.lock

Some generated files are not rendered by default. Learn more about how customized files appear on GitHub.

2 changes: 1 addition & 1 deletion rosy/Cargo.toml
Original file line number Diff line number Diff line change
@@ -1,6 +1,6 @@
[package]
name = "rosy"
version = "1.3.1"
version = "1.4.0"
edition = "2024"

[lib]
Expand Down
7 changes: 3 additions & 4 deletions rosy/src/program/statements/da/daprv/rosy_output.txt
Original file line number Diff line number Diff line change
@@ -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
------------------------------------------------------------------------------
16 changes: 11 additions & 5 deletions rosy/src/rosy_lib/core/da_ops.rs
Original file line number Diff line number Diff line change
Expand Up @@ -203,6 +203,10 @@ pub fn rosy_dafset(template_da: Vec<DA>) -> Result<()> {
/// is used as a universal mask).
pub fn rosy_dafilt(input: &Vec<DA>, result: &mut Vec<DA>) -> Result<()> {
let filter = get_filter_da()?;
let epsilon = {
let rt = get_runtime().context("DAFILT requires DA initialized")?;
rt.config.epsilon
};

match filter {
None => {
Expand All @@ -217,7 +221,7 @@ pub fn rosy_dafilt(input: &Vec<DA>, result: &mut Vec<DA>) -> 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);
}
Expand All @@ -234,9 +238,9 @@ pub fn rosy_dafilt(input: &Vec<DA>, result: &mut Vec<DA>) -> 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<DA>, 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);
Expand All @@ -246,8 +250,10 @@ pub fn rosy_daran(da: &mut Vec<DA>, 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);
}
}
}
}
Expand Down
67 changes: 32 additions & 35 deletions rosy/src/rosy_lib/core/dapew.rs
Original file line number Diff line number Diff line change
Expand Up @@ -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;

// ============================================================================
Expand Down Expand Up @@ -50,49 +50,42 @@ fn decode_transport_id(id: u64) -> [u8; MAX_VARS] {
pub fn rosy_darea(unit: u64, da: &mut Vec<DA>, 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<String> = 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);
Expand Down Expand Up @@ -122,13 +115,19 @@ pub fn rosy_dapew(unit: u64, da: &Vec<DA>, 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()
Expand All @@ -142,7 +141,7 @@ pub fn rosy_dapew(unit: u64, da: &Vec<DA>, 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)
Expand All @@ -151,24 +150,22 @@ pub fn rosy_dapew(unit: u64, da: &Vec<DA>, var_i: usize, order_n: u32) -> Result

for (idx, (monomial, coeff)) in terms.iter().enumerate() {
let order = monomial.total_order;
let exp_parts: Vec<String> = (0..MAX_VARS)
let exp_parts: Vec<String> = (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,
exp_str
));
}

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 {
Expand Down
Loading