Add compressed binary accessibility caches (.agz)

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.*
@@ -23,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
@@ -57,7 +61,43 @@
 ################################################################################  
 ################################################################################
 
+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
+   + 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
  * configure.ac, IntaRNA/Matrix.h, IntaRNA/intarna_config.h.in:
    + prefer native mdspan, falling back to bundled Kokkos with a public
      INTARNA_USE_STD_MDSPAN define; allow explicit --with-mdspan selection

README.md

@@ -190,6 +190,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, ..)
@@ -2112,6 +2114,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
@@ -2135,6 +2138,47 @@ 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.
+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
 
 If you have precomputed data, e.g. the file `plfold_lunp` with unpaired probabilities

configure.ac

@@ -381,6 +381,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 <sstream>
+#include <boost/archive/binary_iarchive.hpp>
+#include <boost/archive/binary_oarchive.hpp>]],
+		[[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([], [])
 
@@ -404,7 +427,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"
 ])
 
 

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

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

doc/benchmarks/accessibility.cpp

@@ -0,0 +1,144 @@
+// Compare generic/direct matrix export and raw/gzip binary I/O on identical data.
+// Usage: accessibility-benchmark SEQUENCE_LENGTH OUTPUT_DIRECTORY [FASTA_FILE]
+#include <IntaRNA/AccessibilityVrna.h>
+#include <IntaRNA/AccessibilityFromStream.h>
+#include <boost/iostreams/device/file_descriptor.hpp>
+#include <boost/iostreams/filter/gzip.hpp>
+#include <boost/iostreams/filtering_stream.hpp>
+#include <algorithm>
+#include <array>
+#include <chrono>
+#include <cstdint>
+#include <filesystem>
+#include <fstream>
+#include <iostream>
+#include <memory>
+#include <stdexcept>
+#include <string>
+
+INITIALIZE_EASYLOGGINGPP
+
+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<AccessibilityFromStream> 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<AccessibilityFromStream>(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<const Accessibility &>(*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<double>(wrote-start).count() << ','
+				<< std::chrono::duration<double>(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;
+	}
+}