From 2abfe0fdb3d61fb1241f56f0bffdcc643d68de25 Mon Sep 17 00:00:00 2001 From: "Peter J. Jones" Date: Fri, 17 Jul 2026 14:16:11 +0200 Subject: [PATCH 01/31] Don't project unneeded columns --- include/mzpeak/data/array_index.h | 4 ++++ src/data/array_index.cpp | 23 +++++++++++++++++++++++ src/data/signals.cpp | 2 ++ 3 files changed, 29 insertions(+) diff --git a/include/mzpeak/data/array_index.h b/include/mzpeak/data/array_index.h index a1ca2e6..8d800a6 100644 --- a/include/mzpeak/data/array_index.h +++ b/include/mzpeak/data/array_index.h @@ -88,6 +88,10 @@ class ArrayIndex final { /// denoted as a CURIE from the PSI-MS controlled vocabulary. Some /// values are only usable with the chunked layout. std::optional transform = {}; + + /// Return `true` if this entry needs to be projected in a query + /// in order to properly decode the dimension it represents. + bool needed_for_decoding() const; }; /** diff --git a/src/data/array_index.cpp b/src/data/array_index.cpp index 6324057..ea501a7 100644 --- a/src/data/array_index.cpp +++ b/src/data/array_index.cpp @@ -16,6 +16,29 @@ directory of this repository. namespace MzPeak::Data { +/******************************************************************************/ +bool ArrayIndex::Entry::needed_for_decoding() const +{ + switch (buffer_format) { + case MzPeak::Schema::BufferFormat::Point: + return true; + case MzPeak::Schema::BufferFormat::ChunkStart: + return false; + case MzPeak::Schema::BufferFormat::ChunkEnd: + return false; + case MzPeak::Schema::BufferFormat::ChunkValues: + return true; + case MzPeak::Schema::BufferFormat::ChunkEncoding: + return true; + case MzPeak::Schema::BufferFormat::ChunkSecondary: + return true; + case MzPeak::Schema::BufferFormat::ChunkTransform: + return true; + } + + std::unreachable(); +} + /******************************************************************************/ bool ArrayIndex::Dimension::needs_delta_model() const { diff --git a/src/data/signals.cpp b/src/data/signals.cpp index 6144bc1..325413a 100644 --- a/src/data/signals.cpp +++ b/src/data/signals.cpp @@ -136,6 +136,8 @@ Signals::select(const std::vector& projection, for (const auto& dim : projection) { for (const auto& entry : dim.entries) { + if (!entry.needed_for_decoding()) continue; + auto field = impl_->array_index_->entry_column(*impl_->parquet_->groups(), entry); if (!field.has_value()) { From e995d31b1406b6b9fe1c4d82f0ef0c9947a74b13 Mon Sep 17 00:00:00 2001 From: "Peter J. Jones" Date: Fri, 17 Jul 2026 15:27:03 +0200 Subject: [PATCH 02/31] Additional types and functions needed for chunked layout decoding --- include/mzpeak/data/array_index.h | 38 +++++++++++ include/mzpeak/data/encoding.h | 54 +++++++++------- include/mzpeak/exception.h | 32 +++++++++ src/data/array_index.cpp | 104 ++++++++++++++++++++++++++++++ 4 files changed, 205 insertions(+), 23 deletions(-) diff --git a/include/mzpeak/data/array_index.h b/include/mzpeak/data/array_index.h index 8d800a6..cceed96 100644 --- a/include/mzpeak/data/array_index.h +++ b/include/mzpeak/data/array_index.h @@ -30,6 +30,23 @@ using namespace MzPeak::Schema; */ class ArrayIndex final { public: + /** + * How dimensions in the Parquet file are encoded. + */ + enum class Layout { + /// Point layout uses one column per dimension. + Point, + + /// Chunked layout differentiates between the main axis and + /// secondary axes. Finding the correct column requires the + /// `buffer_format` member of the `Entry`. + Chunked, + + /// The Parquet file uses an unknown layout and must be decoded + /// manually. + Unknown, + }; + /** * A type to describe each entry in the index. */ @@ -92,6 +109,10 @@ class ArrayIndex final { /// Return `true` if this entry needs to be projected in a query /// in order to properly decode the dimension it represents. bool needed_for_decoding() const; + + /// Return `true` if this entry stores values for the associated + /// dimension. + bool is_value_entry() const; }; /** @@ -113,11 +134,20 @@ class ArrayIndex final { /// The index entries that make up this dimension. std::vector entries; + /// Is this dimension on the main axis? + bool is_main_axis() const; + /// Does this dimension need a delta model for decoding? bool needs_delta_model() const; /// Return an associated Util::Type or throw an exception. Util::Type type_or_throw() const; + + /// Return the entry that holds the (possibly encoded) values for + /// this dimension. Throws an exception of the dimension is + /// malformed and thus doesn't include any of the expected + /// entries. + const Entry& values_entry() const; }; /// Default constructor. @@ -139,6 +169,11 @@ class ArrayIndex final { */ const std::string& prefix() const; + /** + * Return the file layout. + */ + Layout layout() const; + /** * Get a list of entry definitions. */ @@ -171,6 +206,9 @@ class ArrayIndex final { // Root node. std::string prefix_ = "point"; + // Layout. + Layout layout_; + // Entries; std::vector entries_; diff --git a/include/mzpeak/data/encoding.h b/include/mzpeak/data/encoding.h index bf914e8..942ba1b 100644 --- a/include/mzpeak/data/encoding.h +++ b/include/mzpeak/data/encoding.h @@ -16,7 +16,6 @@ top-level directory of this repository. #include "mzpeak/data/null_marking.h" #include "mzpeak/data/signals.h" #include "mzpeak/exception.h" -#include "mzpeak/schema/group.h" #include "mzpeak/schema/psi/data_type.h" #include "mzpeak/util/slice.h" #include "mzpeak/util/types.h" @@ -57,7 +56,8 @@ template class Decoder { void decode(const ArrayIndex::Dimension&, std::vector&) const; template - void point(const Schema::Column&, const N& null_decoder, std::vector&) const; + void + point(const ArrayIndex::Dimension&, const N& null_decoder, std::vector&) const; template void remap(const ArrayIndex::Dimension& dim, std::vector& v) const; @@ -129,42 +129,50 @@ template template void Decoder::decode(const ArrayIndex::Dimension& dim, std::vector& v) const { - const auto& entries = dim.entries; - - if (entries.empty()) { - std::string msg("unable to decode dimension, wrong encoding: "); - throw ParquetError(msg + dim.name); - } else if (entries.size() == 1 && - entries[0].buffer_format == Schema::BufferFormat::Point) { - - auto field = - signals_->array_index()->entry_column(*signals_->groups(), entries[0]); - - if (!field.has_value()) { - throw ParquetError("unable to decode dimension, not in schema: " + dim.name); - } - + switch (signals_->array_index()->layout()) { + case ArrayIndex::Layout::Point: if (dim.needs_delta_model()) { using N = NullMarking::Decoder; - point(field.value(), N{delta_estimator_}, v); + point(dim, N{delta_estimator_}, v); } else { using N = Util::Decoders::NullToZero; - point(field.value(), N{}, v); + point(dim, N{}, v); } - } else { + break; + + case ArrayIndex::Layout::Chunked: throw("not implemented"); - // return decode_chunked(arrays); + break; + + case ArrayIndex::Layout::Unknown: + throw UnknownLayoutError("cannot decode dimension: " + dim.name); } } /******************************************************************************/ template template -void Decoder::point(const Schema::Column& col, +void Decoder::point(const ArrayIndex::Dimension& dim, const N& null_decoder, std::vector& v) const { - slice_->array(col, v, Util::Decoders::Scalar, N>(null_decoder)); + const auto& primary_entry = dim.values_entry(); + + if (primary_entry.buffer_format != Schema::BufferFormat::Point) { + std::string msg("file uses point layout, but " + dim.name); + msg += " is not using the point buffer_format"; + throw InvalidFormatError(msg); + } + + auto col = + signals_->array_index()->entry_column(*signals_->groups(), primary_entry); + + if (!col.has_value()) { + throw ParquetError("unable to decode dimension, not in schema: " + dim.name); + } + + slice_->array(col.value(), v, + Util::Decoders::Scalar, N>(null_decoder)); } } // namespace MzPeak::Data::Encoding diff --git a/include/mzpeak/exception.h b/include/mzpeak/exception.h index 1ba4df9..7cad95c 100644 --- a/include/mzpeak/exception.h +++ b/include/mzpeak/exception.h @@ -57,6 +57,38 @@ class ParquetError final : public Exception { ~ParquetError() = default; }; +/** + * The mzPeak file is malformed and the error could not be recovered + * from. + */ +class InvalidFormatError final : public Exception { +public: + /// Constructor. + InvalidFormatError(const std::string& msg) + : Exception(msg) + { + } + + /// Destructor. + ~InvalidFormatError() = default; +}; + +/** + * Automatic decoding of signal data is only supported for standard + * file layouts such as point and chunked. + */ +class UnknownLayoutError final : public Exception { +public: + /// Constructor. + UnknownLayoutError(const std::string& msg) + : Exception(msg) + { + } + + /// Destructor. + ~UnknownLayoutError() = default; +}; + /** * Attempt to access an invalid iterator. */ diff --git a/src/data/array_index.cpp b/src/data/array_index.cpp index ea501a7..340a199 100644 --- a/src/data/array_index.cpp +++ b/src/data/array_index.cpp @@ -16,6 +16,18 @@ directory of this repository. namespace MzPeak::Data { +/******************************************************************************/ +ArrayIndex::Layout group_name_to_layout(const std::string& name) +{ + if (name == "point") { + return ArrayIndex::Layout::Point; + } else if (name == "chunk") { + return ArrayIndex::Layout::Chunked; + } else { + return ArrayIndex::Layout::Unknown; + } +} + /******************************************************************************/ bool ArrayIndex::Entry::needed_for_decoding() const { @@ -39,6 +51,54 @@ bool ArrayIndex::Entry::needed_for_decoding() const std::unreachable(); } +/******************************************************************************/ +bool ArrayIndex::Entry::is_value_entry() const +{ + switch (buffer_format) { + case MzPeak::Schema::BufferFormat::Point: + return true; + case MzPeak::Schema::BufferFormat::ChunkStart: + return false; + case MzPeak::Schema::BufferFormat::ChunkEnd: + return false; + case MzPeak::Schema::BufferFormat::ChunkValues: + return true; + case MzPeak::Schema::BufferFormat::ChunkEncoding: + return false; + case MzPeak::Schema::BufferFormat::ChunkSecondary: + return true; + case MzPeak::Schema::BufferFormat::ChunkTransform: + return true; + } + + std::unreachable(); +} + +/******************************************************************************/ +bool ArrayIndex::Dimension::is_main_axis() const +{ + for (const auto& entry : entries) { + switch (entry.buffer_format) { + case MzPeak::Schema::BufferFormat::Point: + return true; + case MzPeak::Schema::BufferFormat::ChunkStart: + return true; + case MzPeak::Schema::BufferFormat::ChunkEnd: + return true; + case MzPeak::Schema::BufferFormat::ChunkValues: + return true; + case MzPeak::Schema::BufferFormat::ChunkEncoding: + return true; + case MzPeak::Schema::BufferFormat::ChunkSecondary: + return false; + case MzPeak::Schema::BufferFormat::ChunkTransform: + continue; // Could be main or secondary. + } + } + + return false; +} + /******************************************************************************/ bool ArrayIndex::Dimension::needs_delta_model() const { @@ -56,10 +116,51 @@ Util::Type ArrayIndex::Dimension::type_or_throw() const throw TypeError("dimension " + name + " does not have a type set!"); } +/******************************************************************************/ +const ArrayIndex::Entry& ArrayIndex::Dimension::values_entry() const +{ + // There are a few buffer formats that indicate that an entry is + // definitely the column that stores dimension values. However, + // some of them (i.e. `ChunkTransform`) are ambitious so we need to + // consider them after all other entries have been considered. + // + // We don't assume the `entries` vector is in an particular order + // here. + Entry const* chunk_transform = nullptr; + + for (const auto& entry : entries) { + switch (entry.buffer_format) { + case MzPeak::Schema::BufferFormat::Point: + return entry; + case MzPeak::Schema::BufferFormat::ChunkStart: + continue; + case MzPeak::Schema::BufferFormat::ChunkEnd: + continue; + case MzPeak::Schema::BufferFormat::ChunkValues: + return entry; + case MzPeak::Schema::BufferFormat::ChunkEncoding: + continue; + case MzPeak::Schema::BufferFormat::ChunkSecondary: + return entry; + case MzPeak::Schema::BufferFormat::ChunkTransform: + chunk_transform = &entry; + continue; + } + } + + if (chunk_transform != nullptr && !is_main_axis()) { + return *chunk_transform; + } else { + std::string msg("dimension " + name + " lacks a data values column"); + throw InvalidFormatError(msg); + } +} + /******************************************************************************/ ArrayIndex::ArrayIndex(EntityType entity_type, const json::object& obj) : entity_type_(entity_type) , prefix_(obj.at("prefix").as_string()) + , layout_(group_name_to_layout(prefix_)) , entries_() , num_entities_() { @@ -139,6 +240,9 @@ EntityType ArrayIndex::entity_type() const { return entity_type_; } /******************************************************************************/ const std::string& ArrayIndex::prefix() const { return prefix_; } +/******************************************************************************/ +ArrayIndex::Layout ArrayIndex::layout() const { return layout_; } + /******************************************************************************/ const std::vector& ArrayIndex::entries() const { From 456b8beadf96a66eef6faaea89a3cb7492fd3a52 Mon Sep 17 00:00:00 2001 From: "Peter J. Jones" Date: Wed, 29 Jul 2026 10:55:16 +0200 Subject: [PATCH 03/31] decoding: Add some safety guards around list decoding --- include/mzpeak/util/decoders.h | 12 ++++++++++++ meson.build | 1 + src/util/decoders.cpp | 23 +++++++++++++++++++++++ 3 files changed, 36 insertions(+) create mode 100644 src/util/decoders.cpp diff --git a/include/mzpeak/util/decoders.h b/include/mzpeak/util/decoders.h index 368b2d4..58fbcc4 100644 --- a/include/mzpeak/util/decoders.h +++ b/include/mzpeak/util/decoders.h @@ -46,6 +46,12 @@ concept from_arrow_array = { t.decode(a, r) } -> std::same_as; }; +/******************************************************************************/ +/** + * Return `true` if the given Arrow array is a `ListArray`. + */ +bool is_list_array(const std::shared_ptr&); + /******************************************************************************/ /** * A NULL decoder that always skips NULL values. @@ -168,6 +174,12 @@ class List final : Helper> { /// Decoding function. void decode(const std::shared_ptr& src, C& dst) { + if (!is_list_array(src)) { + std::string msg("expected an arrow list array but found: "); + msg += src->type()->name(); + throw TypeError(msg); + } + std::shared_ptr casted = std::static_pointer_cast(src); diff --git a/meson.build b/meson.build index 62a13f7..97d2d0f 100644 --- a/meson.build +++ b/meson.build @@ -42,6 +42,7 @@ lib_sources = [ 'src/spectra.cpp', 'src/spectrum.cpp', 'src/util/arrow.cpp', + 'src/util/decoders.cpp', 'src/util/executor.cpp', 'src/util/manager.cpp', 'src/util/parquet.cpp', diff --git a/src/util/decoders.cpp b/src/util/decoders.cpp new file mode 100644 index 0000000..c6b38df --- /dev/null +++ b/src/util/decoders.cpp @@ -0,0 +1,23 @@ +/* + +This file is part of the mzpeak project. It is subject to the license +specified in the LICENSE file which can be found in the top-level +directory of this repository. + +*/ + +#include "mzpeak/util/decoders.h" + +namespace MzPeak::Util::Decoders { + +/******************************************************************************/ +bool is_list_array(const std::shared_ptr& ary) +{ + auto t = ary->type_id(); + + return t == arrow::Type::LIST || t == arrow::Type::FIXED_SIZE_LIST || + t == arrow::Type::LARGE_LIST || t == arrow::Type::LIST_VIEW || + t == arrow::Type::LARGE_LIST_VIEW; +} + +} // namespace MzPeak::Util::Decoders From 2ccd05aa3f74b4cba0ad2f40f73518f07d293ce8 Mon Sep 17 00:00:00 2001 From: "Peter J. Jones" Date: Wed, 29 Jul 2026 10:56:00 +0200 Subject: [PATCH 04/31] docs: Add documentation for template arguments --- include/mzpeak/util/slice.h | 12 ++++++++++++ 1 file changed, 12 insertions(+) diff --git a/include/mzpeak/util/slice.h b/include/mzpeak/util/slice.h index 5128dda..676f498 100644 --- a/include/mzpeak/util/slice.h +++ b/include/mzpeak/util/slice.h @@ -55,6 +55,12 @@ class Slice final { /** * Decode the first non-null value. + * + * Template Parameters: + * + * - T: The Decoder class to use (see MzPeak::Util::Decoders) + * + * - R: The destination object to update with the decoded value */ template > void singleton(const Column&, R&, T&& = {}) const; @@ -63,6 +69,12 @@ class Slice final { * Exact and decode an array. * * Use one of the decoders defined in `decoders.h`, or write your own. + * + * Template Parameters: + * + * - T: The Decoder class to use (see MzPeak::Util::Decoders) + * + * - V: The destination object to fill with decoded values */ template > requires Decoders::from_arrow_array From dba49e3723cb1d9122cd2b5139d6c0a518dc906a Mon Sep 17 00:00:00 2001 From: "Peter J. Jones" Date: Wed, 29 Jul 2026 10:56:45 +0200 Subject: [PATCH 05/31] decoding: Count each element of list arrays When array elements are lists, account for the elements in the list when reserving space in the destination vector. --- include/mzpeak/util/slice.h | 6 +++++- 1 file changed, 5 insertions(+), 1 deletion(-) diff --git a/include/mzpeak/util/slice.h b/include/mzpeak/util/slice.h index 676f498..dc7ea37 100644 --- a/include/mzpeak/util/slice.h +++ b/include/mzpeak/util/slice.h @@ -134,7 +134,11 @@ void Slice::array(const Column& field, V& v, T&& t) const std::size_t size{}; for (const auto& chunk : *chunks) { - size += chunk->length(); + if (Decoders::is_list_array(chunk)) { + size += std::static_pointer_cast(chunk)->length(); + } else { + size += chunk->length(); + } } v.reserve(v.size() + size); From 834d44fdce168581a2a70f25d12f7e92d8fa45e7 Mon Sep 17 00:00:00 2001 From: "Peter J. Jones" Date: Wed, 29 Jul 2026 10:58:11 +0200 Subject: [PATCH 06/31] decoding: The list decoder should pass a null decoder to the scalar class This will be needed to decoding null values in the chunk encoding. --- include/mzpeak/util/decoders.h | 25 +++++++++++++++++++++---- 1 file changed, 21 insertions(+), 4 deletions(-) diff --git a/include/mzpeak/util/decoders.h b/include/mzpeak/util/decoders.h index 58fbcc4..17dcdb6 100644 --- a/include/mzpeak/util/decoders.h +++ b/include/mzpeak/util/decoders.h @@ -155,19 +155,33 @@ class Scalar final : Helper> { /** * A decoder where array elements are lists. * - * NULL `ListArray` elements, and NULL elements inside the - * `ListArray` are skipped. + * NULL `ListArray` elements are skipped. + * + * NULL elements inside the `ListArray` elements are decoded using the + * (optionally) provided null decoder. */ -template > +template , typename N = NullSkip> requires Decoders::scalar_or_container_of class List final : Helper> { public: /// Decodes vectors of type T. using value_type = std::vector; + /// The range or scalar type. + using range_type = C; + + /// The null decoder type. + using null_decoder_type = N; + /// Constructor. List() {} + // Constructor where you can pass a null decoder to the scalar decoder. + List(const null_decoder_type& null_decoder) + : scalar_decoder_(null_decoder) + { + } + /// Destructor. ~List() = default; @@ -188,11 +202,14 @@ class List final : Helper> { std::shared_ptr values(casted->value_slice(i)); value_type res; res.reserve(values->length()); - Scalar().decode(values, res); + scalar_decoder_.decode(values, res); this->push(dst, res); } } } + +private: + Scalar scalar_decoder_; }; } // namespace MzPeak::Util::Decoders From 91418b50ffeb0f105715309ce38baa3da732f943 Mon Sep 17 00:00:00 2001 From: "Peter J. Jones" Date: Wed, 29 Jul 2026 14:41:29 +0200 Subject: [PATCH 07/31] decoding: Add the ability to transform and flatten arrays --- include/mzpeak/data/encoding.h | 37 +++++++-------- include/mzpeak/util/decoders.h | 86 ++++++++++++++++++++++++++++++---- include/mzpeak/util/slice.h | 13 +++++ 3 files changed, 107 insertions(+), 29 deletions(-) diff --git a/include/mzpeak/data/encoding.h b/include/mzpeak/data/encoding.h index 942ba1b..878428e 100644 --- a/include/mzpeak/data/encoding.h +++ b/include/mzpeak/data/encoding.h @@ -56,8 +56,9 @@ template class Decoder { void decode(const ArrayIndex::Dimension&, std::vector&) const; template - void - point(const ArrayIndex::Dimension&, const N& null_decoder, std::vector&) const; + void decode_with_nulls(const ArrayIndex::Dimension&, + const N& null_decoder, + std::vector&) const; template void remap(const ArrayIndex::Dimension& dim, std::vector& v) const; @@ -131,19 +132,16 @@ void Decoder::decode(const ArrayIndex::Dimension& dim, std::vector& v) con { switch (signals_->array_index()->layout()) { case ArrayIndex::Layout::Point: + case ArrayIndex::Layout::Chunked: if (dim.needs_delta_model()) { using N = NullMarking::Decoder; - point(dim, N{delta_estimator_}, v); + decode_with_nulls(dim, N{delta_estimator_}, v); } else { using N = Util::Decoders::NullToZero; - point(dim, N{}, v); + decode_with_nulls(dim, N{}, v); } break; - case ArrayIndex::Layout::Chunked: - throw("not implemented"); - break; - case ArrayIndex::Layout::Unknown: throw UnknownLayoutError("cannot decode dimension: " + dim.name); } @@ -152,18 +150,12 @@ void Decoder::decode(const ArrayIndex::Dimension& dim, std::vector& v) con /******************************************************************************/ template template -void Decoder::point(const ArrayIndex::Dimension& dim, - const N& null_decoder, - std::vector& v) const +void Decoder::decode_with_nulls(const ArrayIndex::Dimension& dim, + const N& null_decoder, + std::vector& v) const { const auto& primary_entry = dim.values_entry(); - if (primary_entry.buffer_format != Schema::BufferFormat::Point) { - std::string msg("file uses point layout, but " + dim.name); - msg += " is not using the point buffer_format"; - throw InvalidFormatError(msg); - } - auto col = signals_->array_index()->entry_column(*signals_->groups(), primary_entry); @@ -171,8 +163,15 @@ void Decoder::point(const ArrayIndex::Dimension& dim, throw ParquetError("unable to decode dimension, not in schema: " + dim.name); } - slice_->array(col.value(), v, - Util::Decoders::Scalar, N>(null_decoder)); + if (primary_entry.buffer_format == Schema::BufferFormat::Point) { + auto decoder = Util::Decoders::Scalar, N>(null_decoder); + slice_->array(col.value(), v, decoder); + } else { + auto decoder = Util::Decoders::Flattened, N>(null_decoder); + slice_->array(col.value(), v, decoder); + } + + // FIXME: Apply necessary transformations on the decoded array. } } // namespace MzPeak::Data::Encoding diff --git a/include/mzpeak/util/decoders.h b/include/mzpeak/util/decoders.h index 17dcdb6..b5ab8b0 100644 --- a/include/mzpeak/util/decoders.h +++ b/include/mzpeak/util/decoders.h @@ -33,10 +33,6 @@ concept scalar_or_container_of = * `T` is a type that has a `decode` function that can decode values * from an `arrow::Array` and place the result in `R`. The `R` type * can be a container or scalar value. - * - * The `decode` function should return `true` to indicate it can - * continue to decode values. If it returns `false` the chunk - * decoding will stop. */ template concept from_arrow_array = @@ -122,9 +118,6 @@ class Scalar final : Helper> { { } - /// Destructor. - ~Scalar() = default; - /// Decoding function. void decode(const std::shared_ptr& src, C& dst) { @@ -182,9 +175,6 @@ class List final : Helper> { { } - /// Destructor. - ~List() = default; - /// Decoding function. void decode(const std::shared_ptr& src, C& dst) { @@ -212,4 +202,80 @@ class List final : Helper> { Scalar scalar_decoder_; }; +/******************************************************************************/ +/** + * An array transformer that returns its argument unchanged. + */ +struct IdentityTransform { + std::shared_ptr&& operator()(int64_t, + std::shared_ptr&& a) const + { + return a; + } +}; + +/******************************************************************************/ +/** + * A decoder that handles array elements that are lists, and the + * destination object is a list of scalar values. + * + * This class can decode null values in the lists, and also transform + * the lists using a helper object. Once use of the transformer + * object is to decode delta encoding prior to null decoding. + * + * Transformers are called with two arguments: + * + * - The index of the array element currently being decoded. + * + * - The array element itself, as an Arrow Array. + * + * The transformer should return a `std::shared_ptr` + * which contains the transformed values that can be decoded. + */ +template , + typename NullDecoder = NullSkip, + typename Transformer = IdentityTransform> + requires Decoders::scalar_or_container_of +class Flattened final : Helper> { +public: + /// The types of values this decoder can decode. + using value_type = Value; + + // Constructor where you can pass a null decoder to the scalar decoder. + Flattened(const NullDecoder& null_decoder, Transformer transformer = {}) + : scalar_decoder_(null_decoder) + , transformer_(transformer) + { + } + + /// Decoding function. + void decode(const std::shared_ptr& src, Container& dst) + { + if (!is_list_array(src)) { + std::string msg("expected an arrow list array but found: "); + msg += src->type()->name(); + throw TypeError(msg); + } + + std::shared_ptr casted = + std::static_pointer_cast(src); + + for (int64_t i : std::views::iota(0, casted->length())) { + if (casted->IsValid(i)) { + std::shared_ptr values( + transformer_(index_, casted->value_slice(i))); + scalar_decoder_.decode(values, dst); + } + } + + ++index_; + } + +private: + Scalar scalar_decoder_; + Transformer transformer_; + int64_t index_ = 0; +}; + } // namespace MzPeak::Util::Decoders diff --git a/include/mzpeak/util/slice.h b/include/mzpeak/util/slice.h index dc7ea37..5322690 100644 --- a/include/mzpeak/util/slice.h +++ b/include/mzpeak/util/slice.h @@ -80,6 +80,11 @@ class Slice final { requires Decoders::from_arrow_array void array(const Column&, V&, T&& = {}) const; + /****************************************************************************/ + template > + requires Decoders::from_arrow_array + void array(const Column&, V&, T&) const; + private: friend class MzPeak::Util::Executor; @@ -126,6 +131,14 @@ void Slice::singleton(const Column& field, R& dst, T&& t) const template requires Decoders::from_arrow_array void Slice::array(const Column& field, V& v, T&& t) const +{ + array(field, v, t); +} + +/******************************************************************************/ +template + requires Decoders::from_arrow_array +void Slice::array(const Column& field, V& v, T& t) const { std::shared_ptr chunks = raw(field); if (chunks == nullptr) return; From 7ec340809d07b27518d19a7d0a7487983f6eb758 Mon Sep 17 00:00:00 2001 From: "Peter J. Jones" Date: Wed, 29 Jul 2026 15:26:14 +0200 Subject: [PATCH 08/31] decoding: New class to track the chunk encoding method --- include/mzpeak/schema/psi/chunk_encoding.h | 65 +++++++++++++++++++++ meson.build | 2 + src/schema/psi/chunk_encoding.cpp | 66 ++++++++++++++++++++++ 3 files changed, 133 insertions(+) create mode 100644 include/mzpeak/schema/psi/chunk_encoding.h create mode 100644 src/schema/psi/chunk_encoding.cpp diff --git a/include/mzpeak/schema/psi/chunk_encoding.h b/include/mzpeak/schema/psi/chunk_encoding.h new file mode 100644 index 0000000..eed958a --- /dev/null +++ b/include/mzpeak/schema/psi/chunk_encoding.h @@ -0,0 +1,65 @@ +/* + +This file is part of the mzpeak project. It is subject to the license +specified in the LICENSE file which can be found in the top-level +directory of this repository. + +*/ + +#pragma once + +#include + +#include "mzpeak/schema/cv.h" + +namespace MzPeak::Schema::PSI { + +/** + * Encoding methods used in the `chunk_encoding` column. + */ +class ChunkEncoding final { +public: + enum Type { + /// MS:1000576 + /// + /// Values are not encoded nor compressed. + NoCompression, + + /// MS:1003089 + /// + /// Data array compression using mantissa bit truncation, delta + /// prediction and zlib compression. + /// + /// NOTE: Parquet takes care of everything except the delta + /// encoding. + Delta, + + /// MS:1002312 + /// + /// Compression using MS-Numpress linear prediction compression. + NumpressLinear, + }; + + // Internal storage type. + using value_type = std::variant; + + /** + * Construct an encoding type from a CV. + */ + ChunkEncoding(const CV&); + + /** + * Convert an encoding type to a CV. + */ + CV to_cv() const; + + /** + * Return the encoding type if it is known. + */ + std::optional type() const; + +private: + value_type val_; +}; + +} // namespace MzPeak::Schema::PSI diff --git a/meson.build b/meson.build index 97d2d0f..90c7145 100644 --- a/meson.build +++ b/meson.build @@ -37,6 +37,7 @@ lib_sources = [ 'src/schema/file.cpp', 'src/schema/group.cpp', 'src/schema/psi/array_type.cpp', + 'src/schema/psi/chunk_encoding.cpp', 'src/schema/psi/data_type.cpp', 'src/schema/psi/transform.cpp', 'src/spectra.cpp', @@ -76,6 +77,7 @@ install_headers([ 'include/mzpeak/schema/file.h', 'include/mzpeak/schema/group.h', 'include/mzpeak/schema/psi/array_type.h', + 'include/mzpeak/schema/psi/chunk_encoding.h', 'include/mzpeak/schema/psi/data_type.h', 'include/mzpeak/schema/psi/transform.h', 'include/mzpeak/spectra.h', diff --git a/src/schema/psi/chunk_encoding.cpp b/src/schema/psi/chunk_encoding.cpp new file mode 100644 index 0000000..662dd4f --- /dev/null +++ b/src/schema/psi/chunk_encoding.cpp @@ -0,0 +1,66 @@ +/* + +This file is part of the mzpeak project. It is subject to the license +specified in the LICENSE file which can be found in the top-level +directory of this repository. + +*/ + +#include + +#include "mzpeak/schema/psi/chunk_encoding.h" + +namespace MzPeak::Schema::PSI { + +/******************************************************************************/ +ChunkEncoding::ChunkEncoding(const CV& cv) + : val_(cv) +{ + if (cv.code() == "MS") { + const auto& accession = cv.accession(); + + if (accession == "1000576") { + val_ = NoCompression; + } else if (accession == "1003089") { + val_ = Delta; + } else if (accession == "1002312") { + val_ = NumpressLinear; + } + } +} + +/******************************************************************************/ +CV ChunkEncoding::to_cv() const +{ + return std::visit( + [](auto&& v) -> CV { + using T = std::decay_t; + if constexpr (std::is_same_v) { + return v; + } else if constexpr (std::is_same_v) { + switch (v) { + case NoCompression: + return CV("MS", "1000576"); + case Delta: + return CV("MS", "1003089"); + case NumpressLinear: + return CV("MS", "1002312"); + } + + std::unreachable(); + } + }, + val_); +} + +/******************************************************************************/ +std::optional ChunkEncoding::type() const +{ + if (std::holds_alternative(val_)) { + return std::get(val_); + } else { + return {}; + } +} + +} // namespace MzPeak::Schema::PSI From eec9c4784f2292cea1977394f6876c9378a54f58 Mon Sep 17 00:00:00 2001 From: "Peter J. Jones" Date: Wed, 29 Jul 2026 17:11:29 +0200 Subject: [PATCH 09/31] algorithm: Add function for decoding "delta encoding will nulls" --- include/mzpeak/exception.h | 15 +++++++ include/mzpeak/util/algorithm.h | 75 +++++++++++++++++++++++++++++++++ include/mzpeak/util/types.h | 9 ++++ meson.build | 1 + test/algorithm_test.cpp | 51 ++++++++++++++++++++++ 5 files changed, 151 insertions(+) create mode 100644 test/algorithm_test.cpp diff --git a/include/mzpeak/exception.h b/include/mzpeak/exception.h index 7cad95c..f8b9d9f 100644 --- a/include/mzpeak/exception.h +++ b/include/mzpeak/exception.h @@ -119,4 +119,19 @@ class TypeError final : public Exception { ~TypeError() = default; }; +/** + * Failed to allocate memory. + */ +class AllocationError final : public Exception { +public: + /// Constructor. + AllocationError(const std::string& msg) + : Exception(msg) + { + } + + /// Destructor. + ~AllocationError() = default; +}; + } // namespace MzPeak diff --git a/include/mzpeak/util/algorithm.h b/include/mzpeak/util/algorithm.h index 99388ba..159597a 100644 --- a/include/mzpeak/util/algorithm.h +++ b/include/mzpeak/util/algorithm.h @@ -8,12 +8,15 @@ top-level directory of this repository. #pragma once +#include #include #include #include #include #include +#include "mzpeak/util/types.h" + namespace MzPeak::Util::Algorithm { /** @@ -78,4 +81,76 @@ template T median_delta(const std::vector& values, T or_else) return boost::math::statistics::median(ds.begin(), ds.end()); } +/** + * Decode an Arrow array that was encoded with delta encoding. + * + * Parameters: + * + * - start: The excluded starting value. + * + * - src: The array to decode. + * + * Returns a newly allocated array containing the decoded values. + * Null values transferred from the source array to the returned array + * unchanged. + */ +template +std::shared_ptr +null_delta_decode(typename type_traits::value_type start, + const std::shared_ptr& src) +{ + using ValueType = typename type_traits::value_type; + ValueType zero = {}; + + std::optional last(start); + int64_t length = src->length(); + + using Builder = type_traits::builder_type; + Builder builder; + + // N.B.: "The start point is *excluded* from the chunk-values array." + if (!builder.Reserve(length + 1).ok()) { + throw AllocationError("unable to decode delta encoded array"); + } + + using ArrayType = typename type_traits::array_type; + std::shared_ptr casted = std::static_pointer_cast(src); + + auto append = [&](const std::optional& v) -> void { + arrow::Status status; + + if (v.has_value()) { + status = builder.Append(v.value()); + } else { + status = builder.AppendNull(); + } + + if (!status.ok()) { + throw AllocationError("unable to decode delta encoded element"); + } + }; + + // N.B.: "The start point is *excluded* from the chunk-values array." + append({start}); + + for (int64_t index : std::views::iota(0, length)) { + if (casted->IsValid(index)) { + ValueType delta = casted->Value(index); + last = last.value_or(zero) + delta; + append(last); + } else { + last = {}; + append(last); + } + } + + std::shared_ptr result; + + if (!builder.Finish(&result).ok()) { + throw AllocationError("failed to decode delta encoded array"); + } + + return result; +} + } // namespace MzPeak::Util::Algorithm diff --git a/include/mzpeak/util/types.h b/include/mzpeak/util/types.h index 72ebaf0..47a643d 100644 --- a/include/mzpeak/util/types.h +++ b/include/mzpeak/util/types.h @@ -103,6 +103,7 @@ template <> struct type_traits { using value_type = int8_t; using parquet_type = parquet::Int32Type; using array_type = arrow::Int8Array; + using builder_type = arrow::Int8Builder; }; template <> struct type_traits { @@ -110,6 +111,7 @@ template <> struct type_traits { using value_type = uint8_t; using parquet_type = parquet::Int32Type; using array_type = arrow::UInt8Array; + using builder_type = arrow::UInt8Builder; }; template <> struct type_traits { @@ -117,6 +119,7 @@ template <> struct type_traits { using value_type = int32_t; using parquet_type = parquet::Int32Type; using array_type = arrow::Int32Array; + using builder_type = arrow::Int32Builder; }; template <> struct type_traits { @@ -124,6 +127,7 @@ template <> struct type_traits { using value_type = uint32_t; using parquet_type = parquet::Int32Type; using array_type = arrow::UInt32Array; + using builder_type = arrow::UInt32Builder; }; template <> struct type_traits { @@ -131,6 +135,7 @@ template <> struct type_traits { using value_type = int64_t; using parquet_type = parquet::Int64Type; using array_type = arrow::Int64Array; + using builder_type = arrow::Int64Builder; }; template <> struct type_traits { @@ -138,6 +143,7 @@ template <> struct type_traits { using value_type = uint64_t; using parquet_type = parquet::Int64Type; using array_type = arrow::UInt64Array; + using builder_type = arrow::UInt64Builder; }; template <> struct type_traits { @@ -145,6 +151,7 @@ template <> struct type_traits { using value_type = float; using parquet_type = parquet::FloatType; using array_type = arrow::FloatArray; + using builder_type = arrow::FloatBuilder; }; template <> struct type_traits { @@ -152,6 +159,7 @@ template <> struct type_traits { using value_type = double; using parquet_type = parquet::DoubleType; using array_type = arrow::DoubleArray; + using builder_type = arrow::DoubleBuilder; }; template <> struct type_traits { @@ -159,6 +167,7 @@ template <> struct type_traits { using value_type = std::string_view; using parquet_type = parquet::ByteArrayType; using array_type = arrow::StringArray; + using builder_type = arrow::StringBuilder; }; /******************************************************************************/ diff --git a/meson.build b/meson.build index 90c7145..c6983e8 100644 --- a/meson.build +++ b/meson.build @@ -136,6 +136,7 @@ libmzpeak_a = static_library( # Testing test_names = [ + 'algorithm', 'array_index', 'arrow', 'directory', diff --git a/test/algorithm_test.cpp b/test/algorithm_test.cpp new file mode 100644 index 0000000..438f5ee --- /dev/null +++ b/test/algorithm_test.cpp @@ -0,0 +1,51 @@ +/* + +This file is part of the mzpeak project. It is subject to the license +specified in the LICENSE file which can be found in the top-level +directory of this repository. + +*/ + +#define BOOST_TEST_MODULE Algorithm +#include + +#include +#include + +#include "mzpeak/util/algorithm.h" + +/******************************************************************************/ +BOOST_AUTO_TEST_CASE(null_delta_decode_with_no_nulls) +{ + using namespace MzPeak::Util; + + arrow::Status status; + arrow::Int64Builder builder; + int64_t vals[] = {1, 2, 3}; + int64_t start = 2; + + status = builder.AppendValues(vals, sizeof(vals) / sizeof(int64_t)); + BOOST_TEST(status.ok()); + + std::shared_ptr input; + status = builder.Finish(&input); + BOOST_TEST(status.ok()); + + std::shared_ptr output = + Algorithm::null_delta_decode(start, input); + + std::shared_ptr casted = + std::static_pointer_cast(output); + + std::vector decoded; + decoded.reserve(casted->length()); + + for (int64_t index : std::views::iota(0, casted->length())) { + decoded.push_back(casted->Value(index)); + } + + BOOST_TEST(decoded.size() == casted->length()); + + std::vector expected = {2, 3, 5, 8}; + BOOST_TEST(decoded == expected); +} From a970cbc70d2242fa4090c3b5a8fd889edc587e8f Mon Sep 17 00:00:00 2001 From: "Peter J. Jones" Date: Thu, 30 Jul 2026 15:19:44 +0200 Subject: [PATCH 10/31] dimension: Find array index entries by their buffer format --- include/mzpeak/data/array_index.h | 3 +++ src/data/array_index.cpp | 13 +++++++++++++ 2 files changed, 16 insertions(+) diff --git a/include/mzpeak/data/array_index.h b/include/mzpeak/data/array_index.h index cceed96..ed0268c 100644 --- a/include/mzpeak/data/array_index.h +++ b/include/mzpeak/data/array_index.h @@ -148,6 +148,9 @@ class ArrayIndex final { /// malformed and thus doesn't include any of the expected /// entries. const Entry& values_entry() const; + + /// Find the first entry with the given buffer format. + std::optional entry_with(BufferFormat) const; }; /// Default constructor. diff --git a/src/data/array_index.cpp b/src/data/array_index.cpp index 340a199..41bb189 100644 --- a/src/data/array_index.cpp +++ b/src/data/array_index.cpp @@ -156,6 +156,19 @@ const ArrayIndex::Entry& ArrayIndex::Dimension::values_entry() const } } +/******************************************************************************/ +std::optional +ArrayIndex::Dimension::entry_with(BufferFormat format) const +{ + auto it = std::ranges::find(entries, format, &ArrayIndex::Entry::buffer_format); + + if (it == entries.end()) { + return {}; + } else { + return *it; + } +} + /******************************************************************************/ ArrayIndex::ArrayIndex(EntityType entity_type, const json::object& obj) : entity_type_(entity_type) From 730cdf0be37161e42265722f9a96aded0bc671b0 Mon Sep 17 00:00:00 2001 From: "Peter J. Jones" Date: Thu, 30 Jul 2026 15:21:26 +0200 Subject: [PATCH 11/31] chunk encoding: The chunk_start column is needed for decoding When using delta encoding the start value is necessary to reconstruct the remaining values. --- src/data/array_index.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/data/array_index.cpp b/src/data/array_index.cpp index 41bb189..267e265 100644 --- a/src/data/array_index.cpp +++ b/src/data/array_index.cpp @@ -35,7 +35,7 @@ bool ArrayIndex::Entry::needed_for_decoding() const case MzPeak::Schema::BufferFormat::Point: return true; case MzPeak::Schema::BufferFormat::ChunkStart: - return false; + return true; case MzPeak::Schema::BufferFormat::ChunkEnd: return false; case MzPeak::Schema::BufferFormat::ChunkValues: From 83a740c248d711828a46a511434c9e6ea462eb76 Mon Sep 17 00:00:00 2001 From: "Peter J. Jones" Date: Thu, 30 Jul 2026 15:25:27 +0200 Subject: [PATCH 12/31] signals: Rename (field -> column) and overload function The name `column` is more appropriate, and the overloads reduce code bloat in other compilation units. --- include/mzpeak/data/signals.h | 18 ++++++++++++++++-- src/data/signals.cpp | 18 ++++++++++++++++-- test/array_index_test.cpp | 4 ++-- 3 files changed, 34 insertions(+), 6 deletions(-) diff --git a/include/mzpeak/data/signals.h b/include/mzpeak/data/signals.h index b62800a..7cef70c 100644 --- a/include/mzpeak/data/signals.h +++ b/include/mzpeak/data/signals.h @@ -52,11 +52,25 @@ class Signals { const Util::Query&); /** - * Low-level interface for accessing a group field given its name. + * Low-level interface for accessing a column given its name. * * Useful if you need to manually construct queries. */ - std::optional field(const std::string_view&) const; + std::optional column(const std::string_view&) const; + + /** + * Low-level interface for accessing a column given an array index entry. + */ + std::optional column(const ArrayIndex::Entry&) const; + + /** + * Low-level interface for accessing a column given a dimension and + * a buffer format. + * + * Returns the first matching column. + */ + std::optional column(const ArrayIndex::Dimension&, + Schema::BufferFormat) const; /** * Low-level interface for accessing the schema encoded as a map of diff --git a/src/data/signals.cpp b/src/data/signals.cpp index 325413a..9e5b660 100644 --- a/src/data/signals.cpp +++ b/src/data/signals.cpp @@ -102,11 +102,25 @@ std::size_t Signals::record_count() const } /******************************************************************************/ -std::optional Signals::field(const std::string_view& name) const +std::optional Signals::column(const std::string_view& name) const { return impl_->parquet_->field(impl_->array_index_->prefix(), name); } +/******************************************************************************/ +std::optional Signals::column(const ArrayIndex::Entry& entry) const +{ + return impl_->array_index_->entry_column(*impl_->parquet_->groups(), entry); +} + +/******************************************************************************/ +std::optional Signals::column(const ArrayIndex::Dimension& dim, + Schema::BufferFormat format) const +{ + return dim.entry_with(format).and_then( + [this](const auto& entry) { return column(entry); }); +} + /******************************************************************************/ const std::shared_ptr& Signals::groups() const { @@ -118,7 +132,7 @@ Util::Query::Builder Signals::index() const { auto entity_type = impl_->array_index_->entity_type(); auto field_name = Schema::entity_type_to_string(entity_type) + "_index"; - auto index_field = field(field_name); + auto index_field = column(field_name); if (!index_field.has_value()) { throw ParquetError("parquet file is missing the index column: " + field_name); diff --git a/test/array_index_test.cpp b/test/array_index_test.cpp index 1fda385..11746d7 100644 --- a/test/array_index_test.cpp +++ b/test/array_index_test.cpp @@ -78,8 +78,8 @@ BOOST_AUTO_TEST_CASE(can_read_mz_array) Data::Signals data(std::move(parquet)); - auto index_field = data.field("spectrum_index"); - auto mz_column = data.field("mz"); + auto index_field = data.column("spectrum_index"); + auto mz_column = data.column("mz"); BOOST_TEST(index_field.has_value()); BOOST_TEST(mz_column.has_value()); From ec6e25c05ec7cce6c8063ec493f8d5d5748e39de Mon Sep 17 00:00:00 2001 From: "Peter J. Jones" Date: Thu, 30 Jul 2026 15:26:54 +0200 Subject: [PATCH 13/31] algo: Resolve "undetermined type" issues related to arrow:Builder --- include/mzpeak/util/algorithm.h | 1 + 1 file changed, 1 insertion(+) diff --git a/include/mzpeak/util/algorithm.h b/include/mzpeak/util/algorithm.h index 159597a..256d210 100644 --- a/include/mzpeak/util/algorithm.h +++ b/include/mzpeak/util/algorithm.h @@ -9,6 +9,7 @@ top-level directory of this repository. #pragma once #include +#include #include #include #include From 08adb7097cadcbd42003a7666f3b0a4f159e5702 Mon Sep 17 00:00:00 2001 From: "Peter J. Jones" Date: Thu, 30 Jul 2026 16:40:03 +0200 Subject: [PATCH 14/31] decoding: Allow transformers to directly decode an array This is necessary for numpress decoders that need to turn an array of bytes into an array of doubles. --- include/mzpeak/util/decoders.h | 31 ++++++++++++++++++++++++++----- 1 file changed, 26 insertions(+), 5 deletions(-) diff --git a/include/mzpeak/util/decoders.h b/include/mzpeak/util/decoders.h index b5ab8b0..90fb694 100644 --- a/include/mzpeak/util/decoders.h +++ b/include/mzpeak/util/decoders.h @@ -229,8 +229,13 @@ struct IdentityTransform { * * - The array element itself, as an Arrow Array. * - * The transformer should return a `std::shared_ptr` - * which contains the transformed values that can be decoded. + * The transformer can return one of two types: + * + * - `std::shared_ptr` which contains the transformed + * values that can then be decoded. + * + * - A container of decoded values when can be inserted into the + * destination vector and further decoding can be skipped. */ template , @@ -242,6 +247,10 @@ class Flattened final : Helper> { /// The types of values this decoder can decode. using value_type = Value; + /// The type of return value allowed from transformers. + using transform_result_type = + std::variant, std::shared_ptr>; + // Constructor where you can pass a null decoder to the scalar decoder. Flattened(const NullDecoder& null_decoder, Transformer transformer = {}) : scalar_decoder_(null_decoder) @@ -263,9 +272,21 @@ class Flattened final : Helper> { for (int64_t i : std::views::iota(0, casted->length())) { if (casted->IsValid(i)) { - std::shared_ptr values( - transformer_(index_, casted->value_slice(i))); - scalar_decoder_.decode(values, dst); + transform_result_type values(transformer_(index_, casted->value_slice(i))); + + std::visit( + [&](auto&& v) -> void { + using U = std::decay_t; + + if constexpr (std::is_same_v>) { + scalar_decoder_.decode(v, dst); + } else if constexpr (std::is_same_v>) { + dst.insert(dst.end(), v->begin(), v->end()); + } else { + static_assert(false_type, "invalid transform result"); + } + }, + values); } } From e62107fdd560faa3c15a8f5e8abc9a231817d72d Mon Sep 17 00:00:00 2001 From: "Peter J. Jones" Date: Thu, 30 Jul 2026 16:41:07 +0200 Subject: [PATCH 15/31] schema: Prefer chunk_transform to chunk_values If the `chunk_transform` column is present then it is used instead of the `chunk_values` column, which will also be present but NULL. --- src/data/array_index.cpp | 8 ++++++-- 1 file changed, 6 insertions(+), 2 deletions(-) diff --git a/src/data/array_index.cpp b/src/data/array_index.cpp index 267e265..c19b200 100644 --- a/src/data/array_index.cpp +++ b/src/data/array_index.cpp @@ -127,6 +127,7 @@ const ArrayIndex::Entry& ArrayIndex::Dimension::values_entry() const // We don't assume the `entries` vector is in an particular order // here. Entry const* chunk_transform = nullptr; + Entry const* chunk_values = nullptr; for (const auto& entry : entries) { switch (entry.buffer_format) { @@ -137,7 +138,8 @@ const ArrayIndex::Entry& ArrayIndex::Dimension::values_entry() const case MzPeak::Schema::BufferFormat::ChunkEnd: continue; case MzPeak::Schema::BufferFormat::ChunkValues: - return entry; + chunk_values = &entry; + continue; case MzPeak::Schema::BufferFormat::ChunkEncoding: continue; case MzPeak::Schema::BufferFormat::ChunkSecondary: @@ -148,8 +150,10 @@ const ArrayIndex::Entry& ArrayIndex::Dimension::values_entry() const } } - if (chunk_transform != nullptr && !is_main_axis()) { + if (chunk_transform != nullptr) { return *chunk_transform; + } else if (chunk_values != nullptr) { + return *chunk_values; } else { std::string msg("dimension " + name + " lacks a data values column"); throw InvalidFormatError(msg); From 8481ea810820c1dfd42a81bb7f587590d5c78545 Mon Sep 17 00:00:00 2001 From: "Peter J. Jones" Date: Thu, 30 Jul 2026 17:06:20 +0200 Subject: [PATCH 16/31] decoding: Accommodate encoding schemes that use a starting value The destination vector needs more space to avoid re-allocations and transformers need a way to return one extra value. --- include/mzpeak/util/decoders.h | 14 ++++++++++++-- include/mzpeak/util/slice.h | 14 +++++++++++++- 2 files changed, 25 insertions(+), 3 deletions(-) diff --git a/include/mzpeak/util/decoders.h b/include/mzpeak/util/decoders.h index 90fb694..6cdae2a 100644 --- a/include/mzpeak/util/decoders.h +++ b/include/mzpeak/util/decoders.h @@ -229,11 +229,15 @@ struct IdentityTransform { * * - The array element itself, as an Arrow Array. * - * The transformer can return one of two types: + * The transformer can return one of the following types: * * - `std::shared_ptr` which contains the transformed * values that can then be decoded. * + * - A pair where the first element is a starting value that should + * be inserted into the destination and the second element is an + * Arrow array to decode. + * * - A container of decoded values when can be inserted into the * destination vector and further decoding can be skipped. */ @@ -249,7 +253,9 @@ class Flattened final : Helper> { /// The type of return value allowed from transformers. using transform_result_type = - std::variant, std::shared_ptr>; + std::variant, + std::pair>, + std::shared_ptr>; // Constructor where you can pass a null decoder to the scalar decoder. Flattened(const NullDecoder& null_decoder, Transformer transformer = {}) @@ -277,9 +283,13 @@ class Flattened final : Helper> { std::visit( [&](auto&& v) -> void { using U = std::decay_t; + using P = std::pair>; if constexpr (std::is_same_v>) { scalar_decoder_.decode(v, dst); + } else if constexpr (std::is_same_v) { + dst.push_back(v.first); + scalar_decoder_.decode(v.second, dst); } else if constexpr (std::is_same_v>) { dst.insert(dst.end(), v->begin(), v->end()); } else { diff --git a/include/mzpeak/util/slice.h b/include/mzpeak/util/slice.h index 5322690..9d043ff 100644 --- a/include/mzpeak/util/slice.h +++ b/include/mzpeak/util/slice.h @@ -148,7 +148,19 @@ void Slice::array(const Column& field, V& v, T& t) const for (const auto& chunk : *chunks) { if (Decoders::is_list_array(chunk)) { - size += std::static_pointer_cast(chunk)->length(); + // N.B. When decoding lists they often are accompanied by a + // starting value (i.e. chunk_start) that must also be + // accommodated. + // + // For example, both delta encoding and "basic" encoding + // result in an array that is one element longer than the one + // stored in Parquet. + // + // In the numpress case we'll over-allocate by a small amount. + // Perhaps this code needs to be made smarter or moved + // somewhere else so it has enough context to decide how many + // bytes to allocate. + size += std::static_pointer_cast(chunk)->length() + 1; } else { size += chunk->length(); } From bb83d796bbd095fd74b8b8feece7e838920336ce Mon Sep 17 00:00:00 2001 From: "Peter J. Jones" Date: Thu, 6 Aug 2026 11:51:20 +0200 Subject: [PATCH 17/31] Squashed 'subprojects/msnumpress/' content from commit 94090cf git-subtree-dir: subprojects/msnumpress git-subtree-split: 94090cf72457fe91ff709a0827ff175e53b83cd6 --- README.md | 146 ++++ meson.build | 38 + pom.xml | 33 + src/main/R/RMSNumpress/DESCRIPTION | 14 + src/main/R/RMSNumpress/LICENSE | 3 + src/main/R/RMSNumpress/NAMESPACE | 3 + src/main/R/RMSNumpress/R/RcppExports.R | 189 +++++ src/main/R/RMSNumpress/inst/LICENSE.md | 31 + .../R/RMSNumpress/man/RMSNumpress-package.Rd | 131 +++ src/main/R/RMSNumpress/man/decodeLinear.Rd | 43 + src/main/R/RMSNumpress/man/decodePic.Rd | 27 + src/main/R/RMSNumpress/man/decodeSlof.Rd | 38 + src/main/R/RMSNumpress/man/encodeLinear.Rd | 42 + src/main/R/RMSNumpress/man/encodePic.Rd | 26 + src/main/R/RMSNumpress/man/encodeSlof.Rd | 28 + .../man/optimalLinearFixedPoint.Rd | 17 + .../man/optimalLinearFixedPointMass.Rd | 25 + .../RMSNumpress/man/optimalSlofFixedPoint.Rd | 17 + src/main/R/RMSNumpress/src/MSNumpress.cpp | 802 ++++++++++++++++++ src/main/R/RMSNumpress/src/RMSNumpress.cpp | 239 ++++++ src/main/R/RMSNumpress/src/RcppExports.cpp | 127 +++ .../R/RMSNumpress/src/include/MSNumpress.hpp | 331 ++++++++ src/main/R/RMSNumpress/tests/testthat.R | 4 + .../tests/testthat/test_RMSNumpress.R | 91 ++ src/main/cpp/MSNumpress.cpp | 782 +++++++++++++++++ src/main/cpp/MSNumpress.hpp | 330 +++++++ src/main/cpp/MSNumpressTest.cpp | 786 +++++++++++++++++ src/main/csharp/MSNumpress.cs | 606 +++++++++++++ src/main/csharp/MSNumpressTest.cs | 333 ++++++++ src/main/java/ms/numpress/IntDecoder.java | 76 ++ src/main/java/ms/numpress/MSNumpress.java | 572 +++++++++++++ src/main/python/MSNumpress.pxd | 61 ++ src/main/python/PyMSNumpress.pyx | 260 ++++++ src/main/python/setup.py | 59 ++ src/main/python/test_pymsnumpress.py | 144 ++++ src/test/java/ms/numpress/MSNumpressTest.java | 375 ++++++++ 36 files changed, 6829 insertions(+) create mode 100644 README.md create mode 100644 meson.build create mode 100755 pom.xml create mode 100644 src/main/R/RMSNumpress/DESCRIPTION create mode 100644 src/main/R/RMSNumpress/LICENSE create mode 100644 src/main/R/RMSNumpress/NAMESPACE create mode 100644 src/main/R/RMSNumpress/R/RcppExports.R create mode 100644 src/main/R/RMSNumpress/inst/LICENSE.md create mode 100644 src/main/R/RMSNumpress/man/RMSNumpress-package.Rd create mode 100644 src/main/R/RMSNumpress/man/decodeLinear.Rd create mode 100644 src/main/R/RMSNumpress/man/decodePic.Rd create mode 100644 src/main/R/RMSNumpress/man/decodeSlof.Rd create mode 100644 src/main/R/RMSNumpress/man/encodeLinear.Rd create mode 100644 src/main/R/RMSNumpress/man/encodePic.Rd create mode 100644 src/main/R/RMSNumpress/man/encodeSlof.Rd create mode 100644 src/main/R/RMSNumpress/man/optimalLinearFixedPoint.Rd create mode 100644 src/main/R/RMSNumpress/man/optimalLinearFixedPointMass.Rd create mode 100644 src/main/R/RMSNumpress/man/optimalSlofFixedPoint.Rd create mode 100644 src/main/R/RMSNumpress/src/MSNumpress.cpp create mode 100644 src/main/R/RMSNumpress/src/RMSNumpress.cpp create mode 100644 src/main/R/RMSNumpress/src/RcppExports.cpp create mode 100644 src/main/R/RMSNumpress/src/include/MSNumpress.hpp create mode 100644 src/main/R/RMSNumpress/tests/testthat.R create mode 100644 src/main/R/RMSNumpress/tests/testthat/test_RMSNumpress.R create mode 100644 src/main/cpp/MSNumpress.cpp create mode 100644 src/main/cpp/MSNumpress.hpp create mode 100644 src/main/cpp/MSNumpressTest.cpp create mode 100644 src/main/csharp/MSNumpress.cs create mode 100644 src/main/csharp/MSNumpressTest.cs create mode 100644 src/main/java/ms/numpress/IntDecoder.java create mode 100644 src/main/java/ms/numpress/MSNumpress.java create mode 100644 src/main/python/MSNumpress.pxd create mode 100644 src/main/python/PyMSNumpress.pyx create mode 100644 src/main/python/setup.py create mode 100644 src/main/python/test_pymsnumpress.py create mode 100644 src/test/java/ms/numpress/MSNumpressTest.java diff --git a/README.md b/README.md new file mode 100644 index 0000000..9c2d39d --- /dev/null +++ b/README.md @@ -0,0 +1,146 @@ +MS Numpress +=========== + +Implementations of two compression schemes for numeric data from mass spectrometers. + +The library provides implementations of 3 different algorithms, +1 designed to compress first order smooth data like retention +time or M/Z arrays, and 2 for compressing non smooth data with +lower requirements on precision like ion count arrays. + +Native implementations and unit test are provided in C++, Java, and C#. We provide Python bindings through [PyMSNumpress](https://pypi.org/project/PyMSNumpress/) which can be installed as follows: + + pip install PyMSNumpress + +We provide R bindings through [RMSNumpress](https://CRAN.R-project.org/package=RMSNumpress) which can be installed as follows: +``` +install.packages("RMSNumpress") +``` + +If you use R via Anaconda, it is appropriate to install `RMSNumpress` through **conda** using the **conda-forge** channel. The feedstock can be found [here](https://github.com/conda-forge/r-rmsnumpress-feedstock). + +First make sure you have the conda-forge channel +``` +conda config --add channels conda-forge +``` + +Then you can install RMSNumpress via conda +``` +conda install r-rmsnumpress +``` + + +### C++ library tests + +For C++, move to `src/main/cpp` and compile and run tests (on LINUX) with + + g++ MSNumpress.cpp MSNumpressTest.cpp -o test && ./test + +### Java (maven) library tests + +Ensure that maven (2.2+) is installed. Then, in this directory, run + + mvn test + +### Python library tests + +Ensure that Cython and the Python headers are installed on your system. Then +move to `src/main/python` and compile and run tests (on LINUX) with + + python setup.py build_ext --inplace + nosetests test_pymsnumpress.py + +### C# library tests + +Ensure that a version of Visual Studio is installed on your system. Then open a Visual Studio Cross Tools Command Prompt, +move to `src\main\csharp` and compile\run tests (on WINDOWS) with + + csc /target:library MSNumpress.cs MSNumpressTest.cs /reference:"C:\Program Files (x86)\Microsoft Visual Studio 14.0\Common7\IDE\PublicAssemblies\Microsoft.VisualStudio.QualityTools.UnitTestFramework.dll" + MSTest /testcontainer:MSNumpress.dll + +NOTE: The example above is for Visual Studio Community 2015 (v14.0). If you use a different version, your path to the unit test reference DLL will be slightly different. + +### R library tests + +Ensure that `Rcpp` and `devtools` is installed. Then move to `src/main/R/RMSNumpress/` and compile and run tests (on LINUX) with + +``` +R -e "Rcpp::compileAttributes(); devtools::test()" +``` + +Numpress Pic +------------ +### MS Numpress positive integer compression + +Intended for ion count data, this compression simply rounds values +to the nearest integer, and stores these integers in a truncated +form which is effective for values relatively close to zero. + + +Numpress Slof +------------- +### MS Numpress short logged float compression + +Also targeting ion count data, this compression takes the natural +logarithm of values, multiplies by a scaling factor and rounds to +the nearest integer. For typical ion count dynamic range these values +fits into two byte integers, so only the two least significant bytes +of the integer are stored. + +The scaling factor can be chosen manually, but the library also contains +a function for retrieving the optimal Slof scaling factor for a given data array. +Since the scaling factor is variable, it is stored as a regular double +precision float first in the encoding, and automatically parsed during decoding. + +Numpress Lin +------------ +### MS Numpress linear prediction compression + +This compression uses a fixed point representation, achieve by +multiplication by a scaling factor and rounding to the nearest integer. +To exploit the assumed linearity of the data, linear prediction is +then used in the following way. + +The first two values are stored without compression as 4 byte integers. +For each following value a linear prediction is made from the two previous +values: + + Xpred = (X(n) - X(n-1)) + X(n) + Xres = Xpred - X(n+1) + +The residual `Xres` is then stored, using the same truncated integer +representation as in Numpress Pic. + +The scaling factor can be chosen manually, but the library also contains +a function for retrieving the optimal Lin scaling factor for a given data array. +Since the scaling factor is variable, it is stored as a regular double +precision float first in the encoding, and automatically parsed during decoding. + +Truncated integer representation +--------------------------------- + +This encoding works on a 4 byte integer, by truncating initial zeros or ones. +If the initial (most significant) half byte is 0x0 or 0xf, the number of such +halfbytes starting from the most significant is stored in a halfbyte. This initial +count is then followed by the rest of the ints halfbytes, in little-endian order. +A count halfbyte c of + + 0 <= c <= 8 is interpreted as an initial c 0x0 halfbytes + 9 <= c <= 15 is interpreted as an initial (c-8) 0xf halfbytes + +Examples: + + int c rest + 0 => 0x8 + -1 => 0xf 0xf + 23 => 0x6 0x7 0x1 + + + +License +------- + +This code is open source. It is dual licenced under the Apache 2.0 license as +well as the 3-clause BSD licence. See the LICENCE-BSD and the LICENCE-APACHE +file for the licences. + diff --git a/meson.build b/meson.build new file mode 100644 index 0000000..441d045 --- /dev/null +++ b/meson.build @@ -0,0 +1,38 @@ +# https://mesonbuild.com/ +project('msnumpres', 'cpp', + version : '0.2.3') + +################################################################################ +# C++ Library and Test +cpp_sources = [ 'src/main/cpp/MSNumpress.cpp' ] +cpp_include = include_directories('src/main/cpp') + +install_headers([ + 'src/main/cpp/MSNumpress.hpp' +], subdir: 'MSNumpress') + +cpp_lib_so = library( + 'numpress', cpp_sources, + include_directories : cpp_include, + install : true) + +cpp_lib_a = static_library( + 'numpress', cpp_sources, + include_directories : cpp_include, + install : true) + + +test('MSNumpressTest', + executable('MSNumpressTest', + 'src/main/cpp/MSNumpressTest.cpp', + include_directories : cpp_include, + link_with : cpp_lib_a), + verbose : false) + +msnumpres_cpp_lib_a = declare_dependency( + include_directories : cpp_include, + link_with : cpp_lib_a) + +msnumpres_cpp_lib_so = declare_dependency( + include_directories : cpp_include, + link_with : cpp_lib_so) diff --git a/pom.xml b/pom.xml new file mode 100755 index 0000000..aad52df --- /dev/null +++ b/pom.xml @@ -0,0 +1,33 @@ + + 4.0.0 + se.lth.immun + MsNumpress + 0.1.21 + 2013 + + + + + junit + junit + 4.9 + test + + + + + + src/main/java + src/test/java + + + maven-compiler-plugin + 2.3.2 + + 1.6 + 1.6 + + + + + diff --git a/src/main/R/RMSNumpress/DESCRIPTION b/src/main/R/RMSNumpress/DESCRIPTION new file mode 100644 index 0000000..7626a17 --- /dev/null +++ b/src/main/R/RMSNumpress/DESCRIPTION @@ -0,0 +1,14 @@ +Package: RMSNumpress +Type: Package +Title: 'Rcpp' Bindings to Native C++ Implementation of MS Numpress +Version: 1.0.1 +Date: 2021-02-04 +Authors@R: c(person("Justin", "Sing", email = "justincsing@gmail.com", role = c("cre","aut")), + person("Johan", "Teleman", email = "johan.teleman@immun.lth.se", role = "aut")) +Description: 'Rcpp' bindings to the native C++ implementation of MS Numpress, that provides two compression schemes for numeric data from mass spectrometers. The library provides implementations of 3 different algorithms, 1 designed to compress first order smooth data like retention time or M/Z arrays, and 2 for compressing non smooth data with lower requirements on precision like ion count arrays. Refer to the publication (Teleman et al., (2014) ) for more details. +License: BSD_3_clause + file LICENSE +Imports: Rcpp (>= 1.0.3) +LinkingTo: Rcpp +Suggests: + testthat +RoxygenNote: 7.0.2 diff --git a/src/main/R/RMSNumpress/LICENSE b/src/main/R/RMSNumpress/LICENSE new file mode 100644 index 0000000..024b33e --- /dev/null +++ b/src/main/R/RMSNumpress/LICENSE @@ -0,0 +1,3 @@ +YEAR: 2020 +COPYRIGHT HOLDER: Justin Sing +ORGANIZATION: University of Toronto diff --git a/src/main/R/RMSNumpress/NAMESPACE b/src/main/R/RMSNumpress/NAMESPACE new file mode 100644 index 0000000..8ad0d37 --- /dev/null +++ b/src/main/R/RMSNumpress/NAMESPACE @@ -0,0 +1,3 @@ +useDynLib(RMSNumpress, .registration=TRUE) +exportPattern("^[[:alpha:]]+") +importFrom(Rcpp, evalCpp) diff --git a/src/main/R/RMSNumpress/R/RcppExports.R b/src/main/R/RMSNumpress/R/RcppExports.R new file mode 100644 index 0000000..164cd39 --- /dev/null +++ b/src/main/R/RMSNumpress/R/RcppExports.R @@ -0,0 +1,189 @@ +# Generated by using Rcpp::compileAttributes() -> do not edit by hand +# Generator token: 10BE3573-1514-4C36-9D1C-5A225CD40393 + +#' optimalLinearFixedPointMass +#' +#' Compute the optimal linear fixed point with a desired m/z accuracy. +#' +#' @note If the desired accuracy cannot be reached without overflowing 64 +#' bit integers, then a negative value is returned. You need to check for +#' this and in that case abandon numpress or use optimalLinearFixedPoint +#' which returns the largest safe value. +#' +#' @param data pointer to array of double to be encoded (need memorycont. repr.) +#' @param mass_acc desired m/z accuracy in Th +#' @return the linear fixed point that satisfies the accuracy requirement (or -1 in case of failure). +optimalLinearFixedPointMass <- function(data, mass_acc) { + .Call(`_RMSNumpress_optimalLinearFixedPointMass`, data, mass_acc) +} + +#' optimalLinearFixedPoint +#' +#' Compute the maximal linear fixed point that prevents integer overflow. +#' +#' @param data pointer to array of double to be encoded (need memorycont. repr.) +#' @return the linear fixed point safe to use +#' +optimalLinearFixedPoint <- function(data) { + .Call(`_RMSNumpress_optimalLinearFixedPoint`, data) +} + +#' optimalSlofFixedPoint +#' +#' Compute the maximal natural logarithm fixed point that prevents integer overflow. +#' +#' @param data pointer to array of double to be encoded (need memorycont. repr.) +#' @return the slof fixed point safe to use +#' +optimalSlofFixedPoint <- function(data) { + .Call(`_RMSNumpress_optimalSlofFixedPoint`, data) +} + +#' encodeLinear +#' +#' Encodes the doubles in data by first using a \cr +#' - lossy conversion to a 4 byte 5 decimal fixed point representation \cr +#' - storing the residuals from a linear prediction after first two values \cr +#' - encoding by encodeInt (see above) \cr +#' +#' The resulting binary is maximally 8 + dataSize * 5 bytes, but much less if the +#' data is reasonably smooth on the first order. +#' +#' This encoding is suitable for typical m/z or retention time binary arrays. +#' On a test set, the encoding was empirically show to be accurate to at least 0.002 ppm. +#' +#' @param data pointer to array of double to be encoded (need memorycont. repr.) +#' @param fixedPoint the scaling factor used for getting the fixed point repr. +#' This is stored in the binary and automatically extracted +#' on decoding (see optimalLinearFixedPoint or optimalLinearFixedPointMass) +#' @return the number of encoded bytes +#' +#' @seealso [\code{\link{decodeLinear}}] +#' +#' @examples +#' \dontrun{ +#' ## Retention time array +#' rt_array <- c(4313.0, 4316.4, 4319.8, 4323.2, 4326.6, 4330.1) +#' ## encode retention time array +#' rt_encoded <- encodeLinear(rt_array, 500) +#' #> [1] 40 7f 40 00 00 00 00 00 d4 e7 20 00 78 ee 20 00 88 86 23 +#' } +encodeLinear <- function(data, fixedPoint) { + .Call(`_RMSNumpress_encodeLinear`, data, fixedPoint) +} + +#' decodeLinear +#' +#' Decodes data encoded by encodeLinear. +#' +#' result vector guaranteed to be shorter or equal to (|data| - 8) * 2 +#' +#' Note that this method may throw a const char* if it deems the input data to be corrupt, i.e. +#' that the last encoded int does not use the last byte in the data. In addition the last encoded +#' int need to use either the last halfbyte, or the second last followed by a 0x0 halfbyte. +#' +#' @param data pointer to array of bytes to be decoded (need memorycont. repr.) +#' @return the number of decoded doubles, or -1 if dataSize < 4 or 4 < dataSize < 8 +#' +#' @seealso [\code{\link{encodeLinear}}] +#' +#' @examples +#' \dontrun{ +#' ## Retention time data that is encoded with encodeLinear and is zlib compressed +#' ### NOTE: For the sake of this example, I have broken the raw vector into several parts +#' ### to avoid Rd line widths (>100 characters) issues with CRAN build checks. +#' rt_raw1 <- c("78", "9c", "73", "50", "61", "00", "83", "aa", "15", "0c", "0c", "73", "80") +#' rt_raw2 <- c("b8", "a3", "5d", "fe", "47", "07", "84", "28", "fc", "8f", "c4", "40", "e5") +#' rt_raw3 <- c("61", "51", "84", "a9", "85", "08", "e1", "06", "00", "06", "be", "41", "cf") +#' ## Add all character representation of raw data back together and convert back to hex raw vector +#' rt_blob <- as.raw(as.hexmode(c(rt_raw1, rt_raw2, rt_raw3 ))) +#' ## Decompress blob +#' rt_blob_uncompressed <- memDecompress(rt_blob, type = "gzip", asChar = FALSE) +#' ## Decode to rentention time double values +#' rt_array <- decodeLinear(rt_blob_uncompressed) +#' } +decodeLinear <- function(data) { + .Call(`_RMSNumpress_decodeLinear`, data) +} + +#' encodeSlof +#' +#' Encodes ion counts by taking the natural logarithm, and storing a +#' fixed point representation of this. This is calculated as +#' +#' unsigned short fp = log(d + 1) * fixedPoint + 0.5 +#' +#' the result vector is exactly |data| * 2 + 8 bytes long +#' +#' @param data pointer to array of double to be encoded (need memorycont. repr.) +#' @param fixedPoint fixed point to use for encoding (see optimalSlofFixedPoint) +#' @return the number of encoded bytes +#' +#' @seealso [\code{\link{decodeSlof}}] +encodeSlof <- function(data, fixedPoint) { + .Call(`_RMSNumpress_encodeSlof`, data, fixedPoint) +} + +#' decodeSlof +#' +#' Decodes data encoded by encodeSlof +#' +#' The return will include exactly (|data| - 8) / 2 doubles. +#' +#' Note that this method may throw a const char* if it deems the input data to be corrupt. +#' +#' @param data pointer to array of bytes to be decoded (need memorycont. repr.) +#' @return the number of decoded doubles +#' +#' @seealso [\code{\link{encodeSlof}}] +#' @examples +#' \dontrun{ +#' ## Intensity array to encode +#' ### NOTE: For the sake of this example, I have broken the intensity vector into several parts +#' ### to avoid Rd line widths (>100 characters) issues with CRAN build checks. +#' int_array1 <- c(0.71773432, 0.43443741, 1.71883610, 0.13220307, 0.90664242) +#' int_array2 <- c(0.00000000, 0.00000000, 0.64213755, 0.43443741, 0.47221479) +#' ## Comcatenate into one intensity array +#' int_array <- c(int_array1, int_array2) +#' ## Encode intensity array using encodeSlof +#' int_encode <- encodeSlof( int_array, 16 ) +#' } +decodeSlof <- function(data) { + .Call(`_RMSNumpress_decodeSlof`, data) +} + +#' encodePic +#' +#' Encodes ion counts by simply rounding to the nearest 4 byte integer, +#' and compressing each integer with encodeInt. +#' +#' The handleable range is therefore 0 -> 4294967294. +#' The resulting binary is maximally dataSize * 5 bytes, but much less if the +#' data is close to 0 on average. +#' +#' @param data pointer to array of double to be encoded (need memorycont. repr.) +#' @return the number of encoded bytes +#' +#' @seealso [\code{\link{decodePic}}] +encodePic <- function(data) { + .Call(`_RMSNumpress_encodePic`, data) +} + +#' decodePic +#' +#' Decodes data encoded by encodePic +#' +#' result vector guaranteed to be shorter of equal to |data| * 2 +#' +#' Note that this method may throw a const char* if it deems the input data to be corrupt, i.e. +#' that the last encoded int does not use the last byte in the data. In addition the last encoded +#' int need to use either the last halfbyte, or the second last followed by a 0x0 halfbyte. +#' +#' @param data pointer to array of bytes to be decoded (need memorycont. repr.) +#' @return the number of decoded doubles +#' +#' @seealso [\code{\link{encodePic}}] +decodePic <- function(data) { + .Call(`_RMSNumpress_decodePic`, data) +} + diff --git a/src/main/R/RMSNumpress/inst/LICENSE.md b/src/main/R/RMSNumpress/inst/LICENSE.md new file mode 100644 index 0000000..8cc3f91 --- /dev/null +++ b/src/main/R/RMSNumpress/inst/LICENSE.md @@ -0,0 +1,31 @@ +# BSD_3_clause License + +Copyright (c) 2020 Justin Sing + +Redistribution and use in source and binary forms, with or without +modification, are permitted provided that the following conditions are +met: + + Redistributions of source code must retain the above copyright + notice, this list of conditions and the following disclaimer. + + Redistributions in binary form must reproduce the above copyright + notice, this list of conditions and the following disclaimer in + the documentation and/or other materials provided with the + distribution. + + Neither the name of the nor the names of its + contributors may be used to endorse or promote products derived + from this software without specific prior written permission. + +THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS +"AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT +LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR +A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT +HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, +SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT +LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, +DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY +THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT +(INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE +OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. diff --git a/src/main/R/RMSNumpress/man/RMSNumpress-package.Rd b/src/main/R/RMSNumpress/man/RMSNumpress-package.Rd new file mode 100644 index 0000000..8d9f482 --- /dev/null +++ b/src/main/R/RMSNumpress/man/RMSNumpress-package.Rd @@ -0,0 +1,131 @@ +\name{RMSNumpress-package} +\alias{RMSNumpress-package} +\alias{RMSNumpress} +\docType{package} +\title{ + Rcpp bindings to native C++ implementation of MS Numpress +} +\description{ +MS Numpress + +=========== + + +Implementations of two compression schemes for numeric data from mass spectrometers. + +The library provides implementations of 3 different algorithms, +1 designed to compress first order smooth data like retention +time or M/Z arrays, and 2 for compressing non smooth data with +lower requirements on precision like ion count arrays. + +Numpress Pic + +=========== + +MS Numpress positive integer compression + +Intended for ion count data, this compression simply rounds values to the nearest integer, and stores these integers in a truncated form which is effective for values relatively close to zero. + +Numpress Slof + +=========== + +MS Numpress short logged float compression + +Also targeting ion count data, this compression takes the natural logarithm of values, multiplies by a scaling factor and rounds to the nearest integer. For typical ion count dynamic range these values fits into two byte integers, so only the two least significant bytes of the integer are stored. + +The scaling factor can be chosen manually, but the library also contains a function for retrieving the optimal Slof scaling factor for a given data array. Since the scaling factor is variable, it is stored as a regular double precision float first in the encoding, and automatically parsed during decoding. + +Numpress Lin + +=========== + +MS Numpress linear prediction compression + +This compression uses a fixed point representation, achieve by multiplication by a scaling factor and rounding to the nearest integer. To exploit the assumed linearity of the data, linear prediction is then used in the following way. + +The first two values are stored without compression as 4 byte integers. For each following value a linear prediction is made from the two previous values: + + + +Xpred = (X(n) - X(n-1)) + X(n) + +Xres = Xpred - X(n+1) + + + + +The residual Xres is then stored, using the same truncated integer representation as in Numpress Pic. + +The scaling factor can be chosen manually, but the library also contains a function for retrieving the optimal Lin scaling factor for a given data array. Since the scaling factor is variable, it is stored as a regular double precision float first in the encoding, and automatically parsed during decoding. + +Truncated integer representation + +=========== + +This encoding works on a 4 byte integer, by truncating initial zeros or ones. If the initial (most significant) half byte is 0x0 or 0xf, the number of such halfbytes starting from the most significant is stored in a halfbyte. This initial count is then followed by the rest of the ints halfbytes, in little-endian order. A count halfbyte c of + + + +0 <= c <= 8 is interpreted as an initial c 0x0 halfbytes + +9 <= c <= 15 is interpreted as an initial (c-8) 0xf halfbytes + + + + +Examples: + +int c rest + +0 => 0x8 + +-1 => 0xf 0xf + +23 => 0x6 0x7 0x1 + + + +} +\author{ +Maintainer: Justin Sing +} +\references{ + See: https://github.com/ms-numpress/ms-numpress +} +\keyword{ package } +\seealso{ + \code{\link{encodeLinear}}, + \code{\link{decodeLinear}}, + \code{\link{encodeSlof}}, + \code{\link{decodeSlof}}, + \code{\link{encodePic}}, + \code{\link{decodePic}}, + \code{\link{optimalLinearFixedPoint}}, + \code{\link{optimalSlofFixedPoint}}, + \code{\link{optimalLinearFixedPointMass}}, +} +\examples{ + \dontrun{ + # Encode Numpress Linear + ## Retention time array + rt_array <- c(4313.0, 4316.4, 4319.8, 4323.2, 4326.6, 4330.1) + ## encode retention time array + rt_encoded <- encodeLinear(rt_array, 500) + #> [1] 40 7f 40 00 00 00 00 00 d4 e7 20 00 78 ee 20 00 88 86 23 + + # Decode Numpress Linear + ## Retention time data that is encoded with encodeLinear and is zlib compressed + ### NOTE: For the sake of this example, I have broken the raw vector into several parts + ### to avoid Rd line widths (>100 characters) issues with CRAN build checks. + rt_raw1 <- c("78", "9c", "73", "50", "61", "00", "83", "aa", "15", "0c", "0c", "73", "80") + rt_raw2 <- c("b8", "a3", "5d", "fe", "47", "07", "84", "28", "fc", "8f", "c4", "40", "e5") + rt_raw3 <- c("61", "51", "84", "a9", "85", "08", "e1", "06", "00", "06", "be", "41", "cf") + ## Add all character representation of raw data back together and convert back to hex raw vector + rt_blob <- as.raw(as.hexmode(c(rt_raw1, rt_raw2, rt_raw3 ))) + ## Decompress blob + rt_blob_uncompressed <- as.raw(Rcompression::uncompress( rt_blob, asText = FALSE )) + ## Decode to rentention time double values + rt_array <- decodeLinear(rt_blob_uncompressed) + } +} diff --git a/src/main/R/RMSNumpress/man/decodeLinear.Rd b/src/main/R/RMSNumpress/man/decodeLinear.Rd new file mode 100644 index 0000000..e407b5e --- /dev/null +++ b/src/main/R/RMSNumpress/man/decodeLinear.Rd @@ -0,0 +1,43 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/RcppExports.R +\name{decodeLinear} +\alias{decodeLinear} +\title{decodeLinear} +\usage{ +decodeLinear(data) +} +\arguments{ +\item{data}{pointer to array of bytes to be decoded (need memorycont. repr.)} +} +\value{ +the number of decoded doubles, or -1 if dataSize < 4 or 4 < dataSize < 8 +} +\description{ +Decodes data encoded by encodeLinear. +} +\details{ +result vector guaranteed to be shorter or equal to (|data| - 8) * 2 + +Note that this method may throw a const char* if it deems the input data to be corrupt, i.e. +that the last encoded int does not use the last byte in the data. In addition the last encoded +int need to use either the last halfbyte, or the second last followed by a 0x0 halfbyte. +} +\examples{ +\dontrun{ +## Retention time data that is encoded with encodeLinear and is zlib compressed +### NOTE: For the sake of this example, I have broken the raw vector into several parts +### to avoid Rd line widths (>100 characters) issues with CRAN build checks. +rt_raw1 <- c("78", "9c", "73", "50", "61", "00", "83", "aa", "15", "0c", "0c", "73", "80") +rt_raw2 <- c("b8", "a3", "5d", "fe", "47", "07", "84", "28", "fc", "8f", "c4", "40", "e5") +rt_raw3 <- c("61", "51", "84", "a9", "85", "08", "e1", "06", "00", "06", "be", "41", "cf") +## Add all character representation of raw data back together and convert back to hex raw vector +rt_blob <- as.raw(as.hexmode(c(rt_raw1, rt_raw2, rt_raw3 ))) +## Decompress blob +rt_blob_uncompressed <- as.raw(Rcompression::uncompress( rt_blob, asText = FALSE )) +## Decode to rentention time double values +rt_array <- decodeLinear(rt_blob_uncompressed) +} +} +\seealso{ +[\code{\link{encodeLinear}}] +} diff --git a/src/main/R/RMSNumpress/man/decodePic.Rd b/src/main/R/RMSNumpress/man/decodePic.Rd new file mode 100644 index 0000000..5b3c095 --- /dev/null +++ b/src/main/R/RMSNumpress/man/decodePic.Rd @@ -0,0 +1,27 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/RcppExports.R +\name{decodePic} +\alias{decodePic} +\title{decodePic} +\usage{ +decodePic(data) +} +\arguments{ +\item{data}{pointer to array of bytes to be decoded (need memorycont. repr.)} +} +\value{ +the number of decoded doubles +} +\description{ +Decodes data encoded by encodePic + + result vector guaranteed to be shorter of equal to |data| * 2 +} +\details{ +Note that this method may throw a const char* if it deems the input data to be corrupt, i.e. +that the last encoded int does not use the last byte in the data. In addition the last encoded + int need to use either the last halfbyte, or the second last followed by a 0x0 halfbyte. +} +\seealso{ +[\code{\link{encodePic}}] +} diff --git a/src/main/R/RMSNumpress/man/decodeSlof.Rd b/src/main/R/RMSNumpress/man/decodeSlof.Rd new file mode 100644 index 0000000..32c561d --- /dev/null +++ b/src/main/R/RMSNumpress/man/decodeSlof.Rd @@ -0,0 +1,38 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/RcppExports.R +\name{decodeSlof} +\alias{decodeSlof} +\title{decodeSlof} +\usage{ +decodeSlof(data) +} +\arguments{ +\item{data}{pointer to array of bytes to be decoded (need memorycont. repr.)} +} +\value{ +the number of decoded doubles +} +\description{ +Decodes data encoded by encodeSlof + + The return will include exactly (|data| - 8) / 2 doubles. +} +\details{ +Note that this method may throw a const char* if it deems the input data to be corrupt. +} +\examples{ +\dontrun{ +## Intensity array to encode +### NOTE: For the sake of this example, I have broken the intensity vector into several parts +### to avoid Rd line widths (>100 characters) issues with CRAN build checks. +int_array1 <- c(0.71773432, 0.43443741, 1.71883610, 0.13220307, 0.90664242) +int_array2 <- c(0.00000000, 0.00000000, 0.64213755, 0.43443741, 0.47221479) +## Comcatenate into one intensity array +int_array <- c(int_array1, int_array2) +## Encode intensity array using encodeSlof +int_encode <- encodeSlof( int_array, 16 ) +} +} +\seealso{ +[\code{\link{encodeSlof}}] +} diff --git a/src/main/R/RMSNumpress/man/encodeLinear.Rd b/src/main/R/RMSNumpress/man/encodeLinear.Rd new file mode 100644 index 0000000..1e95127 --- /dev/null +++ b/src/main/R/RMSNumpress/man/encodeLinear.Rd @@ -0,0 +1,42 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/RcppExports.R +\name{encodeLinear} +\alias{encodeLinear} +\title{encodeLinear} +\usage{ +encodeLinear(data, fixedPoint) +} +\arguments{ +\item{data}{pointer to array of double to be encoded (need memorycont. repr.)} + +\item{fixedPoint}{the scaling factor used for getting the fixed point repr. +This is stored in the binary and automatically extracted +on decoding (see optimalLinearFixedPoint or optimalLinearFixedPointMass)} +} +\value{ +the number of encoded bytes +} +\description{ +Encodes the doubles in data by first using a \cr + - lossy conversion to a 4 byte 5 decimal fixed point representation \cr + - storing the residuals from a linear prediction after first two values \cr + - encoding by encodeInt (see above) \cr + + The resulting binary is maximally 8 + dataSize * 5 bytes, but much less if the + data is reasonably smooth on the first order. + + This encoding is suitable for typical m/z or retention time binary arrays. + On a test set, the encoding was empirically show to be accurate to at least 0.002 ppm. +} +\examples{ +\dontrun{ +## Retention time array +rt_array <- c(4313.0, 4316.4, 4319.8, 4323.2, 4326.6, 4330.1) +## encode retention time array +rt_encoded <- encodeLinear(rt_array, 500) +#> [1] 40 7f 40 00 00 00 00 00 d4 e7 20 00 78 ee 20 00 88 86 23 +} +} +\seealso{ +[\code{\link{decodeLinear}}] +} diff --git a/src/main/R/RMSNumpress/man/encodePic.Rd b/src/main/R/RMSNumpress/man/encodePic.Rd new file mode 100644 index 0000000..8aadb6a --- /dev/null +++ b/src/main/R/RMSNumpress/man/encodePic.Rd @@ -0,0 +1,26 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/RcppExports.R +\name{encodePic} +\alias{encodePic} +\title{encodePic} +\usage{ +encodePic(data) +} +\arguments{ +\item{data}{pointer to array of double to be encoded (need memorycont. repr.)} +} +\value{ +the number of encoded bytes +} +\description{ +Encodes ion counts by simply rounding to the nearest 4 byte integer, +and compressing each integer with encodeInt. +} +\details{ +The handleable range is therefore 0 -> 4294967294. +The resulting binary is maximally dataSize * 5 bytes, but much less if the + data is close to 0 on average. +} +\seealso{ +[\code{\link{decodePic}}] +} diff --git a/src/main/R/RMSNumpress/man/encodeSlof.Rd b/src/main/R/RMSNumpress/man/encodeSlof.Rd new file mode 100644 index 0000000..156c16a --- /dev/null +++ b/src/main/R/RMSNumpress/man/encodeSlof.Rd @@ -0,0 +1,28 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/RcppExports.R +\name{encodeSlof} +\alias{encodeSlof} +\title{encodeSlof} +\usage{ +encodeSlof(data, fixedPoint) +} +\arguments{ +\item{data}{pointer to array of double to be encoded (need memorycont. repr.)} + +\item{fixedPoint}{fixed point to use for encoding (see optimalSlofFixedPoint)} +} +\value{ +the number of encoded bytes +} +\description{ +Encodes ion counts by taking the natural logarithm, and storing a + fixed point representation of this. This is calculated as + + unsigned short fp = log(d + 1) * fixedPoint + 0.5 +} +\details{ +the result vector is exactly |data| * 2 + 8 bytes long +} +\seealso{ +[\code{\link{decodeSlof}}] +} diff --git a/src/main/R/RMSNumpress/man/optimalLinearFixedPoint.Rd b/src/main/R/RMSNumpress/man/optimalLinearFixedPoint.Rd new file mode 100644 index 0000000..3905c36 --- /dev/null +++ b/src/main/R/RMSNumpress/man/optimalLinearFixedPoint.Rd @@ -0,0 +1,17 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/RcppExports.R +\name{optimalLinearFixedPoint} +\alias{optimalLinearFixedPoint} +\title{optimalLinearFixedPoint} +\usage{ +optimalLinearFixedPoint(data) +} +\arguments{ +\item{data}{pointer to array of double to be encoded (need memorycont. repr.)} +} +\value{ +the linear fixed point safe to use +} +\description{ +Compute the maximal linear fixed point that prevents integer overflow. +} diff --git a/src/main/R/RMSNumpress/man/optimalLinearFixedPointMass.Rd b/src/main/R/RMSNumpress/man/optimalLinearFixedPointMass.Rd new file mode 100644 index 0000000..f0f25f4 --- /dev/null +++ b/src/main/R/RMSNumpress/man/optimalLinearFixedPointMass.Rd @@ -0,0 +1,25 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/RcppExports.R +\name{optimalLinearFixedPointMass} +\alias{optimalLinearFixedPointMass} +\title{optimalLinearFixedPointMass} +\usage{ +optimalLinearFixedPointMass(data, mass_acc) +} +\arguments{ +\item{data}{pointer to array of double to be encoded (need memorycont. repr.)} + +\item{mass_acc}{desired m/z accuracy in Th} +} +\value{ +the linear fixed point that satisfies the accuracy requirement (or -1 in case of failure). +} +\description{ +Compute the optimal linear fixed point with a desired m/z accuracy. +} +\note{ +If the desired accuracy cannot be reached without overflowing 64 +bit integers, then a negative value is returned. You need to check for + this and in that case abandon numpress or use optimalLinearFixedPoint + which returns the largest safe value. +} diff --git a/src/main/R/RMSNumpress/man/optimalSlofFixedPoint.Rd b/src/main/R/RMSNumpress/man/optimalSlofFixedPoint.Rd new file mode 100644 index 0000000..1500997 --- /dev/null +++ b/src/main/R/RMSNumpress/man/optimalSlofFixedPoint.Rd @@ -0,0 +1,17 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/RcppExports.R +\name{optimalSlofFixedPoint} +\alias{optimalSlofFixedPoint} +\title{optimalSlofFixedPoint} +\usage{ +optimalSlofFixedPoint(data) +} +\arguments{ +\item{data}{pointer to array of double to be encoded (need memorycont. repr.)} +} +\value{ +the slof fixed point safe to use +} +\description{ +Compute the maximal natural logarithm fixed point that prevents integer overflow. +} diff --git a/src/main/R/RMSNumpress/src/MSNumpress.cpp b/src/main/R/RMSNumpress/src/MSNumpress.cpp new file mode 100644 index 0000000..e3b6e17 --- /dev/null +++ b/src/main/R/RMSNumpress/src/MSNumpress.cpp @@ -0,0 +1,802 @@ +/* + MSNumpress.cpp + johan.teleman@immun.lth.se + + Copyright 2013 Johan Teleman + + Licensed under the Apache License, Version 2.0 (the "License"); + you may not use this file except in compliance with the License. + You may obtain a copy of the License at + + http://www.apache.org/licenses/LICENSE-2.0 + + Unless required by applicable law or agreed to in writing, software + distributed under the License is distributed on an "AS IS" BASIS, + WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + See the License for the specific language governing permissions and + limitations under the License. + */ + +#include +using namespace Rcpp; +using namespace std; + +#include +#include +#include +#include +#include "include/MSNumpress.hpp" +#include + +namespace ms { +namespace numpress { +namespace MSNumpress { + +using std::cout; +using std::cerr; +using std::endl; +using std::min; +using std::max; +using std::abs; + +// This is only valid on systems were ints use more bytes than chars... + +const int ONE = 1; +static bool is_little_endian() { + return *((char*)&(ONE)) == 1; +} +bool IS_LITTLE_ENDIAN = is_little_endian(); + + + +///////////////////////////////////////////////////////////// + +static void encodeFixedPoint( + double fixedPoint, + unsigned char *result +) { + int i; + unsigned char *fp = (unsigned char*)&fixedPoint; + for (i=0; i<8; i++) { + result[i] = fp[IS_LITTLE_ENDIAN ? (7-i) : i]; + } +} + + + +static double decodeFixedPoint( + const unsigned char *data +) { + int i; + double fixedPoint; + unsigned char *fp = (unsigned char*)&fixedPoint; + + for (i=0; i<8; i++) { + fp[i] = data[IS_LITTLE_ENDIAN ? (7-i) : i]; + } + + return fixedPoint; +} + +///////////////////////////////////////////////////////////// + +/** + * Encodes the int x as a number of halfbytes in res. + * res_length is incremented by the number of halfbytes, + * which will be 1 <= n <= 9 + */ +static void encodeInt( + const unsigned int x, + unsigned char* res, + size_t *res_length +) { + // get the bit pattern of a signed int x_inp + unsigned int m; + unsigned char i, l; // numbers between 0 and 9 + + unsigned int mask = 0xf0000000; + unsigned int init = x & mask; + + if (init == 0) { + l = 8; + for (i=0; i<8; i++) { + m = mask >> (4*i); + if ((x & m) != 0) { + l = i; + break; + } + } + res[0] = l; + for (i=l; i<8; i++) { + res[1+i-l] = static_cast( x >> (4*(i-l)) ); + } + *res_length += 1+8-l; + + } else if (init == mask) { + l = 7; + for (i=0; i<8; i++) { + m = mask >> (4*i); + if ((x & m) != m) { + l = i; + break; + } + } + res[0] = l + 8; + for (i=l; i<8; i++) { + res[1+i-l] = static_cast( x >> (4*(i-l)) ); + } + *res_length += 1+8-l; + + } else { + res[0] = 0; + for (i=0; i<8; i++) { + res[1+i] = static_cast( x >> (4*i) ); + } + *res_length += 9; + + } +} + + + +/** + * Decodes an int from the half bytes in bp. Lossless reverse of encodeInt + * + * @param data ptr to the char data to decode + * @param di position in the char data array to start decoding (will be advanced) + * @param max_di size of data array + * @param half helper variable (do not change between multiple calls) + * @param res result (a 32 bit integer) + * + * @note the helper variable indicates whether we look at the first half byte + * or second half byte of the current data (thus whether to interpret the first + * half byte of data[*di] or the second half byte). + * + */ +static void decodeInt( + const unsigned char *data, + size_t *di, + size_t max_di, + size_t *half, + unsigned int *res +) { + size_t n, i; + unsigned int mask, m; + unsigned char head; + unsigned char hb; + + // Extract the first half byte, specifying the number of leading zero half + // bytes of the final integer. + // If half is zero, we look at the first half byte, otherwise we look at + // the second (lower) half byte and advance the counter to the next char. + if (*half == 0) { + head = data[*di] >> 4; + } else { + head = data[*di] & 0xf; + (*di)++; + } + + *half = 1-(*half); // switch to other half byte + *res = 0; + + if (head <= 8) { + n = head; + } else { // we have n leading ones, fill n half bytes in res with 0xf + n = head - 8; + mask = 0xf0000000; + for (i=0; i> (4*i); + *res = *res | m; + } + } + + if (n == 8) { + return; + } + + if (*di + ((8 - n) - (1 - *half)) / 2 >= max_di) { + throw "[MSNumpress::decodeInt] Corrupt input data! "; + } + + for (i=n; i<8; i++) { + if (*half == 0) { + hb = data[*di] >> 4; + } else { + hb = data[*di] & 0xf; + (*di)++; + } + *res = *res | ( static_cast(hb) << ((i-n)*4)); + *half = 1 - (*half); + } +} + + + + +///////////////////////////////////////////////////////////// + +double optimalLinearFixedPointMass( + const double *data, + size_t dataSize, + double mass_acc +) { + if (dataSize < 3) return 0; // we just encode the first two points as floats + + // We calculate the maximal fixedPoint we need to achieve a specific mass + // accuracy. Note that the maximal error we will make by encoding as int is + // 0.5 due to rounding errors. + double maxFP = 0.5 / mass_acc; + + // There is a maximal value for the FP given by the int length (32bit) + // which means we cannot choose a value higher than that. In case we cannot + // achieve the desired accuracy, return failure (-1). + double maxFP_overflow = optimalLinearFixedPoint(data, dataSize); + if (maxFP > maxFP_overflow) return -1; + + return maxFP; +} + +double optimalLinearFixedPoint( + const double *data, + size_t dataSize +) { + /* + * safer impl - apparently not needed though + * + if (dataSize == 0) return 0; + + double maxDouble = 0; + double x; + + for (size_t i=0; i(data[0] * fixedPoint + 0.5); + for (i=0; i<4; i++) { + result[8+i] = (ints[1] >> (i*8)) & 0xff; + } + + if (dataSize == 1) return 12; + + ints[2] = static_cast(data[1] * fixedPoint + 0.5); + for (i=0; i<4; i++) { + result[12+i] = (ints[2] >> (i*8)) & 0xff; + } + + halfByteCount = 0; + ri = 16; + + for (i=2; i LLONG_MAX ) { + throw "[MSNumpress::encodeLinear] Next number overflows LLONG_MAX."; + } + + ints[2] = static_cast(data[i] * fixedPoint + 0.5); + extrapol = ints[1] + (ints[1] - ints[0]); + + if (THROW_ON_OVERFLOW && + ( ints[2] - extrapol > INT_MAX + || ints[2] - extrapol < INT_MIN )) { + throw "[MSNumpress::encodeLinear] Cannot encode a number that exceeds the bounds of [-INT_MAX, INT_MAX]."; + } + + diff = static_cast(ints[2] - extrapol); + //printf("%lu %lu %lu, extrapol: %ld diff: %d \n", ints[0], ints[1], ints[2], extrapol, diff); + encodeInt( + static_cast(diff), + &halfBytes[halfByteCount], + &halfByteCount + ); + /* + printf("%d (%d): ", diff, (int)halfByteCount); + for (size_t j=0; j( + (halfBytes[hbi-1] << 4) | (halfBytes[hbi] & 0xf) + ); + //printf("%x \n", result[ri]); + ri++; + } + if (halfByteCount % 2 != 0) { + halfBytes[0] = halfBytes[halfByteCount-1]; + halfByteCount = 1; + } else { + halfByteCount = 0; + } + } + if (halfByteCount == 1) { + result[ri] = static_cast(halfBytes[0] << 4); + ri++; + } + return ri; +} + +size_t decodeLinear( + const unsigned char *data, + const size_t dataSize, + double *result +) { + size_t i; + size_t ri = 0; + unsigned int init, buff; + int diff; + long long ints[3]; + //double d; + size_t di; + size_t half; + long long extrapol; + long long y; + double fixedPoint; + + //printf("Decoding %d bytes with fixed point %f\n", (int)dataSize, fixedPoint); + + if (dataSize == 8) return 0; + + if (dataSize < 8) + throw "[MSNumpress::decodeLinear] Corrupt input data: not enough bytes to read fixed point! "; + + fixedPoint = decodeFixedPoint(data); + + + if (dataSize < 12) + throw "[MSNumpress::decodeLinear] Corrupt input data: not enough bytes to read first value! "; + + ints[1] = 0; + for (i=0; i<4; i++) { + ints[1] = ints[1] | ((0xff & (init = data[8+i])) << (i*8)); + } + result[0] = ints[1] / fixedPoint; + + if (dataSize == 12) return 1; + if (dataSize < 16) + throw "[MSNumpress::decodeLinear] Corrupt input data: not enough bytes to read second value! "; + + ints[2] = 0; + for (i=0; i<4; i++) { + ints[2] = ints[2] | ((0xff & (init = data[12+i])) << (i*8)); + } + result[1] = ints[2] / fixedPoint; + + half = 0; + ri = 2; + di = 16; + + //printf(" di ri half int[0] int[1] extrapol diff\n"); + + while (di < dataSize) { + if (di == (dataSize - 1) && half == 1) { + if ((data[di] & 0xf) == 0x0) { + break; + } + } + //printf("%7d %7d %7d %lu %lu %ld", di, ri, half, ints[0], ints[1], extrapol); + + ints[0] = ints[1]; + ints[1] = ints[2]; + decodeInt(data, &di, dataSize, &half, &buff); + diff = static_cast(buff); + + extrapol = ints[1] + (ints[1] - ints[0]); + y = extrapol + diff; + //printf(" %d \n", diff); + result[ri++] = y / fixedPoint; + ints[2] = y; + } + + return ri; +} + + +void encodeLinear( + const std::vector &data, + std::vector &result, + double fixedPoint +) { + size_t dataSize = data.size(); + result.resize(dataSize * 5 + 8); + size_t encodedLength = encodeLinear(&data[0], dataSize, &result[0], fixedPoint); + result.resize(encodedLength); +} + +void decodeLinear( + const std::vector &data, + std::vector &result +) { + size_t dataSize = data.size(); + result.resize((dataSize - 8) * 2); + size_t decodedLength = decodeLinear(&data[0], dataSize, &result[0]); + result.resize(decodedLength); +} + +///////////////////////////////////////////////////////////// + + +size_t encodeSafe( + const double *data, + const size_t dataSize, + unsigned char *result +) { + size_t i, j, ri = 0; + double latest[3]; + double extrapol, diff; + const unsigned char *fp; + + //printf("d0 d1 d2 extrapol diff\n"); + + if (dataSize == 0) return ri; + + latest[1] = data[0]; + fp = (unsigned char*)data; + for (i=0; i<8; i++) { + result[ri++] = fp[IS_LITTLE_ENDIAN ? (7-i) : i]; + } + + if (dataSize == 1) return ri; + + latest[2] = data[1]; + fp = (unsigned char*)&(data[1]); + for (i=0; i<8; i++) { + result[ri++] = fp[IS_LITTLE_ENDIAN ? (7-i) : i]; + } + + fp = (unsigned char*)&diff; + for (i=2; i INT_MAX || data[i] < -0.5) ){ + throw "[MSNumpress::encodePic] Cannot use Pic to encode a number larger than INT_MAX or smaller than 0."; + } + x = static_cast(data[i] + 0.5); + //printf("%d %d %d, extrapol: %d diff: %d \n", ints[0], ints[1], ints[2], extrapol, diff); + encodeInt(x, &halfBytes[halfByteCount], &halfByteCount); + + for (hbi=1; hbi < halfByteCount; hbi+=2) { + result[ri] = static_cast( + (halfBytes[hbi-1] << 4) | (halfBytes[hbi] & 0xf) + ); + //printf("%x \n", result[ri]); + ri++; + } + if (halfByteCount % 2 != 0) { + halfBytes[0] = halfBytes[halfByteCount-1]; + halfByteCount = 1; + } else { + halfByteCount = 0; + } + } + if (halfByteCount == 1) { + result[ri] = static_cast(halfBytes[0] << 4); + ri++; + } + return ri; +} + + + +size_t decodePic( + const unsigned char *data, + const size_t dataSize, + double *result +) { + size_t ri; + unsigned int x; + size_t di; + size_t half; + + //printf("ri di half dSize count\n"); + + half = 0; + ri = 0; + di = 0; + + while (di < dataSize) { + if (di == (dataSize - 1) && half == 1) { + if ((data[di] & 0xf) == 0x0) { + break; + } + } + + decodeInt(&data[0], &di, dataSize, &half, &x); + + //printf("%7d %7d %7d %7d %7d\n", ri, di, half, dataSize, count); + + //printf("count: %d \n", count); + result[ri++] = static_cast(x); + } + + return ri; +} + + + +void encodePic( + const std::vector &data, + std::vector &result +) { + size_t dataSize = data.size(); + result.resize(dataSize * 5); + size_t encodedLength = encodePic(&data[0], dataSize, &result[0]); + result.resize(encodedLength); +} + + + +void decodePic( + const std::vector &data, + std::vector &result +) { + size_t dataSize = data.size(); + result.resize(dataSize * 2); + size_t decodedLength = decodePic(&data[0], dataSize, &result[0]); + result.resize(decodedLength); +} + + +///////////////////////////////////////////////////////////// + + +double optimalSlofFixedPoint( + const double *data, + size_t dataSize +) { + if (dataSize == 0) return 0; + + double maxDouble = 1; + double x; + double fp; + + for (size_t i=0; i USHRT_MAX + */ + if ( maxDouble*fp > USHRT_MAX ){ + fp = fp - 1; + } + //cout << " max val: " << maxDouble << endl; + //cout << "fixed point: " << fp << endl; + + return fp; +} + + + +size_t encodeSlof( + const double *data, + size_t dataSize, + unsigned char *result, + double fixedPoint +) { + size_t i, ri; + double temp; + unsigned short x; + encodeFixedPoint(fixedPoint, result); + + ri = 8; + for (i=0; i USHRT_MAX ) { + // cout << std::setprecision(50) << "[MSNumpress::encodeSlof] Warning!! issue with data point: " << data[i] << ", log(data[i]+1) * fixedPoint ~= USHRT_MAX." << endl; + // cout << std::setprecision(50) << "[MSNumpress::encodeSlof] Warning!! fixedPoint: " << fixedPoint << endl; + // cout << std::setprecision(50) << "[MSNumpress::encodeSlof] Warning!! temp: " << temp << endl; + throw "[MSNumpress::encodeSlof] Cannot encode a number that overflows USHRT_MAX."; + } + + x = static_cast(temp + 0.5); + result[ri++] = x & 0xff; + result[ri++] = (x >> 8) & 0xff; + } + return ri; +} + + + +size_t decodeSlof( + const unsigned char *data, + const size_t dataSize, + double *result +) { + size_t i, ri; + unsigned short x; + double fixedPoint; + + if (dataSize < 8) + throw "[MSNumpress::decodeSlof] Corrupt input data: not enough bytes to read fixed point! "; + + ri = 0; + fixedPoint = decodeFixedPoint(data); + + for (i=8; i(data[i] | (data[i+1] << 8)); + result[ri++] = exp(x / fixedPoint) - 1; + } + return ri; +} + + + +void encodeSlof( + const std::vector &data, + std::vector &result, + double fixedPoint +) { + size_t dataSize = data.size(); + result.resize(dataSize * 2 + 8); + size_t encodedLength = encodeSlof(&data[0], dataSize, &result[0], fixedPoint); + result.resize(encodedLength); +} + + + +void decodeSlof( + const std::vector &data, + std::vector &result +) { + size_t dataSize = data.size(); + result.resize((dataSize - 8) / 2); + size_t decodedLength = decodeSlof(&data[0], dataSize, &result[0]); + result.resize(decodedLength); +} + +} +} // namespace numpress +} // namespace ms diff --git a/src/main/R/RMSNumpress/src/RMSNumpress.cpp b/src/main/R/RMSNumpress/src/RMSNumpress.cpp new file mode 100644 index 0000000..57520ac --- /dev/null +++ b/src/main/R/RMSNumpress/src/RMSNumpress.cpp @@ -0,0 +1,239 @@ +/* + RMSNumpress.cpp + justincsing@gmail.com + + Copyright 2020 Justin Sing + + Licensed under the Apache License, Version 2.0 (the "License"); +you may not use this file except in compliance with the License. +You may obtain a copy of the License at + + http://www.apache.org/licenses/LICENSE-2.0 + + Unless required by applicable law or agreed to in writing, software + distributed under the License is distributed on an "AS IS" BASIS, + WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +See the License for the specific language governing permissions and + limitations under the License. +*/ +#include +using namespace Rcpp; +#include "include/MSNumpress.hpp" + +//' optimalLinearFixedPointMass +//' +//' Compute the optimal linear fixed point with a desired m/z accuracy. +//' +//' @note If the desired accuracy cannot be reached without overflowing 64 +//' bit integers, then a negative value is returned. You need to check for +//' this and in that case abandon numpress or use optimalLinearFixedPoint +//' which returns the largest safe value. +//' +//' @param data pointer to array of double to be encoded (need memorycont. repr.) +//' @param mass_acc desired m/z accuracy in Th +//' @return the linear fixed point that satisfies the accuracy requirement (or -1 in case of failure). +// [[Rcpp::export]] +double optimalLinearFixedPointMass( + const std::vector &data, + double mass_acc +) { + size_t dataSize = data.size(); + double result = ms::numpress::MSNumpress::optimalLinearFixedPointMass(&data[0], dataSize, mass_acc); + return result; +} + +//' optimalLinearFixedPoint +//' +//' Compute the maximal linear fixed point that prevents integer overflow. +//' +//' @param data pointer to array of double to be encoded (need memorycont. repr.) +//' @return the linear fixed point safe to use +//' +// [[Rcpp::export]] +double optimalLinearFixedPoint(const std::vector &data){ + size_t dataSize = data.size(); + double result = ms::numpress::MSNumpress::optimalLinearFixedPoint( &data[0], dataSize ); + return result; +} + +//' optimalSlofFixedPoint +//' +//' Compute the maximal natural logarithm fixed point that prevents integer overflow. +//' +//' @param data pointer to array of double to be encoded (need memorycont. repr.) +//' @return the slof fixed point safe to use +//' +// [[Rcpp::export]] +double optimalSlofFixedPoint(const std::vector &data){ + size_t dataSize = data.size(); + double result = ms::numpress::MSNumpress::optimalSlofFixedPoint( &data[0], dataSize ); + return result; +} + +//' encodeLinear +//' +//' Encodes the doubles in data by first using a \cr +//' - lossy conversion to a 4 byte 5 decimal fixed point representation \cr +//' - storing the residuals from a linear prediction after first two values \cr +//' - encoding by encodeInt (see above) \cr +//' +//' The resulting binary is maximally 8 + dataSize * 5 bytes, but much less if the +//' data is reasonably smooth on the first order. +//' +//' This encoding is suitable for typical m/z or retention time binary arrays. +//' On a test set, the encoding was empirically show to be accurate to at least 0.002 ppm. +//' +//' @param data pointer to array of double to be encoded (need memorycont. repr.) +//' @param fixedPoint the scaling factor used for getting the fixed point repr. +//' This is stored in the binary and automatically extracted +//' on decoding (see optimalLinearFixedPoint or optimalLinearFixedPointMass) +//' @return the number of encoded bytes +//' +//' @seealso [\code{\link{decodeLinear}}] +//' +//' @examples +//' \dontrun{ +//' ## Retention time array +//' rt_array <- c(4313.0, 4316.4, 4319.8, 4323.2, 4326.6, 4330.1) +//' ## encode retention time array +//' rt_encoded <- encodeLinear(rt_array, 500) +//' #> [1] 40 7f 40 00 00 00 00 00 d4 e7 20 00 78 ee 20 00 88 86 23 +//' } +// [[Rcpp::export]] +std::vector encodeLinear(const std::vector &data, + double fixedPoint) { + std::vector result; + ms::numpress::MSNumpress::encodeLinear(data, result, fixedPoint); + return result; + } + +//' decodeLinear +//' +//' Decodes data encoded by encodeLinear. +//' +//' result vector guaranteed to be shorter or equal to (|data| - 8) * 2 +//' +//' Note that this method may throw a const char* if it deems the input data to be corrupt, i.e. +//' that the last encoded int does not use the last byte in the data. In addition the last encoded +//' int need to use either the last halfbyte, or the second last followed by a 0x0 halfbyte. +//' +//' @param data pointer to array of bytes to be decoded (need memorycont. repr.) +//' @return the number of decoded doubles, or -1 if dataSize < 4 or 4 < dataSize < 8 +//' +//' @seealso [\code{\link{encodeLinear}}] +//' +//' @examples +//' \dontrun{ +//' ## Retention time data that is encoded with encodeLinear and is zlib compressed +//' ### NOTE: For the sake of this example, I have broken the raw vector into several parts +//' ### to avoid Rd line widths (>100 characters) issues with CRAN build checks. +//' rt_raw1 <- c("78", "9c", "73", "50", "61", "00", "83", "aa", "15", "0c", "0c", "73", "80") +//' rt_raw2 <- c("b8", "a3", "5d", "fe", "47", "07", "84", "28", "fc", "8f", "c4", "40", "e5") +//' rt_raw3 <- c("61", "51", "84", "a9", "85", "08", "e1", "06", "00", "06", "be", "41", "cf") +//' ## Add all character representation of raw data back together and convert back to hex raw vector +//' rt_blob <- as.raw(as.hexmode(c(rt_raw1, rt_raw2, rt_raw3 ))) +//' ## Decompress blob +//' rt_blob_uncompressed <- memDecompress(rt_blob, type = "gzip", asChar = FALSE) +//' ## Decode to rentention time double values +//' rt_array <- decodeLinear(rt_blob_uncompressed) +//' } +// [[Rcpp::export]] +std::vector decodeLinear(const std::vector &data) { + std::vector result; + ms::numpress::MSNumpress::decodeLinear(data, result); + return result; +} + +//' encodeSlof +//' +//' Encodes ion counts by taking the natural logarithm, and storing a +//' fixed point representation of this. This is calculated as +//' +//' unsigned short fp = log(d + 1) * fixedPoint + 0.5 +//' +//' the result vector is exactly |data| * 2 + 8 bytes long +//' +//' @param data pointer to array of double to be encoded (need memorycont. repr.) +//' @param fixedPoint fixed point to use for encoding (see optimalSlofFixedPoint) +//' @return the number of encoded bytes +//' +//' @seealso [\code{\link{decodeSlof}}] +// [[Rcpp::export]] +std::vector encodeSlof(const std::vector &data, + double fixedPoint) { + std::vector result; + ms::numpress::MSNumpress::encodeSlof(data, result, fixedPoint); + return result; +} + +//' decodeSlof +//' +//' Decodes data encoded by encodeSlof +//' +//' The return will include exactly (|data| - 8) / 2 doubles. +//' +//' Note that this method may throw a const char* if it deems the input data to be corrupt. +//' +//' @param data pointer to array of bytes to be decoded (need memorycont. repr.) +//' @return the number of decoded doubles +//' +//' @seealso [\code{\link{encodeSlof}}] +//' @examples +//' \dontrun{ +//' ## Intensity array to encode +//' ### NOTE: For the sake of this example, I have broken the intensity vector into several parts +//' ### to avoid Rd line widths (>100 characters) issues with CRAN build checks. +//' int_array1 <- c(0.71773432, 0.43443741, 1.71883610, 0.13220307, 0.90664242) +//' int_array2 <- c(0.00000000, 0.00000000, 0.64213755, 0.43443741, 0.47221479) +//' ## Comcatenate into one intensity array +//' int_array <- c(int_array1, int_array2) +//' ## Encode intensity array using encodeSlof +//' int_encode <- encodeSlof( int_array, 16 ) +//' } +// [[Rcpp::export]] +std::vector decodeSlof(const std::vector &data) { + std::vector result; + ms::numpress::MSNumpress::decodeSlof(data, result); + return result; +} + +//' encodePic +//' +//' Encodes ion counts by simply rounding to the nearest 4 byte integer, +//' and compressing each integer with encodeInt. +//' +//' The handleable range is therefore 0 -> 4294967294. +//' The resulting binary is maximally dataSize * 5 bytes, but much less if the +//' data is close to 0 on average. +//' +//' @param data pointer to array of double to be encoded (need memorycont. repr.) +//' @return the number of encoded bytes +//' +//' @seealso [\code{\link{decodePic}}] +// [[Rcpp::export]] +std::vector encodePic(const std::vector &data) { + std::vector result; + ms::numpress::MSNumpress::encodePic(data, result); + return result; +} + +//' decodePic +//' +//' Decodes data encoded by encodePic +//' +//' result vector guaranteed to be shorter of equal to |data| * 2 +//' +//' Note that this method may throw a const char* if it deems the input data to be corrupt, i.e. +//' that the last encoded int does not use the last byte in the data. In addition the last encoded +//' int need to use either the last halfbyte, or the second last followed by a 0x0 halfbyte. +//' +//' @param data pointer to array of bytes to be decoded (need memorycont. repr.) +//' @return the number of decoded doubles +//' +//' @seealso [\code{\link{encodePic}}] +// [[Rcpp::export]] +std::vector decodePic(const std::vector &data) { + std::vector result; + ms::numpress::MSNumpress::decodePic(data, result); + return result; +} diff --git a/src/main/R/RMSNumpress/src/RcppExports.cpp b/src/main/R/RMSNumpress/src/RcppExports.cpp new file mode 100644 index 0000000..e77e8d4 --- /dev/null +++ b/src/main/R/RMSNumpress/src/RcppExports.cpp @@ -0,0 +1,127 @@ +// Generated by using Rcpp::compileAttributes() -> do not edit by hand +// Generator token: 10BE3573-1514-4C36-9D1C-5A225CD40393 + +#include + +using namespace Rcpp; + +// optimalLinearFixedPointMass +double optimalLinearFixedPointMass(const std::vector& data, double mass_acc); +RcppExport SEXP _RMSNumpress_optimalLinearFixedPointMass(SEXP dataSEXP, SEXP mass_accSEXP) { +BEGIN_RCPP + Rcpp::RObject rcpp_result_gen; + Rcpp::RNGScope rcpp_rngScope_gen; + Rcpp::traits::input_parameter< const std::vector& >::type data(dataSEXP); + Rcpp::traits::input_parameter< double >::type mass_acc(mass_accSEXP); + rcpp_result_gen = Rcpp::wrap(optimalLinearFixedPointMass(data, mass_acc)); + return rcpp_result_gen; +END_RCPP +} +// optimalLinearFixedPoint +double optimalLinearFixedPoint(const std::vector& data); +RcppExport SEXP _RMSNumpress_optimalLinearFixedPoint(SEXP dataSEXP) { +BEGIN_RCPP + Rcpp::RObject rcpp_result_gen; + Rcpp::RNGScope rcpp_rngScope_gen; + Rcpp::traits::input_parameter< const std::vector& >::type data(dataSEXP); + rcpp_result_gen = Rcpp::wrap(optimalLinearFixedPoint(data)); + return rcpp_result_gen; +END_RCPP +} +// optimalSlofFixedPoint +double optimalSlofFixedPoint(const std::vector& data); +RcppExport SEXP _RMSNumpress_optimalSlofFixedPoint(SEXP dataSEXP) { +BEGIN_RCPP + Rcpp::RObject rcpp_result_gen; + Rcpp::RNGScope rcpp_rngScope_gen; + Rcpp::traits::input_parameter< const std::vector& >::type data(dataSEXP); + rcpp_result_gen = Rcpp::wrap(optimalSlofFixedPoint(data)); + return rcpp_result_gen; +END_RCPP +} +// encodeLinear +std::vector encodeLinear(const std::vector& data, double fixedPoint); +RcppExport SEXP _RMSNumpress_encodeLinear(SEXP dataSEXP, SEXP fixedPointSEXP) { +BEGIN_RCPP + Rcpp::RObject rcpp_result_gen; + Rcpp::RNGScope rcpp_rngScope_gen; + Rcpp::traits::input_parameter< const std::vector& >::type data(dataSEXP); + Rcpp::traits::input_parameter< double >::type fixedPoint(fixedPointSEXP); + rcpp_result_gen = Rcpp::wrap(encodeLinear(data, fixedPoint)); + return rcpp_result_gen; +END_RCPP +} +// decodeLinear +std::vector decodeLinear(const std::vector& data); +RcppExport SEXP _RMSNumpress_decodeLinear(SEXP dataSEXP) { +BEGIN_RCPP + Rcpp::RObject rcpp_result_gen; + Rcpp::RNGScope rcpp_rngScope_gen; + Rcpp::traits::input_parameter< const std::vector& >::type data(dataSEXP); + rcpp_result_gen = Rcpp::wrap(decodeLinear(data)); + return rcpp_result_gen; +END_RCPP +} +// encodeSlof +std::vector encodeSlof(const std::vector& data, double fixedPoint); +RcppExport SEXP _RMSNumpress_encodeSlof(SEXP dataSEXP, SEXP fixedPointSEXP) { +BEGIN_RCPP + Rcpp::RObject rcpp_result_gen; + Rcpp::RNGScope rcpp_rngScope_gen; + Rcpp::traits::input_parameter< const std::vector& >::type data(dataSEXP); + Rcpp::traits::input_parameter< double >::type fixedPoint(fixedPointSEXP); + rcpp_result_gen = Rcpp::wrap(encodeSlof(data, fixedPoint)); + return rcpp_result_gen; +END_RCPP +} +// decodeSlof +std::vector decodeSlof(const std::vector& data); +RcppExport SEXP _RMSNumpress_decodeSlof(SEXP dataSEXP) { +BEGIN_RCPP + Rcpp::RObject rcpp_result_gen; + Rcpp::RNGScope rcpp_rngScope_gen; + Rcpp::traits::input_parameter< const std::vector& >::type data(dataSEXP); + rcpp_result_gen = Rcpp::wrap(decodeSlof(data)); + return rcpp_result_gen; +END_RCPP +} +// encodePic +std::vector encodePic(const std::vector& data); +RcppExport SEXP _RMSNumpress_encodePic(SEXP dataSEXP) { +BEGIN_RCPP + Rcpp::RObject rcpp_result_gen; + Rcpp::RNGScope rcpp_rngScope_gen; + Rcpp::traits::input_parameter< const std::vector& >::type data(dataSEXP); + rcpp_result_gen = Rcpp::wrap(encodePic(data)); + return rcpp_result_gen; +END_RCPP +} +// decodePic +std::vector decodePic(const std::vector& data); +RcppExport SEXP _RMSNumpress_decodePic(SEXP dataSEXP) { +BEGIN_RCPP + Rcpp::RObject rcpp_result_gen; + Rcpp::RNGScope rcpp_rngScope_gen; + Rcpp::traits::input_parameter< const std::vector& >::type data(dataSEXP); + rcpp_result_gen = Rcpp::wrap(decodePic(data)); + return rcpp_result_gen; +END_RCPP +} + +static const R_CallMethodDef CallEntries[] = { + {"_RMSNumpress_optimalLinearFixedPointMass", (DL_FUNC) &_RMSNumpress_optimalLinearFixedPointMass, 2}, + {"_RMSNumpress_optimalLinearFixedPoint", (DL_FUNC) &_RMSNumpress_optimalLinearFixedPoint, 1}, + {"_RMSNumpress_optimalSlofFixedPoint", (DL_FUNC) &_RMSNumpress_optimalSlofFixedPoint, 1}, + {"_RMSNumpress_encodeLinear", (DL_FUNC) &_RMSNumpress_encodeLinear, 2}, + {"_RMSNumpress_decodeLinear", (DL_FUNC) &_RMSNumpress_decodeLinear, 1}, + {"_RMSNumpress_encodeSlof", (DL_FUNC) &_RMSNumpress_encodeSlof, 2}, + {"_RMSNumpress_decodeSlof", (DL_FUNC) &_RMSNumpress_decodeSlof, 1}, + {"_RMSNumpress_encodePic", (DL_FUNC) &_RMSNumpress_encodePic, 1}, + {"_RMSNumpress_decodePic", (DL_FUNC) &_RMSNumpress_decodePic, 1}, + {NULL, NULL, 0} +}; + +RcppExport void R_init_RMSNumpress(DllInfo *dll) { + R_registerRoutines(dll, NULL, CallEntries, NULL, NULL); + R_useDynamicSymbols(dll, FALSE); +} diff --git a/src/main/R/RMSNumpress/src/include/MSNumpress.hpp b/src/main/R/RMSNumpress/src/include/MSNumpress.hpp new file mode 100644 index 0000000..525104d --- /dev/null +++ b/src/main/R/RMSNumpress/src/include/MSNumpress.hpp @@ -0,0 +1,331 @@ +/* + MSNumpress.hpp + johan.teleman@immun.lth.se + + Copyright 2013 Johan Teleman + + Licensed under the Apache License, Version 2.0 (the "License"); + you may not use this file except in compliance with the License. + You may obtain a copy of the License at + + http://www.apache.org/licenses/LICENSE-2.0 + + Unless required by applicable law or agreed to in writing, software + distributed under the License is distributed on an "AS IS" BASIS, + WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + See the License for the specific language governing permissions and + limitations under the License. + */ +/* + ==================== encodeInt ==================== + Some of the encodings described below use a integer compression referred to simply as + + encodeInt() + + This encoding works on a 4 byte integer, by truncating initial zeros or ones. + If the initial (most significant) half byte is 0x0 or 0xf, the number of such + halfbytes starting from the most significant is stored in a halfbyte. This initial + count is then followed by the rest of the ints halfbytes, in little-endian order. + A count halfbyte c of + + 0 <= c <= 8 is interpreted as an initial c 0x0 halfbytes + 9 <= c <= 15 is interpreted as an initial (c-8) 0xf halfbytes + + Ex: + int c rest + 0 => 0x8 + -1 => 0xf 0xf + 23 => 0x6 0x7 0x1 + */ + +#ifndef _MSNUMPRESS_HPP_ +#define _MSNUMPRESS_HPP_ + +#include +#include +#include + +// defines whether to throw an exception when a number cannot be encoded safely +// with the given parameters +#ifndef THROW_ON_OVERFLOW +#define THROW_ON_OVERFLOW true +#endif + +namespace ms { +namespace numpress { + +namespace MSNumpress { + + /** + * Compute the maximal linear fixed point that prevents integer overflow. + * + * @data pointer to array of double to be encoded (need memorycont. repr.) + * @dataSize number of doubles from *data to encode + * + * @return the linear fixed point safe to use + */ + double optimalLinearFixedPoint( + const double *data, + size_t dataSize); + + /** + * Compute the optimal linear fixed point with a desired m/z accuracy. + * + * @note If the desired accuracy cannot be reached without overflowing 64 + * bit integers, then a negative value is returned. You need to check for + * this and in that case abandon numpress or use optimalLinearFixedPoint + * which returns the largest safe value. + * + * @data pointer to array of double to be encoded (need memorycont. repr.) + * @dataSize number of doubles from *data to encode + * @mass_acc desired m/z accuracy in Th + * + * @return the linear fixed point that satisfies the accuracy requirement (or -1 in case of failure). + */ + double optimalLinearFixedPointMass( + const double *data, + size_t dataSize, + double mass_acc); + + /** + * Encodes the doubles in data by first using a + * - lossy conversion to a 4 byte 5 decimal fixed point representation + * - storing the residuals from a linear prediction after first two values + * - encoding by encodeInt (see above) + * + * The resulting binary is maximally 8 + dataSize * 5 bytes, but much less if the + * data is reasonably smooth on the first order. + * + * This encoding is suitable for typical m/z or retention time binary arrays. + * On a test set, the encoding was empirically show to be accurate to at least 0.002 ppm. + * + * @data pointer to array of double to be encoded (need memorycont. repr.) + * @dataSize number of doubles from *data to encode + * @result pointer to where resulting bytes should be stored + * @fixedPoint the scaling factor used for getting the fixed point repr. + * This is stored in the binary and automatically extracted + * on decoding. + * @return the number of encoded bytes + */ + size_t encodeLinear( + const double *data, + const size_t dataSize, + unsigned char *result, + double fixedPoint); + + /** + * Calls lower level encodeLinear while handling vector sizes appropriately + * + * @data vector of doubles to be encoded + * @result vector of resulting bytes (will be resized to the number of bytes) + */ + void encodeLinear( + const std::vector &data, + std::vector &result, + double fixedPoint); + + /** + * Decodes data encoded by encodeLinear. + * + * result vector guaranteed to be shorter or equal to (|data| - 8) * 2 + * + * Note that this method may throw a const char* if it deems the input data to be corrupt, i.e. + * that the last encoded int does not use the last byte in the data. In addition the last encoded + * int need to use either the last halfbyte, or the second last followed by a 0x0 halfbyte. + * + * @data pointer to array of bytes to be decoded (need memorycont. repr.) + * @dataSize number of bytes from *data to decode + * @result pointer to were resulting doubles should be stored + * @return the number of decoded doubles, or -1 if dataSize < 4 or 4 < dataSize < 8 + */ + size_t decodeLinear( + const unsigned char *data, + const size_t dataSize, + double *result); + + /** + * Calls lower level decodeLinear while handling vector sizes appropriately + * + * Note that this method may throw a const char* if it deems the input data to be corrupt, i.e.. + * that the last encoded int does not use the last byte in the data. In addition the last encoded + * int need to use either the last halfbyte, or the second last followed by a 0x0 halfbyte. + * + * @data vector of bytes to be decoded + * @result vector of resulting double (will be resized to the number of doubles) + */ + void decodeLinear( + const std::vector &data, + std::vector &result); + +///////////////////////////////////////////////////////////// + + + /** + * Encodes the doubles in data by storing the residuals from a linear prediction after first two values. + * + * The resulting binary is the same size as the input data. + * + * This encoding is suitable for typical m/z or retention time binary arrays, and is + * intended to be used before zlib compression to improve compression. + * + * @data pointer to array of doubles to be encoded (need memorycont. repr.) + * @dataSize number of doubles from *data to encode + * @result pointer to were resulting bytes should be stored + */ + size_t encodeSafe( + const double *data, + const size_t dataSize, + unsigned char *result); + + + /** + * Decodes data encoded by encodeSafe. + * + * result vector is the same size as the input data. + * + * Might throw const char* is something goes wrong during decoding. + * + * @data pointer to array of bytes to be decoded (need memorycont. repr.) + * @dataSize number of bytes from *data to decode + * @result pointer to were resulting doubles should be stored + * @return the number of decoded bytes + */ + size_t decodeSafe( + const unsigned char *data, + const size_t dataSize, + double *result); + +///////////////////////////////////////////////////////////// + + /** + * Encodes ion counts by simply rounding to the nearest 4 byte integer, + * and compressing each integer with encodeInt. + * + * The handleable range is therefore 0 -> 4294967294. + * The resulting binary is maximally dataSize * 5 bytes, but much less if the + * data is close to 0 on average. + * + * @data pointer to array of double to be encoded (need memorycont. repr.) + * @dataSize number of doubles from *data to encode + * @result pointer to were resulting bytes should be stored + * @return the number of encoded bytes + */ + size_t encodePic( + const double *data, + const size_t dataSize, + unsigned char *result); + + /** + * Calls lower level encodePic while handling vector sizes appropriately + * + * @data vector of doubles to be encoded + * @result vector of resulting bytes (will be resized to the number of bytes) + */ + void encodePic( + const std::vector &data, + std::vector &result); + + /** + * Decodes data encoded by encodePic + * + * result vector guaranteed to be shorter of equal to |data| * 2 + * + * Note that this method may throw a const char* if it deems the input data to be corrupt, i.e. + * that the last encoded int does not use the last byte in the data. In addition the last encoded + * int need to use either the last halfbyte, or the second last followed by a 0x0 halfbyte. + * + * @data pointer to array of bytes to be decoded (need memorycont. repr.) + * @dataSize number of bytes from *data to decode + * @result pointer to were resulting doubles should be stored + * @return the number of decoded doubles + */ + size_t decodePic( + const unsigned char *data, + const size_t dataSize, + double *result); + + /** + * Calls lower level decodePic while handling vector sizes appropriately + * + * Note that this method may throw a const char* if it deems the input data to be corrupt, i.e. + * that the last encoded int does not use the last byte in the data. In addition the last encoded + * int need to use either the last halfbyte, or the second last followed by a 0x0 halfbyte. + * + * @data vector of bytes to be decoded + * @result vector of resulting double (will be resized to the number of doubles) + */ + void decodePic( + const std::vector &data, + std::vector &result); + +///////////////////////////////////////////////////////////// + + + double optimalSlofFixedPoint( + const double *data, + size_t dataSize); + + /** + * Encodes ion counts by taking the natural logarithm, and storing a + * fixed point representation of this. This is calculated as + * + * unsigned short fp = log(d + 1) * fixedPoint + 0.5 + * + * the result vector is exactly |data| * 2 + 8 bytes long + * + * @data pointer to array of double to be encoded (need memorycont. repr.) + * @dataSize number of doubles from *data to encode + * @result pointer to were resulting bytes should be stored + * @return the number of encoded bytes + */ + size_t encodeSlof( + const double *data, + const size_t dataSize, + unsigned char *result, + double fixedPoint); + + /** + * Calls lower level encodeSlof while handling vector sizes appropriately + * + * @data vector of doubles to be encoded + * @result vector of resulting bytes (will be resized to the number of bytes) + */ + void encodeSlof( + const std::vector &data, + std::vector &result, + double fixedPoint); + + /** + * Decodes data encoded by encodeSlof + * + * The return will include exactly (|data| - 8) / 2 doubles. + * + * Note that this method may throw a const char* if it deems the input data to be corrupt. + * + * @data pointer to array of bytes to be decoded (need memorycont. repr.) + * @dataSize number of bytes from *data to decode + * @result pointer to were resulting doubles should be stored + * @return the number of decoded doubles + */ + size_t decodeSlof( + const unsigned char *data, + const size_t dataSize, + double *result); + + /** + * Calls lower level decodeSlof while handling vector sizes appropriately + * + * Note that this method may throw a const char* if it deems the input data to be corrupt. + * + * @data vector of bytes to be decoded + * @result vector of resulting double (will be resized to the number of doubles) + */ + void decodeSlof( + const std::vector &data, + std::vector &result); + +} // namespace MSNumpress +} // namespace msdata +} // namespace pwiz + +#endif // _MSNUMPRESS_HPP_ diff --git a/src/main/R/RMSNumpress/tests/testthat.R b/src/main/R/RMSNumpress/tests/testthat.R new file mode 100644 index 0000000..895d5b2 --- /dev/null +++ b/src/main/R/RMSNumpress/tests/testthat.R @@ -0,0 +1,4 @@ +library(testthat) +library(RMSNumpress) + +test_check("RMSNumpress") diff --git a/src/main/R/RMSNumpress/tests/testthat/test_RMSNumpress.R b/src/main/R/RMSNumpress/tests/testthat/test_RMSNumpress.R new file mode 100644 index 0000000..9658a4a --- /dev/null +++ b/src/main/R/RMSNumpress/tests/testthat/test_RMSNumpress.R @@ -0,0 +1,91 @@ +#!/bin/usr/Rscript +library(RMSNumpress) + +data = c(100, 101, 102, 103) +data_long = c(100.0, + 200.0, + 300.00005, + 400.00010, + 450.00010, + 455.00010, + 700.00010) +data_slof = c(100.0, + 200.0, + 300.00005, + 400.00010) +fp_slof = 10000 +linear_result = (as.raw(as.hexmode(c("40", "f8", "6a", "00", "00", "00", "00", "00", "80", "96", "98", "00", "20", "1d", "9a", "00", "88")))) + +test_that("encodeLinear encodes doubles in data", { + expect_equal( encodeLinear(data, 100000.0), linear_result ) + expect_equal( length(encodeLinear(data, 100000.0)), 17) + expect_equal( encodeLinear(data, 100000.0)[1], as.raw(as.hexmode("40"))) + + expect_equal( length(encodeLinear(data_long, 5.0)), 22) + + expect_equal( length(encodeLinear(data_long, 500.0)), 25) + + expect_equal( length(encodeLinear(data_long, 5e4)), 29) + + expect_equal( length(encodeLinear(data_long, 5e5)), 30) + + expect_equal( length(encodeLinear(data_long, 5e6)), 31) + + # accurate to 3 sign digits + result = encodeLinear(data_long, 500.0) + decoded = decodeLinear(result) + + expect_equal(decoded[1], 100, tolerance=3) + expect_equal(decoded[2], 200, tolerance=3) + expect_equal(decoded[3], 300, tolerance=3) + expect_equal(decoded[4], 400.00010, tolerance=3) + expect_equal(decoded[5], 450.00010, tolerance=3) + expect_equal(decoded[6], 455.00010, tolerance=3) + expect_equal(decoded[7], 700.00010, tolerance=3) +}) + +test_that("decodeLinear decode data encoded with encodeLinear", { + expect_equal( decodeLinear(as.raw(linear_result)), data ) + expect_equal( length(decodeLinear(as.raw(linear_result))), 4 ) + expect_equal( decodeLinear(as.raw(linear_result))[1], 100 ) +}) + +test_that("encodePic encodes ion counts", { + encode = encodePic(data) + expect_equal( length(encode), 6 ) + + result = decodePic(encode) + expect_equal( length(result), 4 ) + expect_equal( result[1], 100 ) + expect_equal( result, data ) +}) + +test_that("encodeSlof encodes ion counts", { + encode = encodeSlof(data_slof, fp_slof) + expect_equal( length(encode), 16 ) + + result = decodeSlof(encode) + expect_equal( length(result), 4 ) + + expect_true( abs(result[1] - 100) < 1) + expect_true( abs(result[2] - 200) < 1) + expect_true( abs(result[3] - 300) < 1) + expect_true( abs(result[4] - 400) < 1) +}) + +test_that("Compute the maximal linear fixed point that prevents integer overflow", { + pt = optimalLinearFixedPoint(data) + expect_equal(pt, 21262214.0) +}) + +test_that("Compute the maximal slof fixed point that prevents integer overflow", { + pt = optimalSlofFixedPoint(data) + expect_equal(pt, 14110.0) +}) + +test_that("optimal linear fixed point with a desired m/z accuracy", { + pt = optimalLinearFixedPointMass(data, 0.001) + expect_equal(pt, 500.0) + pt = optimalLinearFixedPointMass(data, 1e-10) + expect_equal(pt, -1) +}) \ No newline at end of file diff --git a/src/main/cpp/MSNumpress.cpp b/src/main/cpp/MSNumpress.cpp new file mode 100644 index 0000000..480c743 --- /dev/null +++ b/src/main/cpp/MSNumpress.cpp @@ -0,0 +1,782 @@ +/* + MSNumpress.cpp + johan.teleman@immun.lth.se + + Copyright 2013 Johan Teleman + + Licensed under the Apache License, Version 2.0 (the "License"); + you may not use this file except in compliance with the License. + You may obtain a copy of the License at + + http://www.apache.org/licenses/LICENSE-2.0 + + Unless required by applicable law or agreed to in writing, software + distributed under the License is distributed on an "AS IS" BASIS, + WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + See the License for the specific language governing permissions and + limitations under the License. + */ + +#include +#include +#include +#include +#include "MSNumpress.hpp" + +namespace ms { +namespace numpress { +namespace MSNumpress { + +using std::cout; +using std::cerr; +using std::endl; +using std::min; +using std::max; +using std::abs; + +// This is only valid on systems were ints use more bytes than chars... + +const int ONE = 1; +static bool is_little_endian() { + return *((char*)&(ONE)) == 1; +} +bool IS_LITTLE_ENDIAN = is_little_endian(); + + + +///////////////////////////////////////////////////////////// + +static void encodeFixedPoint( + double fixedPoint, + unsigned char *result +) { + int i; + unsigned char *fp = (unsigned char*)&fixedPoint; + for (i=0; i<8; i++) { + result[i] = fp[IS_LITTLE_ENDIAN ? (7-i) : i]; + } +} + + + +static double decodeFixedPoint( + const unsigned char *data +) { + int i; + double fixedPoint; + unsigned char *fp = (unsigned char*)&fixedPoint; + + for (i=0; i<8; i++) { + fp[i] = data[IS_LITTLE_ENDIAN ? (7-i) : i]; + } + + return fixedPoint; +} + +///////////////////////////////////////////////////////////// + +/** + * Encodes the int x as a number of halfbytes in res. + * res_length is incremented by the number of halfbytes, + * which will be 1 <= n <= 9 + */ +static void encodeInt( + const unsigned int x, + unsigned char* res, + size_t *res_length +) { + // get the bit pattern of a signed int x_inp + unsigned int m; + unsigned char i, l; // numbers between 0 and 9 + + unsigned int mask = 0xf0000000; + unsigned int init = x & mask; + + if (init == 0) { + l = 8; + for (i=0; i<8; i++) { + m = mask >> (4*i); + if ((x & m) != 0) { + l = i; + break; + } + } + res[0] = l; + for (i=l; i<8; i++) { + res[1+i-l] = static_cast( x >> (4*(i-l)) ); + } + *res_length += 1+8-l; + + } else if (init == mask) { + l = 7; + for (i=0; i<8; i++) { + m = mask >> (4*i); + if ((x & m) != m) { + l = i; + break; + } + } + res[0] = l + 8; + for (i=l; i<8; i++) { + res[1+i-l] = static_cast( x >> (4*(i-l)) ); + } + *res_length += 1+8-l; + + } else { + res[0] = 0; + for (i=0; i<8; i++) { + res[1+i] = static_cast( x >> (4*i) ); + } + *res_length += 9; + + } +} + + + +/** + * Decodes an int from the half bytes in bp. Lossless reverse of encodeInt + * + * @param data ptr to the char data to decode + * @param di position in the char data array to start decoding (will be advanced) + * @param max_di size of data array + * @param half helper variable (do not change between multiple calls) + * @param res result (a 32 bit integer) + * + * @note the helper variable indicates whether we look at the first half byte + * or second half byte of the current data (thus whether to interpret the first + * half byte of data[*di] or the second half byte). + * + */ +static void decodeInt( + const unsigned char *data, + size_t *di, + size_t max_di, + size_t *half, + unsigned int *res +) { + size_t n, i; + unsigned int mask, m; + unsigned char head; + unsigned char hb; + + // Extract the first half byte, specifying the number of leading zero half + // bytes of the final integer. + // If half is zero, we look at the first half byte, otherwise we look at + // the second (lower) half byte and advance the counter to the next char. + if (*half == 0) { + head = data[*di] >> 4; + } else { + head = data[*di] & 0xf; + (*di)++; + } + + *half = 1-(*half); // switch to other half byte + *res = 0; + + if (head <= 8) { + n = head; + } else { // we have n leading ones, fill n half bytes in res with 0xf + n = head - 8; + mask = 0xf0000000; + for (i=0; i> (4*i); + *res = *res | m; + } + } + + if (n == 8) { + return; + } + + if (*di + ((8 - n) - (1 - *half)) / 2 >= max_di) { + throw "[MSNumpress::decodeInt] Corrupt input data! "; + } + + for (i=n; i<8; i++) { + if (*half == 0) { + hb = data[*di] >> 4; + } else { + hb = data[*di] & 0xf; + (*di)++; + } + *res = *res | ( static_cast(hb) << ((i-n)*4)); + *half = 1 - (*half); + } +} + + + + +///////////////////////////////////////////////////////////// + +double optimalLinearFixedPointMass( + const double *data, + size_t dataSize, + double mass_acc +) { + if (dataSize < 3) return 0; // we just encode the first two points as floats + + // We calculate the maximal fixedPoint we need to achieve a specific mass + // accuracy. Note that the maximal error we will make by encoding as int is + // 0.5 due to rounding errors. + double maxFP = 0.5 / mass_acc; + + // There is a maximal value for the FP given by the int length (32bit) + // which means we cannot choose a value higher than that. In case we cannot + // achieve the desired accuracy, return failure (-1). + double maxFP_overflow = optimalLinearFixedPoint(data, dataSize); + if (maxFP > maxFP_overflow) return -1; + + return maxFP; +} + +double optimalLinearFixedPoint( + const double *data, + size_t dataSize +) { + /* + * safer impl - apparently not needed though + * + if (dataSize == 0) return 0; + + double maxDouble = 0; + double x; + + for (size_t i=0; i(data[0] * fixedPoint + 0.5); + for (i=0; i<4; i++) { + result[8+i] = (ints[1] >> (i*8)) & 0xff; + } + + if (dataSize == 1) return 12; + + ints[2] = static_cast(data[1] * fixedPoint + 0.5); + for (i=0; i<4; i++) { + result[12+i] = (ints[2] >> (i*8)) & 0xff; + } + + halfByteCount = 0; + ri = 16; + + for (i=2; i LLONG_MAX ) { + throw "[MSNumpress::encodeLinear] Next number overflows LLONG_MAX."; + } + + ints[2] = static_cast(data[i] * fixedPoint + 0.5); + extrapol = ints[1] + (ints[1] - ints[0]); + + if (THROW_ON_OVERFLOW && + ( ints[2] - extrapol > INT_MAX + || ints[2] - extrapol < INT_MIN )) { + throw "[MSNumpress::encodeLinear] Cannot encode a number that exceeds the bounds of [-INT_MAX, INT_MAX]."; + } + + diff = static_cast(ints[2] - extrapol); + //printf("%lu %lu %lu, extrapol: %ld diff: %d \n", ints[0], ints[1], ints[2], extrapol, diff); + encodeInt( + static_cast(diff), + &halfBytes[halfByteCount], + &halfByteCount + ); + /* + printf("%d (%d): ", diff, (int)halfByteCount); + for (size_t j=0; j( + (halfBytes[hbi-1] << 4) | (halfBytes[hbi] & 0xf) + ); + //printf("%x \n", result[ri]); + ri++; + } + if (halfByteCount % 2 != 0) { + halfBytes[0] = halfBytes[halfByteCount-1]; + halfByteCount = 1; + } else { + halfByteCount = 0; + } + } + if (halfByteCount == 1) { + result[ri] = static_cast(halfBytes[0] << 4); + ri++; + } + return ri; +} + + + +size_t decodeLinear( + const unsigned char *data, + const size_t dataSize, + double *result +) { + size_t i; + size_t ri = 0; + unsigned int init, buff; + int diff; + long long ints[3]; + //double d; + size_t di; + size_t half; + long long extrapol; + long long y; + double fixedPoint; + + //printf("Decoding %d bytes with fixed point %f\n", (int)dataSize, fixedPoint); + + if (dataSize == 8) return 0; + + if (dataSize < 8) + throw "[MSNumpress::decodeLinear] Corrupt input data: not enough bytes to read fixed point! "; + + fixedPoint = decodeFixedPoint(data); + + + if (dataSize < 12) + throw "[MSNumpress::decodeLinear] Corrupt input data: not enough bytes to read first value! "; + + ints[1] = 0; + for (i=0; i<4; i++) { + ints[1] = ints[1] | ((0xff & (init = data[8+i])) << (i*8)); + } + result[0] = ints[1] / fixedPoint; + + if (dataSize == 12) return 1; + if (dataSize < 16) + throw "[MSNumpress::decodeLinear] Corrupt input data: not enough bytes to read second value! "; + + ints[2] = 0; + for (i=0; i<4; i++) { + ints[2] = ints[2] | ((0xff & (init = data[12+i])) << (i*8)); + } + result[1] = ints[2] / fixedPoint; + + half = 0; + ri = 2; + di = 16; + + //printf(" di ri half int[0] int[1] extrapol diff\n"); + + while (di < dataSize) { + if (di == (dataSize - 1) && half == 1) { + if ((data[di] & 0xf) == 0x0) { + break; + } + } + //printf("%7d %7d %7d %lu %lu %ld", di, ri, half, ints[0], ints[1], extrapol); + + ints[0] = ints[1]; + ints[1] = ints[2]; + decodeInt(data, &di, dataSize, &half, &buff); + diff = static_cast(buff); + + extrapol = ints[1] + (ints[1] - ints[0]); + y = extrapol + diff; + //printf(" %d \n", diff); + result[ri++] = y / fixedPoint; + ints[2] = y; + } + + return ri; +} + + + +void encodeLinear( + const std::vector &data, + std::vector &result, + double fixedPoint +) { + size_t dataSize = data.size(); + result.resize(dataSize * 5 + 8); + size_t encodedLength = encodeLinear(&data[0], dataSize, &result[0], fixedPoint); + result.resize(encodedLength); +} + + + +void decodeLinear( + const std::vector &data, + std::vector &result +) { + size_t dataSize = data.size(); + result.resize((dataSize - 8) * 2); + size_t decodedLength = decodeLinear(&data[0], dataSize, &result[0]); + result.resize(decodedLength); +} + +///////////////////////////////////////////////////////////// + + +size_t encodeSafe( + const double *data, + const size_t dataSize, + unsigned char *result +) { + size_t i, j, ri = 0; + double latest[3]; + double extrapol, diff; + const unsigned char *fp; + + //printf("d0 d1 d2 extrapol diff\n"); + + if (dataSize == 0) return ri; + + latest[1] = data[0]; + fp = (unsigned char*)data; + for (i=0; i<8; i++) { + result[ri++] = fp[IS_LITTLE_ENDIAN ? (7-i) : i]; + } + + if (dataSize == 1) return ri; + + latest[2] = data[1]; + fp = (unsigned char*)&(data[1]); + for (i=0; i<8; i++) { + result[ri++] = fp[IS_LITTLE_ENDIAN ? (7-i) : i]; + } + + fp = (unsigned char*)&diff; + for (i=2; i INT_MAX || data[i] < -0.5) ){ + throw "[MSNumpress::encodePic] Cannot use Pic to encode a number larger than INT_MAX or smaller than 0."; + } + x = static_cast(data[i] + 0.5); + //printf("%d %d %d, extrapol: %d diff: %d \n", ints[0], ints[1], ints[2], extrapol, diff); + encodeInt(x, &halfBytes[halfByteCount], &halfByteCount); + + for (hbi=1; hbi < halfByteCount; hbi+=2) { + result[ri] = static_cast( + (halfBytes[hbi-1] << 4) | (halfBytes[hbi] & 0xf) + ); + //printf("%x \n", result[ri]); + ri++; + } + if (halfByteCount % 2 != 0) { + halfBytes[0] = halfBytes[halfByteCount-1]; + halfByteCount = 1; + } else { + halfByteCount = 0; + } + } + if (halfByteCount == 1) { + result[ri] = static_cast(halfBytes[0] << 4); + ri++; + } + return ri; +} + + + +size_t decodePic( + const unsigned char *data, + const size_t dataSize, + double *result +) { + size_t ri; + unsigned int x; + size_t di; + size_t half; + + //printf("ri di half dSize count\n"); + + half = 0; + ri = 0; + di = 0; + + while (di < dataSize) { + if (di == (dataSize - 1) && half == 1) { + if ((data[di] & 0xf) == 0x0) { + break; + } + } + + decodeInt(&data[0], &di, dataSize, &half, &x); + + //printf("%7d %7d %7d %7d %7d\n", ri, di, half, dataSize, count); + + //printf("count: %d \n", count); + result[ri++] = static_cast(x); + } + + return ri; +} + + + +void encodePic( + const std::vector &data, + std::vector &result +) { + size_t dataSize = data.size(); + result.resize(dataSize * 5); + size_t encodedLength = encodePic(&data[0], dataSize, &result[0]); + result.resize(encodedLength); +} + + + +void decodePic( + const std::vector &data, + std::vector &result +) { + size_t dataSize = data.size(); + result.resize(dataSize * 2); + size_t decodedLength = decodePic(&data[0], dataSize, &result[0]); + result.resize(decodedLength); +} + + +///////////////////////////////////////////////////////////// + + +double optimalSlofFixedPoint( + const double *data, + size_t dataSize +) { + if (dataSize == 0) return 0; + + double maxDouble = 1; + double x; + double fp; + + for (size_t i=0; i USHRT_MAX ) { + throw "[MSNumpress::encodeSlof] Cannot encode a number that overflows USHRT_MAX."; + } + + x = static_cast(temp + 0.5); + result[ri++] = x & 0xff; + result[ri++] = (x >> 8) & 0xff; + } + return ri; +} + + + +size_t decodeSlof( + const unsigned char *data, + const size_t dataSize, + double *result +) { + size_t i, ri; + unsigned short x; + double fixedPoint; + + if (dataSize < 8) + throw "[MSNumpress::decodeSlof] Corrupt input data: not enough bytes to read fixed point! "; + + ri = 0; + fixedPoint = decodeFixedPoint(data); + + for (i=8; i(data[i] | (data[i+1] << 8)); + result[ri++] = exp(x / fixedPoint) - 1; + } + return ri; +} + + + +void encodeSlof( + const std::vector &data, + std::vector &result, + double fixedPoint +) { + size_t dataSize = data.size(); + result.resize(dataSize * 2 + 8); + size_t encodedLength = encodeSlof(&data[0], dataSize, &result[0], fixedPoint); + result.resize(encodedLength); +} + + + +void decodeSlof( + const std::vector &data, + std::vector &result +) { + size_t dataSize = data.size(); + result.resize((dataSize - 8) / 2); + size_t decodedLength = decodeSlof(&data[0], dataSize, &result[0]); + result.resize(decodedLength); +} + +} +} // namespace numpress +} // namespace ms diff --git a/src/main/cpp/MSNumpress.hpp b/src/main/cpp/MSNumpress.hpp new file mode 100644 index 0000000..985820f --- /dev/null +++ b/src/main/cpp/MSNumpress.hpp @@ -0,0 +1,330 @@ +/* + MSNumpress.hpp + johan.teleman@immun.lth.se + + Copyright 2013 Johan Teleman + + Licensed under the Apache License, Version 2.0 (the "License"); + you may not use this file except in compliance with the License. + You may obtain a copy of the License at + + http://www.apache.org/licenses/LICENSE-2.0 + + Unless required by applicable law or agreed to in writing, software + distributed under the License is distributed on an "AS IS" BASIS, + WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + See the License for the specific language governing permissions and + limitations under the License. + */ +/* + ==================== encodeInt ==================== + Some of the encodings described below use a integer compression referred to simply as + + encodeInt() + + This encoding works on a 4 byte integer, by truncating initial zeros or ones. + If the initial (most significant) half byte is 0x0 or 0xf, the number of such + halfbytes starting from the most significant is stored in a halfbyte. This initial + count is then followed by the rest of the ints halfbytes, in little-endian order. + A count halfbyte c of + + 0 <= c <= 8 is interpreted as an initial c 0x0 halfbytes + 9 <= c <= 15 is interpreted as an initial (c-8) 0xf halfbytes + + Ex: + int c rest + 0 => 0x8 + -1 => 0xf 0xf + 23 => 0x6 0x7 0x1 + */ + +#ifndef _MSNUMPRESS_HPP_ +#define _MSNUMPRESS_HPP_ + +#include +#include + +// defines whether to throw an exception when a number cannot be encoded safely +// with the given parameters +#ifndef THROW_ON_OVERFLOW +#define THROW_ON_OVERFLOW true +#endif + +namespace ms { +namespace numpress { + +namespace MSNumpress { + + /** + * Compute the maximal linear fixed point that prevents integer overflow. + * + * @data pointer to array of double to be encoded (need memorycont. repr.) + * @dataSize number of doubles from *data to encode + * + * @return the linear fixed point safe to use + */ + double optimalLinearFixedPoint( + const double *data, + size_t dataSize); + + /** + * Compute the optimal linear fixed point with a desired m/z accuracy. + * + * @note If the desired accuracy cannot be reached without overflowing 64 + * bit integers, then a negative value is returned. You need to check for + * this and in that case abandon numpress or use optimalLinearFixedPoint + * which returns the largest safe value. + * + * @data pointer to array of double to be encoded (need memorycont. repr.) + * @dataSize number of doubles from *data to encode + * @mass_acc desired m/z accuracy in Th + * + * @return the linear fixed point that satisfies the accuracy requirement (or -1 in case of failure). + */ + double optimalLinearFixedPointMass( + const double *data, + size_t dataSize, + double mass_acc); + + /** + * Encodes the doubles in data by first using a + * - lossy conversion to a 4 byte 5 decimal fixed point representation + * - storing the residuals from a linear prediction after first two values + * - encoding by encodeInt (see above) + * + * The resulting binary is maximally 8 + dataSize * 5 bytes, but much less if the + * data is reasonably smooth on the first order. + * + * This encoding is suitable for typical m/z or retention time binary arrays. + * On a test set, the encoding was empirically show to be accurate to at least 0.002 ppm. + * + * @data pointer to array of double to be encoded (need memorycont. repr.) + * @dataSize number of doubles from *data to encode + * @result pointer to where resulting bytes should be stored + * @fixedPoint the scaling factor used for getting the fixed point repr. + * This is stored in the binary and automatically extracted + * on decoding. + * @return the number of encoded bytes + */ + size_t encodeLinear( + const double *data, + const size_t dataSize, + unsigned char *result, + double fixedPoint); + + /** + * Calls lower level encodeLinear while handling vector sizes appropriately + * + * @data vector of doubles to be encoded + * @result vector of resulting bytes (will be resized to the number of bytes) + */ + void encodeLinear( + const std::vector &data, + std::vector &result, + double fixedPoint); + + /** + * Decodes data encoded by encodeLinear. + * + * result vector guaranteed to be shorter or equal to (|data| - 8) * 2 + * + * Note that this method may throw a const char* if it deems the input data to be corrupt, i.e. + * that the last encoded int does not use the last byte in the data. In addition the last encoded + * int need to use either the last halfbyte, or the second last followed by a 0x0 halfbyte. + * + * @data pointer to array of bytes to be decoded (need memorycont. repr.) + * @dataSize number of bytes from *data to decode + * @result pointer to were resulting doubles should be stored + * @return the number of decoded doubles, or -1 if dataSize < 4 or 4 < dataSize < 8 + */ + size_t decodeLinear( + const unsigned char *data, + const size_t dataSize, + double *result); + + /** + * Calls lower level decodeLinear while handling vector sizes appropriately + * + * Note that this method may throw a const char* if it deems the input data to be corrupt, i.e.. + * that the last encoded int does not use the last byte in the data. In addition the last encoded + * int need to use either the last halfbyte, or the second last followed by a 0x0 halfbyte. + * + * @data vector of bytes to be decoded + * @result vector of resulting double (will be resized to the number of doubles) + */ + void decodeLinear( + const std::vector &data, + std::vector &result); + +///////////////////////////////////////////////////////////// + + + /** + * Encodes the doubles in data by storing the residuals from a linear prediction after first two values. + * + * The resulting binary is the same size as the input data. + * + * This encoding is suitable for typical m/z or retention time binary arrays, and is + * intended to be used before zlib compression to improve compression. + * + * @data pointer to array of doubles to be encoded (need memorycont. repr.) + * @dataSize number of doubles from *data to encode + * @result pointer to were resulting bytes should be stored + */ + size_t encodeSafe( + const double *data, + const size_t dataSize, + unsigned char *result); + + + /** + * Decodes data encoded by encodeSafe. + * + * result vector is the same size as the input data. + * + * Might throw const char* is something goes wrong during decoding. + * + * @data pointer to array of bytes to be decoded (need memorycont. repr.) + * @dataSize number of bytes from *data to decode + * @result pointer to were resulting doubles should be stored + * @return the number of decoded bytes + */ + size_t decodeSafe( + const unsigned char *data, + const size_t dataSize, + double *result); + +///////////////////////////////////////////////////////////// + + /** + * Encodes ion counts by simply rounding to the nearest 4 byte integer, + * and compressing each integer with encodeInt. + * + * The handleable range is therefore 0 -> 4294967294. + * The resulting binary is maximally dataSize * 5 bytes, but much less if the + * data is close to 0 on average. + * + * @data pointer to array of double to be encoded (need memorycont. repr.) + * @dataSize number of doubles from *data to encode + * @result pointer to were resulting bytes should be stored + * @return the number of encoded bytes + */ + size_t encodePic( + const double *data, + const size_t dataSize, + unsigned char *result); + + /** + * Calls lower level encodePic while handling vector sizes appropriately + * + * @data vector of doubles to be encoded + * @result vector of resulting bytes (will be resized to the number of bytes) + */ + void encodePic( + const std::vector &data, + std::vector &result); + + /** + * Decodes data encoded by encodePic + * + * result vector guaranteed to be shorter of equal to |data| * 2 + * + * Note that this method may throw a const char* if it deems the input data to be corrupt, i.e. + * that the last encoded int does not use the last byte in the data. In addition the last encoded + * int need to use either the last halfbyte, or the second last followed by a 0x0 halfbyte. + * + * @data pointer to array of bytes to be decoded (need memorycont. repr.) + * @dataSize number of bytes from *data to decode + * @result pointer to were resulting doubles should be stored + * @return the number of decoded doubles + */ + size_t decodePic( + const unsigned char *data, + const size_t dataSize, + double *result); + + /** + * Calls lower level decodePic while handling vector sizes appropriately + * + * Note that this method may throw a const char* if it deems the input data to be corrupt, i.e. + * that the last encoded int does not use the last byte in the data. In addition the last encoded + * int need to use either the last halfbyte, or the second last followed by a 0x0 halfbyte. + * + * @data vector of bytes to be decoded + * @result vector of resulting double (will be resized to the number of doubles) + */ + void decodePic( + const std::vector &data, + std::vector &result); + +///////////////////////////////////////////////////////////// + + + double optimalSlofFixedPoint( + const double *data, + size_t dataSize); + + /** + * Encodes ion counts by taking the natural logarithm, and storing a + * fixed point representation of this. This is calculated as + * + * unsigned short fp = log(d + 1) * fixedPoint + 0.5 + * + * the result vector is exactly |data| * 2 + 8 bytes long + * + * @data pointer to array of double to be encoded (need memorycont. repr.) + * @dataSize number of doubles from *data to encode + * @result pointer to were resulting bytes should be stored + * @return the number of encoded bytes + */ + size_t encodeSlof( + const double *data, + const size_t dataSize, + unsigned char *result, + double fixedPoint); + + /** + * Calls lower level encodeSlof while handling vector sizes appropriately + * + * @data vector of doubles to be encoded + * @result vector of resulting bytes (will be resized to the number of bytes) + */ + void encodeSlof( + const std::vector &data, + std::vector &result, + double fixedPoint); + + /** + * Decodes data encoded by encodeSlof + * + * The return will include exactly (|data| - 8) / 2 doubles. + * + * Note that this method may throw a const char* if it deems the input data to be corrupt. + * + * @data pointer to array of bytes to be decoded (need memorycont. repr.) + * @dataSize number of bytes from *data to decode + * @result pointer to were resulting doubles should be stored + * @return the number of decoded doubles + */ + size_t decodeSlof( + const unsigned char *data, + const size_t dataSize, + double *result); + + /** + * Calls lower level decodeSlof while handling vector sizes appropriately + * + * Note that this method may throw a const char* if it deems the input data to be corrupt. + * + * @data vector of bytes to be decoded + * @result vector of resulting double (will be resized to the number of doubles) + */ + void decodeSlof( + const std::vector &data, + std::vector &result); + +} // namespace MSNumpress +} // namespace msdata +} // namespace pwiz + +#endif // _MSNUMPRESS_HPP_ diff --git a/src/main/cpp/MSNumpressTest.cpp b/src/main/cpp/MSNumpressTest.cpp new file mode 100644 index 0000000..a56543a --- /dev/null +++ b/src/main/cpp/MSNumpressTest.cpp @@ -0,0 +1,786 @@ +/* + MSNumpressTest.cpp + johan.teleman@immun.lth.se + + Copyright 2013 Johan Teleman + + Licensed under the Apache License, Version 2.0 (the "License"); + you may not use this file except in compliance with the License. + You may obtain a copy of the License at + + http://www.apache.org/licenses/LICENSE-2.0 + + Unless required by applicable law or agreed to in writing, software + distributed under the License is distributed on an "AS IS" BASIS, + WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + See the License for the specific language governing permissions and + limitations under the License. + */ + +/* + Compile and run tests (on LINUX) with + + > g++ MSNumpress.cpp MSNumpressTest.cpp -o test && ./test + + */ + +#include "MSNumpress.hpp" +#include +#include +#include +#include +#include + +using std::cout; +using std::endl; +using std::abs; +using std::max; + + +double ENC_TWO_BYTE_FIXED_POINT = 3000.0; + + + +void encodeLinear1() { + + double mzs[1]; + + mzs[0] = 100.0; + + size_t nMzs = 1; + unsigned char encoded[12]; + size_t encodedBytes = ms::numpress::MSNumpress::encodeLinear(&mzs[0], nMzs, &encoded[0], 100000.0); + + assert(12 == encodedBytes); + assert(0x80 == encoded[8]); + assert(0x96 == encoded[9]); + assert(0x98 == encoded[10]); + assert(0x00 == encoded[11]); + + cout << "+ pass encodeLinear1 " << endl << endl; +} + +void encodeLinear() { + double mzs[4]; + + mzs[0] = 100.0; + mzs[1] = 200.0; + mzs[2] = 300.00005; + mzs[3] = 400.00010; + + size_t nMzs = 4; + unsigned char encoded[20]; + size_t encodedBytes = ms::numpress::MSNumpress::encodeLinear(&mzs[0], nMzs, &encoded[0], 100000.0); + + assert(18 == encodedBytes); + assert(0x80 == encoded[8]); + assert(0x96 == encoded[9]); + assert(0x98 == encoded[10]); + assert(0x00 == encoded[11]); + assert(0x75 == encoded[16]); + assert(0x80 == encoded[17]); + + cout << "+ pass encodeLinear " << endl << endl; +} + +void decodeLinearNice() { + + double mzs[4]; + + mzs[0] = 100.0; + mzs[1] = 200.0; + mzs[2] = 300.00005; + mzs[3] = 400.00010; + + size_t nMzs = 4; + unsigned char encoded[28]; + double fixedPoint = ms::numpress::MSNumpress::optimalLinearFixedPoint(&mzs[0], nMzs); + size_t encodedBytes = ms::numpress::MSNumpress::encodeLinear(&mzs[0], nMzs, &encoded[0], fixedPoint); + + double decoded[4]; + size_t numDecoded = ms::numpress::MSNumpress::decodeLinear(&encoded[0], encodedBytes, &decoded[0]); + assert(4 == numDecoded); + assert(abs(100.0 - decoded[0]) < 0.000005); + assert(abs(200.0 - decoded[1]) < 0.000005); + assert(abs(300.00005 - decoded[2]) < 0.000005); + assert(abs(400.00010 - decoded[3]) < 0.000005); + + cout << "+ pass decodeLinearNice " << endl << endl; +} + +void decodeLinearNiceLowFP() { + + double mzs[7]; + + mzs[0] = 100.0; + mzs[1] = 200.0; + mzs[2] = 300.00005; + mzs[3] = 400.00010; + mzs[4] = 450.00010; + mzs[5] = 455.00010; + mzs[6] = 700.00010; + + size_t nMzs = 7; + unsigned char encoded[33]; // max length is 33 bytes + + // check for fixed points + { + double fixedPoint; + fixedPoint = ms::numpress::MSNumpress::optimalLinearFixedPointMass(&mzs[0], nMzs, 0.1); + assert( abs(5 - fixedPoint) < 0.000005); + + fixedPoint = ms::numpress::MSNumpress::optimalLinearFixedPointMass(&mzs[0], nMzs, 1e-3); + assert( abs(500 - fixedPoint) < 0.000005); + + fixedPoint = ms::numpress::MSNumpress::optimalLinearFixedPointMass(&mzs[0], nMzs, 1e-5); + assert( abs(50000 - fixedPoint) < 0.000005); + + fixedPoint = ms::numpress::MSNumpress::optimalLinearFixedPointMass(&mzs[0], nMzs, 1e-7); + assert( abs(5000000 - fixedPoint) < 0.000005); + + // cannot fulfill accuracy of 1e-8 + fixedPoint = ms::numpress::MSNumpress::optimalLinearFixedPointMass(&mzs[0], nMzs, 1e-8); + assert( abs(-1 - fixedPoint) < 0.000005); + } + + { + double fixedPoint = ms::numpress::MSNumpress::optimalLinearFixedPointMass(&mzs[0], nMzs, 0.001); + size_t encodedBytes = ms::numpress::MSNumpress::encodeLinear(&mzs[0], nMzs, &encoded[0], fixedPoint); + + double decoded[7]; + size_t numDecoded = ms::numpress::MSNumpress::decodeLinear(&encoded[0], encodedBytes, &decoded[0]); + assert(25 == encodedBytes); + assert(7 == numDecoded); + + assert(abs(100.0 - decoded[0]) < 0.001); + assert(abs(200.0 - decoded[1]) < 0.001); + assert(abs(300.00005 - decoded[2]) < 0.001); + assert(abs(400.00010 - decoded[3]) < 0.001); + } + + double mz_err[5]; + double encodedLength[5]; + + // for higher accuracy, we get longer encoded lengths + mz_err[0] = 0.1; encodedLength[0] = 22; + mz_err[1] = 1e-3; encodedLength[1] = 25; + mz_err[2] = 1e-5; encodedLength[2] = 29; + mz_err[3] = 1e-6; encodedLength[3] = 30; + mz_err[4] = 1e-7; encodedLength[4] = 31; + + for (int k = 0; k < 4; k++) + { + double fixedPoint = ms::numpress::MSNumpress::optimalLinearFixedPointMass(&mzs[0], nMzs, mz_err[k]); + size_t encodedBytes = ms::numpress::MSNumpress::encodeLinear(&mzs[0], nMzs, &encoded[0], fixedPoint); + + double decoded[7]; + size_t numDecoded = ms::numpress::MSNumpress::decodeLinear(&encoded[0], encodedBytes, &decoded[0]); + assert( encodedLength[k] == encodedBytes); + assert(7 == numDecoded); + assert(abs(100.0 - decoded[0]) < mz_err[k]); + assert(abs(200.0 - decoded[1]) < mz_err[k]); + assert(abs(300.00005 - decoded[2]) < mz_err[k]); + assert(abs(400.00010 - decoded[3]) < mz_err[k]); + assert(abs(450.00010 - decoded[4]) < mz_err[k]); + assert(abs(455.00010 - decoded[5]) < mz_err[k]); + assert(abs(700.00010 - decoded[6]) < mz_err[k]); + } + + cout << "+ pass decodeLinearNiceFP " << endl << endl; +} + +void decodeLinearWierd() { + double mzs[4]; + + mzs[0] = 100.0; + mzs[1] = 200.0; + mzs[2] = 300.00005; + mzs[3] = 0.00010; + + size_t nMzs = 4; + unsigned char encoded[28]; + double fixedPoint = ms::numpress::MSNumpress::optimalLinearFixedPoint(&mzs[0], nMzs); + size_t encodedBytes = ms::numpress::MSNumpress::encodeLinear(&mzs[0], nMzs, &encoded[0], fixedPoint); + + double decoded[4]; + size_t numDecoded = ms::numpress::MSNumpress::decodeLinear(&encoded[0], encodedBytes, &decoded[0]); + assert(4 == numDecoded); + assert(abs(100.0 - decoded[0]) < 0.000005); + assert(abs(200.0 - decoded[1]) < 0.000005); + assert(abs(300.00005 - decoded[2]) < 0.000005); + assert(abs(0.00010 - decoded[3]) < 0.000005); + + cout << "+ pass decodeLinearWierd " << endl << endl; +} + +void decodeLinearWierd_llong_overflow() { + double mzs[4]; + + mzs[0] = 100.0; + mzs[1] = 200.0; + mzs[2] = 30000000.00005; + mzs[3] = 0.000010; + + size_t nMzs = 4; + unsigned char encoded[28]; + + try { + ms::numpress::MSNumpress::encodeLinear(&mzs[0], nMzs, &encoded[0], 1000000000000); + cout << "- fail test decodeLinearWierd_llong_overflow: didn't throw exception for corrupt input " << endl << endl; + assert(0 == 1); + } catch (const char *err) { + assert( std::string(err) == std::string("[MSNumpress::encodeLinear] Next number overflows LLONG_MAX.")); + } + cout << "+ pass decodeLinearWierd_llong_overflow " << endl << endl; +} + +void decodeLinearWierd_int_overflow() { + double mzs[4]; + + mzs[0] = 100.0; + mzs[1] = 200.0; + mzs[2] = 30000000.00005; + mzs[3] = 0.00006; + + size_t nMzs = 4; + unsigned char encoded[28]; + + try { + ms::numpress::MSNumpress::encodeLinear(&mzs[0], nMzs, &encoded[0], 1000000000); + cout << "- fail test decodeLinearWierd_int_overflow: didn't throw exception for corrupt input " << endl << endl; + assert(0 == 1); + } catch (const char *err) { + assert( std::string(err) == std::string("[MSNumpress::encodeLinear] Cannot encode a number that exceeds the bounds of [-INT_MAX, INT_MAX].")); + } + cout << "+ pass decodeLinearWierd3 " << endl << endl; +} + +void decodeLinearWierd_int_underflow() { + double mzs[4]; + + mzs[0] = 30000000.00005; + mzs[1] = 60000000.00005; + mzs[2] = 30000000.00005; + mzs[3] = 0.00006; + + size_t nMzs = 4; + unsigned char encoded[28]; + + try { + ms::numpress::MSNumpress::encodeLinear(&mzs[0], nMzs, &encoded[0], 1000000000); + cout << "- fail test decodeLinearWierd_int_underflow: didn't throw exception for corrupt input " << endl << endl; + assert(0 == 1); + } catch (const char *err) { + assert( std::string(err) == std::string("[MSNumpress::encodeLinear] Cannot encode a number that exceeds the bounds of [-INT_MAX, INT_MAX].")); + } + cout << "+ pass decodeLinearWierd3 " << endl << endl; +} + +void decodeLinearCorrupt1() { + unsigned char encoded[20] = {0}; + double decoded[4]; + + try { + ms::numpress::MSNumpress::decodeLinear(&encoded[0], 11, &decoded[0]); + cout << "- fail decodeLinearCorrupt1: didn't throw exception for corrupt input " << endl << endl; + assert(0 == 1); + } catch (const char *err) { + + } + + try { + ms::numpress::MSNumpress::decodeLinear(&encoded[0], 14, &decoded[0]); + cout << "- fail decodeLinearCorrupt1: didn't throw exception for corrupt input " << endl << endl; + assert(0 == 1); + } catch (const char *err) { + + } + + cout << "+ pass decodeLinearCorrupt 1 " << endl << endl; +} + + +void decodeLinearCorrupt2() { + + double mzs[4]; + + mzs[0] = 100.0; + mzs[1] = 200.0; + mzs[2] = 300.00005; + mzs[3] = 0.00010; + + size_t nMzs = 4; + unsigned char encoded[28]; + double fixedPoint = ms::numpress::MSNumpress::optimalLinearFixedPoint(&mzs[0], nMzs); + size_t encodedBytes = ms::numpress::MSNumpress::encodeLinear(&mzs[0], nMzs, &encoded[0], fixedPoint); + + double decoded[4]; + try { + ms::numpress::MSNumpress::decodeLinear(&encoded[0], encodedBytes-1, &decoded[0]); + cout << "- fail decodeLinearCorrupt2: didn't throw exception for corrupt input " << endl << endl; + assert(0 == 1); + } catch (const char *err) { + + } + + cout << "+ pass decodeLinearCorrupt 2 " << endl << endl; +} + + + +void optimalLinearFixedPoint() { + srand(123459); + + size_t n = 1000; + double mzs[1000]; + mzs[0] = 300 + (rand() % 1000) / 1000.0; + for (size_t i=1; i= mLim) { + cout << "error " << error << " above limit " << mLim << endl; + assert(error < mLim); + } + } + cout << "+ size compressed: " << encodedBytes / double(n*8) * 100 << "% " << endl; + cout << "+ max error: " << m << " limit: " << mLim << endl; + cout << "+ pass encodeDecodeLinearStraight " << endl << endl; +} + + + +void encodeDecodeSafeStraight() { + double error; + double eLim = 1.0e-300; + size_t n = 15; + double mzs[15]; + for (size_t i=0; i= eLim) { + cout << "error " << error << " is non-zero ( >= " << eLim << " )" << endl; + assert(error == 0); + } + } + cout << "+ pass encodeDecodeSafeStraight " << endl << endl; +} + + + +void encodeDecodeSafe() { + srand(123459); + + double error; + double eLim = 1.0e-300; + size_t n = 1000; + double mzs[1000]; + mzs[0] = 300 + rand() / double(RAND_MAX); + for (size_t i=1; i= eLim) { + cout << "error " << error << " is non-zero ( >= " << eLim << " )" << endl; + assert(error == 0); + } + } + cout << "+ pass encodeDecodeSafe " << endl << endl; +} + + + +void encodeDecodeLinear() { + srand(123459); + + size_t n = 1000; + double mzs[1000]; + mzs[0] = 300 + rand() / double(RAND_MAX); + for (size_t i=1; i= mLim) { + cout << "error " << error << " above limit " << mLim << endl; + assert(error < mLim); + } + } + cout << "+ size compressed: " << encodedBytes / double(n*8) * 100 << "% " << endl; + cout << "+ max error: " << m << " limit: " << mLim << endl; + cout << "+ pass encodeDecodeLinear " << endl << endl; +} + + + +void encodeDecodeLinear5() { + srand(123662); + + size_t n = 1000; + double mzs[1000]; + mzs[0] = 100 + (rand() % 1000) / 1000.0; + for (size_t i=1; i= mLim) { + cout << endl << ics[i] << " " << decoded[i] << endl; + assert(error < mLim); + } + } else { + error = abs((ics[i] - decoded[i]) / ((ics[i] + decoded[i])/2)); + rm = max(rm, error); + if (error >= rmLim) { + cout << endl << ics[i] << " " << decoded[i] << endl; + assert(error < rmLim); + } + } + cout << "+ max error: " << m << " limit: " << mLim << endl; + cout << "+ max rel error: " << rm << " limit: " << rmLim << endl; + cout << "+ pass encodeDecodeSlof " << endl << endl; +} + + + +void encodeDecodeSlof5() { + srand(123459); + + size_t n = 1000; + double ics[1000]; + for (size_t i=0; i result; + + // set data to [ 100, 102, 140, 92, 33, 80, 145 ]; // Base64 is "ZGaMXCFQkQ==" + std::vector data; + data.resize(32); + data[0] = 100; + data[1] = 102; + data[2] = 140; + data[3] = 92; + data[4] = 33; + data[5] = 80; + data[6] = 145; + + try { + ms::numpress::MSNumpress::decodePic(data, result); + cout << "- fail testErroneousDecodePic: didn't throw exception for corrupt input " << endl << endl; + assert(0 == 1); + } catch (const char *err) { + + } + + cout << "+ pass testErroneousDecodePic " << endl << endl; +} + + +int main() { + optimalLinearFixedPoint(); + optimalLinearFixedPointMass(); + encodeLinear1(); + encodeLinear(); + decodeLinearNice(); + decodeLinearNiceLowFP(); + decodeLinearWierd(); + decodeLinearCorrupt1(); + decodeLinearCorrupt2(); + encodeDecodeLinearStraight(); + encodeDecodeLinear(); + encodeDecodePic(); + encodeDecodeSafeStraight(); + encodeDecodeSafe(); + optimalSlofFixedPoint(); + encodeDecodeSlof(); + encodeDecodeLinear5(); + encodeDecodePic5(); + encodeDecodeSlof5(); + testErroneousDecodePic(); + + cout << "=== all tests succeeded! ===" << endl; + return 0; +} diff --git a/src/main/csharp/MSNumpress.cs b/src/main/csharp/MSNumpress.cs new file mode 100644 index 0000000..4b511c5 --- /dev/null +++ b/src/main/csharp/MSNumpress.cs @@ -0,0 +1,606 @@ +/* + MSNumpress.cs + rfellers@gmail.com + Copyright 2017 Ryan Fellers + + Based on: + MSNumpress.java and IntDecoder.java + johan.teleman@immun.lth.se + + Copyright 2013 Johan Teleman + Licensed under the Apache License, Version 2.0 (the "License"); + you may not use this file except in compliance with the License. + You may obtain a copy of the License at + http://www.apache.org/licenses/LICENSE-2.0 + Unless required by applicable law or agreed to in writing, software + distributed under the License is distributed on an "AS IS" BASIS, + WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + See the License for the specific language governing permissions and + limitations under the License. + */ + +using System; +using System.Diagnostics; + +/// +/// Implementations of two compression schemes for numeric data from mass spectrometers. +/// +public class MSNumpress +{ + /// + /// MS Numpress linear prediction compression + /// + public const string ACC_NUMPRESS_LINEAR = "MS:1002312"; + + /// + /// MS Numpress positive integer compression + /// + public const string ACC_NUMPRESS_PIC = "MS:1002313"; + + /// + /// MS Numpress short logged float compression + /// + public const string ACC_NUMPRESS_SLOF = "MS:1002314"; + + /// + /// Convenience function for decoding binary data encoded by MSNumpress. + /// + /// The PSI-MS obo CV accession of the encoded data. + /// array of double to be encoded. + /// number of doubles from data to encode. + /// The decoded doubles + /// + /// Cannot decode numLin data, need at least 8 initial bytes for fixed point. + /// or + /// Corrupt numLin data! + /// or + /// Cannot decode numPic data, need at least 8 initial bytes for fixed point. + /// or + /// Corrupt numPic data! + /// or + /// '" + cvAccession + "' is not a numpress compression term + /// + /// + /// If the passed cvAccession is one of + /// + /// ACC_NUMPRESS_LINEAR = "MS:1002312" + /// ACC_NUMPRESS_PIC = "MS:1002313" + /// ACC_NUMPRESS_SLOF = "MS:1002314" + /// + /// the corresponding decode function will be called. + /// + public static double[] decode(string cvAccession, byte[] data, int dataSize) + { + if (cvAccession == ACC_NUMPRESS_LINEAR) + { + if (dataSize < 8 || data.Length < 8) + throw new ArgumentException("Cannot decode numLin data, need at least 8 initial bytes for fixed point."); + + double[] buffer = new double[dataSize * 2]; + int nbrOfDoubles = MSNumpress.decodeLinear(data, dataSize, buffer); + if (nbrOfDoubles < 0) + throw new ArgumentException("Corrupt numLin data!"); + + double[] result = new double[nbrOfDoubles]; + Array.Copy(buffer, 0, result, 0, nbrOfDoubles); + + return result; + } + + if (cvAccession == ACC_NUMPRESS_SLOF) + { + double[] result = new double[(dataSize - 8) / 2]; + MSNumpress.decodeSlof(data, dataSize, result); + + return result; + } + + if (cvAccession == ACC_NUMPRESS_PIC) + { + if (dataSize < 8 || data.Length < 8) + throw new ArgumentException("Cannot decode numPic data, need at least 8 initial bytes for fixed point."); + + double[] buffer = new double[dataSize * 2]; + int nbrOfDoubles = MSNumpress.decodePic(data, dataSize, buffer); + if (nbrOfDoubles < 0) + throw new ArgumentException("Corrupt numPic data!"); + + double[] result = new double[nbrOfDoubles]; + Array.Copy(buffer, 0, result, 0, nbrOfDoubles); + return result; + + } + + throw new ArgumentException("'" + cvAccession + "' is not a numpress compression term"); + } + + /// + /// This encoding works on a 4 byte integer, by truncating initial zeros or ones. + /// + /// the int to be encoded + /// the byte array were halfbytes are stored + /// position in res were halfbytes are written + /// the number of resulting halfbytes + /// + /// If the initial (most significant) half byte is 0x0 or 0xf, the number of such + /// halfbytes starting from the most significant is stored in a halfbyte. This initial + /// count is then followed by the rest of the ints halfbytes, in little-endian order. + /// A count halfbyte c of + /// + /// 0 <= c <= 8 is interpreted as an initial c 0x0 halfbytes + /// 9 <= c <= 15 is interpreted as an initial (c-8) 0xf halfbytes + /// + /// Ex: + /// int c rest + /// 0 => 0x8 + /// -1 => 0xf 0xf + /// 23 => 0x6 0x7 0x1 + /// + /// @x the int to be encoded + /// @res the byte array were halfbytes are stored + /// @resOffset position in res were halfbytes are written + /// @return the number of resulting halfbytes + /// + public static int encodeInt(long x, byte[] res, int resOffset) + { + byte i, l; + long m; + long mask = 0xf0000000; + long init = x & mask; + + if (init == 0) + { + l = 8; + for (i = 0; i < 8; i++) + { + m = mask >> (4 * i); + if ((x & m) != 0) + { + l = i; + break; + } + } + res[resOffset] = l; + for (i = l; i < 8; i++) + res[resOffset + 1 + i - l] = (byte)(0xf & (x >> (4 * (i - l)))); + + return 1 + 8 - l; + + } + else if (init == mask) + { + l = 7; + for (i = 0; i < 8; i++) + { + m = mask >> (4 * i); + if ((x & m) != m) + { + l = i; + break; + } + } + res[resOffset] = (byte)(l | 8); + for (i = l; i < 8; i++) + res[resOffset + 1 + i - l] = (byte)(0xf & (x >> (4 * (i - l)))); + + return 1 + 8 - l; + + } + else + { + res[resOffset] = 0; + for (i = 0; i < 8; i++) + res[resOffset + 1 + i] = (byte)(0xf & (x >> (4 * i))); + + return 9; + + } + } + + public static void encodeFixedPoint(double fixedPoint, byte[] result) + { + //long fp = double.doubleToLongBits(fixedPoint); + long fp = BitConverter.DoubleToInt64Bits(fixedPoint); // RTF + + for (int i = 0; i < 8; i++) + { + result[7 - i] = (byte)((fp >> (8 * i)) & 0xff); + } + } + + public static double decodeFixedPoint(byte[] data) + { + long fp = 0; + for (int i = 0; i < 8; i++) + { + fp = fp | ((0xFFL & data[7 - i]) << (8 * i)); + } + + //return double.longBitsToDouble(fp); + return BitConverter.Int64BitsToDouble(fp); + } + + ///////////////////////////////////////////////////////////////////////////////// + + public static double optimalLinearFixedPoint(double[] data, int dataSize) + { + if (dataSize == 0) return 0; + if (dataSize == 1) return Math.Floor(0xFFFFFFFFL / data[0]); + double maxDouble = Math.Max(data[0], data[1]); + + for (int i = 2; i < dataSize; i++) + { + double extrapol = data[i - 1] + (data[i - 1] - data[i - 2]); + double diff = data[i] - extrapol; + maxDouble = Math.Max(maxDouble, Math.Ceiling(Math.Abs(diff) + 1)); + } + + return Math.Floor(0x7FFFFFFFL / maxDouble); + } + + /// + /// Encodes data using MS Numpress linear prediction compression. + /// + /// array of doubles to be encoded + /// number of doubles from data to encode + /// array were resulting bytes should be stored + /// the scaling factor used for getting the fixed point repr. This is stored in the binary and automatically extracted on decoding. + /// the number of encoded bytes + /// + /// Encodes the doubles in data by first using a + /// - lossy conversion to a 4 byte 5 decimal fixed point repressentation + /// - storing the residuals from a linear prediction after first two values + /// - encoding by encodeInt (see above) + /// + /// The resulting binary is maximally 8 + dataSize * 5 bytes, but much less if the + /// data is reasonably smooth on the first order. + /// + /// This encoding is suitable for typical m/z or retention time binary arrays. + /// On a test set, the encoding was empirically show to be accurate to at least 0.002 ppm. + /// + public static int encodeLinear(double[] data, int dataSize, byte[] result, double fixedPoint) + { + long[] ints = new long[3]; + int i; + int ri = 16; + byte[] halfBytes = new byte[10]; + int halfByteCount = 0; + int hbi; + long extrapol; + long diff; + + encodeFixedPoint(fixedPoint, result); + + if (dataSize == 0) return 8; + + ints[1] = (long)(data[0] * fixedPoint + 0.5); + for (i = 0; i < 4; i++) + { + result[8 + i] = (byte)((ints[1] >> (i * 8)) & 0xff); + } + + if (dataSize == 1) return 12; + + ints[2] = (long)(data[1] * fixedPoint + 0.5); + for (i = 0; i < 4; i++) + { + result[12 + i] = (byte)((ints[2] >> (i * 8)) & 0xff); + } + + halfByteCount = 0; + ri = 16; + + for (i = 2; i < dataSize; i++) + { + ints[0] = ints[1]; + ints[1] = ints[2]; + ints[2] = (long)(data[i] * fixedPoint + 0.5); + extrapol = ints[1] + (ints[1] - ints[0]); + diff = ints[2] - extrapol; + halfByteCount += encodeInt(diff, halfBytes, halfByteCount); + + for (hbi = 1; hbi < halfByteCount; hbi += 2) + result[ri++] = (byte)((halfBytes[hbi - 1] << 4) | (halfBytes[hbi] & 0xf)); + + if (halfByteCount % 2 != 0) + { + halfBytes[0] = halfBytes[halfByteCount - 1]; + halfByteCount = 1; + } + else + halfByteCount = 0; + } + + if (halfByteCount == 1) + result[ri++] = (byte)(halfBytes[0] << 4); + + return ri; + } + + /// + /// Decodes data using MS Numpress linear prediction compression. + /// + /// array of bytes to be decoded + /// number of bytes from data to decode + /// array were resulting doubles should be stored + /// the number of decoded doubles, or -1 if dataSize < 4 or 4 < dataSize < 8 + /// + /// Result vector guaranteed to be shorter or equal to (|data| - 8) * 2 + /// + /// Note that this method may throw a ArrayIndexOutOfBoundsException if it deems the input data to + /// be corrupt, i.e. that the last encoded int does not use the last byte in the data. In addition + /// the last encoded int need to use either the last halfbyte, or the second last followed by a + /// 0x0 halfbyte. + /// + public static int decodeLinear(byte[] data, int dataSize, double[] result) + { + int ri = 2; + long[] ints = new long[3]; + long extrapol = 0; + long y; + var dec = new IntDecoder(data, 16); + + if (dataSize == 8) return 0; + if (dataSize < 8) return -1; + double fixedPoint = decodeFixedPoint(data); + if (dataSize < 12) return -1; + + ints[1] = 0; + for (int i = 0; i < 4; i++) + { + ints[1] = ints[1] | ((0xFFL & data[8 + i]) << (i * 8)); + } + result[0] = ints[1] / fixedPoint; + + if (dataSize == 12) return 1; + if (dataSize < 16) return -1; + + ints[2] = 0; + for (int i = 0; i < 4; i++) + { + ints[2] = ints[2] | ((0xFFL & data[12 + i]) << (i * 8)); + } + result[1] = ints[2] / fixedPoint; + + while (dec.pos < dataSize) + { + if (dec.pos == (dataSize - 1) && dec.half) + if ((data[dec.pos] & 0xf) != 0x8) + break; + + ints[0] = ints[1]; + ints[1] = ints[2]; + ints[2] = dec.next(); + + extrapol = ints[1] + (ints[1] - ints[0]); + y = extrapol + ints[2]; + result[ri++] = y / fixedPoint; + ints[2] = y; + } + + return ri; + } + + ///////////////////////////////////////////////////////////////////////////////// + + /// + /// Encodes ion counts by simply rounding to the nearest 4 byte integer, and compressing each integer with encodeInt. + /// + /// array of doubles to be encoded + /// number of doubles from data to encode + /// array were resulting bytes should be stored + /// the number of encoded bytes + /// + /// The handleable range is therefore 0 -> 4294967294. + /// The resulting binary is maximally dataSize * 5 bytes, but much less if the + /// data is close to 0 on average. + /// + public static int encodePic(double[] data, int dataSize, byte[] result) + { + long count; + int ri = 0; + int hbi = 0; + byte[] halfBytes = new byte[10]; + int halfByteCount = 0; + + for (int i = 0; i < dataSize; i++) + { + count = (long)(data[i] + 0.5); + halfByteCount += encodeInt(count, halfBytes, halfByteCount); + + for (hbi = 1; hbi < halfByteCount; hbi += 2) + result[ri++] = (byte)((halfBytes[hbi - 1] << 4) | (halfBytes[hbi] & 0xf)); + + if (halfByteCount % 2 != 0) + { + halfBytes[0] = halfBytes[halfByteCount - 1]; + halfByteCount = 1; + } + else + halfByteCount = 0; + + } + if (halfByteCount == 1) + result[ri++] = (byte)(halfBytes[0] << 4); + + return ri; + } + + /// + /// Decodes data encoded by encodePic. + /// + /// array of bytes to be decoded (need memorycont. repr.) + /// number of bytes from data to decode + /// array were resulting doubles should be stored + /// the number of decoded doubles + /// + /// Result vector guaranteed to be shorter of equal to |data| * 2 + /// + /// Note that this method may throw a ArrayIndexOutOfBoundsException if it deems the input data to + /// be corrupt, i.e. that the last encoded int does not use the last byte in the data. In addition + /// the last encoded int need to use either the last halfbyte, or the second last followed by a + /// 0x0 halfbyte. + /// + public static int decodePic(byte[] data, int dataSize, double[] result) + { + int ri = 0; + long count; + IntDecoder dec = new IntDecoder(data, 0); + + while (dec.pos < dataSize) + { + if (dec.pos == (dataSize - 1) && dec.half) + if ((data[dec.pos] & 0xf) != 0x8) + break; + + count = dec.next(); + result[ri++] = count; + } + return ri; + } + + ///////////////////////////////////////////////////////////////////////////////// + + public static double optimalSlofFixedPoint(double[] data, int dataSize) + { + if (dataSize == 0) return 0; + + double maxDouble = 1; + double x; + double fp; + + for (int i = 0; i < dataSize; i++) + { + x = Math.Log(data[i] + 1); + maxDouble = Math.Max(maxDouble, x); + } + + fp = Math.Floor(0xFFFF / maxDouble); + + return fp; + } + + /// + /// Encodes ion counts by taking the natural logarithm, and storing a fixed point representation of this. + /// + /// array of doubles to be encoded + /// number of doubles from data to encode + /// array were resulting bytes should be stored + /// the scaling factor used for getting the fixed point repr. This is stored in the binary and automatically extracted on decoding. + /// the number of encoded bytes + /// + /// Encodes ion counts by taking the natural logarithm, and storing a + /// fixed point representation of this. This is calculated as + /// + /// unsigned short fp = log(d+1) * fixedPoint + 0.5 + /// + /// the result vector is exactly |data| * 2 + 8 bytes long + /// + public static int encodeSlof(double[] data, int dataSize, byte[] result, double fixedPoint) + { + int x; + int ri = 8; + + encodeFixedPoint(fixedPoint, result); + + for (int i = 0; i < dataSize; i++) + { + x = (int)(Math.Log(data[i] + 1) * fixedPoint + 0.5); + + result[ri++] = (byte)(0xff & x); + result[ri++] = (byte)(x >> 8); + } + return ri; + } + + /// + /// Decodes data encoded by encodeSlof. + /// + /// array of bytes to be decoded (need memorycont. repr.) + /// number of bytes from data to decode + /// array were resulting doubles should be stored + /// the number of decoded doubles + /// + /// The result vector will be exactly (|data| - 8) / 2 doubles. + /// returns the number of doubles read, or -1 is there is a problem decoding. + /// + public static int decodeSlof(byte[] data, int dataSize, double[] result) + { + int x; + int ri = 0; + + if (dataSize < 8) return -1; + double fixedPoint = decodeFixedPoint(data); + + if (dataSize % 2 != 0) return -1; + + for (int i = 8; i < dataSize; i += 2) + { + x = (0xff & data[i]) | ((0xff & data[i + 1]) << 8); + result[ri++] = Math.Exp((0xffff & x) / fixedPoint) - 1; + } + return ri; + } + + /// + /// Decodes ints from the half bytes in bytes. Lossless reverse of encodeInt, although not symmetrical in input arguments. + /// + public class IntDecoder + { + public int pos = 0; + public bool half = false; + public byte[] bytes; + + public IntDecoder(byte[] _bytes, int _pos) + { + bytes = _bytes; + pos = _pos; + } + + public long next() + { + int head; + int i, n; + long res = 0; + long mask, m; + int hb; + + if (!half) + head = (0xff & bytes[pos]) >> 4; + else + head = 0xf & bytes[pos++]; + + half = !half; + + if (head <= 8) + n = head; + else + { + // leading ones, fill in res + n = head - 8; + mask = unchecked((int)0xF0000000); + + for (i = 0; i < n; i++) + { + m = mask >> (4 * i); + res = res | m; + } + } + + if (n == 8) return 0; + + for (i = n; i < 8; i++) + { + if (!half) + hb = (0xff & bytes[pos]) >> 4; + else + hb = 0xf & bytes[pos++]; + + res = (int)res | (hb << ((i - n) * 4)); + half = !half; + } + + return res; + } + } +} \ No newline at end of file diff --git a/src/main/csharp/MSNumpressTest.cs b/src/main/csharp/MSNumpressTest.cs new file mode 100644 index 0000000..9bdc71d --- /dev/null +++ b/src/main/csharp/MSNumpressTest.cs @@ -0,0 +1,333 @@ +/* + MSNumpressTest.cs + rfellers@gmail.com + Copyright 2017 Ryan Fellers + + Based on: + + MSNumpressTest.java + johan.teleman@immun.lth.se + + Copyright 2013 Johan Teleman + + Licensed under the Apache License, Version 2.0 (the "License"); + you may not use this file except in compliance with the License. + You may obtain a copy of the License at + + http://www.apache.org/licenses/LICENSE-2.0 + + Unless required by applicable law or agreed to in writing, software + distributed under the License is distributed on an "AS IS" BASIS, + WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + See the License for the specific language governing permissions and + limitations under the License. + */ + +using System; +using Microsoft.VisualStudio.TestTools.UnitTesting; +using System.Linq; + +[TestClass] +public class MSNumpressTest +{ + [TestMethod] + public void encodeInt() + { + byte[] res = new byte[10]; + int l; + l = MSNumpress.encodeInt(0L, res, 0); + Assert.AreEqual(1, l); + Assert.AreEqual(8, res[0] & 0xf); + + l = MSNumpress.encodeInt(-1L, res, 0); + Assert.AreEqual(2, l); + Assert.AreEqual(0xf, res[0] & 0xf); + Assert.AreEqual(0xf, res[1] & 0xf); + + l = MSNumpress.encodeInt(35L, res, 0); + Assert.AreEqual(3, l); + Assert.AreEqual(6, res[0] & 0xf); + Assert.AreEqual(0x3, res[1] & 0xf); + Assert.AreEqual(0x2, res[2] & 0xf); + + l = MSNumpress.encodeInt(370000000L, res, 0); + Assert.AreEqual(9, l); + Assert.AreEqual(0, res[0] & 0xf); + Assert.AreEqual(0x0, res[1] & 0xf); + Assert.AreEqual(0x8, res[2] & 0xf); + Assert.AreEqual(0x0, res[3] & 0xf); + Assert.AreEqual(0xc, res[4] & 0xf); + Assert.AreEqual(0xd, res[5] & 0xf); + Assert.AreEqual(0x0, res[6] & 0xf); + Assert.AreEqual(0x6, res[7] & 0xf); + Assert.AreEqual(0x1, res[8] & 0xf); + } + + [TestMethod] + public void decodeInt() + { + byte[] res = new byte[10]; + res[0] = 0x75; + res[1] = 0x87; + res[2] = 0x10; + res[3] = 0x08; + res[4] = 0x0c; + res[5] = 0xd0; + res[6] = 0x61; + var dec = new MSNumpress.IntDecoder(res, 0); + + long l; + l = dec.next(); + Assert.AreEqual(5, l); + + l = dec.next(); + Assert.AreEqual(0, l); + + l = dec.next(); + Assert.AreEqual(1, l); + + l = dec.next(); + Assert.AreEqual(370000000L, l); + } + + [TestMethod] + public void encodeFixedPoint() + { + byte[] encoded = new byte[8]; + MSNumpress.encodeFixedPoint(1.00, encoded); + Assert.AreEqual(0x3f, 0xff & encoded[0]); + Assert.AreEqual(0xf0, 0xff & encoded[1]); + Assert.AreEqual(0x0, 0xff & encoded[2]); + Assert.AreEqual(0x0, 0xff & encoded[3]); + Assert.AreEqual(0x0, 0xff & encoded[4]); + Assert.AreEqual(0x0, 0xff & encoded[5]); + Assert.AreEqual(0x0, 0xff & encoded[6]); + Assert.AreEqual(0x0, 0xff & encoded[7]); + } + + [TestMethod] + public void encodeDecodeFixedPoint() + { + double fp = 300.21941382293625; + byte[] encoded = new byte[8]; + MSNumpress.encodeFixedPoint(fp, encoded); + double decoded = MSNumpress.decodeFixedPoint(encoded); + Assert.AreEqual(fp, decoded, 0); + } + + [TestMethod] + public void encodeLinear() + { + double[] mzs = { 100.0, 200.0, 300.00005, 400.00010 }; + byte[] encoded = new byte[40]; + int encodedBytes = MSNumpress.encodeLinear(mzs, 4, encoded, 100000.0); + Assert.AreEqual(18, encodedBytes); + Assert.AreEqual(0x80, 0xff & encoded[8]); + Assert.AreEqual(0x96, 0xff & encoded[9]); + Assert.AreEqual(0x98, 0xff & encoded[10]); + Assert.AreEqual(0x00, 0xff & encoded[11]); + Assert.AreEqual(0x75, 0xff & encoded[16]); + Assert.AreEqual(0x80, 0xf0 & encoded[17]); + } + + [TestMethod] + public void encodeDecodeLinearEmpty() + { + byte[] encoded = new byte[8]; + int encodedBytes = MSNumpress.encodeLinear(new double[0], 0, encoded, 100000.0); + double[] decoded = new double[0]; + int decodedBytes = MSNumpress.decodeLinear(encoded, 8, decoded); + Assert.AreEqual(0, decodedBytes); + } + + [TestMethod] + public void decodeLinearNice() + { + double[] mzs = { 100.0, 200.0, 300.00005, 400.00010 }; + byte[] encoded = new byte[28]; + int encodedBytes = MSNumpress.encodeLinear(mzs, 4, encoded, 100000.0); + double[] decoded = new double[4]; + int numDecoded = MSNumpress.decodeLinear(encoded, encodedBytes, decoded); + Assert.AreEqual(4, numDecoded); + Assert.AreEqual(100.0, decoded[0], 0.000005); + Assert.AreEqual(200.0, decoded[1], 0.000005); + Assert.AreEqual(300.00005, decoded[2], 0.000005); + Assert.AreEqual(400.00010, decoded[3], 0.000005); + } + + [TestMethod] + public void decodeLinearWierd() + { + double[] mzs = { 100.0, 200.0, 4000.00005, 0.00010 }; + byte[] encoded = new byte[28]; + double fixedPoint = MSNumpress.optimalLinearFixedPoint(mzs, 4); + int encodedBytes = MSNumpress.encodeLinear(mzs, 4, encoded, fixedPoint); + double[] decoded = new double[4]; + int numDecoded = MSNumpress.decodeLinear(encoded, encodedBytes, decoded); + Assert.AreEqual(100.0, decoded[0], 0.000005); + Assert.AreEqual(200.0, decoded[1], 0.000005); + Assert.AreEqual(4000.00005, decoded[2], 0.000005); + Assert.AreEqual(0.00010, decoded[3], 0.000005); + } + + [TestMethod] + public void encodeDecodeLinear() + { + Random random = new Random(); + int n = 1000; + double[] mzs = new double[n]; + mzs[0] = 300 + random.NextDouble(); + for (int i = 1; i < n; i++) + mzs[i] = mzs[i - 1] + random.NextDouble(); + + byte[] encoded = new byte[n * 5]; + double fixedPoint = MSNumpress.optimalLinearFixedPoint(mzs, n); + int encodedBytes = MSNumpress.encodeLinear(mzs, n, encoded, fixedPoint); + double[] decoded = new double[n]; + int decodedDoubles = MSNumpress.decodeLinear(encoded, encodedBytes, decoded); + + Assert.AreEqual(n, decodedDoubles); + + var list = Enumerable.Range(0, 1000).Select(i => decoded[i] / mzs[i]).ToList(); + + for (int i = 0; i < n; i++) + Assert.AreEqual(mzs[i], decoded[i], 0.000005); + } + + [TestMethod] + public void encodeDecodePic() + { + Random random = new Random(); + int n = 1000; + double[] ics = new double[n]; + for (int i = 0; i < n; i++) + ics[i] = Math.Pow(10, 6 * random.NextDouble()); + + byte[] encoded = new byte[n * 5]; + int encodedBytes = MSNumpress.encodePic(ics, n, encoded); + double[] decoded = new double[n]; + int decodedDoubles = MSNumpress.decodePic(encoded, encodedBytes, decoded); + + Assert.AreEqual(n, decodedDoubles); + + for (int i = 0; i < n; i++) + Assert.AreEqual(ics[i], decoded[i], 0.5); + } + + [TestMethod] + public void encodeDecodeSlof() + { + Random random = new Random(); + int n = 1000; + double[] ics = new double[n]; + for (int i = 0; i < n; i++) + ics[i] = Math.Pow(10, 6 * random.NextDouble()); + + byte[] encoded = new byte[n * 2 + 8]; + double fixedPoint = MSNumpress.optimalSlofFixedPoint(ics, n); + int encodedBytes = MSNumpress.encodeSlof(ics, n, encoded, fixedPoint); + double[] decoded = new double[n]; + int decodedDoubles = MSNumpress.decodeSlof(encoded, encodedBytes, decoded); + + Assert.AreEqual(n, decodedDoubles); + + for (int i = 0; i < n; i++) + Assert.AreEqual(0.0, (ics[i] - decoded[i]) / ((ics[i] + decoded[i]) / 2), 0.0005); + } + + [TestMethod] + public void encodeDecodeLinear5() + { + Random random = new Random(); + int n = 1000; + double[] mzs = new double[n]; + mzs[0] = 300 + random.NextDouble(); + for (int i = 1; i < n; i++) + mzs[i] = mzs[i - 1] + random.NextDouble(); + + byte[] encoded = new byte[n * 5]; + double[] decoded = new double[n]; + double[] firstDecoded = new double[n]; + double fixedPoint = MSNumpress.optimalLinearFixedPoint(mzs, n); + + int encodedBytes = MSNumpress.encodeLinear(mzs, n, encoded, fixedPoint); + int decodedDoubles = MSNumpress.decodeLinear(encoded, encodedBytes, decoded); + + for (int i = 0; i < n; i++) + firstDecoded[i] = decoded[i]; + + for (int i = 0; i < 5; i++) + { + MSNumpress.encodeLinear(decoded, n, encoded, fixedPoint); + MSNumpress.decodeLinear(encoded, encodedBytes, decoded); + } + + Assert.AreEqual(n, decodedDoubles); + + for (int i = 0; i < n; i++) + Assert.AreEqual(firstDecoded[i], decoded[i], double.Epsilon); + } + + [TestMethod] + public void encodeDecodePic5() + { + Random random = new Random(); + int n = 1000; + double[] ics = new double[n]; + for (int i = 0; i < n; i++) + ics[i] = Math.Pow(10, 6 * random.NextDouble()); + + byte[] encoded = new byte[n * 5]; + double[] decoded = new double[n]; + double[] firstDecoded = new double[n]; + + int encodedBytes = MSNumpress.encodePic(ics, n, encoded); + int decodedDoubles = MSNumpress.decodePic(encoded, encodedBytes, decoded); + + for (int i = 0; i < n; i++) + firstDecoded[i] = decoded[i]; + + for (int i = 0; i < 5; i++) + { + MSNumpress.encodePic(decoded, n, encoded); + MSNumpress.decodePic(encoded, encodedBytes, decoded); + } + + Assert.AreEqual(n, decodedDoubles); + + for (int i = 0; i < n; i++) + Assert.AreEqual(firstDecoded[i], decoded[i], double.Epsilon); + } + + [TestMethod] + public void encodeDecodeSlof5() + { + Random random = new Random(); + int n = 1000; + double[] ics = new double[n]; + for (int i = 0; i < n; i++) + ics[i] = Math.Pow(10, 6 * random.NextDouble()); + + byte[] encoded = new byte[n * 2 + 8]; + double[] decoded = new double[n]; + double[] firstDecoded = new double[n]; + double fixedPoint = MSNumpress.optimalSlofFixedPoint(ics, n); + + int encodedBytes = MSNumpress.encodeSlof(ics, n, encoded, fixedPoint); + int decodedDoubles = MSNumpress.decodeSlof(encoded, encodedBytes, decoded); + + for (int i = 0; i < n; i++) + firstDecoded[i] = decoded[i]; + + for (int i = 0; i < 5; i++) + { + MSNumpress.encodeSlof(decoded, n, encoded, fixedPoint); + MSNumpress.decodeSlof(encoded, encodedBytes, decoded); + } + + Assert.AreEqual(n, decodedDoubles); + + for (int i = 0; i < n; i++) + Assert.AreEqual(firstDecoded[i], decoded[i], double.Epsilon); + } +} \ No newline at end of file diff --git a/src/main/java/ms/numpress/IntDecoder.java b/src/main/java/ms/numpress/IntDecoder.java new file mode 100644 index 0000000..1d3e857 --- /dev/null +++ b/src/main/java/ms/numpress/IntDecoder.java @@ -0,0 +1,76 @@ +/* + IntDecoder.java + johan.teleman@immun.lth.se + + Copyright 2013 Johan Teleman + + Licensed under the Apache License, Version 2.0 (the "License"); + you may not use this file except in compliance with the License. + You may obtain a copy of the License at + + http://www.apache.org/licenses/LICENSE-2.0 + + Unless required by applicable law or agreed to in writing, software + distributed under the License is distributed on an "AS IS" BASIS, + WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + See the License for the specific language governing permissions and + limitations under the License. + */ + package ms.numpress; + +/** + * Decodes ints from the half bytes in bytes. Lossless reverse of encodeInt, + * although not symmetrical in input arguments. + */ +class IntDecoder { + + int pos = 0; + boolean half = false; + byte[] bytes; + + public IntDecoder(byte[] _bytes, int _pos) { + bytes = _bytes; + pos = _pos; + } + + public long next() { + int head; + int i, n; + long res = 0; + long mask, m; + int hb; + + if (!half) + head = (0xff & bytes[pos]) >> 4; + else + head = 0xf & bytes[pos++]; + + half = !half; + + if (head <= 8) + n = head; + else { + // leading ones, fill in res + n = head - 8; + mask = 0xf0000000; + for (i=0; i> (4*i); + res = res | m; + } + } + + if (n == 8) return 0; + + for (i=n; i<8; i++) { + if (!half) + hb = (0xff & bytes[pos]) >> 4; + else + hb = 0xf & bytes[pos++]; + + res = res | (hb << ((i-n)*4)); + half = !half; + } + + return res; + } +} diff --git a/src/main/java/ms/numpress/MSNumpress.java b/src/main/java/ms/numpress/MSNumpress.java new file mode 100644 index 0000000..97ffdf8 --- /dev/null +++ b/src/main/java/ms/numpress/MSNumpress.java @@ -0,0 +1,572 @@ +/* + MSNumpress.java + johan.teleman@immun.lth.se + + Copyright 2013 Johan Teleman + + Licensed under the Apache License, Version 2.0 (the "License"); + you may not use this file except in compliance with the License. + You may obtain a copy of the License at + + http://www.apache.org/licenses/LICENSE-2.0 + + Unless required by applicable law or agreed to in writing, software + distributed under the License is distributed on an "AS IS" BASIS, + WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + See the License for the specific language governing permissions and + limitations under the License. + */ +package ms.numpress; + +public class MSNumpress { + + ///PSI-MS obo accession numbers. + public static final String ACC_NUMPRESS_LINEAR = "MS:1002312"; + public static final String ACC_NUMPRESS_PIC = "MS:1002313"; + public static final String ACC_NUMPRESS_SLOF = "MS:1002314"; + + + /** + * Convenience function for decoding binary data encoded by MSNumpress. If + * the passed cvAccession is one of + * + * ACC_NUMPRESS_LINEAR = "MS:1002312" + * ACC_NUMPRESS_PIC = "MS:1002313" + * ACC_NUMPRESS_SLOF = "MS:1002314" + * + * the corresponding decode function will be called. + * + * @cvAccession The PSI-MS obo CV accession of the encoded data. + * @data array of double to be encoded + * @dataSize number of doubles from data to encode + * @return The decoded doubles + */ + public static double[] decode( + String cvAccession, + byte[] data, + int dataSize + ) { + if (cvAccession.equals(ACC_NUMPRESS_LINEAR)) { + if (dataSize < 8 || data.length < 8) + throw new IllegalArgumentException("Cannot decode numLin data, need at least 8 initial bytes for fixed point."); + + double[] buffer = new double[dataSize * 2]; + int nbrOfDoubles = MSNumpress.decodeLinear(data, dataSize, buffer); + if (nbrOfDoubles < 0) + throw new IllegalArgumentException("Corrupt numLin data!"); + + double[] result = new double[nbrOfDoubles]; + System.arraycopy(buffer, 0, result, 0, nbrOfDoubles); + return result; + + } else if (cvAccession.equals(ACC_NUMPRESS_SLOF)) { + double[] result = new double[(dataSize-8) / 2]; + MSNumpress.decodeSlof(data, dataSize, result); + return result; + + } else if (cvAccession.equals(ACC_NUMPRESS_PIC)) { + if (dataSize < 8 || data.length < 8) + throw new IllegalArgumentException("Cannot decode numPic data, need at least 8 initial bytes for fixed point."); + + double[] buffer = new double[dataSize * 2]; + int nbrOfDoubles = MSNumpress.decodePic(data, dataSize, buffer); + if (nbrOfDoubles < 0) + throw new IllegalArgumentException("Corrupt numPic data!"); + + double[] result = new double[nbrOfDoubles]; + System.arraycopy(buffer, 0, result, 0, nbrOfDoubles); + return result; + + } + + throw new IllegalArgumentException("'"+cvAccession+"' is not a numpress compression term"); + } + + + /** + * This encoding works on a 4 byte integer, by truncating initial zeros or ones. + * If the initial (most significant) half byte is 0x0 or 0xf, the number of such + * halfbytes starting from the most significant is stored in a halfbyte. This initial + * count is then followed by the rest of the ints halfbytes, in little-endian order. + * A count halfbyte c of + * + * 0 <= c <= 8 is interpreted as an initial c 0x0 halfbytes + * 9 <= c <= 15 is interpreted as an initial (c-8) 0xf halfbytes + * + * Ex: + * int c rest + * 0 => 0x8 + * -1 => 0xf 0xf + * 23 => 0x6 0x7 0x1 + * + * @x the int to be encoded + * @res the byte array were halfbytes are stored + * @resOffset position in res were halfbytes are written + * @return the number of resulting halfbytes + */ + protected static int encodeInt( + long x, + byte[] res, + int resOffset + ) { + byte i, l; + long m; + long mask = 0xf0000000; + long init = x & mask; + + if (init == 0) { + l = 8; + for (i=0; i<8; i++) { + m = mask >> (4*i); + if ((x & m) != 0) { + l = i; + break; + } + } + res[resOffset] = l; + for (i=l; i<8; i++) + res[resOffset+1+i-l] = (byte)(0xf & (x >> (4*(i-l)))); + + return 1+8-l; + + } else if (init == mask) { + l = 7; + for (i=0; i<8; i++) { + m = mask >> (4*i); + if ((x & m) != m) { + l = i; + break; + } + } + res[resOffset] = (byte)(l | 8); + for (i=l; i<8; i++) + res[resOffset+1+i-l] = (byte)(0xf & (x >> (4*(i-l)))); + + return 1+8-l; + + } else { + res[resOffset] = 0; + for (i=0; i<8; i++) + res[resOffset+1+i] = (byte)(0xf & (x >> (4*i))); + + return 9; + + } + } + + + + public static void encodeFixedPoint( + double fixedPoint, + byte[] result + ) { + long fp = Double.doubleToLongBits(fixedPoint); + for (int i=0; i<8; i++) { + result[7-i] = (byte)((fp >> (8*i)) & 0xff); + } + } + + + + public static double decodeFixedPoint( + byte[] data + ) { + long fp = 0; + for (int i=0; i<8; i++) { + fp = fp | ((0xFFl & data[7-i]) << (8*i)); + } + return Double.longBitsToDouble(fp); + } + + + + + ///////////////////////////////////////////////////////////////////////////////// + + /** + * Compute the maximal linear fixed point that prevents integer overflow. + * + * @data array of doubles to be encoded + * @dataSize number of doubles from data to encode + * + * @return the largest linear fixed point safe to use + */ + public static double optimalLinearFixedPoint( + double[] data, + int dataSize + ) { + if (dataSize == 0) return 0; + if (dataSize == 1) return Math.floor(0xFFFFFFFFl / data[0]); + double maxDouble = Math.max(data[0], data[1]); + + for (int i=2; i maxFpOverflow) return -1; + + return maxFp; + } + + + /** + * Encodes the doubles in data by first using a + * - lossy conversion to a 4 byte 5 decimal fixed point repressentation + * - storing the residuals from a linear prediction after first two values + * - encoding by encodeInt (see above) + * + * The resulting binary is maximally 8 + dataSize * 5 bytes, but much less if the + * data is reasonably smooth on the first order. + * + * This encoding is suitable for typical m/z or retention time binary arrays. + * On a test set, the encoding was empirically show to be accurate to at least 0.002 ppm. + * + * @data array of doubles to be encoded + * @dataSize number of doubles from data to encode + * @result array were resulting bytes should be stored + * @fixedPoint the scaling factor used for getting the fixed point repr. + * This is stored in the binary and automatically extracted + * on decoding. + * @return the number of encoded bytes + */ + public static int encodeLinear( + double[] data, + int dataSize, + byte[] result, + double fixedPoint + ) { + long[] ints = new long[3]; + int i; + int ri = 16; + byte halfBytes[] = new byte[10]; + int halfByteCount = 0; + int hbi; + long extrapol; + long diff; + + encodeFixedPoint(fixedPoint, result); + + if (dataSize == 0) return 8; + + ints[1] = (long)(data[0] * fixedPoint + 0.5); + for (i=0; i<4; i++) { + result[8+i] = (byte)((ints[1] >> (i*8)) & 0xff); + } + + if (dataSize == 1) return 12; + + ints[2] = (long)(data[1] * fixedPoint + 0.5); + for (i=0; i<4; i++) { + result[12+i] = (byte)((ints[2] >> (i*8)) & 0xff); + } + + halfByteCount = 0; + ri = 16; + + for (i=2; i 4294967294. + * The resulting binary is maximally dataSize * 5 bytes, but much less if the + * data is close to 0 on average. + * + * @data array of doubles to be encoded + * @dataSize number of doubles from data to encode + * @result array were resulting bytes should be stored + * @return the number of encoded bytes + */ + public static int encodePic( + double[] data, + int dataSize, + byte[] result + ) { + long count; + int ri = 0; + int hbi = 0; + byte halfBytes[] = new byte[10]; + int halfByteCount = 0; + + //printf("Encoding %d doubles\n", (int)dataSize); + + for (int i=0; i> 8); + } + return ri; + } + + + /** + * Decodes data encoded by encodeSlof + * + * The result vector will be exactly (|data| - 8) / 2 doubles. + * returns the number of doubles read, or -1 is there is a problem decoding. + * + * @data array of bytes to be decoded (need memorycont. repr.) + * @dataSize number of bytes from data to decode + * @result array were resulting doubles should be stored + * @return the number of decoded doubles + */ + public static int decodeSlof( + byte[] data, + int dataSize, + double[] result + ) { + int x; + int ri = 0; + + if (dataSize < 8) return -1; + double fixedPoint = decodeFixedPoint(data); + + if (dataSize % 2 != 0) return -1; + + for (int i=8; i>> data = [100, 101, 102, 103] + >>> encoded = []; decoded = [] + >>> PyMSNumpress.encodeLinear(data, encoded, 500.0) + >>> encoded + [64, 127, 64, 0, 0, 0, 0, 0, 80, 195, 0, 0, 68, 197, 0, 0, 136] + >>> PyMSNumpress.decodeLinear(encoded, decoded) + >>> decoded + [100.0, 101.0, 102.0, 103.0] + +The interface expects Python lists of ordinal numbers, these can be converted +to byte strings with "ord" and "chr" if desired: + + >>> bstr = "".join([chr(e) for e in encoded]) + >>> blist = [ord(b) for b in bstr] + + PyMSNumpress.pyx + roest@imsb.biol.ethz.ch + + Copyright 2013 Hannes Roest + + Licensed under the Apache License, Version 2.0 (the "License"); + you may not use this file except in compliance with the License. + You may obtain a copy of the License at + + http://www.apache.org/licenses/LICENSE-2.0 + + Unless required by applicable law or agreed to in writing, software + distributed under the License is distributed on an "AS IS" BASIS, + WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + See the License for the specific language governing permissions and + limitations under the License. +""" + + +from libcpp.vector cimport vector as libcpp_vector +from cython.operator cimport dereference as deref, preincrement as inc, address as address +from MSNumpress cimport encodeLinear as _encodeLinear +from MSNumpress cimport decodeLinear as _decodeLinear +from MSNumpress cimport optimalLinearFixedPoint as _optimalLinearFixedPoint +from MSNumpress cimport optimalLinearFixedPointMass as _optimalLinearFixedPointMass +from MSNumpress cimport encodeSlof as _encodeSlof +from MSNumpress cimport decodeSlof as _decodeSlof +from MSNumpress cimport optimalSlofFixedPoint as _optimalSlofFixedPoint +from MSNumpress cimport encodePic as _encodePic +from MSNumpress cimport decodePic as _decodePic + +def optimalLinearFixedPointMass(data, mz): + """ + + Compute the optimal linear fixed point with a desired m/z accuracy. + + @note If the desired accuracy cannot be reached without overflowing 64 + bit integers, then a negative value is returned. You need to check for + this and in that case abandon numpress or use optimalLinearFixedPoint + which returns the largest safe value. + + @data pointer to array of double to be encoded (need memorycont. repr.) + @mass_acc desired m/z accuracy in Th + @return the linear fixed point that satisfies the accuracy requirement (or -1 in case of failure). + + """ + dataSize = len(data) + cdef libcpp_vector[double] c_data = data + + cdef double result = _optimalLinearFixedPointMass( &c_data[0], dataSize, mz) + + return result + +def optimalLinearFixedPoint(data): + """ + + Compute the maximal linear fixed point that prevents integer overflow. + + @data pointer to array of double to be encoded (need memorycont. repr.) + @return the linear fixed point safe to use + + """ + dataSize = len(data) + cdef libcpp_vector[double] c_data = data + + cdef double result = _optimalLinearFixedPoint( &c_data[0], dataSize) + + return result + +def optimalSlofFixedPoint(data): + dataSize = len(data) + cdef libcpp_vector[double] c_data = data + + cdef double result = _optimalSlofFixedPoint( &c_data[0], dataSize) + + return result + +def decodeLinear(data, result): + """ + + Decodes data encoded by encodeLinear. + + result vector guaranteed to be shorter or equal to (|data| - 8) * 2 + + Note that this method may throw a const char* if it deems the input data to be corrupt, i.e. + that the last encoded int does not use the last byte in the data. In addition the last encoded + int need to use either the last halfbyte, or the second last followed by a 0x0 halfbyte. + + @data pointer to array of bytes to be decoded (need memorycont. repr.) + @result pointer to were resulting doubles should be stored + @return the number of decoded doubles, or -1 if dataSize < 4 or 4 < dataSize < 8 + + """ + cdef libcpp_vector[unsigned char] c_data = data + cdef libcpp_vector[double] c_result + + _decodeLinear(c_data, c_result) + + cdef libcpp_vector[double].iterator it_result = c_result.begin() + while it_result != c_result.end(): + result.append( deref(it_result) ) + inc(it_result) + +def encodeLinear(data, result, fixedPoint): + """ + + Encodes the doubles in data by first using a + - lossy conversion to a 4 byte 5 decimal fixed point representation + - storing the residuals from a linear prediction after first two values + - encoding by encodeInt (see above) + + The resulting binary is maximally 8 + dataSize * 5 bytes, but much less if the + data is reasonably smooth on the first order. + + This encoding is suitable for typical m/z or retention time binary arrays. + On a test set, the encoding was empirically show to be accurate to at least 0.002 ppm. + + @data pointer to array of double to be encoded (need memorycont. repr.) + @result pointer to where resulting bytes should be stored + @fixedPoint the scaling factor used for getting the fixed point repr. + This is stored in the binary and automatically extracted + on decoding (see optimalLinearFixedPoint or optimalLinearFixedPointMass) + @return the number of encoded bytes + + """ + cdef double c_fixedPoint = fixedPoint + cdef libcpp_vector[double] c_data = data + cdef libcpp_vector[unsigned char] c_result + + _encodeLinear(c_data, c_result, c_fixedPoint) + + cdef libcpp_vector[unsigned char].iterator it_result = c_result.begin() + while it_result != c_result.end(): + result.append( deref(it_result) ) + inc(it_result) + +def decodeSlof(data, result): + """ + + Decodes data encoded by encodeSlof + + The return will include exactly (|data| - 8) / 2 doubles. + + Note that this method may throw a const char* if it deems the input data to be corrupt. + + @data pointer to array of bytes to be decoded (need memorycont. repr.) + @result pointer to were resulting doubles should be stored + @return the number of decoded doubles + """ + cdef libcpp_vector[unsigned char] c_data = data + cdef libcpp_vector[double] c_result + + _decodeSlof(c_data, c_result) + + cdef libcpp_vector[double].iterator it_result = c_result.begin() + while it_result != c_result.end(): + result.append( deref(it_result) ) + inc(it_result) + +def encodeSlof(data, result, fixedPoint): + """ + + Encodes ion counts by taking the natural logarithm, and storing a + fixed point representation of this. This is calculated as + + unsigned short fp = log(d + 1) * fixedPoint + 0.5 + + the result vector is exactly |data| * 2 + 8 bytes long + + @data pointer to array of double to be encoded (need memorycont. repr.) + @result pointer to were resulting bytes should be stored + @fixedPoint fixed point to use for encoding (see optimalSlofFixedPoint) + @return the number of encoded bytes + """ + cdef double c_fixedPoint = fixedPoint + cdef libcpp_vector[double] c_data = data + cdef libcpp_vector[unsigned char] c_result + + _encodeSlof(c_data, c_result, c_fixedPoint) + + cdef libcpp_vector[unsigned char].iterator it_result = c_result.begin() + while it_result != c_result.end(): + result.append( deref(it_result) ) + inc(it_result) + +def decodePic(data, result): + """ + + Decodes data encoded by encodePic + + result vector guaranteed to be shorter of equal to |data| * 2 + + Note that this method may throw a const char* if it deems the input data to be corrupt, i.e. + that the last encoded int does not use the last byte in the data. In addition the last encoded + int need to use either the last halfbyte, or the second last followed by a 0x0 halfbyte. + + @data pointer to array of bytes to be decoded (need memorycont. repr.) + @result pointer to were resulting doubles should be stored + @return the number of decoded doubles + """ + cdef libcpp_vector[unsigned char] c_data = data + cdef libcpp_vector[double] c_result + + _decodePic(c_data, c_result) + + cdef libcpp_vector[double].iterator it_result = c_result.begin() + while it_result != c_result.end(): + result.append( deref(it_result) ) + inc(it_result) + +def encodePic(data, result): + """ + + Encodes ion counts by simply rounding to the nearest 4 byte integer, + and compressing each integer with encodeInt. + + The handleable range is therefore 0 -> 4294967294. + The resulting binary is maximally dataSize * 5 bytes, but much less if the + data is close to 0 on average. + + @data pointer to array of double to be encoded (need memorycont. repr.) + @result pointer to were resulting bytes should be stored + @return the number of encoded bytes + + """ + cdef libcpp_vector[double] c_data = data + cdef libcpp_vector[unsigned char] c_result + + _encodePic(c_data, c_result) + + cdef libcpp_vector[unsigned char].iterator it_result = c_result.begin() + while it_result != c_result.end(): + result.append( deref(it_result) ) + inc(it_result) + diff --git a/src/main/python/setup.py b/src/main/python/setup.py new file mode 100644 index 0000000..a8de752 --- /dev/null +++ b/src/main/python/setup.py @@ -0,0 +1,59 @@ +""" + setup.py + roest@imsb.biol.ethz.ch + + Copyright 2013 Hannes Roest + + Licensed under the Apache License, Version 2.0 (the "License"); + you may not use this file except in compliance with the License. + You may obtain a copy of the License at + + http://www.apache.org/licenses/LICENSE-2.0 + + Unless required by applicable law or agreed to in writing, software + distributed under the License is distributed on an "AS IS" BASIS, + WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + See the License for the specific language governing permissions and + limitations under the License. +""" + +""" +These Python bindings rely on the original C++ code for the actual function +calls. To compile the bindings, you will need Cython and the Python headers +installed on your system. + +To build and test, run + +$ python setup.py build_ext --inplace +$ nosetests test_pymsnumpress.py + +""" + +import os, shutil +from setuptools import setup, Extension +from Cython.Distutils import build_ext + +# copy C++ files +try: + numpress_file = os.path.join( "..", "cpp", "MSNumpress.cpp") + shutil.copy(numpress_file, ".") + numpress_file = os.path.join( "..", "cpp", "MSNumpress.hpp") + shutil.copy(numpress_file, ".") +except IOError: + pass + +ext_modules = [Extension("PyMSNumpress", + ["PyMSNumpress.pyx", "MSNumpress.cpp"], + language='c++', + )] + +setup( + ext_modules = ext_modules, + cmdclass = {'build_ext': build_ext}, + + name="PyMSNumpress", + + version="0.2.3" + +) + diff --git a/src/main/python/test_pymsnumpress.py b/src/main/python/test_pymsnumpress.py new file mode 100644 index 0000000..89bab64 --- /dev/null +++ b/src/main/python/test_pymsnumpress.py @@ -0,0 +1,144 @@ +#!/usr/bin/python +# -*- coding: utf-8 -*- +""" + setup.py + roest@imsb.biol.ethz.ch + + Copyright 2013 Hannes Roest + + Licensed under the Apache License, Version 2.0 (the "License"); + you may not use this file except in compliance with the License. + You may obtain a copy of the License at + + http://www.apache.org/licenses/LICENSE-2.0 + + Unless required by applicable law or agreed to in writing, software + distributed under the License is distributed on an "AS IS" BASIS, + WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + See the License for the specific language governing permissions and + limitations under the License. +""" + +import unittest +import subprocess as sub +import os + +import PyMSNumpress + +class TestMSNumpress(unittest.TestCase): + + def setUp(self): + self.data = [100, 101, 102, 103] + self.data_slof = [ + 100.0, + 200.0, + 300.00005, + 400.00010, + ] + self.data_long = [ + 100.0, + 200.0, + 300.00005, + 400.00010, + 450.00010, + 455.00010, + 700.00010 + ] + self.fp_slof = 10000 + self.linear_result = [64, 248, 106, 0, 0, 0, 0, 0, 128, 150, 152, 0, 32, 29, 154, 0, 136] + + def test_encodeLinear(self): + + result = [] + PyMSNumpress.encodeLinear(self.data,result,100000.0) + + self.assertEqual(len(result), 17) + self.assertEqual(result[0], 64) + self.assertEqual(result, self.linear_result) + + result = [] + PyMSNumpress.encodeLinear(self.data_long,result,5.0) + self.assertEqual(len(result), 22) + + result = [] + PyMSNumpress.encodeLinear(self.data_long,result,500.0) + self.assertEqual(len(result), 25) + + result = [] + PyMSNumpress.encodeLinear(self.data_long,result,5e4) + self.assertEqual(len(result), 29) + + result = [] + PyMSNumpress.encodeLinear(self.data_long,result,5e5) + self.assertEqual(len(result), 30) + + result = [] + PyMSNumpress.encodeLinear(self.data_long,result,5e6) + self.assertEqual(len(result), 31) + + # accurate to 3 sign digits + result = [] + decoded = [] + PyMSNumpress.encodeLinear(self.data_long,result,500.0) + PyMSNumpress.decodeLinear(result, decoded) + + self.assertAlmostEqual(decoded[0], 100, 3) + self.assertAlmostEqual(decoded[1], 200, 3) + self.assertAlmostEqual(decoded[2], 300, 3) + self.assertAlmostEqual(decoded[3], 400.00010, 3) + self.assertAlmostEqual(decoded[4], 450.00010, 3) + self.assertAlmostEqual(decoded[5], 455.00010, 3) + self.assertAlmostEqual(decoded[6], 700.00010, 3) + + + def test_decodeLinear(self): + + result = [] + PyMSNumpress.decodeLinear(self.linear_result, result) + + self.assertEqual(len(result), 4) + self.assertAlmostEqual(result[0], 100) + self.assertAlmostEqual(result, self.data) + + def test_encodePic(self): + + result = [] + encoded = [] + PyMSNumpress.encodePic(self.data,encoded) + self.assertEqual(len(encoded), 6) + PyMSNumpress.decodePic(encoded, result) + + self.assertEqual(len(result), 4) + self.assertAlmostEqual(result[0], 100) + self.assertAlmostEqual(result, self.data) + + def test_optimalSlofFixedPoint(self): + pt = PyMSNumpress.optimalSlofFixedPoint(self.data) + self.assertAlmostEqual(pt, 14110.0) + + def test_optimalLinearFixedPoint(self): + pt = PyMSNumpress.optimalLinearFixedPoint(self.data) + self.assertAlmostEqual(pt, 21262214.0) + + def test_optimalLinearFixedPointMass(self): + pt = PyMSNumpress.optimalLinearFixedPointMass(self.data, 0.001) + self.assertAlmostEqual(pt, 500.0) + pt = PyMSNumpress.optimalLinearFixedPointMass(self.data, 1e-10) + self.assertAlmostEqual(pt, -1) + + def test_encodeSlof(self): + + result = [] + encoded = [] + PyMSNumpress.encodeSlof(self.data_slof, encoded, self.fp_slof) + self.assertEqual(len(encoded), 16) + PyMSNumpress.decodeSlof(encoded, result) + + self.assertTrue( abs(result[0] - 100) < 1) + self.assertTrue( abs(result[1] - 200) < 1) + self.assertTrue( abs(result[2] - 300) < 1) + self.assertTrue( abs(result[3] - 400) < 1) + + +if __name__ == '__main__': + unittest.main() diff --git a/src/test/java/ms/numpress/MSNumpressTest.java b/src/test/java/ms/numpress/MSNumpressTest.java new file mode 100644 index 0000000..45459fc --- /dev/null +++ b/src/test/java/ms/numpress/MSNumpressTest.java @@ -0,0 +1,375 @@ +/* + MSNumpressTest.java + johan.teleman@immun.lth.se + + Copyright 2013 Johan Teleman + + Licensed under the Apache License, Version 2.0 (the "License"); + you may not use this file except in compliance with the License. + You may obtain a copy of the License at + + http://www.apache.org/licenses/LICENSE-2.0 + + Unless required by applicable law or agreed to in writing, software + distributed under the License is distributed on an "AS IS" BASIS, + WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + See the License for the specific language governing permissions and + limitations under the License. + */ +package ms.numpress; + +import static org.junit.Assert.*; + +import org.junit.Test; + +public class MSNumpressTest { + + @Test + public void encodeInt() { + byte[] res = new byte[10]; + int l; + l = MSNumpress.encodeInt(0l, res, 0); + assertEquals(1, l); + assertEquals(8, res[0] & 0xf); + + l = MSNumpress.encodeInt(-1l, res, 0); + assertEquals(2, l); + assertEquals(0xf, res[0] & 0xf); + assertEquals(0xf, res[1] & 0xf); + + l = MSNumpress.encodeInt(35l, res, 0); + assertEquals(3, l); + assertEquals(6, res[0] & 0xf); + assertEquals(0x3, res[1] & 0xf); + assertEquals(0x2, res[2] & 0xf); + + l = MSNumpress.encodeInt(370000000l, res, 0); + assertEquals(9, l); + assertEquals(0, res[0] & 0xf); + assertEquals(0x0, res[1] & 0xf); + assertEquals(0x8, res[2] & 0xf); + assertEquals(0x0, res[3] & 0xf); + assertEquals(0xc, res[4] & 0xf); + assertEquals(0xd, res[5] & 0xf); + assertEquals(0x0, res[6] & 0xf); + assertEquals(0x6, res[7] & 0xf); + assertEquals(0x1, res[8] & 0xf); + } + + + @Test + public void decodeInt() { + byte[] res = new byte[10]; + res[0] = (byte)0x75; + res[1] = (byte)0x87; + res[2] = (byte)0x10; + res[3] = (byte)0x08; + res[4] = (byte)0x0c; + res[5] = (byte)0xd0; + res[6] = (byte)0x61; + IntDecoder dec = new IntDecoder(res, 0); + + long l; + l = dec.next(); + assertEquals(5, l); + + l = dec.next(); + assertEquals(0, l); + + l = dec.next(); + assertEquals(1, l); + + l = dec.next(); + assertEquals(370000000L, l); + } + + + + @Test + public void encodeFixedPoint() { + byte[] encoded = new byte[8]; + MSNumpress.encodeFixedPoint(1.00, encoded); + assertEquals(0x3f, 0xff & encoded[0]); + assertEquals(0xf0, 0xff & encoded[1]); + assertEquals(0x0, 0xff & encoded[2]); + assertEquals(0x0, 0xff & encoded[3]); + assertEquals(0x0, 0xff & encoded[4]); + assertEquals(0x0, 0xff & encoded[5]); + assertEquals(0x0, 0xff & encoded[6]); + assertEquals(0x0, 0xff & encoded[7]); + } + + + @Test + public void encodeDecodeFixedPoint() { + double fp = 300.21941382293625; + byte[] encoded = new byte[8]; + MSNumpress.encodeFixedPoint(fp, encoded); + double decoded = MSNumpress.decodeFixedPoint(encoded); + assertEquals(fp, decoded, 0); + } + + + @Test + public void computeLinearFixedPointAtPPM() { + double[] mzs = {100.0, 200.0, 300.00005, 400.00010, 450.00010, 455.00010, 700.00010}; + + double fp; + fp = MSNumpress.optimalLinearFixedPointMass( mzs, mzs.length, 0.1 ); + assertEquals( 5, fp, 0.000005 ); + + fp = MSNumpress.optimalLinearFixedPointMass( mzs, mzs.length, 1e-3 ); + assertEquals( 500, fp, 0.000005 ); + + fp = MSNumpress.optimalLinearFixedPointMass( mzs, mzs.length, 1e-5 ); + assertEquals( 50000, fp, 0.000005 ); + + fp = MSNumpress.optimalLinearFixedPointMass( mzs, mzs.length, 1e-7 ); + assertEquals( 5000000, fp, 0.000005 ); + + // cannot fulfill accuracy of 1e-8 + fp = MSNumpress.optimalLinearFixedPointMass( mzs, mzs.length, 1e-8 ); + assertEquals( -1.0, fp, 0 ); + } + + + @Test + public void encodeLinear() { + double[] mzs = {100.0, 200.0, 300.00005, 400.00010}; + byte[] encoded = new byte[40]; + int encodedBytes = MSNumpress.encodeLinear(mzs, 4, encoded, 100000.0); + assertEquals(18, encodedBytes); + assertEquals(0x80, 0xff & encoded[8]); + assertEquals(0x96, 0xff & encoded[9]); + assertEquals(0x98, 0xff & encoded[10]); + assertEquals(0x00, 0xff & encoded[11]); + assertEquals(0x75, 0xff & encoded[16]); + assertEquals(0x80, 0xf0 & encoded[17]); + } + + + @Test + public void encodeLinearWithFixedPointAtPPM() { + double[] mzs = {100.0, 200.0, 300.00005, 400.00010, 450.00010, 455.00010, 700.00010}; + byte[] buffer = new byte[100]; + double[] decodeBuffer = new double[100]; + + double[] mzErrs = {0.1, 1e-3, 1e-5, 1e-6, 1e-7}; + int[] expectedLengths = {22, 25, 29, 30, 31}; + + for ( int i = 0; i < mzErrs.length; i++ ) + { + double mzErr = mzErrs[i]; + int expectedLength = expectedLengths[i]; + + double fp = MSNumpress.optimalLinearFixedPointMass( mzs, mzs.length, mzErr ); + int encodedLength = MSNumpress.encodeLinear( mzs, mzs.length, buffer, fp ); + int decodedLength = MSNumpress.decodeLinear( buffer, encodedLength, decodeBuffer ); + + assertEquals( expectedLength, encodedLength ); + assertEquals( mzs.length, decodedLength ); + + for ( int j = 0; j < mzs.length; j++ ) + { + assertEquals( mzs[j], decodeBuffer[j], mzErr ); + } + } + } + + + @Test + public void encodeDecodeLinearEmpty() { + byte[] encoded = new byte[8]; + int encodedBytes = MSNumpress.encodeLinear(new double[0], 0, encoded, 100000.0); + double[] decoded = new double[0]; + int decodedBytes = MSNumpress.decodeLinear(encoded, 8, decoded); + assertEquals(0, decodedBytes); + } + + + @Test + public void decodeLinearNice() { + double[] mzs = {100.0, 200.0, 300.00005, 400.00010}; + byte[] encoded = new byte[28]; + int encodedBytes = MSNumpress.encodeLinear(mzs, 4, encoded, 100000.0); + double[] decoded = new double[4]; + int numDecoded = MSNumpress.decodeLinear(encoded, encodedBytes, decoded); + assertEquals(4, numDecoded); + assertEquals(100.0, decoded[0], 0.000005); + assertEquals(200.0, decoded[1], 0.000005); + assertEquals(300.00005, decoded[2], 0.000005); + assertEquals(400.00010, decoded[3], 0.000005); + } + + + @Test + public void decodeLinearWierd() { + double[] mzs = {100.0, 200.0, 4000.00005, 0.00010}; + byte[] encoded = new byte[28]; + double fixedPoint = MSNumpress.optimalLinearFixedPoint(mzs, 4); + int encodedBytes = MSNumpress.encodeLinear(mzs, 4, encoded, fixedPoint); + double[] decoded = new double[4]; + int numDecoded = MSNumpress.decodeLinear(encoded, encodedBytes, decoded); + assertEquals(4, numDecoded); + assertEquals(100.0, decoded[0], 0.000005); + assertEquals(200.0, decoded[1], 0.000005); + assertEquals(4000.00005, decoded[2], 0.000005); + assertEquals(0.00010, decoded[3], 0.000005); + } + + + @Test + public void encodeDecodeLinear() { + + int n = 1000; + double[] mzs = new double[n]; + mzs[0] = 300 + Math.random(); + for (int i=1; i Date: Fri, 31 Jul 2026 11:25:22 +0200 Subject: [PATCH 18/31] dependency: Add MS-Numpress --- .gitignore | 1 + meson.build | 5 +++++ subprojects/README.md | 19 +++++++++++++++++++ 3 files changed, 25 insertions(+) create mode 100644 subprojects/README.md diff --git a/.gitignore b/.gitignore index 4dbe29f..1fb35d9 100644 --- a/.gitignore +++ b/.gitignore @@ -1 +1,2 @@ /.cache/ +/subprojects/.wraplock diff --git a/meson.build b/meson.build index c6983e8..39eb679 100644 --- a/meson.build +++ b/meson.build @@ -98,6 +98,9 @@ install_headers([ 'include/mzpeak/util/types.h', ], subdir : 'mzpeak') +# Sub-projects: +numpress = subproject('msnumpress') + # Library and executable targets: include = include_directories('include') @@ -120,6 +123,8 @@ deps = [ dependency('parquet', # Part of Arrow above. version : arrow_ver_str, required : true), + + numpress.get_variable('msnumpres_cpp_lib_so'), ] libmzpeak_so = library( diff --git a/subprojects/README.md b/subprojects/README.md new file mode 100644 index 0000000..bc41dc7 --- /dev/null +++ b/subprojects/README.md @@ -0,0 +1,19 @@ +# Sub-projects + +Dependencies that are not packaged separately. + +Every sub-project in this directory is a Git subtree and can be +updated to match upstream. + +## MS-Numpress + + - Original Repository: https://github.com/ms-numpress/ms-numpress + - Peter's Maintained Fork: https://github.com/pjones/ms-numpress + - License: http://www.apache.org/licenses/LICENSE-2.0 + +To update: + +``` +git subtree pull --squash --prefix=subprojects/msnumpress \ + https://github.com/pjones/ms-numpress.git pjones +``` From 400a8131cd41730f6320449d8110f7f7b26d8a47 Mon Sep 17 00:00:00 2001 From: "Peter J. Jones" Date: Fri, 31 Jul 2026 13:43:29 +0200 Subject: [PATCH 19/31] numpress: Implement a wrapper around MSNumpress --- include/mzpeak/util/numpress.h | 54 ++++++++++++++++++++++++++++++++++ meson.build | 2 ++ src/util/numpress.cpp | 48 ++++++++++++++++++++++++++++++ 3 files changed, 104 insertions(+) create mode 100644 include/mzpeak/util/numpress.h create mode 100644 src/util/numpress.cpp diff --git a/include/mzpeak/util/numpress.h b/include/mzpeak/util/numpress.h new file mode 100644 index 0000000..3845928 --- /dev/null +++ b/include/mzpeak/util/numpress.h @@ -0,0 +1,54 @@ +/* + +This file is part of the mzpeak project. It is subject to the license +specified in the LICENSE file which can be found in the top-level +directory of this repository. + +*/ + +#pragma once + +#include +#include +#include +#include + +namespace MzPeak::Util::Numpress { + +/** + * Decode a vector of bytes into a vector of doubles. + * + * The bytes need to be encoded using the MS-Numpress Linear encoding. + */ +void decode_linear(const std::vector&, std::vector&); + +/** + * Decode an arrow array of `uint8_t` values. + */ +std::shared_ptr> +decode_linear(const std::shared_ptr&); + +/** + * Decode and perform type conversion if necessary. + */ +template +std::shared_ptr> +decode_linear_convert(const std::shared_ptr& src) +{ + std::shared_ptr> doubles = decode_linear(src); + + if constexpr (std::is_same_v) { + return doubles; + } else { + std::shared_ptr> result = std::make_shared>(); + result->reserve(doubles->size()); + + for (const auto& d : *doubles) { + result->push_back(static_cast(d)); + } + + return result; + } +}; + +} // namespace MzPeak::Util::Numpress diff --git a/meson.build b/meson.build index 39eb679..9ebdb05 100644 --- a/meson.build +++ b/meson.build @@ -46,6 +46,7 @@ lib_sources = [ 'src/util/decoders.cpp', 'src/util/executor.cpp', 'src/util/manager.cpp', + 'src/util/numpress.cpp', 'src/util/parquet.cpp', 'src/util/planner.cpp', 'src/util/projection.cpp', @@ -90,6 +91,7 @@ install_headers([ 'include/mzpeak/util/enumerable_proxy.h', 'include/mzpeak/util/executor.h', 'include/mzpeak/util/manager.h', + 'include/mzpeak/util/numpress.h', 'include/mzpeak/util/parquet.h', 'include/mzpeak/util/planner.h', 'include/mzpeak/util/projection.h', diff --git a/src/util/numpress.cpp b/src/util/numpress.cpp new file mode 100644 index 0000000..56561f4 --- /dev/null +++ b/src/util/numpress.cpp @@ -0,0 +1,48 @@ +/* + +This file is part of the mzpeak project. It is subject to the license +specified in the LICENSE file which can be found in the top-level +directory of this repository. + +*/ + +#include + +#include "mzpeak/exception.h" +#include "mzpeak/util/decoders.h" +#include "mzpeak/util/numpress.h" + +namespace MzPeak::Util::Numpress { + +/******************************************************************************/ +void decode_linear(const std::vector& input, std::vector& output) +{ + try { + ms::numpress::MSNumpress::decodeLinear(input, output); + } catch (const char* msg) { + throw InvalidFormatError(msg); + } +} + +/******************************************************************************/ +std::shared_ptr> +decode_linear(const std::shared_ptr& src) +{ + if (src->type_id() != arrow::Type::UINT8) { + std::string msg("numpress decoding requested but source array is not uint8"); + throw InvalidFormatError(msg); + } + + std::vector bytes; + bytes.reserve(src->length()); + + Decoders::Scalar decoder; + decoder.decode(src, bytes); + + auto values = std::make_shared>(); + decode_linear(bytes, *values); + + return values; +} + +} // namespace MzPeak::Util::Numpress From 4a1784eddde2508bd2686956a6dad02fb73cb5d4 Mon Sep 17 00:00:00 2001 From: "Peter J. Jones" Date: Fri, 31 Jul 2026 13:49:08 +0200 Subject: [PATCH 20/31] chunk: Implement primary axis chunk decoding --- include/mzpeak/data/encoding.h | 24 ++-- include/mzpeak/data/transformer/primary.h | 137 ++++++++++++++++++++++ meson.build | 1 + 3 files changed, 154 insertions(+), 8 deletions(-) create mode 100644 include/mzpeak/data/transformer/primary.h diff --git a/include/mzpeak/data/encoding.h b/include/mzpeak/data/encoding.h index 878428e..86680ee 100644 --- a/include/mzpeak/data/encoding.h +++ b/include/mzpeak/data/encoding.h @@ -15,6 +15,7 @@ top-level directory of this repository. #include "mzpeak/data/array_index.h" #include "mzpeak/data/null_marking.h" #include "mzpeak/data/signals.h" +#include "mzpeak/data/transformer/primary.h" #include "mzpeak/exception.h" #include "mzpeak/schema/psi/data_type.h" #include "mzpeak/util/slice.h" @@ -155,23 +156,30 @@ void Decoder::decode_with_nulls(const ArrayIndex::Dimension& dim, std::vector& v) const { const auto& primary_entry = dim.values_entry(); - - auto col = - signals_->array_index()->entry_column(*signals_->groups(), primary_entry); + auto col = signals_->column(primary_entry); if (!col.has_value()) { throw ParquetError("unable to decode dimension, not in schema: " + dim.name); } + auto go = [&](auto&& decoder) -> void { slice_->array(col.value(), v, decoder); }; + if (primary_entry.buffer_format == Schema::BufferFormat::Point) { auto decoder = Util::Decoders::Scalar, N>(null_decoder); - slice_->array(col.value(), v, decoder); + go(decoder); } else { - auto decoder = Util::Decoders::Flattened, N>(null_decoder); - slice_->array(col.value(), v, decoder); + if (dim.is_main_axis()) { + using Transformer = Transformer::Primary::Decoder; + Transformer transformer(signals_, slice_, dim); + auto decoder = Util::Decoders::Flattened, N, Transformer>( + null_decoder, std::move(transformer)); + go(decoder); + } else { + // FIXME: Apply necessary transformations on the decoded array. + auto decoder = Util::Decoders::Flattened, N>(null_decoder); + go(decoder); + } } - - // FIXME: Apply necessary transformations on the decoded array. } } // namespace MzPeak::Data::Encoding diff --git a/include/mzpeak/data/transformer/primary.h b/include/mzpeak/data/transformer/primary.h new file mode 100644 index 0000000..872a0ce --- /dev/null +++ b/include/mzpeak/data/transformer/primary.h @@ -0,0 +1,137 @@ +/* + +This file is part of the mzpeak project. It is subject to the license +specified in the LICENSE file which can be found in the top-level +directory of this repository. + +*/ + +#pragma once + +#include "mzpeak/data/array_index.h" +#include "mzpeak/data/signals.h" +#include "mzpeak/schema/psi/chunk_encoding.h" +#include "mzpeak/util/algorithm.h" +#include "mzpeak/util/numpress.h" +#include "mzpeak/util/slice.h" +#include "mzpeak/util/types.h" + +namespace MzPeak::Data::Transformer::Primary { + +/** + * Transform Arrow arrays that have been encoded with one of the chunk + * encoding schemes on the primary axis. + */ +template class Decoder final { +public: + /// The type of result this function object returns. + using result_type = std::variant, + std::pair>, + std::shared_ptr>>; + + /// Constructor. + Decoder(std::shared_ptr, + std::shared_ptr, + const ArrayIndex::Dimension&); + + /** + * Transform the given array. + * + * This array could come from the `chunk_values` column, or from the + * `chunk_transform` column. + */ + result_type operator()(int64_t, const std::shared_ptr&) const; + +private: + std::string dim_name_; + std::vector chunk_encoding_; + std::vector chunk_start_; +}; + +/******************************************************************************/ +template +Decoder::Decoder(std::shared_ptr signals, + std::shared_ptr slice, + const ArrayIndex::Dimension& dim) + : dim_name_(dim.name) + , chunk_encoding_() + , chunk_start_() +{ + auto decode = [&](Schema::BufferFormat format, + std::vector& dest) -> void { + std::optional column = signals->column(dim, format); + + if (!column.has_value()) { + std::string msg("while decoding " + dim.name); + msg += " a needed column with buffer format "; + msg += Schema::buffer_format_to_string(format); + msg += " was not found"; + throw InvalidFormatError(msg); + } + + slice->array(*column, dest, Util::Decoders::Scalar()); + }; + + // Decode the `chunk_encoding` column. + std::vector encodings; + decode(Schema::BufferFormat::ChunkEncoding, encodings); + chunk_encoding_.reserve(encodings.size()); + + for (const auto& s : encodings) { + std::optional cv = Schema::CV::from_string(s); + + if (!cv.has_value()) { + throw InvalidFormatError("invalid chunk encoding CV: " + std::string(s)); + } + + chunk_encoding_.emplace_back(*cv); + } + + // Decode the `chunk_start` column. + decode(Schema::BufferFormat::ChunkStart, chunk_start_); + + // Sanity checks: + if (chunk_encoding_.size() != chunk_start_.size()) { + std::string msg("while decoding " + dim.name); + msg += " the chunk encoding and chunk start columns"; + msg += " have different lengths"; + throw InvalidFormatError(msg); + } +} + +/******************************************************************************/ +template +Decoder::result_type +Decoder::operator()(int64_t index, const std::shared_ptr& src) const +{ + if (index < 0 || static_cast(index) >= chunk_encoding_.size()) { + std::string msg("while decoding " + dim_name_); + msg += " the chunk_values/chunk_transform column is out of bounds "; + msg += std::to_string(index) + " >= " + std::to_string(chunk_encoding_.size()); + throw InvalidFormatError(msg); + } + + std::optional type = + chunk_encoding_[index].type(); + + if (!type.has_value()) { + std::string msg("while decoding " + dim_name_); + msg += " unknown chunk encoding method: "; + msg += chunk_encoding_[index].to_cv().to_string(); + throw InvalidFormatError(msg); + } + + switch (type.value()) { + case Schema::PSI::ChunkEncoding::Type::NoCompression: + return std::make_pair(chunk_start_[index], src); + case Schema::PSI::ChunkEncoding::Type::Delta: + return Util::Algorithm::null_delta_decode>( + chunk_start_[index], src); + case Schema::PSI::ChunkEncoding::Type::NumpressLinear: + return Util::Numpress::decode_linear_convert(src); + } + + std::unreachable(); +} + +} // namespace MzPeak::Data::Transformer::Primary diff --git a/meson.build b/meson.build index 9ebdb05..d36da6e 100644 --- a/meson.build +++ b/meson.build @@ -62,6 +62,7 @@ install_headers([ 'include/mzpeak/data/encoding.h', 'include/mzpeak/data/null_marking.h', 'include/mzpeak/data/signals.h', + 'include/mzpeak/data/transformer/primary.h', 'include/mzpeak/exception.h', 'include/mzpeak/index.h', 'include/mzpeak/io/archive.h', From b40ebb86550921cb4d252c7f04eee7d0036963dd Mon Sep 17 00:00:00 2001 From: "Peter J. Jones" Date: Thu, 6 Aug 2026 17:23:50 +0200 Subject: [PATCH 21/31] numpress: Simple way to track numpress compression type --- include/mzpeak/schema/group.h | 7 +++++++ include/mzpeak/util/numpress.h | 7 +++++++ src/schema/group.cpp | 10 ++++++++++ 3 files changed, 24 insertions(+) diff --git a/include/mzpeak/schema/group.h b/include/mzpeak/schema/group.h index 9b4550e..4e041d7 100644 --- a/include/mzpeak/schema/group.h +++ b/include/mzpeak/schema/group.h @@ -15,6 +15,7 @@ top-level directory of this repository. #include "mzpeak/schema/cv.h" #include "mzpeak/schema/file.h" +#include "mzpeak/util/numpress.h" #include "mzpeak/util/types.h" // Forward declarations. @@ -119,6 +120,12 @@ class Group final { */ void type(Util::Type); + /** + * Return the numpress method type if the column name indicates + * this is a numpress compressed column of `uint8_t`. + */ + std::optional possibly_numpress() const; + private: friend class Group; diff --git a/include/mzpeak/util/numpress.h b/include/mzpeak/util/numpress.h index 3845928..bb0b214 100644 --- a/include/mzpeak/util/numpress.h +++ b/include/mzpeak/util/numpress.h @@ -15,6 +15,13 @@ directory of this repository. namespace MzPeak::Util::Numpress { +/** + * The kind of Numpress compression. + */ +enum Type { + Linear, +}; + /** * Decode a vector of bytes into a vector of doubles. * diff --git a/src/schema/group.cpp b/src/schema/group.cpp index 78e1702..8650ddd 100644 --- a/src/schema/group.cpp +++ b/src/schema/group.cpp @@ -122,6 +122,16 @@ const std::optional& Group::Field::type() const { return type_; } /******************************************************************************/ void Group::Field::type(Util::Type type) { type_ = type; } +/******************************************************************************/ +std::optional Group::Field::possibly_numpress() const +{ + if (schema_name_.contains("numpress_linear")) { + return Util::Numpress::Linear; + } else { + return {}; + } +} + /******************************************************************************/ Group::Group(const parquet::schema::GroupNode& node, const Schema::File& file) : name_("root") From 6f5906b40d17d21f1ce9153133d27b8b24caac64 Mon Sep 17 00:00:00 2001 From: "Peter J. Jones" Date: Thu, 6 Aug 2026 13:15:40 +0200 Subject: [PATCH 22/31] decoding: Better length guessing, more list type handling - Improve length guessing when numpress is used - Deal with all of the list types --- include/mzpeak/util/decoders.h | 129 +++++++++++++++++++-------------- include/mzpeak/util/slice.h | 18 +---- src/util/decoders.cpp | 49 +++++++++++-- 3 files changed, 118 insertions(+), 78 deletions(-) diff --git a/include/mzpeak/util/decoders.h b/include/mzpeak/util/decoders.h index 6cdae2a..34f4153 100644 --- a/include/mzpeak/util/decoders.h +++ b/include/mzpeak/util/decoders.h @@ -13,6 +13,7 @@ top-level directory of this repository. #include #include +#include "mzpeak/schema/group.h" #include "mzpeak/util/compat.h" // IWYU pragma: keep #include "mzpeak/util/types.h" @@ -44,9 +45,50 @@ concept from_arrow_array = /******************************************************************************/ /** - * Return `true` if the given Arrow array is a `ListArray`. + * If the given array is a "list of lists" then visit each element of + * the outer list. The given function is called on non-null elements + * and given the index to the list element. + * + * Returns the length of the outer list. */ -bool is_list_array(const std::shared_ptr&); +template int64_t visit(const std::shared_ptr& ary, F f) +{ + auto go = [&f](const std::shared_ptr& list) -> int64_t { + for (int64_t index : std::views::iota(0, list->length())) { + if (list->IsValid(index)) { + std::invoke(f, index, list->value_slice(index)); + } + } + + return list->length(); + }; + + auto type = ary->type_id(); + + if (type == arrow::Type::LIST) { + return go(std::static_pointer_cast(ary)); + } else if (type == arrow::Type::FIXED_SIZE_LIST) { + return go(std::static_pointer_cast(ary)); + } else if (type == arrow::Type::LARGE_LIST) { + return go(std::static_pointer_cast(ary)); + } else if (type == arrow::Type::LIST_VIEW) { + return go(std::static_pointer_cast(ary)); + } else if (type == arrow::Type::LARGE_LIST_VIEW) { + return go(std::static_pointer_cast(ary)); + } else { + std::string msg("expected an arrow list array but found: "); + msg += ary->type()->name(); + throw TypeError(msg); + } +} + +/******************************************************************************/ +/** + * Try to figure out how many bytes should be reserved to decode the + * given array. + */ +std::size_t guess_array_length(const Schema::Column&, + const std::shared_ptr&); /******************************************************************************/ /** @@ -178,24 +220,12 @@ class List final : Helper> { /// Decoding function. void decode(const std::shared_ptr& src, C& dst) { - if (!is_list_array(src)) { - std::string msg("expected an arrow list array but found: "); - msg += src->type()->name(); - throw TypeError(msg); - } - - std::shared_ptr casted = - std::static_pointer_cast(src); - - for (int64_t i : std::views::iota(0, casted->length())) { - if (!casted->IsNull(i)) { - std::shared_ptr values(casted->value_slice(i)); - value_type res; - res.reserve(values->length()); - scalar_decoder_.decode(values, res); - this->push(dst, res); - } - } + visit(src, [&](int64_t, const std::shared_ptr& values) { + value_type res; + res.reserve(values->length()); + scalar_decoder_.decode(values, res); + this->push(dst, res); + }); } private: @@ -207,8 +237,8 @@ class List final : Helper> { * An array transformer that returns its argument unchanged. */ struct IdentityTransform { - std::shared_ptr&& operator()(int64_t, - std::shared_ptr&& a) const + const std::shared_ptr& + operator()(int64_t, const std::shared_ptr& a) const { return a; } @@ -267,40 +297,27 @@ class Flattened final : Helper> { /// Decoding function. void decode(const std::shared_ptr& src, Container& dst) { - if (!is_list_array(src)) { - std::string msg("expected an arrow list array but found: "); - msg += src->type()->name(); - throw TypeError(msg); - } - - std::shared_ptr casted = - std::static_pointer_cast(src); - - for (int64_t i : std::views::iota(0, casted->length())) { - if (casted->IsValid(i)) { - transform_result_type values(transformer_(index_, casted->value_slice(i))); - - std::visit( - [&](auto&& v) -> void { - using U = std::decay_t; - using P = std::pair>; - - if constexpr (std::is_same_v>) { - scalar_decoder_.decode(v, dst); - } else if constexpr (std::is_same_v) { - dst.push_back(v.first); - scalar_decoder_.decode(v.second, dst); - } else if constexpr (std::is_same_v>) { - dst.insert(dst.end(), v->begin(), v->end()); - } else { - static_assert(false_type, "invalid transform result"); - } - }, - values); - } - } - - ++index_; + index_ += visit(src, [&](int64_t i, const std::shared_ptr& elm) { + transform_result_type values(transformer_(index_ + i, elm)); + + std::visit( + [&](auto&& v) -> void { + using U = std::decay_t; + using P = std::pair>; + + if constexpr (std::is_same_v>) { + scalar_decoder_.decode(v, dst); + } else if constexpr (std::is_same_v) { + dst.push_back(v.first); + scalar_decoder_.decode(v.second, dst); + } else if constexpr (std::is_same_v>) { + dst.insert(dst.end(), v->begin(), v->end()); + } else { + static_assert(false_type, "invalid transform result"); + } + }, + values); + }); } private: diff --git a/include/mzpeak/util/slice.h b/include/mzpeak/util/slice.h index 9d043ff..d29c1e9 100644 --- a/include/mzpeak/util/slice.h +++ b/include/mzpeak/util/slice.h @@ -147,23 +147,7 @@ void Slice::array(const Column& field, V& v, T& t) const std::size_t size{}; for (const auto& chunk : *chunks) { - if (Decoders::is_list_array(chunk)) { - // N.B. When decoding lists they often are accompanied by a - // starting value (i.e. chunk_start) that must also be - // accommodated. - // - // For example, both delta encoding and "basic" encoding - // result in an array that is one element longer than the one - // stored in Parquet. - // - // In the numpress case we'll over-allocate by a small amount. - // Perhaps this code needs to be made smarter or moved - // somewhere else so it has enough context to decide how many - // bytes to allocate. - size += std::static_pointer_cast(chunk)->length() + 1; - } else { - size += chunk->length(); - } + size += Decoders::guess_array_length(field, chunk); } v.reserve(v.size() + size); diff --git a/src/util/decoders.cpp b/src/util/decoders.cpp index c6b38df..40d3328 100644 --- a/src/util/decoders.cpp +++ b/src/util/decoders.cpp @@ -11,13 +11,52 @@ directory of this repository. namespace MzPeak::Util::Decoders { /******************************************************************************/ -bool is_list_array(const std::shared_ptr& ary) +std::size_t guess_array_length(const Schema::Column& column, + const std::shared_ptr& ary) { - auto t = ary->type_id(); + std::optional numpress = column.second->possibly_numpress(); - return t == arrow::Type::LIST || t == arrow::Type::FIXED_SIZE_LIST || - t == arrow::Type::LARGE_LIST || t == arrow::Type::LIST_VIEW || - t == arrow::Type::LARGE_LIST_VIEW; + auto count = + [&numpress](const std::shared_ptr& nums) -> std::size_t { + if (numpress.has_value()) { + switch (numpress.value()) { + case Numpress::Linear: + return (nums->length() - 8) * 2; + } + + std::unreachable(); + } else { + return static_cast(nums->length()); + } + }; + + auto for_list = [&](const std::shared_ptr& list) -> std::size_t { + std::size_t size = {}; + + for (int64_t index : std::views::iota(0, list->length())) { + if (list->IsValid(index)) { + size += count(list->value_slice(index)); + } + } + + return size; + }; + + auto type = ary->type_id(); + + if (type == arrow::Type::LIST) { + return for_list(std::static_pointer_cast(ary)); + } else if (type == arrow::Type::FIXED_SIZE_LIST) { + return for_list(std::static_pointer_cast(ary)); + } else if (type == arrow::Type::LARGE_LIST) { + return for_list(std::static_pointer_cast(ary)); + } else if (type == arrow::Type::LIST_VIEW) { + return for_list(std::static_pointer_cast(ary)); + } else if (type == arrow::Type::LARGE_LIST_VIEW) { + return for_list(std::static_pointer_cast(ary)); + } else { + return count(ary); + } } } // namespace MzPeak::Util::Decoders From e011dde2080e9edce9d1fd343987b7722a77091a Mon Sep 17 00:00:00 2001 From: "Peter J. Jones" Date: Tue, 11 Aug 2026 14:40:11 +0200 Subject: [PATCH 23/31] delta: Account for initial null values --- include/mzpeak/util/algorithm.h | 18 +++++++++++++++++- 1 file changed, 17 insertions(+), 1 deletion(-) diff --git a/include/mzpeak/util/algorithm.h b/include/mzpeak/util/algorithm.h index 256d210..e0d1f35 100644 --- a/include/mzpeak/util/algorithm.h +++ b/include/mzpeak/util/algorithm.h @@ -132,7 +132,23 @@ null_delta_decode(typename type_traits::value_type start, }; // N.B.: "The start point is *excluded* from the chunk-values array." - append({start}); + // + // If the first value in the list is NULL then the initial starting + // value is ignored and the delta starts at 0. However, if the + // first *two* values are NULL then we still treat the starting + // point as 0 but add the starting value as the first element in the + // results. + // + // This is due to delta encoding coming after null marking. + if (length > 0 && casted->IsNull(0)) { + if (length > 1 && casted->IsNull(1)) { + append({start}); + } + + last = {}; + } else { + append({start}); + } for (int64_t index : std::views::iota(0, length)) { if (casted->IsValid(index)) { From 3cdb3482238777da7faa45b80d8ba62369cf4c4f Mon Sep 17 00:00:00 2001 From: "Peter J. Jones" Date: Tue, 11 Aug 2026 14:41:14 +0200 Subject: [PATCH 24/31] numpress: Implement wrappers for the remaining algorithms --- include/mzpeak/util/numpress.h | 93 ++++++++++++++++++++++++++++++---- src/util/decoders.cpp | 8 +-- src/util/numpress.cpp | 80 ++++++++++++++++++++++++++--- 3 files changed, 157 insertions(+), 24 deletions(-) diff --git a/include/mzpeak/util/numpress.h b/include/mzpeak/util/numpress.h index bb0b214..3db3885 100644 --- a/include/mzpeak/util/numpress.h +++ b/include/mzpeak/util/numpress.h @@ -11,6 +11,7 @@ directory of this repository. #include #include #include +#include #include namespace MzPeak::Util::Numpress { @@ -20,8 +21,40 @@ namespace MzPeak::Util::Numpress { */ enum Type { Linear, + SLOF, + PIC, }; +/** + * Return the number of elements that should be reserved in order to + * decode `n` bytes encoding with `Type` `t`. + */ +std::size_t decoding_space_needed(std::size_t n, Type t); + +/** + * Possibly create a new vector and copy all of the elements from + * `doubles` with a static cast to `T`. + * + * If `T` is `double` then return the input vector unchanged. + */ +template +std::shared_ptr> +cast(const std::shared_ptr>& doubles) +{ + if constexpr (std::is_same_v) { + return doubles; + } else { + std::shared_ptr> result = std::make_shared>(); + result->reserve(doubles->size()); + + for (const auto& d : *doubles) { + result->push_back(static_cast(d)); + } + + return result; + } +} + /** * Decode a vector of bytes into a vector of doubles. * @@ -43,19 +76,57 @@ std::shared_ptr> decode_linear_convert(const std::shared_ptr& src) { std::shared_ptr> doubles = decode_linear(src); + return cast(doubles); +} - if constexpr (std::is_same_v) { - return doubles; - } else { - std::shared_ptr> result = std::make_shared>(); - result->reserve(doubles->size()); +/** + * Decode a vector of bytes into a vector of doubles. + * + * The bytes need to be encoded using the MS-Numpress short logged + * float compression encoding. + */ +void decode_slof(const std::vector&, std::vector&); - for (const auto& d : *doubles) { - result->push_back(static_cast(d)); - } +/** + * Decode an arrow array of `uint8_t` values. + */ +std::shared_ptr> +decode_slof(const std::shared_ptr&); - return result; - } -}; +/** + * Decode and perform type conversion if necessary. + */ +template +std::shared_ptr> +decode_slof_convert(const std::shared_ptr& src) +{ + std::shared_ptr> doubles = decode_slof(src); + return cast(doubles); +} + +/** + * Decode a vector of bytes into a vector of doubles. + * + * The bytes need to be encoding using MS-Numpress positive integer + * compression encoding. + */ +void decode_pic(const std::vector&, std::vector&); + +/** + * Decode an arrow array of `uint8_t` values. + */ +std::shared_ptr> +decode_pic(const std::shared_ptr&); + +/** + * Decode and perform type conversion if necessary. + */ +template +std::shared_ptr> +decode_pic_convert(const std::shared_ptr& src) +{ + std::shared_ptr> doubles = decode_pic(src); + return cast(doubles); +} } // namespace MzPeak::Util::Numpress diff --git a/src/util/decoders.cpp b/src/util/decoders.cpp index 40d3328..a262df7 100644 --- a/src/util/decoders.cpp +++ b/src/util/decoders.cpp @@ -7,6 +7,7 @@ directory of this repository. */ #include "mzpeak/util/decoders.h" +#include "mzpeak/util/numpress.h" namespace MzPeak::Util::Decoders { @@ -19,12 +20,7 @@ std::size_t guess_array_length(const Schema::Column& column, auto count = [&numpress](const std::shared_ptr& nums) -> std::size_t { if (numpress.has_value()) { - switch (numpress.value()) { - case Numpress::Linear: - return (nums->length() - 8) * 2; - } - - std::unreachable(); + return Numpress::decoding_space_needed(nums->length(), numpress.value()); } else { return static_cast(nums->length()); } diff --git a/src/util/numpress.cpp b/src/util/numpress.cpp index 56561f4..39e9b57 100644 --- a/src/util/numpress.cpp +++ b/src/util/numpress.cpp @@ -9,24 +9,39 @@ directory of this repository. #include #include "mzpeak/exception.h" +#include "mzpeak/util/compat.h" // IWYU pragma: keep #include "mzpeak/util/decoders.h" #include "mzpeak/util/numpress.h" namespace MzPeak::Util::Numpress { /******************************************************************************/ -void decode_linear(const std::vector& input, std::vector& output) +using decoder = + std::move_only_function&, std::vector&)>; + +/******************************************************************************/ +std::size_t decoding_space_needed(std::size_t n, Type t) { - try { - ms::numpress::MSNumpress::decodeLinear(input, output); - } catch (const char* msg) { - throw InvalidFormatError(msg); + switch (t) { + case Numpress::Linear: + // Need C++26 for saturating_sub :( + if (n <= 8) return 0; + return (n - 8) * 2; + case Numpress::SLOF: + // FIXME: Could this be a typo in the numpress lib? + // Need C++26 for saturating_sub :( + if (n <= 8) return 0; + return (n - 8) / 2; + case Numpress::PIC: + return n * 2; } + + std::unreachable(); } /******************************************************************************/ std::shared_ptr> -decode_linear(const std::shared_ptr& src) +from_arrow(const std::shared_ptr& src, decoder f) { if (src->type_id() != arrow::Type::UINT8) { std::string msg("numpress decoding requested but source array is not uint8"); @@ -40,9 +55,60 @@ decode_linear(const std::shared_ptr& src) decoder.decode(src, bytes); auto values = std::make_shared>(); - decode_linear(bytes, *values); + f(bytes, *values); return values; } +/******************************************************************************/ +void decode_linear(const std::vector& input, std::vector& output) +{ + try { + ms::numpress::MSNumpress::decodeLinear(input, output); + } catch (const char* msg) { + throw InvalidFormatError(msg); + } +} + +/******************************************************************************/ +std::shared_ptr> +decode_linear(const std::shared_ptr& src) +{ + return from_arrow(src, [](auto& i, auto& o) -> void { decode_linear(i, o); }); +} + +/******************************************************************************/ +void decode_slof(const std::vector& input, std::vector& output) +{ + try { + ms::numpress::MSNumpress::decodeSlof(input, output); + } catch (const char* msg) { + throw InvalidFormatError(msg); + } +} + +/******************************************************************************/ +std::shared_ptr> +decode_slof(const std::shared_ptr& src) +{ + return from_arrow(src, [](auto& i, auto& o) -> void { decode_slof(i, o); }); +} + +/******************************************************************************/ +void decode_pic(const std::vector& input, std::vector& output) +{ + try { + ms::numpress::MSNumpress::decodePic(input, output); + } catch (const char* msg) { + throw InvalidFormatError(msg); + } +} + +/******************************************************************************/ +std::shared_ptr> +decode_pic(const std::shared_ptr& src) +{ + return from_arrow(src, [](auto& i, auto& o) -> void { decode_pic(i, o); }); +} + } // namespace MzPeak::Util::Numpress From 6ce1b3cf3a96797de035c0b39ce720055f633663 Mon Sep 17 00:00:00 2001 From: "Peter J. Jones" Date: Tue, 11 Aug 2026 14:41:57 +0200 Subject: [PATCH 25/31] transform: Add additional transform types for numpress --- include/mzpeak/schema/psi/transform.h | 14 ++++++++++++++ src/schema/psi/transform.cpp | 13 +++++++++++++ 2 files changed, 27 insertions(+) diff --git a/include/mzpeak/schema/psi/transform.h b/include/mzpeak/schema/psi/transform.h index 2bdcb0f..bf3ac71 100644 --- a/include/mzpeak/schema/psi/transform.h +++ b/include/mzpeak/schema/psi/transform.h @@ -23,16 +23,30 @@ class Transform final { * Known transform types. */ enum Type { + /// MS:1003901 + /// /// Apply an algorithm to remove excess zero intensity value data /// points from a spectrum. Data may be retained for /// interperatbility such as retaining only zeros that flank /// non-zero intensity value data points from a profile spectrum. ZeroIntensityTrim, + /// MS:1003902 + /// /// A zero intensity point trimming algorithm that interpolates /// the m/z coordinate values from the local data or an estimated /// model. ZeroIntensityInterpolation, + + /// MS:1002314 + /// + /// Compression using MS-Numpress short logged float compression. + NumpressSLOF, + + /// MS:1002313 + /// + /// Compression using MS-Numpress positive integer compression. + NumpressPIC, }; /// The transform value can be a type as described by the Type enum, diff --git a/src/schema/psi/transform.cpp b/src/schema/psi/transform.cpp index cec0653..1850374 100644 --- a/src/schema/psi/transform.cpp +++ b/src/schema/psi/transform.cpp @@ -21,6 +21,10 @@ CV type_to_cv(Transform::Type t) return CV("MS", "1003901"); case Transform::ZeroIntensityInterpolation: return CV("MS", "1003902"); + case Transform::NumpressSLOF: + return CV("MS", "1002314"); + case Transform::NumpressPIC: + return CV("MS", "1002313"); } std::unreachable(); @@ -34,8 +38,13 @@ Transform::value_type to_value_type(const CV& cv) return Transform::ZeroIntensityTrim; } else if (cv.accession() == "1003902") { return Transform::ZeroIntensityInterpolation; + } else if (cv.accession() == "1002314") { + return Transform::NumpressSLOF; + } else if (cv.accession() == "1002313") { + return Transform::NumpressPIC; } } + return cv; } @@ -91,6 +100,10 @@ bool Transform::needs_delta_model() const noexcept return false; case ZeroIntensityInterpolation: return true; + case NumpressSLOF: + return false; + case NumpressPIC: + return false; } return {}; From 0f76cd2ea7efe1efacbbe35532db186db143c115 Mon Sep 17 00:00:00 2001 From: "Peter J. Jones" Date: Tue, 11 Aug 2026 14:42:40 +0200 Subject: [PATCH 26/31] chunk: Decode secondary dimensions --- include/mzpeak/data/encoding.h | 7 +- include/mzpeak/data/transformer/secondary.h | 80 +++++++++++++++++++++ meson.build | 1 + 3 files changed, 86 insertions(+), 2 deletions(-) create mode 100644 include/mzpeak/data/transformer/secondary.h diff --git a/include/mzpeak/data/encoding.h b/include/mzpeak/data/encoding.h index 86680ee..88024b1 100644 --- a/include/mzpeak/data/encoding.h +++ b/include/mzpeak/data/encoding.h @@ -16,6 +16,7 @@ top-level directory of this repository. #include "mzpeak/data/null_marking.h" #include "mzpeak/data/signals.h" #include "mzpeak/data/transformer/primary.h" +#include "mzpeak/data/transformer/secondary.h" #include "mzpeak/exception.h" #include "mzpeak/schema/psi/data_type.h" #include "mzpeak/util/slice.h" @@ -175,8 +176,10 @@ void Decoder::decode_with_nulls(const ArrayIndex::Dimension& dim, null_decoder, std::move(transformer)); go(decoder); } else { - // FIXME: Apply necessary transformations on the decoded array. - auto decoder = Util::Decoders::Flattened, N>(null_decoder); + using Transformer = Transformer::Secondary::Decoder; + Transformer transformer(dim); + auto decoder = Util::Decoders::Flattened, N, Transformer>( + null_decoder, std::move(transformer)); go(decoder); } } diff --git a/include/mzpeak/data/transformer/secondary.h b/include/mzpeak/data/transformer/secondary.h new file mode 100644 index 0000000..3d2d7d4 --- /dev/null +++ b/include/mzpeak/data/transformer/secondary.h @@ -0,0 +1,80 @@ +/* + +This file is part of the mzpeak project. It is subject to the license +specified in the LICENSE file which can be found in the top-level +directory of this repository. + +*/ + +#pragma once + +#include "mzpeak/data/array_index.h" +#include "mzpeak/util/types.h" + +namespace MzPeak::Data::Transformer::Secondary { + +/** + * Decoding class for secondary dimensions. + */ +template class Decoder final { +public: + /// The type of result this function object returns. + using result_type = std::variant, + std::pair>, + std::shared_ptr>>; + /// Constructor. + Decoder(const ArrayIndex::Dimension&); + + /** + * Transform the given array. + * + * This array could come from the `chunk_values` column, or from the + * `chunk_transform` column. + */ + result_type operator()(int64_t, const std::shared_ptr&) const; + +private: + std::string dim_name_; + std::optional transform_; +}; + +/******************************************************************************/ +template +Decoder::Decoder(const ArrayIndex::Dimension& dim) + : dim_name_(dim.name) + , transform_(dim.transform) +{ +} + +/******************************************************************************/ +template +Decoder::result_type +Decoder::operator()(int64_t, const std::shared_ptr& src) const +{ + if (!transform_.has_value()) return src; + + std::optional type = transform_.value().type(); + + if (!type.has_value()) { + std::string msg("while decoding " + dim_name_ + ": "); + msg += "unknown transform method: " + transform_.value().to_cv().to_string(); + throw InvalidFormatError(msg); + } + + switch (type.value()) { + case Schema::PSI::Transform::ZeroIntensityTrim: + case Schema::PSI::Transform::ZeroIntensityInterpolation: + // Handled by the null decoding code. + return src; + + case Schema::PSI::Transform::NumpressSLOF: + return Util::Numpress::decode_slof_convert(src); + + case Schema::PSI::Transform::NumpressPIC: + return Util::Numpress::decode_pic_convert(src); + } + + std::unreachable(); +} + +} // namespace MzPeak::Data::Transformer::Secondary diff --git a/meson.build b/meson.build index d36da6e..13d6d8b 100644 --- a/meson.build +++ b/meson.build @@ -63,6 +63,7 @@ install_headers([ 'include/mzpeak/data/null_marking.h', 'include/mzpeak/data/signals.h', 'include/mzpeak/data/transformer/primary.h', + 'include/mzpeak/data/transformer/secondary.h', 'include/mzpeak/exception.h', 'include/mzpeak/index.h', 'include/mzpeak/io/archive.h', From 301241be621f91876f06f6a4016c3a747206158b Mon Sep 17 00:00:00 2001 From: "Peter J. Jones" Date: Tue, 11 Aug 2026 14:43:07 +0200 Subject: [PATCH 27/31] test: Test decoding of all test files --- test/spectra_test.cpp | 67 +++++++++++++++++++++++++++---------------- 1 file changed, 43 insertions(+), 24 deletions(-) diff --git a/test/spectra_test.cpp b/test/spectra_test.cpp index bd3993c..410af51 100644 --- a/test/spectra_test.cpp +++ b/test/spectra_test.cpp @@ -20,35 +20,54 @@ BOOST_AUTO_TEST_CASE(can_read_spectra) { using namespace MzPeak; - auto mzpeak = MzPeak::open("../test/files/small.dir"); - auto spectra = mzpeak.spectra(); + auto go = [](const std::string file_name) { + auto mzpeak = MzPeak::open("../test/files/" + file_name); + auto spectra = mzpeak.spectra(); - BOOST_TEST((spectra.size() == 48)); + BOOST_TEST_CONTEXT("while using the " << file_name << "file") + { + BOOST_TEST((spectra.size() == 48)); - auto spectrum = spectra[0]; - auto mz = spectrum.mz(); + auto spectrum = spectra[0]; + auto mz = spectrum.mz(); - BOOST_TEST(mz.size() == 13589); - BOOST_TEST(mz[0] == 202.607, boost::test_tools::tolerance(0.001)); - BOOST_TEST(mz[mz.size() - 1] == 1999.840, boost::test_tools::tolerance(0.001)); + auto tolerance = boost::test_tools::tolerance(0.001); - // Test some NULL values. - BOOST_TEST(mz[7] == 202.608, boost::test_tools::tolerance(0.001)); - BOOST_TEST(mz[8] == 202.609, boost::test_tools::tolerance(0.001)); - BOOST_TEST(mz[14] == 204.761, boost::test_tools::tolerance(0.001)); - BOOST_TEST(mz[15] == 204.762, boost::test_tools::tolerance(0.001)); + if (file_name == "small.numpress.mzpeak") { + // Some numpress linear values are less precise than their + // matching point or delta values. + tolerance = boost::test_tools::tolerance(0.1); + } - // The m/z values should be monotonically increasing. - for (std::size_t i : std::views::iota(1ul, mz.size())) { - BOOST_TEST(mz[i] > mz[i - 1]); - } + BOOST_TEST(mz.size() == 13589); + BOOST_TEST(mz[0] == 202.607, tolerance); + BOOST_TEST(mz[mz.size() - 1] == 1999.840, tolerance); - auto intensity = spectrum.intensity(); - BOOST_TEST((intensity.size() == mz.size())); - BOOST_TEST(intensity[0] == 0.0, boost::test_tools::tolerance(0.001)); - BOOST_TEST(intensity[1] == 1938.12, boost::test_tools::tolerance(0.001)); - BOOST_TEST(intensity[intensity.size() - 1] == 0.0, - boost::test_tools::tolerance(0.001)); + // Test some NULL values. + BOOST_TEST_REQUIRE(mz[7] == 202.60831, tolerance); + BOOST_TEST_REQUIRE(mz[8] == 202.60856, tolerance); + BOOST_TEST_REQUIRE(mz[14] == 204.761, tolerance); + BOOST_TEST_REQUIRE(mz[15] == 204.762, tolerance); - BOOST_TEST(spectrum.ms_level() == 1u); + // The m/z values should be monotonically increasing. + for (std::size_t i : std::views::iota(1ul, mz.size())) { + BOOST_TEST_REQUIRE(mz[i] > mz[i - 1]); + } + + auto intensity = spectrum.intensity(); + BOOST_TEST_REQUIRE((intensity.size() == mz.size())); + BOOST_TEST_REQUIRE(intensity[0] == 0.0, tolerance); + BOOST_TEST_REQUIRE(intensity[1] == 1938.12, tolerance); + BOOST_TEST_REQUIRE(intensity[7] == 0.0, tolerance); + BOOST_TEST_REQUIRE(intensity[8] == 0.0, tolerance); + BOOST_TEST_REQUIRE(intensity[9] == 1422.17, tolerance); + BOOST_TEST_REQUIRE(intensity[intensity.size() - 1] == 0.0, tolerance); + + BOOST_TEST_REQUIRE(spectrum.ms_level() == 1u); + } + }; + + go("small.mzpeak"); + go("small.chunked.mzpeak"); + go("small.numpress.mzpeak"); } From d9ea73f990fdfc7e6d4faa622b268052a53d352e Mon Sep 17 00:00:00 2001 From: "Peter J. Jones" Date: Tue, 11 Aug 2026 14:56:27 +0200 Subject: [PATCH 28/31] mzp-inspect: Add the ability to dump spectra data --- bin/mzp-inspect.cpp | 22 ++++++++++++++++++++++ 1 file changed, 22 insertions(+) diff --git a/bin/mzp-inspect.cpp b/bin/mzp-inspect.cpp index 5b008c7..ed2dd76 100644 --- a/bin/mzp-inspect.cpp +++ b/bin/mzp-inspect.cpp @@ -146,6 +146,24 @@ int print_fmd_kv(MzPeak::Index& index, return 0; } +/******************************************************************************/ +int dump_spectra(MzPeak::Index& index) +{ + auto spectra = index.spectra(); + + for (std::size_t spectrum_index : std::views::iota(0ul, spectra.size())) { + const auto& spectrum = spectra[spectrum_index]; + const auto& mz = spectrum.mz(); + const auto& intensity = spectrum.intensity(); + + for (std::size_t row : std::views::iota(0ul, mz.size())) { + std::println("{},{:.5f},{:.5f}", spectrum_index, mz[row], intensity[row]); + } + } + + return 0; +} + /******************************************************************************/ int main(int argc, char* argv[]) { @@ -169,6 +187,8 @@ int main(int argc, char* argv[]) desc.add_options()("fmd-key", po::value(), "Used with --fmdkv to print the value of the given key"); + desc.add_options()("spectra", "Print all m/z and intensity values"); + po::positional_options_description pops; pops.add("file", 1); @@ -204,6 +224,8 @@ int main(int argc, char* argv[]) } return print_fmd_kv(index, vmap["fmdkv"].as(), key); + } else if (vmap.count("spectra")) { + dump_spectra(index); } else { std::println("WARN: no command given"); return 1; From 1bf145d2e18febd9eff09ba131925337d3b876ed Mon Sep 17 00:00:00 2001 From: "Peter J. Jones" Date: Tue, 11 Aug 2026 15:08:42 +0200 Subject: [PATCH 29/31] numpress: Gather all numpress-related code into the same translation unit --- include/mzpeak/util/numpress.h | 6 ++++++ src/schema/group.cpp | 6 +----- src/util/numpress.cpp | 14 ++++++++++++++ 3 files changed, 21 insertions(+), 5 deletions(-) diff --git a/include/mzpeak/util/numpress.h b/include/mzpeak/util/numpress.h index 3db3885..1778811 100644 --- a/include/mzpeak/util/numpress.h +++ b/include/mzpeak/util/numpress.h @@ -25,6 +25,12 @@ enum Type { PIC, }; +/** + * Attempt to infer the numpress method type from the name of a + * column. + */ +std::optional type_from_column_name(const std::string&); + /** * Return the number of elements that should be reserved in order to * decode `n` bytes encoding with `Type` `t`. diff --git a/src/schema/group.cpp b/src/schema/group.cpp index 8650ddd..63e1642 100644 --- a/src/schema/group.cpp +++ b/src/schema/group.cpp @@ -125,11 +125,7 @@ void Group::Field::type(Util::Type type) { type_ = type; } /******************************************************************************/ std::optional Group::Field::possibly_numpress() const { - if (schema_name_.contains("numpress_linear")) { - return Util::Numpress::Linear; - } else { - return {}; - } + return Util::Numpress::type_from_column_name(schema_name_); } /******************************************************************************/ diff --git a/src/util/numpress.cpp b/src/util/numpress.cpp index 39e9b57..f0369b8 100644 --- a/src/util/numpress.cpp +++ b/src/util/numpress.cpp @@ -19,6 +19,20 @@ namespace MzPeak::Util::Numpress { using decoder = std::move_only_function&, std::vector&)>; +/******************************************************************************/ +std::optional type_from_column_name(const std::string& name) +{ + if (name.contains("numpress_linear")) { + return Util::Numpress::Linear; + } else if (name.contains("numpress_slof")) { + return Util::Numpress::SLOF; + } else if (name.contains("numpress_pic")) { + return Util::Numpress::PIC; + } else { + return {}; + } +} + /******************************************************************************/ std::size_t decoding_space_needed(std::size_t n, Type t) { From 1caafc1f8f6a72101ee7c6527c721e01cf3574e6 Mon Sep 17 00:00:00 2001 From: "Peter J. Jones" Date: Tue, 11 Aug 2026 15:27:48 +0200 Subject: [PATCH 30/31] Ensure the correct headers are included --- include/mzpeak/data/transformer/primary.h | 2 ++ include/mzpeak/util/types.h | 1 + 2 files changed, 3 insertions(+) diff --git a/include/mzpeak/data/transformer/primary.h b/include/mzpeak/data/transformer/primary.h index 872a0ce..3061d40 100644 --- a/include/mzpeak/data/transformer/primary.h +++ b/include/mzpeak/data/transformer/primary.h @@ -8,6 +8,8 @@ directory of this repository. #pragma once +#include + #include "mzpeak/data/array_index.h" #include "mzpeak/data/signals.h" #include "mzpeak/schema/psi/chunk_encoding.h" diff --git a/include/mzpeak/util/types.h b/include/mzpeak/util/types.h index 47a643d..a3703be 100644 --- a/include/mzpeak/util/types.h +++ b/include/mzpeak/util/types.h @@ -8,6 +8,7 @@ directory of this repository. #pragma once +#include #include #include #include From 6f7230be75909339b006380508d615b081ffdffa Mon Sep 17 00:00:00 2001 From: "Peter J. Jones" Date: Tue, 11 Aug 2026 15:37:35 +0200 Subject: [PATCH 31/31] array_index: Add default values due to default constructor --- include/mzpeak/data/array_index.h | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/include/mzpeak/data/array_index.h b/include/mzpeak/data/array_index.h index ed0268c..74b3e6e 100644 --- a/include/mzpeak/data/array_index.h +++ b/include/mzpeak/data/array_index.h @@ -204,13 +204,13 @@ class ArrayIndex final { private: // The entity type for the entire Parquet file. - EntityType entity_type_; + EntityType entity_type_ = EntityType::Other; // Root node. std::string prefix_ = "point"; // Layout. - Layout layout_; + Layout layout_ = Layout::Unknown; // Entries; std::vector entries_;