From 6b91e0c2853d80c04e8ce4f793ce3c455e698fab Mon Sep 17 00:00:00 2001 From: AlexanderMitrofanov Date: Tue, 29 Sep 2026 19:37:28 +0200 Subject: [PATCH 1/4] Add compressed binary accessibility caches for issue #245 --- ChangeLog | 21 ++++ README.md | 39 +++++++ configure.ac | 25 ++++- doc/Makefile.am | 2 + doc/benchmarks/accessibility.cpp | 57 ++++++++++ doc/benchmarks/accessibility.md | 45 ++++++++ src/IntaRNA/Accessibility.cpp | 11 ++ src/IntaRNA/Accessibility.h | 9 ++ src/IntaRNA/AccessibilityArchive.h | 131 ++++++++++++++++++++++ src/IntaRNA/AccessibilityFromStream.cpp | 24 +++++ src/IntaRNA/AccessibilityFromStream.h | 4 + src/IntaRNA/Makefile.am | 2 + src/IntaRNA/general.cpp | 4 +- src/IntaRNA/general.h | 4 +- src/bin/CommandLineParsing.cpp | 55 ++++++---- tests/AccessibilityBinary_test.cpp | 137 ++++++++++++++++++++++++ tests/Makefile.am | 3 +- tests/runAccessibilityBinary.sh | 56 ++++++++++ 18 files changed, 601 insertions(+), 28 deletions(-) create mode 100644 doc/benchmarks/accessibility.cpp create mode 100644 doc/benchmarks/accessibility.md create mode 100644 src/IntaRNA/AccessibilityArchive.h create mode 100644 tests/AccessibilityBinary_test.cpp create mode 100644 tests/runAccessibilityBinary.sh diff --git a/ChangeLog b/ChangeLog index 786d1b97..14189aab 100644 --- a/ChangeLog +++ b/ChangeLog @@ -14,6 +14,8 @@ ## Interface and handling +- compressed binary .agz accessibility caches for repeated screens (issue #245) + - too short sequences are skipped with warning if they are shorter than the required seed or helix length (were causing an exception and abort before) - default SHAPE method changed to "D", which is the default for ViennaRNA 2.7.* @@ -52,6 +54,25 @@ ################################################################################ ################################################################################ +260929 Alexander Mitrofanov + * IntaRNA/Accessibility, AccessibilityArchive, AccessibilityFromStream : + + versioned Boost binary serialization of exact ED matrices and sequence data + + stream rows with bounded scratch memory, retaining dangling-end intervals + + validate sequence, dimensions, energies, archive format and stream integrity + * bin/CommandLineParsing, IntaRNA/general : + + select compressed binary accessibility input/output via .agz filenames + * support all existing accessibility producers and both P/E input modes + * release accessibility streams when parsing or writing fails + * configure.ac, src/IntaRNA/Makefile.am : + + check and link Boost.Serialization; distribute the internal matrix archive view + * tests/AccessibilityBinary_test.cpp, tests/runAccessibilityBinary.sh, tests/Makefile.am : + + exact round trips, malformed archives, prediction reuse and filename dispatch + * doc/benchmarks/accessibility.cpp, doc/benchmarks/accessibility.md, doc/Makefile.am : + + reproducible compressed text versus binary I/O benchmark with exact ED checks + * README.md and CLI help : + + document binary caches, unchanged text formats and compatibility requirements + * resolves https://github.com/BackofenLab/IntaRNA/issues/245 + 260928 Alexander Mitrofanov * AGENTS.md : + repository-wide AI coding guidance based on the current code and build setup diff --git a/README.md b/README.md index c115f384..c2b0879a 100644 --- a/README.md +++ b/README.md @@ -187,6 +187,8 @@ dependencies: - libboost_program_options - libboost_filesystem - libboost_system + - libboost_iostreams + - libboost_serialization - [Vienna RNA package](http://www.tbi.univie.ac.at/RNA/) version >= 2.4.14 - `pkg-config` for detailed version checks of dependencies - if [cloning from github](#instgithub): GNU autotools (automake, autoconf, ..) @@ -2109,6 +2111,7 @@ formats | RNAplfold-styled ED values | `IntaRNA --out=*Acc:` | | ---- | --- | | .. with gzip-compression | `IntaRNA --out=*:*.gz` | +| IntaRNA binary accessibility (gzip-compressed) | `IntaRNA --out=*Acc:*.agz` or `--out=*Pu:*.agz` | The **RNAplfold** format is a table encoding of a banded upper triangular matrix with band width l. First row contains a header comment on the data starting with @@ -2132,6 +2135,42 @@ example for a sequence of length 5 with a maximal window length of 3. ``` +##### Binary accessibility caches (`.agz`) + +For repeated screens, use a filename ending in `.agz` to store a compressed +Boost binary archive of the accessibility matrix. Both `qAcc:`/`tAcc:` and +`qPu:`/`tPu:` write **exact internal ED values** to this format, including the +extra interval length used for dangling-end probabilities. This avoids decimal +formatting/parsing and probability-to-energy rounding. It works with every +accessibility implementation, including disabled accessibility and loaded data. + +```bash +# Add these options to otherwise identical IntaRNA calls: +IntaRNA [..] --out=qAcc:query.agz --out=tAcc:target.agz +IntaRNA [..] --qAcc=E --qAccFile=query.agz --tAcc=E --tAccFile=target.agz +``` + +The `.agz` extension (case insensitive) selects the binary reader with either +`--qAcc=E`/`--tAcc=E` or `--qAcc=P`/`--tAcc=P`; the stored ED values are used +unchanged in both cases. Sequence content and length must match. The requested +interaction length can be smaller than the stored limit; larger requests are +limited to the stored interaction length. A format version and gzip integrity +checks reject unsupported or damaged files. Existing plain text, `.gz` text, +`STDIN` and `STDOUT` behavior is unchanged; binary CLI I/O requires a filename. +Multi-sequence filenames use the same `-s#` suffix convention described below. + +Reuse caches with the same temperature, energy model, folding window, base-pair +span and folding constraints used to generate them. These settings are not +stored or checked, and changing the temperature does not rescale cached EDs. +Constraints already reflected in ED values are preserved; supplying new +accessibility constraints while reading a cache remains unsupported. + +The binary format uses native Boost.Serialization archives and requires a +compatible architecture and Boost archive version. Use RNAplfold-style text +for portable interchange. For C++ callers, `Accessibility::writeBinary()` and +`AccessibilityFromStream::IntaRNA_Binary` operate on decompressed archive streams; +`newOutputStream()`/`newInputStream()` supply gzip compression for `.agz` files. + ##### Use case examples for read/write accessibilities and unpaired probabilities If you have precomputed data, e.g. the file `plfold_lunp` with unpaired probabilities diff --git a/configure.ac b/configure.ac index ea32bd28..e13234cb 100644 --- a/configure.ac +++ b/configure.ac @@ -327,6 +327,29 @@ AS_IF([test "x$want_boost" != "xno" && test "x$FOUND_BOOST" = "x1"], [ AX_BOOST_PROGRAM_OPTIONS AX_BOOST_REGEX AX_BOOST_IOSTREAMS + AC_MSG_CHECKING([for Boost.Serialization binary archives]) + intarna_saved_CPPFLAGS=$CPPFLAGS + intarna_saved_LDFLAGS=$LDFLAGS + intarna_saved_LIBS=$LIBS + CPPFLAGS="$CPPFLAGS $BOOST_CPPFLAGS" + LDFLAGS="$LDFLAGS $BOOST_LDFLAGS" + LIBS="-lboost_serialization $LIBS" + AC_LINK_IFELSE([AC_LANG_PROGRAM( + [[#include +#include +#include ]], + [[std::stringstream stream; +boost::archive::binary_oarchive output(stream); +int value = 1; +output << value; +boost::archive::binary_iarchive input(stream); +input >> value;]])], + [AC_MSG_RESULT([yes])], + [AC_MSG_RESULT([no]) + AC_MSG_ERROR([Boost.Serialization headers and library are required for .agz accessibility files. Install libboost-serialization-dev or select --with-boost and --with-boost-libdir.])]) + CPPFLAGS=$intarna_saved_CPPFLAGS + LDFLAGS=$intarna_saved_LDFLAGS + LIBS=$intarna_saved_LIBS ]) AX_CHECK_ZLIB([], []) @@ -350,7 +373,7 @@ AS_IF([test $want_boost != "no" && test $FOUND_BOOST != 1], [ ], [ AM_CXXFLAGS="$BOOST_CPPFLAGS $AM_CXXFLAGS" AM_LDFLAGS="$BOOST_LDFLAGS $AM_LDFLAGS" - LIBS="$LIBS -lboost_regex -lboost_program_options -lboost_filesystem -lboost_system -lboost_iostreams -lz" + LIBS="$LIBS -lboost_regex -lboost_program_options -lboost_filesystem -lboost_system -lboost_iostreams -lboost_serialization -lz" ]) diff --git a/doc/Makefile.am b/doc/Makefile.am index 322e9d93..985e6fae 100644 --- a/doc/Makefile.am +++ b/doc/Makefile.am @@ -4,6 +4,8 @@ ################################################################ EXTRA_DIST = \ + benchmarks/accessibility.cpp \ + benchmarks/accessibility.md \ conda.txt \ doxygen.cfg \ latex-deps/adjcalc.sty \ diff --git a/doc/benchmarks/accessibility.cpp b/doc/benchmarks/accessibility.cpp new file mode 100644 index 00000000..c3ecfa34 --- /dev/null +++ b/doc/benchmarks/accessibility.cpp @@ -0,0 +1,57 @@ +// Compare compressed text and binary accessibility I/O after one ViennaRNA fold. +// Usage: accessibility-benchmark SEQUENCE_LENGTH EXISTING_OUTPUT_DIRECTORY +#include +#include +#include +#include +#include +#include +#include +#include + +INITIALIZE_EASYLOGGINGPP + +int main(int argc, char **argv) { + if (argc != 3) return 2; + el::Loggers::reconfigureAllLoggers(el::ConfigurationType::Enabled, "false"); + const size_t n = std::stoull(argv[1]); + std::string sequence; + sequence.reserve(n); + unsigned state = 245; + for (size_t i=0;i> 30]; + } + IntaRNA::RnaSequence rna("benchmark", sequence); + IntaRNA::VrnaHandler vrna(37, "Turner04", false, false); + IntaRNA::AccessibilityVrna source(rna, 100, nullptr, vrna, 150); + using Clock = std::chrono::steady_clock; + for (int run=0;run<3;++run) { + for (bool binary : {false, true}) { + const auto path = std::string(argv[2]) + (binary ? "/cache.agz" : "/cache.txt.gz"); + auto start=Clock::now(); + auto *out=IntaRNA::newOutputStream(path); + if (binary) source.writeBinary(*out); else source.writeRNAplfold_ED_text(*out); + IntaRNA::deleteOutputStream(out); + auto wrote=Clock::now(); + auto *in=IntaRNA::newInputStream(path); + IntaRNA::AccessibilityFromStream read(rna,100,nullptr,*in, + binary ? IntaRNA::AccessibilityFromStream::IntaRNA_Binary : IntaRNA::AccessibilityFromStream::ED_RNAplfold_Text, 1.0); + IntaRNA::deleteInputStream(in); + auto loaded=Clock::now(); + if (read.getMaxLength()!=source.getMaxLength()) return 3; + size_t differences = 0; + IntaRNA::E_type maxDifference = 0; + for(size_t i=0;i(wrote-start).count() << ',' + << std::chrono::duration(loaded-wrote).count() << ',' + << std::filesystem::file_size(path) << ',' << differences << ',' << maxDifference << std::endl; + } + } +} diff --git a/doc/benchmarks/accessibility.md b/doc/benchmarks/accessibility.md new file mode 100644 index 00000000..439a04ea --- /dev/null +++ b/doc/benchmarks/accessibility.md @@ -0,0 +1,45 @@ +# Accessibility I/O benchmark + +This benchmark isolates writing and reading the same accessibility matrix as +RNAplfold-style gzip-compressed ED text and as a binary `.agz` archive. It folds +one deterministic pseudorandom RNA with ViennaRNA (37 C, Turner04, folding window +150, interaction length 100), then alternates formats for three trials. Folding +and verification are outside the timed sections. Every loaded ED cell, including +the dangling-end band, is compared to the computed matrix. Binary input must +match exactly; text conversion differences are reported. + +Build and install IntaRNA, then compile against that installation: + +```bash +export PKG_CONFIG_PATH=/path/to/intarna-install/lib/pkgconfig +c++ -std=c++23 -O3 doc/benchmarks/accessibility.cpp \ + $(pkg-config --cflags --libs IntaRNA) -leasylogging \ + -o accessibility-benchmark +mkdir -p accessibility-benchmark-output +./accessibility-benchmark 100000 accessibility-benchmark-output +``` + +Use the same compiler, dependency include/library paths and runtime library +paths as the IntaRNA build. With dependencies outside system paths, additional +`-L` and `-Wl,-rpath` flags may be necessary. Output columns are +`trial,format,write_seconds,read_seconds,compressed_bytes,different_cells,max_ED_difference`. +ED differences are measured in internal units (hundredths of kcal/mol). Files are overwritten +between trials; the output directory must already exist. These timings include +compression and stream closing, but not an `fsync` to durable storage. Results +depend on sequence, bandwidth, compressor, filesystem and machine load. + +## Example measurement (2026-09-29) + +Linux x86-64, Ryzen 5 7530U, GCC 16.2.0 release build, Boost 1.85.0, +ViennaRNA 2.7.2, 100,000 bases, one folding/I/O thread. Medians of the three +alternating trials above (other build work was running on the machine): + +| Format | Write (s) | Read (s) | Compressed bytes | ED cells differing from source | Maximum ED difference | +| --- | ---: | ---: | ---: | ---: | ---: | +| ED text + gzip | 8.726 | 3.365 | 20,891,786 | 551,861 | 1 (0.01 kcal/mol) | +| Binary `.agz` | 3.582 | 0.349 | 16,620,982 | 0 | 0 | + +In this local run, binary output was about 2.4 times faster to write, 9.6 times +faster to read, and 20% smaller. Text differences arise in existing decimal +conversion; the benchmark does not modify that behavior. These figures describe +I/O of the same computed data, not an end-to-end genome-screen speedup. diff --git a/src/IntaRNA/Accessibility.cpp b/src/IntaRNA/Accessibility.cpp index 9c4d0c09..b5b5ce87 100644 --- a/src/IntaRNA/Accessibility.cpp +++ b/src/IntaRNA/Accessibility.cpp @@ -1,6 +1,8 @@ #include "IntaRNA/Accessibility.h" +#include "IntaRNA/AccessibilityArchive.h" +#include namespace IntaRNA { @@ -8,6 +10,15 @@ namespace IntaRNA { const E_type Accessibility::ED_UPPER_BOUND = (E_type) E_INF; +void +Accessibility::writeBinary( std::ostream & out ) const +{ + boost::archive::binary_oarchive archive(out); + archive << AccessibilityArchive(*this); + out.flush(); + if (!out) throw std::runtime_error("Accessibility::writeBinary: output failure"); +} + //////////////////////////////////////////////////////////////////// std::ostream& diff --git a/src/IntaRNA/Accessibility.h b/src/IntaRNA/Accessibility.h index 7da52c4f..652c7da5 100644 --- a/src/IntaRNA/Accessibility.h +++ b/src/IntaRNA/Accessibility.h @@ -108,6 +108,15 @@ class Accessibility { void writeRNAplfold_ED_text( std::ostream& out ) const; + /** + * Write a native Boost binary archive of exact ED values, including the + * extra interval length needed for dangling ends. Works for every subclass. + * Compression is supplied by the stream (e.g. newOutputStream("file.agz")). + * @param out binary output stream + * @throw std::exception on invalid ED values or output failure + */ + void writeBinary( std::ostream & out ) const; + /** * Prints the accessibility values to stream as upper triangular matrix * @param out the ostream to write to diff --git a/src/IntaRNA/AccessibilityArchive.h b/src/IntaRNA/AccessibilityArchive.h new file mode 100644 index 00000000..a51cbdb4 --- /dev/null +++ b/src/IntaRNA/AccessibilityArchive.h @@ -0,0 +1,131 @@ +#ifndef INTARNA_ACCESSIBILITYARCHIVE_H_ +#define INTARNA_ACCESSIBILITYARCHIVE_H_ + +#include "IntaRNA/Accessibility.h" + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +namespace IntaRNA { + +/** + * Serialization view of the logical accessibility matrix, independent of its + * physical storage. Version 1 stores exact internal ED integers, in rows of + * increasing start position and interval length, including maxLength+1 for + * dangling ends. Only valid cells are stored; scratch memory is O(maxLength). + * + * The input sequence and requested band bound allocations on load. Native + * Boost binary archives require a compatible architecture and Boost version. + * This is an internal file-format helper, not an installed public interface. + */ +class AccessibilityArchive { +public: + typedef boost::numeric::ublas::banded_matrix EdMatrix; + + /** Create a read-only serialization view of any accessibility implementation. */ + explicit AccessibilityArchive( const Accessibility & source ); + + /** Create a loading view; target is resized to the requested available band. */ + AccessibilityArchive( const RnaSequence & sequence, size_t maxLength, EdMatrix & target ); + + /** Return the interaction length retained after loading. */ + size_t getMaxLength() const; + +private: + friend class boost::serialization::access; + const RnaSequence & sequence; + size_t maxLength; + const Accessibility * source; + EdMatrix * target; + + template void save( Archive & archive, unsigned int version ) const; + template void load( Archive & archive, unsigned int version ); + BOOST_SERIALIZATION_SPLIT_MEMBER() +}; + +inline +AccessibilityArchive::AccessibilityArchive( const Accessibility & source ) + : sequence(source.getSequence()), maxLength(source.getMaxLength()), source(&source), target(nullptr) +{} + +inline +AccessibilityArchive::AccessibilityArchive( const RnaSequence & sequence, size_t maxLength, EdMatrix & target ) + : sequence(sequence), maxLength(maxLength), source(nullptr), target(&target) +{} + +inline size_t +AccessibilityArchive::getMaxLength() const +{ + return maxLength; +} + +template void +AccessibilityArchive::save( Archive & archive, unsigned int ) const +{ + const std::uint32_t magic = 0x49414343; // IACC: IntaRNA accessibility + const std::uint32_t formatVersion = 1; + const std::uint64_t length = sequence.size(), band = maxLength; + const E_type infinity = Accessibility::ED_UPPER_BOUND; + archive & magic & formatVersion & length & band & infinity; + archive & boost::serialization::make_array(sequence.asString().data(), sequence.size()); + // Avoid maxLength+1 overflow when the entire sequence is covered. + const size_t width = maxLength < sequence.size() ? maxLength+1 : sequence.size(); + std::vector row(width); + for (size_t i = 0; i < sequence.size(); ++i) { + const size_t count = std::min(width, sequence.size()-i); + for (size_t k = 0; k < count; ++k) { + row[k] = source->getED(i, i+k); + if (row[k] < 0 || row[k] > infinity) + throw std::runtime_error("Accessibility archive: invalid ED value"); + } + archive & boost::serialization::make_array(row.data(), count); + } +} + +template void +AccessibilityArchive::load( Archive & archive, unsigned int ) +{ + std::uint32_t magic, formatVersion; + std::uint64_t length, band; + E_type infinity; + archive & magic & formatVersion & length & band & infinity; + if (magic != 0x49414343 || formatVersion != 1) + throw std::runtime_error("Accessibility archive: invalid signature or unsupported format version"); + if (length == 0 || length != sequence.size() || band > length) + throw std::runtime_error("Accessibility archive: sequence length or matrix band mismatch"); + if (infinity != Accessibility::ED_UPPER_BOUND) + throw std::runtime_error("Accessibility archive: incompatible energy representation"); + std::string storedSequence(sequence.size(), '\0'); + archive & boost::serialization::make_array(storedSequence.data(), storedSequence.size()); + if (storedSequence != sequence.asString()) + throw std::runtime_error("Accessibility archive: sequence mismatch"); + + const size_t storedMaxLength = static_cast(band); + maxLength = std::min(maxLength, storedMaxLength); + const size_t storedWidth = storedMaxLength < sequence.size() ? storedMaxLength+1 : sequence.size(); + const size_t width = maxLength < sequence.size() ? maxLength+1 : sequence.size(); + if (sequence.size() > std::numeric_limits::max() / width / sizeof(E_type)) + throw std::runtime_error("Accessibility archive: matrix size overflow"); + target->resize(sequence.size(), sequence.size(), 0, width-1, false); + std::vector row(storedWidth); + for (size_t i = 0; i < sequence.size(); ++i) { + const size_t count = std::min(storedWidth, sequence.size()-i); + archive & boost::serialization::make_array(row.data(), count); + for (size_t k = 0; k < count; ++k) { + if (row[k] < 0 || row[k] > infinity) + throw std::runtime_error("Accessibility archive: invalid ED value"); + if (k < width) (*target)(i, i+k) = row[k]; + } + } +} + +} // namespace IntaRNA +#endif diff --git a/src/IntaRNA/AccessibilityFromStream.cpp b/src/IntaRNA/AccessibilityFromStream.cpp index d486c189..acc747b0 100644 --- a/src/IntaRNA/AccessibilityFromStream.cpp +++ b/src/IntaRNA/AccessibilityFromStream.cpp @@ -1,5 +1,7 @@ #include "IntaRNA/AccessibilityFromStream.h" +#include "IntaRNA/AccessibilityArchive.h" +#include #include #include @@ -31,6 +33,10 @@ AccessibilityFromStream( } switch( inStreamType ) { + case IntaRNA_Binary : + parseBinary( inStream ); + break; + case Pu_RNAplfold_Text : parsePu_RNAplfold_text( inStream, RT ); break; @@ -52,6 +58,24 @@ AccessibilityFromStream:: ///////////////////////////////////////////////////////////////////////// +void +AccessibilityFromStream::parseBinary( std::istream & inStream ) +{ + try { + boost::archive::binary_iarchive archive(inStream); + AccessibilityArchive matrix(getSequence(), availMaxLength, edValues); + archive >> matrix; + // Read through EOF so gzip trailer/CRC and truncation errors surface. + if (inStream.get() != std::char_traits::eof() || inStream.bad()) + throw std::runtime_error("trailing data or damaged compressed stream"); + availMaxLength = matrix.getMaxLength(); + } catch (const std::exception & error) { + throw std::runtime_error(std::string("AccessibilityFromStream: invalid binary accessibility input: ") + error.what()); + } +} + +///////////////////////////////////////////////////////////////////////// + void AccessibilityFromStream:: parseRNAplfold_text( std::istream & inStream, const Z_type RT, const bool parseProbs ) diff --git a/src/IntaRNA/AccessibilityFromStream.h b/src/IntaRNA/AccessibilityFromStream.h index f71a7353..9d2fdb29 100644 --- a/src/IntaRNA/AccessibilityFromStream.h +++ b/src/IntaRNA/AccessibilityFromStream.h @@ -21,6 +21,7 @@ class AccessibilityFromStream: public Accessibility enum InStreamType { Pu_RNAplfold_Text //! Pu values in RNAplfold text format , ED_RNAplfold_Text //!< ED values in RNAplfold text Pu format + , IntaRNA_Binary //!< native Boost binary ED archive, already decompressed }; public: @@ -90,6 +91,9 @@ class AccessibilityFromStream: public Accessibility //! the ED values for the given sequence EdMatrix edValues; + /** Load and validate a decompressed IntaRNA binary accessibility archive. */ + void parseBinary( std::istream & inStream ); + //! maximal available window size size_t availMaxLength; diff --git a/src/IntaRNA/Makefile.am b/src/IntaRNA/Makefile.am index 56b7ed87..2115c11b 100644 --- a/src/IntaRNA/Makefile.am +++ b/src/IntaRNA/Makefile.am @@ -9,6 +9,8 @@ AUTOMAKE_OPTIONS = std-options subdir-objects AM_DEFAULT_SOURCE_EXT = .cpp +noinst_HEADERS = AccessibilityArchive.h + ############################################################################### # THE INTARNA LIBRARY ############################################################################### diff --git a/src/IntaRNA/general.cpp b/src/IntaRNA/general.cpp index c328c89a..facd5f38 100644 --- a/src/IntaRNA/general.cpp +++ b/src/IntaRNA/general.cpp @@ -44,7 +44,7 @@ std::ostream* newOutputStream(const std::string& out) BOOST_IOS::openmode fopenmode = BOOST_IOS::out; // Gzip-Erkennung - if (boost::iends_with(out, ".gz")) { + if (boost::iends_with(out, ".gz") || boost::iends_with(out, ".agz")) { // Gzip-Kompressor mit explizit vergrößertem Puffer hinzufügen fstream->push(bio::gzip_compressor(), BUFFER_SIZE); fopenmode |= BOOST_IOS::binary; @@ -115,7 +115,7 @@ newInputStream( const std::string & in ) BOOST_IOS::openmode fopenmode = BOOST_IOS::in; // gzipped input file stream - if (in.size()>3 && boost::iequals(in.substr(in.size()-3,3),".gz")) { + if (boost::iends_with(in, ".gz") || boost::iends_with(in, ".agz")) { // gzip compression fstream->push( bio::gzip_decompressor() ); // binary input diff --git a/src/IntaRNA/general.h b/src/IntaRNA/general.h index ab9dfdcc..1a982982 100644 --- a/src/IntaRNA/general.h +++ b/src/IntaRNA/general.h @@ -190,7 +190,7 @@ namespace IntaRNA { * - & std::cerr : if outName == STDERR * - new std::fstream( outName ) : else if outName non-empty * - * If the filename ends in '.gz', gzip compression and binary output is enabled. + * If the filename ends in '.gz' or '.agz', gzip compression and binary output is enabled. * * @param outName the name of the output to open. use STDOUT/STDERR for the * respective output stream or otherwise a filename to be created. @@ -218,7 +218,7 @@ deleteOutputStream( std::ostream *& outStream ); * - & std::cin : if inName == STDIN * - new std::fstream( inName ) : else inName non-empty * - * If the filename ends in '.gz', gzip compression and binary input is enabled. + * If the filename ends in '.gz' or '.agz', gzip compression and binary input is enabled. * * @param inName the name of the input to open. use STDIN for the * respective input stream or otherwise a filename to be created. diff --git a/src/bin/CommandLineParsing.cpp b/src/bin/CommandLineParsing.cpp index 7fa60e4d..22af1cef 100644 --- a/src/bin/CommandLineParsing.cpp +++ b/src/bin/CommandLineParsing.cpp @@ -7,6 +7,7 @@ #include #include #include +#include #if INTARNA_MULITHREADING #include @@ -400,8 +401,8 @@ CommandLineParsing::CommandLineParsing( const Personality personality ) , std::string("accessibility computation :" "\n 'N' no accessibility contributions" "\n 'C' computation of accessibilities" - "\n 'P' unpaired probabilities in RNAplfold format from --qAccFile" - "\n 'E' ED values in RNAplfold Pu-like format from --qAccFile" + "\n 'P' unpaired probabilities in RNAplfold format or binary .agz from --qAccFile" + "\n 'E' ED values in RNAplfold Pu-like format or binary .agz from --qAccFile" ).c_str()) (qAccW.name.c_str() , value(&(qAccW.val)) @@ -434,7 +435,7 @@ CommandLineParsing::CommandLineParsing( const Personality personality ) ).c_str()) ("qAccFile" , value(&(qAccFile)) - , std::string("accessibility computation : the file/stream to be parsed, if --qAcc is to be read from file. Used 'STDIN' if to read from standard input stream.").c_str()) + , std::string("accessibility computation : the file/stream to be parsed, if --qAcc is to be read from file. A .agz suffix selects IntaRNA binary ED data for either P or E; use STDIN for text input.").c_str()) (qIntLenMax.name.c_str() , value(&(qIntLenMax.val)) ->default_value(qIntLenMax.def) @@ -515,8 +516,8 @@ CommandLineParsing::CommandLineParsing( const Personality personality ) , std::string("accessibility computation :" "\n 'N' no accessibility contributions" "\n 'C' computation of accessibilities" - "\n 'P' unpaired probabilities in RNAplfold format from --tAccFile" - "\n 'E' ED values in RNAplfold Pu-like format from --tAccFile" + "\n 'P' unpaired probabilities in RNAplfold format or binary .agz from --tAccFile" + "\n 'E' ED values in RNAplfold Pu-like format or binary .agz from --tAccFile" ).c_str()) (tAccW.name.c_str() , value(&(tAccW.val)) @@ -549,7 +550,7 @@ CommandLineParsing::CommandLineParsing( const Personality personality ) ).c_str()) ("tAccFile" , value(&(tAccFile)) - , std::string("accessibility computation : the file/stream to be parsed, if --tAcc is to be read from file. Used 'STDIN' if to read from standard input stream.").c_str()) + , std::string("accessibility computation : the file/stream to be parsed, if --tAcc is to be read from file. A .agz suffix selects IntaRNA binary ED data for either P or E; use STDIN for text input.").c_str()) (tIntLenMax.name.c_str() , value(&(tIntLenMax.val)) ->default_value(tIntLenMax.def) @@ -917,12 +918,12 @@ CommandLineParsing::CommandLineParsing( const Personality personality ) " ADDITIONAL output:" "\n 'qMinE:' (query) for each position the minimal energy of any interaction covering the position (CSV format)" "\n 'qSpotProb:' (query) for each position the probability that is is covered by an interaction covering (CSV format)" - "\n 'qAcc:' (query) ED accessibility values ('qPu'-like format)." - "\n 'qPu:' (query) unpaired probabilities values (RNAplfold format)." + "\n 'qAcc:' (query) ED accessibility values ('qPu'-like format; .agz for binary)." + "\n 'qPu:' (query) unpaired probability values (RNAplfold format; .agz stores exact ED data)." "\n 'tMinE:' (target) for each position the minimal energy of any interaction covering the position (CSV format)" "\n 'tSpotProb:' (target) for each position the probability that is is covered by an interaction covering (CSV format)" - "\n 'tAcc:' (target) ED accessibility values ('tPu'-like format)." - "\n 'tPu:' (target) unpaired probabilities values (RNAplfold format)." + "\n 'tAcc:' (target) ED accessibility values ('tPu'-like format; .agz for binary)." + "\n 'tPu:' (target) unpaired probability values (RNAplfold format; .agz stores exact ED data)." "\n 'pMinE:' (target+query) for each index pair the minimal energy of any interaction covering the pair (CSV format)" "\n 'spotProb:' (target+query) tracks for a given set of interaction spots their probability to be covered by an interaction. If no spots are provided, probabilities for all index combinations are computed. Spots are encoded by comma-separated 'idxT&idxQ' pairs (target-query). For each spot a probability is provided in concert with the probability that none of the spots (encoded by '0&0') is covered (CSV format). The spot encoding is followed colon-separated by the output stream/file name, eg. '--out=\"spotProb:3&76,59&2:STDERR\"'. NOTE: value has to be quoted due to '&' symbol!" "\nFor each, provide a file name or STDOUT/STDERR to write to the respective output stream." @@ -1938,7 +1939,9 @@ getQueryAccessibility( const size_t sequenceNumber ) const case 'E' : // drop to next handling case 'P' : { // VRNA RNAplfold unpaired probability file output - std::istream * accStream = newInputStream( getFullFilename(qAccFile, NULL, &(seq)) ); + const std::string filename = getFullFilename(qAccFile, NULL, &(seq)); + auto closeInput = [](std::istream * stream) { deleteInputStream(stream); }; + std::unique_ptr accStream(newInputStream(filename), closeInput); if (accStream == NULL) { throw std::runtime_error("accessibility parsing of --qAccFile : could not open file '"+qAccFile+"'"); } @@ -1946,10 +1949,9 @@ getQueryAccessibility( const size_t sequenceNumber ) const , qIntLenMax.val , &accConstraint , *accStream - , (qAcc.val == 'P' ? AccessibilityFromStream::Pu_RNAplfold_Text : AccessibilityFromStream::ED_RNAplfold_Text) + , (boost::iends_with(filename, ".agz") ? AccessibilityFromStream::IntaRNA_Binary + : (qAcc.val == 'P' ? AccessibilityFromStream::Pu_RNAplfold_Text : AccessibilityFromStream::ED_RNAplfold_Text)) , vrnaHandler.getRT() ); - // cleanup - deleteInputStream( accStream ); return acc; } @@ -2013,7 +2015,9 @@ getTargetAccessibility( const size_t sequenceNumber ) const case 'E' : // drop to next handling case 'P' : { // VRNA RNAplfold unpaired probability file output - std::istream * accStream = newInputStream( getFullFilename(tAccFile, &(seq), NULL) ); + const std::string filename = getFullFilename(tAccFile, &(seq), NULL); + auto closeInput = [](std::istream * stream) { deleteInputStream(stream); }; + std::unique_ptr accStream(newInputStream(filename), closeInput); if (accStream == NULL) { throw std::runtime_error("accessibility parsing of --tAccFile : could not open file '"+tAccFile+"'"); } @@ -2022,10 +2026,9 @@ getTargetAccessibility( const size_t sequenceNumber ) const , tIntLenMax.val , &accConstraint , *accStream - , ( tAcc.val == 'P' ? AccessibilityFromStream::Pu_RNAplfold_Text : AccessibilityFromStream::ED_RNAplfold_Text ) + , (boost::iends_with(filename, ".agz") ? AccessibilityFromStream::IntaRNA_Binary + : (tAcc.val == 'P' ? AccessibilityFromStream::Pu_RNAplfold_Text : AccessibilityFromStream::ED_RNAplfold_Text)) , vrnaHandler.getRT() ); - // cleanup - deleteInputStream( accStream ); return acc; } case 'C' : // compute accessibilities @@ -2761,20 +2764,28 @@ writeAccessibility( const Accessibility& acc, const std::string & fileOrStream, return; // setup output stream - std::ostream * out = newOutputStream( fileOrStream ); + auto closeOutput = [](std::ostream * stream) noexcept { + try { deleteOutputStream(stream); } + catch (...) { if (stream != &std::cout && stream != &std::cerr) delete stream; } + }; + std::unique_ptr out(newOutputStream(fileOrStream), closeOutput); if (out == NULL) { throw std::runtime_error("could not open output file '"+fileOrStream +"' for "+(writeED?"accessibility":"unpaired probability")+" output"); } // write data to stream - if (writeED) { + if (boost::iends_with(fileOrStream, ".agz")) { + acc.writeBinary(*out); + } else if (writeED) { acc.writeRNAplfold_ED_text( *out ); } else { acc.writeRNAplfold_Pu_text( *out, vrnaHandler.getRT() ); } - // clean up - deleteOutputStream( out ); + // Explicit close propagates compression/file errors before releasing ownership. + std::ostream * stream = out.get(); + deleteOutputStream(stream); + out.release(); } //////////////////////////////////////////////////////////////////////////// diff --git a/tests/AccessibilityBinary_test.cpp b/tests/AccessibilityBinary_test.cpp new file mode 100644 index 00000000..334d9536 --- /dev/null +++ b/tests/AccessibilityBinary_test.cpp @@ -0,0 +1,137 @@ +#include "catch.hpp" + +#include "IntaRNA/AccessibilityFromStream.h" +#include "IntaRNA/AccessibilityBasePair.h" +#include "IntaRNA/AccessibilityDisabled.h" +#include "IntaRNA/AccessibilityVrna.h" +#include "IntaRNA/ReverseAccessibility.h" +#include +#include +#include + +using namespace IntaRNA; + +namespace { + +void checkBinaryRoundTrip( const Accessibility & source, size_t requestedLength ) +{ + std::stringstream bytes; + source.writeBinary(bytes); + AccessibilityFromStream loaded(source.getSequence(), requestedLength, nullptr, + bytes, AccessibilityFromStream::IntaRNA_Binary, 0.0); + const size_t expectedLength = std::min(source.getMaxLength(), + requestedLength == 0 ? source.getSequence().size() : requestedLength); + REQUIRE(loaded.getMaxLength() == expectedLength); + for (size_t i = 0; i < source.getSequence().size(); ++i) { + for (size_t j = i; j < source.getSequence().size(); ++j) { + const E_type expected = j-i <= expectedLength ? source.getED(i,j) : Accessibility::ED_UPPER_BOUND; + REQUIRE(loaded.getED(i,j) == expected); + } + } + // Re-exporting loaded matrices must retain the dangling-end band too. + std::stringstream second; + loaded.writeBinary(second); + AccessibilityFromStream reloaded(source.getSequence(), 0, nullptr, + second, AccessibilityFromStream::IntaRNA_Binary, 10.0); + REQUIRE(reloaded.getMaxLength() == loaded.getMaxLength()); + for (size_t i = 0; i < source.getSequence().size(); ++i) + for (size_t j = i; j < source.getSequence().size(); ++j) + REQUIRE(reloaded.getED(i,j) == loaded.getED(i,j)); +} + +// Locate the application header without depending on the Boost archive preamble. +size_t headerOffset( const std::string & data ) +{ + const std::uint32_t magic = 0x49414343; + const size_t offset = data.find(std::string(reinterpret_cast(&magic), sizeof(magic))); + REQUIRE(offset != std::string::npos); + return offset; +} + +template std::string replaceField( std::string data, size_t offset, T value ) +{ + std::memcpy(&data.at(offset), &value, sizeof(value)); + return data; +} + +} // namespace + +TEST_CASE("Binary accessibility preserves exact ED matrices", "[AccessibilityBinary]") +{ +#include "testEasyLoggingSetup.icc" + RnaSequence rna("test", "GGGGAAAACCCCUAGC"); + VrnaHandler vrna(37, "Turner04", false, false); + for (size_t length : {size_t(1), size_t(5), rna.size()}) { + AccessibilityVrna folded(rna, length, nullptr, vrna, rna.size()); + AccessibilityBasePair basePair(rna, length, nullptr); + AccessibilityDisabled disabled(rna, length, nullptr); + ReverseAccessibility reversed(folded); + for (const Accessibility * acc : {static_cast(&folded), + static_cast(&basePair), static_cast(&disabled), + static_cast(&reversed)}) { + for (size_t requested : {size_t(0), size_t(1), size_t(3), rna.size()}) + checkBinaryRoundTrip(*acc, requested); + } + } + SECTION("constraints and infinity are retained in the values") { + AccessibilityConstraint constraint(rna, "xp............px", 0, "", "", ""); + AccessibilityVrna folded(rna, 5, &constraint, vrna, rna.size()); + REQUIRE(folded.getED(1,1) == Accessibility::ED_UPPER_BOUND); + checkBinaryRoundTrip(folded, 0); + } + SECTION("text input with only the dangling-end column") { + RnaSequence shortRna("short", "AC"); + std::istringstream text("#unpaired probabilities\n #i$\tl=1\n1\t0.5\n2\t1.0\n"); + AccessibilityFromStream loaded(shortRna, 0, nullptr, text, + AccessibilityFromStream::Pu_RNAplfold_Text, 1.0); + REQUIRE(loaded.getMaxLength() == 0); + checkBinaryRoundTrip(loaded, 0); + } + SECTION("short sequence") { + RnaSequence shortRna("short", "A"); + AccessibilityDisabled disabled(shortRna, 0, nullptr); + checkBinaryRoundTrip(disabled, 0); + } +} + +TEST_CASE("Binary accessibility rejects invalid archives", "[AccessibilityBinary]") +{ +#include "testEasyLoggingSetup.icc" + RnaSequence rna("test", "GGGGAAAACCCC"); + AccessibilityDisabled source(rna, 5, nullptr); + std::stringstream stream; + source.writeBinary(stream); + const std::string valid = stream.str(); + const size_t header = headerOffset(valid); + auto rejects = [&](const std::string & data) { + std::istringstream input(data); + REQUIRE_THROWS_AS(AccessibilityFromStream(rna, 3, nullptr, input, + AccessibilityFromStream::IntaRNA_Binary, 1.0), std::exception); + }; + SECTION("every truncated archive") { + for (size_t size = 0; size < valid.size(); ++size) rejects(valid.substr(0,size)); + } + SECTION("header and matrix values") { + rejects("#unpaired probabilities\n"); + rejects(valid + "trailing data"); + rejects(replaceField(valid, header, std::uint32_t(0))); + rejects(replaceField(valid, header+4, std::uint32_t(2))); + rejects(replaceField(valid, header+8, std::uint64_t(rna.size()+1))); + rejects(replaceField(valid, header+16, std::uint64_t(0))); + rejects(replaceField(valid, header+16, std::uint64_t(-1))); + rejects(replaceField(valid, header+24, E_type(1))); + rejects(replaceField(valid, valid.size()-sizeof(E_type), E_type(-1))); + rejects(replaceField(valid, valid.size()-sizeof(E_type), Accessibility::ED_UPPER_BOUND+1)); + } + SECTION("different sequence of the same length") { + RnaSequence other("other", "AAAAAAAACCCC"); + std::istringstream input(valid); + REQUIRE_THROWS_WITH(AccessibilityFromStream(other, 5, nullptr, input, + AccessibilityFromStream::IntaRNA_Binary, 1.0), Catch::Contains("sequence mismatch")); + } + SECTION("output failure") { + std::ostringstream output; + output.setstate(std::ios::badbit); + REQUIRE_THROWS(source.writeBinary(output)); + } +} diff --git a/tests/Makefile.am b/tests/Makefile.am index 3bc74a69..946632b3 100644 --- a/tests/Makefile.am +++ b/tests/Makefile.am @@ -15,7 +15,7 @@ TEST_EXTENSIONS = $(EXEEXT) .sh SH_LOG_COMPILER = $(SHELL) # the script needed for tests -dist_check_SCRIPTS = runIntaRNA.sh +dist_check_SCRIPTS = runIntaRNA.sh runAccessibilityBinary.sh # the program to build check_PROGRAMS = runApiTests @@ -26,6 +26,7 @@ runApiTests_SOURCES = \ testEasyLoggingSetup.icc \ AccessibilityConstraint_test.cpp \ AccessibilityFromStream_test.cpp \ + AccessibilityBinary_test.cpp \ AccessibilityBasePair_test.cpp \ AccessibilityVrna_test.cpp \ HelixConstraint_test.cpp \ diff --git a/tests/runAccessibilityBinary.sh b/tests/runAccessibilityBinary.sh new file mode 100644 index 00000000..4ef5acae --- /dev/null +++ b/tests/runAccessibilityBinary.sh @@ -0,0 +1,56 @@ +#!/usr/bin/env bash +# Exercise actual filename dispatch, gzip finalization and prediction reuse. +set -euo pipefail +bin="$INTARNABINPATH/src/bin/IntaRNA" +tmp=$(mktemp -d) +trap 'rm -rf "$tmp"' EXIT +query=GGGGAAAACCCCUAGC +target=GCUAGGGGUUUUCCCC +common=(--query="$query" --target="$target" --qIntLenMax=5 --tIntLenMax=5 + --seedBP=3 --outMode=C --outCsvCols=id1,start1,end1,id2,start2,end2,E + --threads=1 --default-log-file=/dev/null) +"$bin" "${common[@]}" --out="qAcc:$tmp/query.agz" --out="tPu:$tmp/target.AGZ" > "$tmp/computed" +test "$(wc -l < "$tmp/computed")" -gt 1 +gzip -t "$tmp/query.agz" "$tmp/target.AGZ" +for mode in E P; do + "$bin" "${common[@]}" --qAcc="$mode" --tAcc="$mode" \ + --qAccFile="$tmp/query.agz" --tAccFile="$tmp/target.AGZ" > "$tmp/reused" + cmp "$tmp/computed" "$tmp/reused" +done +# A smaller requested band can reuse a wider stored matrix. +smaller=("${common[@]}") +smaller[2]=--qIntLenMax=3 +smaller[3]=--tIntLenMax=3 +"$bin" "${smaller[@]}" > "$tmp/smaller" +"$bin" "${smaller[@]}" --qAcc=E --tAcc=E \ + --qAccFile="$tmp/query.agz" --tAccFile="$tmp/target.AGZ" > "$tmp/smaller-reused" +cmp "$tmp/smaller" "$tmp/smaller-reused" +# Preserve the existing compressed text format and dispatch. +"$bin" "${common[@]}" --out="qAcc:$tmp/query.txt.gz" --out="tPu:$tmp/target.txt.gz" > "$tmp/text-computed" +gzip -cd "$tmp/query.txt.gz" > "$tmp/query.txt" +grep -q '^#ensemble' "$tmp/query.txt" +"$bin" "${common[@]}" --qAcc=E --tAcc=P --qAccFile="$tmp/query.txt.gz" \ + --tAccFile="$tmp/target.txt.gz" > "$tmp/text-reused" +cmp "$tmp/text-computed" "$tmp/text-reused" +# Multi-sequence suffixes must keep .agz as the extension. +printf '>q1\n%s\n>q2\n%s\n' "$query" "$target" > "$tmp/queries.fa" +multi=("${common[@]}") +multi[0]="--query=$tmp/queries.fa" +"$bin" "${multi[@]}" --out="qPu:$tmp/multi.agz" > "$tmp/multi-computed" +test -s "$tmp/multi-s1.agz" +test -s "$tmp/multi-s2.agz" +"$bin" "${multi[@]}" --qAcc=P --qAccFile="$tmp/multi.agz" > "$tmp/multi-reused" +cmp "$tmp/multi-computed" "$tmp/multi-reused" +# Damaged trailers must fail even if every ED value has been decompressed. +dd if="$tmp/query.agz" of="$tmp/truncated.agz" bs=1 count="$(( $(wc -c < "$tmp/query.agz") - 4 ))" 2>/dev/null +cp "$tmp/query.agz" "$tmp/checksum.agz" +printf '\000\000\000\000' | dd of="$tmp/checksum.agz" bs=1 \ + seek="$(( $(wc -c < "$tmp/query.agz") - 8 ))" conv=notrunc 2>/dev/null +for damaged in truncated checksum; do + status=0 + "$bin" "${common[@]}" --qAcc=E --qAccFile="$tmp/$damaged.agz" \ + > "$tmp/bad.out" 2> "$tmp/bad.err" || status=$? + # A controlled error, not success or termination by a signal. + test "$status" -eq 255 +done +echo 'Binary accessibility CLI checks passed' From 4c22195c4842021f3b9f9e805b701b4b40ed7ab4 Mon Sep 17 00:00:00 2001 From: Martin Raden Date: Wed, 30 Sep 2026 12:26:31 +0200 Subject: [PATCH 2/4] revised --- src/IntaRNA/AccessibilityArchive.h | 4 ++-- src/bin/CommandLineParsing.cpp | 24 ++++++++++++++---------- 2 files changed, 16 insertions(+), 12 deletions(-) diff --git a/src/IntaRNA/AccessibilityArchive.h b/src/IntaRNA/AccessibilityArchive.h index a51cbdb4..d99b45d0 100644 --- a/src/IntaRNA/AccessibilityArchive.h +++ b/src/IntaRNA/AccessibilityArchive.h @@ -2,8 +2,8 @@ #define INTARNA_ACCESSIBILITYARCHIVE_H_ #include "IntaRNA/Accessibility.h" +#include "IntaRNA/Matrix.h" -#include #include #include #include @@ -28,7 +28,7 @@ namespace IntaRNA { */ class AccessibilityArchive { public: - typedef boost::numeric::ublas::banded_matrix EdMatrix; + typedef UpperBandedMatrix EdMatrix; /** Create a read-only serialization view of any accessibility implementation. */ explicit AccessibilityArchive( const Accessibility & source ); diff --git a/src/bin/CommandLineParsing.cpp b/src/bin/CommandLineParsing.cpp index 22af1cef..e8b0a643 100644 --- a/src/bin/CommandLineParsing.cpp +++ b/src/bin/CommandLineParsing.cpp @@ -401,8 +401,8 @@ CommandLineParsing::CommandLineParsing( const Personality personality ) , std::string("accessibility computation :" "\n 'N' no accessibility contributions" "\n 'C' computation of accessibilities" - "\n 'P' unpaired probabilities in RNAplfold format or binary .agz from --qAccFile" - "\n 'E' ED values in RNAplfold Pu-like format or binary .agz from --qAccFile" + "\n 'P' unpaired probabilities in RNAplfold format or IntaRNA's binary format (.agz file) from --qAccFile" + "\n 'E' ED values in RNAplfold Pu-like format or IntaRNA's binary format (.agz file) from --qAccFile" ).c_str()) (qAccW.name.c_str() , value(&(qAccW.val)) @@ -435,7 +435,9 @@ CommandLineParsing::CommandLineParsing( const Personality personality ) ).c_str()) ("qAccFile" , value(&(qAccFile)) - , std::string("accessibility computation : the file/stream to be parsed, if --qAcc is to be read from file. A .agz suffix selects IntaRNA binary ED data for either P or E; use STDIN for text input.").c_str()) + , std::string("accessibility computation : the file/stream to be parsed, if --qAcc is to be read from file." + " A '.agz' file ending identifies IntaRNA binary ED output format for either P or E, otherwise an RNAplfold-like text-based format is expected." + " Use STDIN for text input from standard input stream.").c_str()) (qIntLenMax.name.c_str() , value(&(qIntLenMax.val)) ->default_value(qIntLenMax.def) @@ -516,8 +518,8 @@ CommandLineParsing::CommandLineParsing( const Personality personality ) , std::string("accessibility computation :" "\n 'N' no accessibility contributions" "\n 'C' computation of accessibilities" - "\n 'P' unpaired probabilities in RNAplfold format or binary .agz from --tAccFile" - "\n 'E' ED values in RNAplfold Pu-like format or binary .agz from --tAccFile" + "\n 'P' unpaired probabilities in RNAplfold format or IntaRNA's binary format (.agz file) from --tAccFile" + "\n 'E' ED values in RNAplfold Pu-like format or IntaRNA's binary format (.agz file) from --tAccFile" ).c_str()) (tAccW.name.c_str() , value(&(tAccW.val)) @@ -550,7 +552,9 @@ CommandLineParsing::CommandLineParsing( const Personality personality ) ).c_str()) ("tAccFile" , value(&(tAccFile)) - , std::string("accessibility computation : the file/stream to be parsed, if --tAcc is to be read from file. A .agz suffix selects IntaRNA binary ED data for either P or E; use STDIN for text input.").c_str()) + , std::string("accessibility computation : the file/stream to be parsed, if --tAcc is to be read from file." + " A '.agz' file ending identifies IntaRNA binary ED output format for either P or E, otherwise an RNAplfold-like text-based format is expected." + " Use STDIN for text input from standard input stream.").c_str()) (tIntLenMax.name.c_str() , value(&(tIntLenMax.val)) ->default_value(tIntLenMax.def) @@ -918,12 +922,12 @@ CommandLineParsing::CommandLineParsing( const Personality personality ) " ADDITIONAL output:" "\n 'qMinE:' (query) for each position the minimal energy of any interaction covering the position (CSV format)" "\n 'qSpotProb:' (query) for each position the probability that is is covered by an interaction covering (CSV format)" - "\n 'qAcc:' (query) ED accessibility values ('qPu'-like format; .agz for binary)." - "\n 'qPu:' (query) unpaired probability values (RNAplfold format; .agz stores exact ED data)." + "\n 'qAcc:' (query) ED accessibility values ('qPu'-like text format OR '.agz' file ending for IntaRNA's binary format)." + "\n 'qPu:' (query) unpaired probability values (RNAplfold format OR '.agz' file ending for IntaRNA's binary format)." "\n 'tMinE:' (target) for each position the minimal energy of any interaction covering the position (CSV format)" "\n 'tSpotProb:' (target) for each position the probability that is is covered by an interaction covering (CSV format)" - "\n 'tAcc:' (target) ED accessibility values ('tPu'-like format; .agz for binary)." - "\n 'tPu:' (target) unpaired probability values (RNAplfold format; .agz stores exact ED data)." + "\n 'tAcc:' (target) ED accessibility values ('tPu'-like format OR '.agz' file ending for IntaRNA's binary format)." + "\n 'tPu:' (target) unpaired probability values (RNAplfold format OR '.agz' file ending for IntaRNA's binary format)." "\n 'pMinE:' (target+query) for each index pair the minimal energy of any interaction covering the pair (CSV format)" "\n 'spotProb:' (target+query) tracks for a given set of interaction spots their probability to be covered by an interaction. If no spots are provided, probabilities for all index combinations are computed. Spots are encoded by comma-separated 'idxT&idxQ' pairs (target-query). For each spot a probability is provided in concert with the probability that none of the spots (encoded by '0&0') is covered (CSV format). The spot encoding is followed colon-separated by the output stream/file name, eg. '--out=\"spotProb:3&76,59&2:STDERR\"'. NOTE: value has to be quoted due to '&' symbol!" "\nFor each, provide a file name or STDOUT/STDERR to write to the respective output stream." From 2e2c5f0a111805822608a180bfa785288bcca501 Mon Sep 17 00:00:00 2001 From: AlexanderMitrofanov Date: Wed, 30 Sep 2026 13:43:26 +0200 Subject: [PATCH 3/4] Serialize accessibility matrices directly and measure gzip trade-offs --- ChangeLog | 21 +++ README.md | 5 + doc/Makefile.am | 2 + doc/benchmarks/accessibility-ecoli.csv | 25 ++++ doc/benchmarks/accessibility-random.csv | 25 ++++ doc/benchmarks/accessibility.cpp | 179 ++++++++++++++++++------ doc/benchmarks/accessibility.md | 133 +++++++++++++----- src/IntaRNA/Accessibility.cpp | 8 +- src/IntaRNA/Accessibility.h | 13 +- src/IntaRNA/AccessibilityArchive.h | 62 +++++--- src/IntaRNA/AccessibilityFromStream.h | 15 ++ src/IntaRNA/AccessibilityVrna.h | 15 ++ src/IntaRNA/Matrix.h | 34 +++++ tests/AccessibilityBinary_test.cpp | 62 ++++++++ tests/Matrix_test.cpp | 21 +++ 15 files changed, 521 insertions(+), 99 deletions(-) create mode 100644 doc/benchmarks/accessibility-ecoli.csv create mode 100644 doc/benchmarks/accessibility-random.csv diff --git a/ChangeLog b/ChangeLog index ca939a5a..2458a32a 100644 --- a/ChangeLog +++ b/ChangeLog @@ -25,6 +25,8 @@ ## Technical changes and Optimizations +- serialize stored accessibility rows directly and load retained rows in place; + benchmark direct/generic export and raw/gzip binary I/O (PR #250) - add repository-specific AI coding guidance in AGENTS.md (issue #247) - replace uBLAS storage matrices with std::vector and std::mdspan, retaining compact upper-band and triangular storage; configure selects native mdspan @@ -59,6 +61,25 @@ ################################################################################ ################################################################################ +260930 Alexander Mitrofanov + * IntaRNA/Matrix.h : + + const/mutable upper-band row views that exclude structural zeros and padding + * IntaRNA/Accessibility, AccessibilityArchive, AccessibilityVrna, AccessibilityFromStream : + * export stored ED rows directly without copies or per-cell getED() calls + * retain generic export for other producers and constrained accessibility data + * deserialize retained rows into final storage; validate discarded band tails + * preserve version 1 archive bytes and existing constraint masking + * tests/AccessibilityBinary_test.cpp, tests/Matrix_test.cpp : + + verify row aliasing, byte-identical direct/generic output, zero interface + lookups on direct export, constrained short sequences and corrupt discarded data + * doc/benchmarks/accessibility.cpp, doc/benchmarks/accessibility.md, doc/Makefile.am : + * compare direct/generic and raw/gzip binary I/O with exact ED verification + + reproducible synthetic RNA and E. coli measurements with per-trial CSV results + * retain gzip: both 100,000-base datasets shrink by 58-59%, despite slower I/O + * README.md : + + document direct matrix I/O and link the compression measurements + * addresses review requests in https://github.com/BackofenLab/IntaRNA/pull/250 + 260929 Alexander Mitrofanov * IntaRNA/Accessibility, AccessibilityArchive, AccessibilityFromStream : + versioned Boost binary serialization of exact ED matrices and sequence data diff --git a/README.md b/README.md index e3f583ba..4db61ce4 100644 --- a/README.md +++ b/README.md @@ -2173,6 +2173,11 @@ compatible architecture and Boost archive version. Use RNAplfold-style text for portable interchange. For C++ callers, `Accessibility::writeBinary()` and `AccessibilityFromStream::IntaRNA_Binary` operate on decompressed archive streams; `newOutputStream()`/`newInputStream()` supply gzip compression for `.agz` files. +Unconstrained `AccessibilityVrna` and `AccessibilityFromStream` export their +stored matrix rows directly; other implementations and constrained data use +`getED()` to preserve their accessibility semantics. Loading fills the retained +matrix rows directly. The [I/O benchmark](doc/benchmarks/accessibility.md) +compares direct and generic export, and quantifies gzip's size/time trade-off. ##### Use case examples for read/write accessibilities and unpaired probabilities diff --git a/doc/Makefile.am b/doc/Makefile.am index 985e6fae..2b004d3d 100644 --- a/doc/Makefile.am +++ b/doc/Makefile.am @@ -6,6 +6,8 @@ EXTRA_DIST = \ benchmarks/accessibility.cpp \ benchmarks/accessibility.md \ + benchmarks/accessibility-random.csv \ + benchmarks/accessibility-ecoli.csv \ conda.txt \ doxygen.cfg \ latex-deps/adjcalc.sty \ diff --git a/doc/benchmarks/accessibility-ecoli.csv b/doc/benchmarks/accessibility-ecoli.csv new file mode 100644 index 00000000..eb485367 --- /dev/null +++ b/doc/benchmarks/accessibility-ecoli.csv @@ -0,0 +1,25 @@ +trial,producer,serialization,format,write_seconds,read_seconds,bytes,different_cells +0,vrna,direct,raw,0.0343495,0.0239433,40479873,0 +0,vrna,direct,gzip,2.9852,0.267999,16850736,0 +0,vrna,generic,raw,0.383549,0.0217923,40479873,0 +0,vrna,generic,gzip,3.01589,0.266227,16850736,0 +0,stream,direct,raw,0.0339346,0.0211565,40479873,0 +0,stream,direct,gzip,2.87353,0.268133,16850736,0 +0,stream,generic,raw,0.161544,0.0225503,40479873,0 +0,stream,generic,gzip,2.93791,0.255312,16850736,0 +1,vrna,generic,gzip,3.00123,0.261814,16850736,0 +1,stream,direct,raw,0.108538,0.0238946,40479873,0 +1,stream,direct,gzip,2.99695,0.260998,16850736,0 +1,stream,generic,raw,0.183074,0.0237418,40479873,0 +1,stream,generic,gzip,3.05163,0.267577,16850736,0 +1,vrna,direct,raw,0.10758,0.0236052,40479873,0 +1,vrna,direct,gzip,3.0188,0.273299,16850736,0 +1,vrna,generic,raw,0.164736,0.023654,40479873,0 +2,stream,generic,raw,0.160865,0.0269454,40479873,0 +2,stream,generic,gzip,3.12204,0.277653,16850736,0 +2,vrna,direct,raw,0.112687,0.0239146,40479873,0 +2,vrna,direct,gzip,3.0273,0.270982,16850736,0 +2,vrna,generic,raw,0.455994,0.0238453,40479873,0 +2,vrna,generic,gzip,3.2056,0.309473,16850736,0 +2,stream,direct,raw,0.452168,0.0293377,40479873,0 +2,stream,direct,gzip,3.2628,0.274283,16850736,0 diff --git a/doc/benchmarks/accessibility-random.csv b/doc/benchmarks/accessibility-random.csv new file mode 100644 index 00000000..ab06fb25 --- /dev/null +++ b/doc/benchmarks/accessibility-random.csv @@ -0,0 +1,25 @@ +trial,producer,serialization,format,write_seconds,read_seconds,bytes,different_cells +0,vrna,direct,raw,0.0313849,0.0215885,40479873,0 +0,vrna,direct,gzip,2.79625,0.241071,16620982,0 +0,vrna,generic,raw,0.164226,0.0219357,40479873,0 +0,vrna,generic,gzip,2.8534,0.241125,16620982,0 +0,stream,direct,raw,0.0320046,0.0212699,40479873,0 +0,stream,direct,gzip,2.82735,0.273768,16620982,0 +0,stream,generic,raw,0.166378,0.0203288,40479873,0 +0,stream,generic,gzip,3.05354,0.246805,16620982,0 +1,vrna,generic,gzip,2.85945,0.245383,16620982,0 +1,stream,direct,raw,0.125027,0.0224652,40479873,0 +1,stream,direct,gzip,2.76234,0.257171,16620982,0 +1,stream,generic,raw,0.146028,0.0204308,40479873,0 +1,stream,generic,gzip,2.83353,0.263496,16620982,0 +1,vrna,direct,raw,0.104695,0.0224409,40479873,0 +1,vrna,direct,gzip,2.94117,0.248038,16620982,0 +1,vrna,generic,raw,0.157238,0.0218882,40479873,0 +2,stream,generic,raw,0.149852,0.0214029,40479873,0 +2,stream,generic,gzip,3.06712,0.262208,16620982,0 +2,vrna,direct,raw,0.108833,0.0224185,40479873,0 +2,vrna,direct,gzip,2.77458,0.259067,16620982,0 +2,vrna,generic,raw,0.153023,0.0218832,40479873,0 +2,vrna,generic,gzip,2.78912,0.241287,16620982,0 +2,stream,direct,raw,0.101867,0.0202582,40479873,0 +2,stream,direct,gzip,2.68447,0.239721,16620982,0 diff --git a/doc/benchmarks/accessibility.cpp b/doc/benchmarks/accessibility.cpp index c3ecfa34..40969316 100644 --- a/doc/benchmarks/accessibility.cpp +++ b/doc/benchmarks/accessibility.cpp @@ -1,57 +1,144 @@ -// Compare compressed text and binary accessibility I/O after one ViennaRNA fold. -// Usage: accessibility-benchmark SEQUENCE_LENGTH EXISTING_OUTPUT_DIRECTORY +// Compare generic/direct matrix export and raw/gzip binary I/O on identical data. +// Usage: accessibility-benchmark SEQUENCE_LENGTH OUTPUT_DIRECTORY [FASTA_FILE] #include #include +#include +#include +#include #include +#include #include -#include +#include #include +#include #include +#include +#include #include INITIALIZE_EASYLOGGINGPP -int main(int argc, char **argv) { - if (argc != 3) return 2; - el::Loggers::reconfigureAllLoggers(el::ConfigurationType::Enabled, "false"); - const size_t n = std::stoull(argv[1]); - std::string sequence; - sequence.reserve(n); - unsigned state = 245; - for (size_t i=0;i> 30]; - } - IntaRNA::RnaSequence rna("benchmark", sequence); - IntaRNA::VrnaHandler vrna(37, "Turner04", false, false); - IntaRNA::AccessibilityVrna source(rna, 100, nullptr, vrna, 150); - using Clock = std::chrono::steady_clock; - for (int run=0;run<3;++run) { - for (bool binary : {false, true}) { - const auto path = std::string(argv[2]) + (binary ? "/cache.agz" : "/cache.txt.gz"); - auto start=Clock::now(); - auto *out=IntaRNA::newOutputStream(path); - if (binary) source.writeBinary(*out); else source.writeRNAplfold_ED_text(*out); - IntaRNA::deleteOutputStream(out); - auto wrote=Clock::now(); - auto *in=IntaRNA::newInputStream(path); - IntaRNA::AccessibilityFromStream read(rna,100,nullptr,*in, - binary ? IntaRNA::AccessibilityFromStream::IntaRNA_Binary : IntaRNA::AccessibilityFromStream::ED_RNAplfold_Text, 1.0); - IntaRNA::deleteInputStream(in); - auto loaded=Clock::now(); - if (read.getMaxLength()!=source.getMaxLength()) return 3; - size_t differences = 0; - IntaRNA::E_type maxDifference = 0; - for(size_t i=0;i(wrote-start).count() << ',' - << std::chrono::duration(loaded-wrote).count() << ',' - << std::filesystem::file_size(path) << ',' << differences << ',' << maxDifference << std::endl; - } - } +namespace { +using namespace IntaRNA; +namespace bio = boost::iostreams; +using Clock = std::chrono::steady_clock; +constexpr size_t interactionLength = 100; +constexpr size_t foldingWindow = 150; + +std::string makeSequence(size_t length, const char * fasta) +{ + std::string sequence; + sequence.reserve(length); + if (fasta) { + std::ifstream input(fasta); + if (!input) throw std::runtime_error("cannot open FASTA input"); + std::string line; + bool inSequence = false; + while (sequence.size() < length && std::getline(input, line)) { + if (!line.empty() && line[0] == '>') { + if (inSequence) break; + inSequence = true; + continue; + } + for (char c : line) { + if (c != ' ' && c != '\t' && c != '\r' && sequence.size() < length) + sequence += c; + } + } + if (sequence.size() != length) throw std::runtime_error("first FASTA sequence is too short"); + } else { + std::uint32_t state = 245; + for (size_t i = 0; i < length; ++i) { + state = state * 1664525u + 1013904223u; + sequence += "ACGU"[state >> 30]; + } + } + return sequence; +} + +void writeCache(const Accessibility & source, const std::string & path, bool compressed, bool generic) +{ + // Match production stream buffering; only the gzip filter differs. + constexpr std::streamsize bufferSize = 512 * 1024; + bio::filtering_ostream output; + if (compressed) output.push(bio::gzip_compressor(), bufferSize); + output.push(bio::file_descriptor_sink(path, std::ios::out | std::ios::binary), bufferSize); + if (generic) source.Accessibility::writeBinary(output); + else source.writeBinary(output); + output.reset(); // include compressor finalization and file close in timing +} + +std::unique_ptr readCache(const RnaSequence & rna, const std::string & path, bool compressed) +{ + bio::filtering_istream input; + if (compressed) input.push(bio::gzip_decompressor()); + input.push(bio::file_descriptor_source(path, std::ios::in | std::ios::binary)); + auto result = std::make_unique(rna, interactionLength, nullptr, + input, AccessibilityFromStream::IntaRNA_Binary, 1.0); + input.reset(); + return result; +} + +size_t countDifferences(const Accessibility & source, const Accessibility & loaded) +{ + if (loaded.getMaxLength() != source.getMaxLength()) throw std::runtime_error("interaction length mismatch"); + size_t differences = 0; + for (size_t i = 0; i < source.getSequence().size(); ++i) + for (size_t j = i; j < source.getSequence().size() && j-i <= source.getMaxLength(); ++j) + if (source.getED(i,j) != loaded.getED(i,j)) ++differences; + return differences; +} + +void benchmark(size_t length, const std::filesystem::path & directory, const char * fasta) +{ + RnaSequence rna("benchmark", makeSequence(length, fasta)); + VrnaHandler vrna(37, "Turner04", false, false); + AccessibilityVrna folded(rna, interactionLength, nullptr, vrna, foldingWindow); + const auto seed = (directory / "seed.bin").string(); + writeCache(folded, seed, false, false); + auto reloaded = readCache(rna, seed, false); + if (countDifferences(folded, *reloaded)) throw std::runtime_error("initial reload differs"); + std::filesystem::remove(seed); + + // Rotate all eight combinations so one format/method is not always first. + std::cout << "trial,producer,serialization,format,write_seconds,read_seconds,bytes,different_cells\n"; + for (size_t trial = 0; trial < 3; ++trial) { + for (size_t step = 0; step < 8; ++step) { + const size_t variant = (step + 3*trial) % 8; + const bool streamSource = variant & 4, generic = variant & 2, compressed = variant & 1; + const Accessibility & source = streamSource ? static_cast(*reloaded) : folded; + const std::string producer = streamSource ? "stream" : "vrna"; + const auto path = (directory / (producer + (compressed ? ".agz" : ".bin"))).string(); + const auto start = Clock::now(); + writeCache(source, path, compressed, generic); + const auto wrote = Clock::now(); + auto loaded = readCache(rna, path, compressed); + const auto read = Clock::now(); + const auto differences = countDifferences(folded, *loaded); + if (differences) throw std::runtime_error("binary ED values differ"); + std::cout << trial << ',' << producer << ',' << (generic ? "generic" : "direct") << ',' + << (compressed ? "gzip" : "raw") << ',' + << std::chrono::duration(wrote-start).count() << ',' + << std::chrono::duration(read-wrote).count() << ',' + << std::filesystem::file_size(path) << ',' << differences << std::endl; + } + } +} +} // namespace + +int main(int argc, char **argv) +{ + if (argc < 3 || argc > 4) { + std::cerr << "Usage: accessibility-benchmark SEQUENCE_LENGTH OUTPUT_DIRECTORY [FASTA_FILE]\n"; + return 2; + } + el::Loggers::reconfigureAllLoggers(el::ConfigurationType::Enabled, "false"); + try { + const size_t length = std::stoull(argv[1]); + if (length < foldingWindow) throw std::runtime_error("use at least 150 bases"); + benchmark(length, argv[2], argc == 4 ? argv[3] : nullptr); + } catch (const std::exception & error) { + std::cerr << error.what() << '\n'; + return 1; + } } diff --git a/doc/benchmarks/accessibility.md b/doc/benchmarks/accessibility.md index 439a04ea..e5504753 100644 --- a/doc/benchmarks/accessibility.md +++ b/doc/benchmarks/accessibility.md @@ -1,12 +1,21 @@ -# Accessibility I/O benchmark +# Accessibility binary I/O benchmark -This benchmark isolates writing and reading the same accessibility matrix as -RNAplfold-style gzip-compressed ED text and as a binary `.agz` archive. It folds -one deterministic pseudorandom RNA with ViennaRNA (37 C, Turner04, folding window -150, interaction length 100), then alternates formats for three trials. Folding -and verification are outside the timed sections. Every loaded ED cell, including -the dangling-end band, is compared to the computed matrix. Binary input must -match exactly; text conversion differences are reported. +This benchmark measures two independent choices on identical ED data: + +- Export stored matrix rows directly, or gather each row through `getED()` using + the generic `Accessibility::writeBinary()` implementation. +- Write/read an uncompressed Boost binary archive, or add the production gzip + filter used for `.agz` files. + +It tests both `AccessibilityVrna` and a previously loaded +`AccessibilityFromStream`. Direct export passes non-owning matrix row views to +Boost without copying ED rows. Constrained accessibility data and other +implementations retain the generic path to preserve their `getED()` semantics. +All reads use the optimized loader, which fills retained rows in place and +validates discarded tails when the requested interaction length is smaller. +The version 1 archive layout is unchanged. + +## Reproduction Build and install IntaRNA, then compile against that installation: @@ -15,31 +24,91 @@ export PKG_CONFIG_PATH=/path/to/intarna-install/lib/pkgconfig c++ -std=c++23 -O3 doc/benchmarks/accessibility.cpp \ $(pkg-config --cflags --libs IntaRNA) -leasylogging \ -o accessibility-benchmark -mkdir -p accessibility-benchmark-output -./accessibility-benchmark 100000 accessibility-benchmark-output +export OMP_NUM_THREADS=1 OPENBLAS_NUM_THREADS=1 +mkdir -p accessibility-benchmark-output/random accessibility-benchmark-output/ecoli +./accessibility-benchmark 100000 accessibility-benchmark-output/random > random.csv +curl -fL --retry 2 \ + 'https://eutils.ncbi.nlm.nih.gov/entrez/eutils/efetch.fcgi?db=nuccore&id=NC_000913.3&rettype=fasta&retmode=text&seq_start=1&seq_stop=100000' \ + -o ecoli-100k.fa +./accessibility-benchmark 100000 accessibility-benchmark-output/ecoli ecoli-100k.fa > ecoli.csv ``` Use the same compiler, dependency include/library paths and runtime library paths as the IntaRNA build. With dependencies outside system paths, additional -`-L` and `-Wl,-rpath` flags may be necessary. Output columns are -`trial,format,write_seconds,read_seconds,compressed_bytes,different_cells,max_ED_difference`. -ED differences are measured in internal units (hundredths of kcal/mol). Files are overwritten -between trials; the output directory must already exist. These timings include -compression and stream closing, but not an `fsync` to durable storage. Results -depend on sequence, bandwidth, compressor, filesystem and machine load. - -## Example measurement (2026-09-29) - -Linux x86-64, Ryzen 5 7530U, GCC 16.2.0 release build, Boost 1.85.0, -ViennaRNA 2.7.2, 100,000 bases, one folding/I/O thread. Medians of the three -alternating trials above (other build work was running on the machine): - -| Format | Write (s) | Read (s) | Compressed bytes | ED cells differing from source | Maximum ED difference | -| --- | ---: | ---: | ---: | ---: | ---: | -| ED text + gzip | 8.726 | 3.365 | 20,891,786 | 551,861 | 1 (0.01 kcal/mol) | -| Binary `.agz` | 3.582 | 0.349 | 16,620,982 | 0 | 0 | - -In this local run, binary output was about 2.4 times faster to write, 9.6 times -faster to read, and 20% smaller. Text differences arise in existing decimal -conversion; the benchmark does not modify that behavior. These figures describe -I/O of the same computed data, not an end-to-end genome-screen speedup. +`-L` and `-Wl,-rpath` flags may be necessary. The output directory must exist; +archive files are overwritten between trials. + +Without a FASTA argument, the sequence is deterministic pseudorandom RNA (the +32-bit generator and seed 245 are in the source). With a FASTA argument, the +benchmark uses the requested number of bases from its first record. The second +measurement uses bases 1-100,000 of +[E. coli K-12 MG1655, NC_000913.3](https://www.ncbi.nlm.nih.gov/nuccore/NC_000913.3). +The downloaded FASTA has SHA-256 +`7ad0766d4eb40e247ad5c2a9de2518738a39c8d261d27efb3bf1218c7525d0bc`. + +Each dataset is folded once with ViennaRNA at 37 C, Turner04, folding window +150 and maximum interaction length 100. Three trials rotate through the eight +producer/method/format combinations. Output columns are +`trial,producer,serialization,format,write_seconds,read_seconds,bytes,different_cells`. +Every reloaded ED cell, including the dangling-end band, is compared exactly to +the original folded matrix; any difference aborts the run. Folding, the initial +stream-source load and verification are outside the timed sections. + +Both formats use Boost file-descriptor streams and the production buffering +(512 KiB for output, default buffering for input); only the gzip filter differs. +Timings include compressor finalization and stream close, without `fsync` or +cache eviction. Reads therefore benefit from the filesystem cache. These are +local I/O wall times, not durable-storage throughput or end-to-end genome-screen +speedups. The short raw writes vary across trials; the CSV files retain that +variation rather than reporting only a best case. + +## Measurements (2026-09-30) + +Linux x86-64, Ryzen 5 7530U, GCC 16.2.0 release build, native `std::mdspan`, Boost +1.85.0 and ViennaRNA 2.7.2. One folding/I/O thread, datasets run sequentially, +with no concurrent IntaRNA builds/tests. Each input has 100,000 bases. The +following values are medians of three trials; all 48 reloads matched every ED +cell exactly. Per-trial results: [synthetic RNA](accessibility-random.csv) and +[E. coli](accessibility-ecoli.csv). + +### Compression after direct matrix export + +| Dataset | Producer | Format | Write (s) | Read (s) | Bytes | +| --- | --- | --- | ---: | ---: | ---: | +| Synthetic | Vrna | raw | 0.105 | 0.022 | 40,479,873 | +| Synthetic | Vrna | gzip | 2.796 | 0.248 | 16,620,982 | +| Synthetic | FromStream | raw | 0.102 | 0.021 | 40,479,873 | +| Synthetic | FromStream | gzip | 2.762 | 0.257 | 16,620,982 | +| E. coli | Vrna | raw | 0.108 | 0.024 | 40,479,873 | +| E. coli | Vrna | gzip | 3.019 | 0.271 | 16,850,736 | +| E. coli | FromStream | raw | 0.109 | 0.024 | 40,479,873 | +| E. coli | FromStream | gzip | 2.997 | 0.268 | 16,850,736 | + +Gzip reduces the archive size by 59% on synthetic RNA and 58% on E. coli; raw +files are about 2.4 times as large. It adds roughly 2.7-2.9 seconds per write and +0.23-0.25 seconds per read in these measurements. Compression is therefore +retained for `.agz`: its CPU cost is substantial, but the uncompressed output is +also substantially larger. This follows the review criterion of removing gzip +only if the size saving is small. C++ callers can already use uncompressed +binary streams; this comparison does not introduce a new CLI format. + +### Direct versus generic export + +All times below are write medians. Both paths generate the same version 1 +payload. Reads do not depend on which writer produced that payload; their +per-trial timings are included in the CSV files. + +| Dataset | Producer | Raw direct (s) | Raw generic (s) | Gzip direct (s) | Gzip generic (s) | +| --- | --- | ---: | ---: | ---: | ---: | +| Synthetic | Vrna | 0.105 | 0.157 | 2.796 | 2.853 | +| Synthetic | FromStream | 0.102 | 0.150 | 2.762 | 3.054 | +| E. coli | Vrna | 0.108 | 0.384 | 3.019 | 3.016 | +| E. coli | FromStream | 0.109 | 0.162 | 2.997 | 3.052 | + +Direct export removes the row scratch buffer and per-cell virtual `getED()` +lookups for unconstrained matrix-backed producers. The raw measurements expose +that benefit; gzip dominates compressed write time, where some differences are +within run-to-run variation. API tests independently check zero `getED()` calls +on direct export, byte-identical generic/direct archives, constraint masking, +and exact round trips. An archive from the original version 1 implementation +was also loaded and re-exported byte-for-byte identically. diff --git a/src/IntaRNA/Accessibility.cpp b/src/IntaRNA/Accessibility.cpp index b5b5ce87..e054879c 100644 --- a/src/IntaRNA/Accessibility.cpp +++ b/src/IntaRNA/Accessibility.cpp @@ -12,9 +12,15 @@ const E_type Accessibility::ED_UPPER_BOUND = (E_type) E_INF; void Accessibility::writeBinary( std::ostream & out ) const +{ + writeBinary(out, nullptr); +} + +void +Accessibility::writeBinary( std::ostream & out, const UpperBandedMatrix * matrix ) const { boost::archive::binary_oarchive archive(out); - archive << AccessibilityArchive(*this); + archive << AccessibilityArchive(*this, matrix); out.flush(); if (!out) throw std::runtime_error("Accessibility::writeBinary: output failure"); } diff --git a/src/IntaRNA/Accessibility.h b/src/IntaRNA/Accessibility.h index 652c7da5..8d4cd478 100644 --- a/src/IntaRNA/Accessibility.h +++ b/src/IntaRNA/Accessibility.h @@ -12,6 +12,8 @@ namespace IntaRNA { +template class UpperBandedMatrix; + /** * Abstract interface that represents accessibility data for a given RNA * sequence. @@ -115,7 +117,7 @@ class Accessibility { * @param out binary output stream * @throw std::exception on invalid ED values or output failure */ - void writeBinary( std::ostream & out ) const; + virtual void writeBinary( std::ostream & out ) const; /** * Prints the accessibility values to stream as upper triangular matrix @@ -164,6 +166,15 @@ class Accessibility { protected: + /** + * Write an archive using stored ED rows when available. Non-empty + * constraints use getED() to preserve any masking applied by the subclass. + * @param out binary output stream + * @param matrix non-owning matrix with the same unconstrained ED values as + * getED(), or nullptr to gather values through getED(); used only during this call + */ + void writeBinary( std::ostream & out, const UpperBandedMatrix * matrix ) const; + //! the RNA sequence the accessibilities correspond to const RnaSequence & seq; diff --git a/src/IntaRNA/AccessibilityArchive.h b/src/IntaRNA/AccessibilityArchive.h index d99b45d0..aea096db 100644 --- a/src/IntaRNA/AccessibilityArchive.h +++ b/src/IntaRNA/AccessibilityArchive.h @@ -20,7 +20,9 @@ namespace IntaRNA { * Serialization view of the logical accessibility matrix, independent of its * physical storage. Version 1 stores exact internal ED integers, in rows of * increasing start position and interval length, including maxLength+1 for - * dangling ends. Only valid cells are stored; scratch memory is O(maxLength). + * dangling ends. Only valid cells are stored. Matrix-backed saves and full-band + * loads use row views directly; generic saves and discarded input tails need + * O(maxLength) scratch memory. The version 1 archive layout is unchanged. * * The input sequence and requested band bound allocations on load. Native * Boost binary archives require a compatible architecture and Boost version. @@ -30,8 +32,13 @@ class AccessibilityArchive { public: typedef UpperBandedMatrix EdMatrix; - /** Create a read-only serialization view of any accessibility implementation. */ - explicit AccessibilityArchive( const Accessibility & source ); + /** + * Create a read-only serialization view. + * @param source accessibility data, alive until serialization finishes + * @param matrix optional non-owning view of the source's unconstrained ED + * values; nullptr or non-empty constraints select the generic getED() path + */ + explicit AccessibilityArchive( const Accessibility & source, const EdMatrix * matrix = nullptr ); /** Create a loading view; target is resized to the requested available band. */ AccessibilityArchive( const RnaSequence & sequence, size_t maxLength, EdMatrix & target ); @@ -44,6 +51,8 @@ class AccessibilityArchive { const RnaSequence & sequence; size_t maxLength; const Accessibility * source; + //! Optional matrix storage for direct output; never owned or modified. + const EdMatrix * sourceMatrix; EdMatrix * target; template void save( Archive & archive, unsigned int version ) const; @@ -52,13 +61,14 @@ class AccessibilityArchive { }; inline -AccessibilityArchive::AccessibilityArchive( const Accessibility & source ) - : sequence(source.getSequence()), maxLength(source.getMaxLength()), source(&source), target(nullptr) +AccessibilityArchive::AccessibilityArchive( const Accessibility & source, const EdMatrix * matrix ) + : sequence(source.getSequence()), maxLength(source.getMaxLength()), source(&source) + , sourceMatrix(source.getAccConstraint().isEmpty() ? matrix : nullptr), target(nullptr) {} inline AccessibilityArchive::AccessibilityArchive( const RnaSequence & sequence, size_t maxLength, EdMatrix & target ) - : sequence(sequence), maxLength(maxLength), source(nullptr), target(&target) + : sequence(sequence), maxLength(maxLength), source(nullptr), sourceMatrix(nullptr), target(&target) {} inline size_t @@ -78,15 +88,25 @@ AccessibilityArchive::save( Archive & archive, unsigned int ) const archive & boost::serialization::make_array(sequence.asString().data(), sequence.size()); // Avoid maxLength+1 overflow when the entire sequence is covered. const size_t width = maxLength < sequence.size() ? maxLength+1 : sequence.size(); - std::vector row(width); + if (sourceMatrix && (sourceMatrix->size1() != sequence.size() || sourceMatrix->size2() != sequence.size())) + throw std::runtime_error("Accessibility archive: source matrix shape mismatch"); + std::vector scratch(sourceMatrix ? 0 : width); for (size_t i = 0; i < sequence.size(); ++i) { const size_t count = std::min(width, sequence.size()-i); - for (size_t k = 0; k < count; ++k) { - row[k] = source->getED(i, i+k); - if (row[k] < 0 || row[k] > infinity) - throw std::runtime_error("Accessibility archive: invalid ED value"); + std::span row; + if (sourceMatrix) { + row = sourceMatrix->row(i); + if (row.size() < count) + throw std::runtime_error("Accessibility archive: source matrix band too narrow"); + row = row.first(count); + } else { + for (size_t k = 0; k < count; ++k) scratch[k] = source->getED(i, i+k); + row = std::span(scratch.data(), count); } - archive & boost::serialization::make_array(row.data(), count); + for (const E_type value : row) + if (value < 0 || value > infinity) + throw std::runtime_error("Accessibility archive: invalid ED value"); + archive & boost::serialization::make_array(row.data(), row.size()); } } @@ -115,15 +135,19 @@ AccessibilityArchive::load( Archive & archive, unsigned int ) if (sequence.size() > std::numeric_limits::max() / width / sizeof(E_type)) throw std::runtime_error("Accessibility archive: matrix size overflow"); target->resize(sequence.size(), sequence.size(), 0, width-1, false); - std::vector row(storedWidth); + // Read retained cells straight into their final storage. Only the discarded + // suffix of a wider input band needs scratch space, and it is still validated. + std::vector discarded(storedWidth-width); for (size_t i = 0; i < sequence.size(); ++i) { const size_t count = std::min(storedWidth, sequence.size()-i); - archive & boost::serialization::make_array(row.data(), count); - for (size_t k = 0; k < count; ++k) { - if (row[k] < 0 || row[k] > infinity) - throw std::runtime_error("Accessibility archive: invalid ED value"); - if (k < width) (*target)(i, i+k) = row[k]; - } + auto row = target->row(i); + archive & boost::serialization::make_array(row.data(), row.size()); + const size_t tailSize = count-row.size(); + if (tailSize) archive & boost::serialization::make_array(discarded.data(), tailSize); + for (const auto values : {std::span(row), std::span(discarded.data(), tailSize)}) + for (const E_type value : values) + if (value < 0 || value > infinity) + throw std::runtime_error("Accessibility archive: invalid ED value"); } } diff --git a/src/IntaRNA/AccessibilityFromStream.h b/src/IntaRNA/AccessibilityFromStream.h index 335eb87d..79493c0e 100644 --- a/src/IntaRNA/AccessibilityFromStream.h +++ b/src/IntaRNA/AccessibilityFromStream.h @@ -83,6 +83,14 @@ class AccessibilityFromStream: public Accessibility getMaxLength() const; + /** + * Write the stored ED matrix directly, without copying rows or calling + * getED() for unconstrained data. Constraints retain the generic getED() + * path. Subclasses changing getED() semantics must also override this method. + * @param out binary output stream; compression is supplied by the stream + */ + void writeBinary( std::ostream & out ) const override; + protected: //! type for the ED value matrix (upper triangular matrix banded by maxLength) @@ -135,6 +143,13 @@ class AccessibilityFromStream: public Accessibility ///////////////////////////////////////////////////////////////////////// +inline +void +AccessibilityFromStream::writeBinary( std::ostream & out ) const +{ + Accessibility::writeBinary(out, &edValues); +} + inline E_type AccessibilityFromStream:: diff --git a/src/IntaRNA/AccessibilityVrna.h b/src/IntaRNA/AccessibilityVrna.h index 35448e21..d828ec3e 100644 --- a/src/IntaRNA/AccessibilityVrna.h +++ b/src/IntaRNA/AccessibilityVrna.h @@ -83,6 +83,14 @@ class AccessibilityVrna : public Accessibility { , const Accessibility & acc ); + /** + * Write the stored ED matrix directly, without copying rows or calling + * getED() for unconstrained data. Constraints retain the generic getED() + * path. Subclasses changing getED() semantics must also override this method. + * @param out binary output stream; compression is supplied by the stream + */ + void writeBinary( std::ostream & out ) const override; + protected: //! type for the ED value matrix (upper triangular matrix banded by maxLength) @@ -141,6 +149,13 @@ class AccessibilityVrna : public Accessibility { ///////////////////////////////////////////////////////////////////////////// ///////////////////////////////////////////////////////////////////////////// +inline +void +AccessibilityVrna::writeBinary( std::ostream & out ) const +{ + Accessibility::writeBinary(out, &edValues); +} + inline E_type AccessibilityVrna:: diff --git a/src/IntaRNA/Matrix.h b/src/IntaRNA/Matrix.h index 94e9472d..c5390f82 100644 --- a/src/IntaRNA/Matrix.h +++ b/src/IntaRNA/Matrix.h @@ -5,6 +5,7 @@ #include #include #include +#include #include #include #include @@ -499,6 +500,21 @@ class UpperBandedMatrix { * @return read-only reference to the cell */ const T &operator()(std::size_t i, std::size_t j) const; + /** + * View the contiguous stored cells (i,i), (i,i+1), ... in one row. + * Excludes structural zeros and row-end padding. The view aliases this + * matrix and must not outlive its storage or be retained across resize, + * move, swap or assignment. + * @param i zero-based row, less than size1() + * @return mutable stored cells; empty when i >= size2() or the band is empty + */ + std::span row(std::size_t i); + /** + * Read-only view of one stored row, with the same bounds and lifetime as row(). + * @param i zero-based row, less than size1() + * @return read-only stored cells, excluding padding and structural zeros + */ + std::span row(std::size_t i) const; /** * Reset stored cells to default values without changing the shape. */ @@ -586,6 +602,24 @@ const T &UpperBandedMatrix::operator()(std::size_t i, std::size_t j) const return i > j || j - i >= band.size2() ? zero : band(i, j - i); } +template +inline +std::span UpperBandedMatrix::row(std::size_t i) +{ + assert(i < size1()); + const auto count = i < columns ? std::min(band.size2(), columns-i) : 0; + return count == 0 ? std::span{} : std::span{&band(i, 0), count}; +} + +template +inline +std::span UpperBandedMatrix::row(std::size_t i) const +{ + assert(i < size1()); + const auto count = i < columns ? std::min(band.size2(), columns-i) : 0; + return count == 0 ? std::span{} : std::span{&band(i, 0), count}; +} + template inline void UpperBandedMatrix::clear() diff --git a/tests/AccessibilityBinary_test.cpp b/tests/AccessibilityBinary_test.cpp index 334d9536..77ebce96 100644 --- a/tests/AccessibilityBinary_test.cpp +++ b/tests/AccessibilityBinary_test.cpp @@ -13,6 +13,21 @@ using namespace IntaRNA; namespace { +// Count interface lookups without changing the underlying ED semantics. +template class CountedAccessibility : public Base { +public: + using Base::Base; + mutable size_t lookups = 0; + E_type getED(size_t from, size_t to) const override; +}; + +template E_type +CountedAccessibility::getED(size_t from, size_t to) const +{ + ++lookups; + return Base::getED(from, to); +} + void checkBinaryRoundTrip( const Accessibility & source, size_t requestedLength ) { std::stringstream bytes; @@ -121,6 +136,8 @@ TEST_CASE("Binary accessibility rejects invalid archives", "[AccessibilityBinary rejects(replaceField(valid, header+16, std::uint64_t(-1))); rejects(replaceField(valid, header+24, E_type(1))); rejects(replaceField(valid, valid.size()-sizeof(E_type), E_type(-1))); + // Invalid data in the discarded suffix of the first row must still fail. + rejects(replaceField(valid, header+28+rna.size()+5*sizeof(E_type), E_type(-1))); rejects(replaceField(valid, valid.size()-sizeof(E_type), Accessibility::ED_UPPER_BOUND+1)); } SECTION("different sequence of the same length") { @@ -135,3 +152,48 @@ TEST_CASE("Binary accessibility rejects invalid archives", "[AccessibilityBinary REQUIRE_THROWS(source.writeBinary(output)); } } + +TEST_CASE("Stored accessibility matrices serialize without interface lookups", "[AccessibilityBinary]") +{ +#include "testEasyLoggingSetup.icc" + RnaSequence rna("test", "GGGGAAAACCCCUAGC"); + VrnaHandler vrna(37, "Turner04", false, false); + for (size_t length : {size_t(1), size_t(5), rna.size()}) { + CountedAccessibility folded(rna, length, nullptr, vrna, rna.size()); + std::stringstream direct, generic; + const Accessibility & polymorphic = folded; + polymorphic.writeBinary(direct); + REQUIRE(folded.lookups == 0); + folded.Accessibility::writeBinary(generic); + REQUIRE(folded.lookups > 0); + // Same version 1 payload, independent of matrix padding and dispatch. + REQUIRE(direct.str() == generic.str()); + for (size_t requested : {size_t(0), size_t(1), size_t(3)}) { + std::istringstream input(direct.str()); + CountedAccessibility loaded(rna, requested, nullptr, + input, AccessibilityFromStream::IntaRNA_Binary, 1.0); + std::stringstream stored, gathered; + static_cast(loaded).writeBinary(stored); + REQUIRE(loaded.lookups == 0); + loaded.Accessibility::writeBinary(gathered); + REQUIRE(loaded.lookups > 0); + REQUIRE(stored.str() == gathered.str()); + } + } +} + +TEST_CASE("Constrained matrix output retains getED masking", "[AccessibilityBinary]") +{ +#include "testEasyLoggingSetup.icc" + // The short-sequence matrix contains zeros even at inaccessible positions. + RnaSequence rna("short", "AC"); + AccessibilityConstraint constraint(rna, "pb", 0, "", "", ""); + VrnaHandler vrna(37, "Turner04", false, false); + CountedAccessibility folded(rna, 1, &constraint, vrna, 0); + std::stringstream direct, generic; + folded.writeBinary(direct); + REQUIRE(folded.lookups > 0); + folded.Accessibility::writeBinary(generic); + REQUIRE(direct.str() == generic.str()); + checkBinaryRoundTrip(folded, 0); +} diff --git a/tests/Matrix_test.cpp b/tests/Matrix_test.cpp index ab9bce12..de0085bb 100644 --- a/tests/Matrix_test.cpp +++ b/tests/Matrix_test.cpp @@ -158,3 +158,24 @@ TEST_CASE("Oversized matrix dimensions fail before allocation", "[Matrix]") { REQUIRE(matrix.size1() == 1); REQUIRE(matrix(0, 0) == 42); } + +TEST_CASE("Upper band row views alias valid cells without padding", "[Matrix]") { + for (std::size_t rows : {0u, 2u, 5u}) { + for (std::size_t columns : {0u, 2u, 5u}) { + for (std::size_t upper : {0u, 1u, 5u}) { + UpperBandedMatrix matrix(rows, columns, 0, upper); + for (std::size_t i = 0; i < rows; ++i) { + auto row = matrix.row(i); + const auto constRow = std::as_const(matrix).row(i); + REQUIRE(row.size() == (i < columns ? std::min(upper+1, columns-i) : 0)); + REQUIRE(constRow.size() == row.size()); + for (std::size_t k = 0; k < row.size(); ++k) { + row[k] = static_cast(100*i+k); + REQUIRE(&constRow[k] == &matrix(i, i+k)); + REQUIRE(constRow[k] == static_cast(100*i+k)); + } + } + } + } + } +} From 793f076e59e3ca16edcdfd8862a3727a2896614d Mon Sep 17 00:00:00 2001 From: Martin Raden Date: Wed, 30 Sep 2026 14:50:21 +0200 Subject: [PATCH 4/4] cleanup Makefiles --- doc/Makefile.am | 4 ---- src/IntaRNA/Makefile.am | 3 +-- 2 files changed, 1 insertion(+), 6 deletions(-) diff --git a/doc/Makefile.am b/doc/Makefile.am index 2b004d3d..322e9d93 100644 --- a/doc/Makefile.am +++ b/doc/Makefile.am @@ -4,10 +4,6 @@ ################################################################ EXTRA_DIST = \ - benchmarks/accessibility.cpp \ - benchmarks/accessibility.md \ - benchmarks/accessibility-random.csv \ - benchmarks/accessibility-ecoli.csv \ conda.txt \ doxygen.cfg \ latex-deps/adjcalc.sty \ diff --git a/src/IntaRNA/Makefile.am b/src/IntaRNA/Makefile.am index e3cb1d5f..22291e61 100644 --- a/src/IntaRNA/Makefile.am +++ b/src/IntaRNA/Makefile.am @@ -9,8 +9,6 @@ AUTOMAKE_OPTIONS = std-options subdir-objects AM_DEFAULT_SOURCE_EXT = .cpp -noinst_HEADERS = AccessibilityArchive.h - ############################################################################### # THE INTARNA LIBRARY ############################################################################### @@ -30,6 +28,7 @@ libIntaRNA_a_HEADERS = \ intarna_config.h \ Matrix.h \ Accessibility.h \ + AccessibilityArchive.h \ AccessibilityConstraint.h \ AccessibilityDisabled.h \ AccessibilityFromStream.h \