diff --git a/.gitignore b/.gitignore index 4dbe29f..1fb35d9 100644 --- a/.gitignore +++ b/.gitignore @@ -1 +1,2 @@ /.cache/ +/subprojects/.wraplock 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; diff --git a/include/mzpeak/data/array_index.h b/include/mzpeak/data/array_index.h index a1ca2e6..74b3e6e 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. */ @@ -88,6 +105,14 @@ 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; + + /// Return `true` if this entry stores values for the associated + /// dimension. + bool is_value_entry() const; }; /** @@ -109,11 +134,23 @@ 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; + + /// Find the first entry with the given buffer format. + std::optional entry_with(BufferFormat) const; }; /// Default constructor. @@ -135,6 +172,11 @@ class ArrayIndex final { */ const std::string& prefix() const; + /** + * Return the file layout. + */ + Layout layout() const; + /** * Get a list of entry definitions. */ @@ -162,11 +204,14 @@ 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::Unknown; + // Entries; std::vector entries_; diff --git a/include/mzpeak/data/encoding.h b/include/mzpeak/data/encoding.h index bf914e8..88024b1 100644 --- a/include/mzpeak/data/encoding.h +++ b/include/mzpeak/data/encoding.h @@ -15,8 +15,9 @@ 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/data/transformer/secondary.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 +58,9 @@ 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 decode_with_nulls(const ArrayIndex::Dimension&, + const N& null_decoder, + std::vector&) const; template void remap(const ArrayIndex::Dimension& dim, std::vector& v) const; @@ -129,42 +132,57 @@ 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: + case ArrayIndex::Layout::Chunked: if (dim.needs_delta_model()) { using N = NullMarking::Decoder; - point(field.value(), N{delta_estimator_}, v); + decode_with_nulls(dim, N{delta_estimator_}, v); } else { using N = Util::Decoders::NullToZero; - point(field.value(), N{}, v); + decode_with_nulls(dim, N{}, v); } - } else { - 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, - 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 { - slice_->array(col, v, Util::Decoders::Scalar, N>(null_decoder)); + const auto& primary_entry = dim.values_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); + go(decoder); + } else { + 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 { + using Transformer = Transformer::Secondary::Decoder; + Transformer transformer(dim); + auto decoder = Util::Decoders::Flattened, N, Transformer>( + null_decoder, std::move(transformer)); + go(decoder); + } + } } } // namespace MzPeak::Data::Encoding 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/include/mzpeak/data/transformer/primary.h b/include/mzpeak/data/transformer/primary.h new file mode 100644 index 0000000..3061d40 --- /dev/null +++ b/include/mzpeak/data/transformer/primary.h @@ -0,0 +1,139 @@ +/* + +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/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/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/include/mzpeak/exception.h b/include/mzpeak/exception.h index 1ba4df9..f8b9d9f 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. */ @@ -87,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/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/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/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/include/mzpeak/util/algorithm.h b/include/mzpeak/util/algorithm.h index 99388ba..e0d1f35 100644 --- a/include/mzpeak/util/algorithm.h +++ b/include/mzpeak/util/algorithm.h @@ -8,12 +8,16 @@ top-level directory of this repository. #pragma once +#include +#include #include #include #include #include #include +#include "mzpeak/util/types.h" + namespace MzPeak::Util::Algorithm { /** @@ -78,4 +82,92 @@ 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." + // + // 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)) { + 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/decoders.h b/include/mzpeak/util/decoders.h index 368b2d4..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" @@ -33,10 +34,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 = @@ -46,6 +43,53 @@ concept from_arrow_array = { t.decode(a, r) } -> std::same_as; }; +/******************************************************************************/ +/** + * 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. + */ +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&); + /******************************************************************************/ /** * A NULL decoder that always skips NULL values. @@ -116,9 +160,6 @@ class Scalar final : Helper> { { } - /// Destructor. - ~Scalar() = default; - /// Decoding function. void decode(const std::shared_ptr& src, C& dst) { @@ -149,38 +190,140 @@ 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() {} - /// Destructor. - ~List() = default; + // Constructor where you can pass a null decoder to the scalar decoder. + List(const null_decoder_type& null_decoder) + : scalar_decoder_(null_decoder) + { + } /// Decoding function. void decode(const std::shared_ptr& src, C& dst) { - std::shared_ptr casted = - std::static_pointer_cast(src); + 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); + }); + } - 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().decode(values, res); - this->push(dst, res); - } - } +private: + Scalar scalar_decoder_; +}; + +/******************************************************************************/ +/** + * An array transformer that returns its argument unchanged. + */ +struct IdentityTransform { + const std::shared_ptr& + operator()(int64_t, const 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 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. + */ +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; + + /// The type of return value allowed from transformers. + using transform_result_type = + 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 = {}) + : scalar_decoder_(null_decoder) + , transformer_(transformer) + { } + + /// Decoding function. + void decode(const std::shared_ptr& src, Container& dst) + { + 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: + Scalar scalar_decoder_; + Transformer transformer_; + int64_t index_ = 0; }; } // namespace MzPeak::Util::Decoders diff --git a/include/mzpeak/util/numpress.h b/include/mzpeak/util/numpress.h new file mode 100644 index 0000000..1778811 --- /dev/null +++ b/include/mzpeak/util/numpress.h @@ -0,0 +1,138 @@ +/* + +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 +#include + +namespace MzPeak::Util::Numpress { + +/** + * The kind of Numpress compression. + */ +enum Type { + Linear, + SLOF, + 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`. + */ +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. + * + * 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); + return cast(doubles); +} + +/** + * 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&); + +/** + * Decode an arrow array of `uint8_t` values. + */ +std::shared_ptr> +decode_slof(const std::shared_ptr&); + +/** + * 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/include/mzpeak/util/slice.h b/include/mzpeak/util/slice.h index 5128dda..d29c1e9 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,11 +69,22 @@ 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 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; @@ -114,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; @@ -122,7 +147,7 @@ void Slice::array(const Column& field, V& v, T&& t) const std::size_t size{}; for (const auto& chunk : *chunks) { - size += chunk->length(); + size += Decoders::guess_array_length(field, chunk); } v.reserve(v.size() + size); diff --git a/include/mzpeak/util/types.h b/include/mzpeak/util/types.h index 72ebaf0..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 @@ -103,6 +104,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 +112,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 +120,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 +128,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 +136,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 +144,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 +152,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 +160,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 +168,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 62a13f7..13d6d8b 100644 --- a/meson.build +++ b/meson.build @@ -37,13 +37,16 @@ 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', 'src/spectrum.cpp', 'src/util/arrow.cpp', + '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', @@ -59,6 +62,8 @@ 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/data/transformer/secondary.h', 'include/mzpeak/exception.h', 'include/mzpeak/index.h', 'include/mzpeak/io/archive.h', @@ -75,6 +80,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', @@ -87,6 +93,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', @@ -95,6 +102,9 @@ install_headers([ 'include/mzpeak/util/types.h', ], subdir : 'mzpeak') +# Sub-projects: +numpress = subproject('msnumpress') + # Library and executable targets: include = include_directories('include') @@ -117,6 +127,8 @@ deps = [ dependency('parquet', # Part of Arrow above. version : arrow_ver_str, required : true), + + numpress.get_variable('msnumpres_cpp_lib_so'), ] libmzpeak_so = library( @@ -133,6 +145,7 @@ libmzpeak_a = static_library( # Testing test_names = [ + 'algorithm', 'array_index', 'arrow', 'directory', diff --git a/src/data/array_index.cpp b/src/data/array_index.cpp index 6324057..c19b200 100644 --- a/src/data/array_index.cpp +++ b/src/data/array_index.cpp @@ -16,6 +16,89 @@ 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 +{ + switch (buffer_format) { + case MzPeak::Schema::BufferFormat::Point: + return true; + case MzPeak::Schema::BufferFormat::ChunkStart: + return true; + 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::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 { @@ -33,10 +116,68 @@ 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; + Entry const* chunk_values = 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: + chunk_values = &entry; + continue; + 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) { + 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); + } +} + +/******************************************************************************/ +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) , prefix_(obj.at("prefix").as_string()) + , layout_(group_name_to_layout(prefix_)) , entries_() , num_entities_() { @@ -116,6 +257,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 { diff --git a/src/data/signals.cpp b/src/data/signals.cpp index 6144bc1..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); @@ -136,6 +150,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()) { diff --git a/src/schema/group.cpp b/src/schema/group.cpp index 78e1702..63e1642 100644 --- a/src/schema/group.cpp +++ b/src/schema/group.cpp @@ -122,6 +122,12 @@ 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 +{ + return Util::Numpress::type_from_column_name(schema_name_); +} + /******************************************************************************/ Group::Group(const parquet::schema::GroupNode& node, const Schema::File& file) : name_("root") 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 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 {}; diff --git a/src/util/decoders.cpp b/src/util/decoders.cpp new file mode 100644 index 0000000..a262df7 --- /dev/null +++ b/src/util/decoders.cpp @@ -0,0 +1,58 @@ +/* + +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" +#include "mzpeak/util/numpress.h" + +namespace MzPeak::Util::Decoders { + +/******************************************************************************/ +std::size_t guess_array_length(const Schema::Column& column, + const std::shared_ptr& ary) +{ + std::optional numpress = column.second->possibly_numpress(); + + auto count = + [&numpress](const std::shared_ptr& nums) -> std::size_t { + if (numpress.has_value()) { + return Numpress::decoding_space_needed(nums->length(), numpress.value()); + } 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 diff --git a/src/util/numpress.cpp b/src/util/numpress.cpp new file mode 100644 index 0000000..f0369b8 --- /dev/null +++ b/src/util/numpress.cpp @@ -0,0 +1,128 @@ +/* + +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/compat.h" // IWYU pragma: keep +#include "mzpeak/util/decoders.h" +#include "mzpeak/util/numpress.h" + +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) +{ + 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> +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"); + throw InvalidFormatError(msg); + } + + std::vector bytes; + bytes.reserve(src->length()); + + Decoders::Scalar decoder; + decoder.decode(src, bytes); + + auto values = std::make_shared>(); + 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 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 +``` diff --git a/subprojects/msnumpress/README.md b/subprojects/msnumpress/README.md new file mode 100644 index 0000000..9c2d39d --- /dev/null +++ b/subprojects/msnumpress/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/subprojects/msnumpress/meson.build b/subprojects/msnumpress/meson.build new file mode 100644 index 0000000..441d045 --- /dev/null +++ b/subprojects/msnumpress/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/subprojects/msnumpress/pom.xml b/subprojects/msnumpress/pom.xml new file mode 100755 index 0000000..aad52df --- /dev/null +++ b/subprojects/msnumpress/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/subprojects/msnumpress/src/main/R/RMSNumpress/DESCRIPTION b/subprojects/msnumpress/src/main/R/RMSNumpress/DESCRIPTION new file mode 100644 index 0000000..7626a17 --- /dev/null +++ b/subprojects/msnumpress/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/subprojects/msnumpress/src/main/R/RMSNumpress/LICENSE b/subprojects/msnumpress/src/main/R/RMSNumpress/LICENSE new file mode 100644 index 0000000..024b33e --- /dev/null +++ b/subprojects/msnumpress/src/main/R/RMSNumpress/LICENSE @@ -0,0 +1,3 @@ +YEAR: 2020 +COPYRIGHT HOLDER: Justin Sing +ORGANIZATION: University of Toronto diff --git a/subprojects/msnumpress/src/main/R/RMSNumpress/NAMESPACE b/subprojects/msnumpress/src/main/R/RMSNumpress/NAMESPACE new file mode 100644 index 0000000..8ad0d37 --- /dev/null +++ b/subprojects/msnumpress/src/main/R/RMSNumpress/NAMESPACE @@ -0,0 +1,3 @@ +useDynLib(RMSNumpress, .registration=TRUE) +exportPattern("^[[:alpha:]]+") +importFrom(Rcpp, evalCpp) diff --git a/subprojects/msnumpress/src/main/R/RMSNumpress/R/RcppExports.R b/subprojects/msnumpress/src/main/R/RMSNumpress/R/RcppExports.R new file mode 100644 index 0000000..164cd39 --- /dev/null +++ b/subprojects/msnumpress/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/subprojects/msnumpress/src/main/R/RMSNumpress/inst/LICENSE.md b/subprojects/msnumpress/src/main/R/RMSNumpress/inst/LICENSE.md new file mode 100644 index 0000000..8cc3f91 --- /dev/null +++ b/subprojects/msnumpress/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/subprojects/msnumpress/src/main/R/RMSNumpress/man/RMSNumpress-package.Rd b/subprojects/msnumpress/src/main/R/RMSNumpress/man/RMSNumpress-package.Rd new file mode 100644 index 0000000..8d9f482 --- /dev/null +++ b/subprojects/msnumpress/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/subprojects/msnumpress/src/main/R/RMSNumpress/man/decodeLinear.Rd b/subprojects/msnumpress/src/main/R/RMSNumpress/man/decodeLinear.Rd new file mode 100644 index 0000000..e407b5e --- /dev/null +++ b/subprojects/msnumpress/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/subprojects/msnumpress/src/main/R/RMSNumpress/man/decodePic.Rd b/subprojects/msnumpress/src/main/R/RMSNumpress/man/decodePic.Rd new file mode 100644 index 0000000..5b3c095 --- /dev/null +++ b/subprojects/msnumpress/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/subprojects/msnumpress/src/main/R/RMSNumpress/man/decodeSlof.Rd b/subprojects/msnumpress/src/main/R/RMSNumpress/man/decodeSlof.Rd new file mode 100644 index 0000000..32c561d --- /dev/null +++ b/subprojects/msnumpress/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/subprojects/msnumpress/src/main/R/RMSNumpress/man/encodeLinear.Rd b/subprojects/msnumpress/src/main/R/RMSNumpress/man/encodeLinear.Rd new file mode 100644 index 0000000..1e95127 --- /dev/null +++ b/subprojects/msnumpress/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/subprojects/msnumpress/src/main/R/RMSNumpress/man/encodePic.Rd b/subprojects/msnumpress/src/main/R/RMSNumpress/man/encodePic.Rd new file mode 100644 index 0000000..8aadb6a --- /dev/null +++ b/subprojects/msnumpress/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/subprojects/msnumpress/src/main/R/RMSNumpress/man/encodeSlof.Rd b/subprojects/msnumpress/src/main/R/RMSNumpress/man/encodeSlof.Rd new file mode 100644 index 0000000..156c16a --- /dev/null +++ b/subprojects/msnumpress/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/subprojects/msnumpress/src/main/R/RMSNumpress/man/optimalLinearFixedPoint.Rd b/subprojects/msnumpress/src/main/R/RMSNumpress/man/optimalLinearFixedPoint.Rd new file mode 100644 index 0000000..3905c36 --- /dev/null +++ b/subprojects/msnumpress/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/subprojects/msnumpress/src/main/R/RMSNumpress/man/optimalLinearFixedPointMass.Rd b/subprojects/msnumpress/src/main/R/RMSNumpress/man/optimalLinearFixedPointMass.Rd new file mode 100644 index 0000000..f0f25f4 --- /dev/null +++ b/subprojects/msnumpress/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/subprojects/msnumpress/src/main/R/RMSNumpress/man/optimalSlofFixedPoint.Rd b/subprojects/msnumpress/src/main/R/RMSNumpress/man/optimalSlofFixedPoint.Rd new file mode 100644 index 0000000..1500997 --- /dev/null +++ b/subprojects/msnumpress/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/subprojects/msnumpress/src/main/R/RMSNumpress/src/MSNumpress.cpp b/subprojects/msnumpress/src/main/R/RMSNumpress/src/MSNumpress.cpp new file mode 100644 index 0000000..e3b6e17 --- /dev/null +++ b/subprojects/msnumpress/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/subprojects/msnumpress/src/main/R/RMSNumpress/src/RMSNumpress.cpp b/subprojects/msnumpress/src/main/R/RMSNumpress/src/RMSNumpress.cpp new file mode 100644 index 0000000..57520ac --- /dev/null +++ b/subprojects/msnumpress/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/subprojects/msnumpress/src/main/R/RMSNumpress/src/RcppExports.cpp b/subprojects/msnumpress/src/main/R/RMSNumpress/src/RcppExports.cpp new file mode 100644 index 0000000..e77e8d4 --- /dev/null +++ b/subprojects/msnumpress/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/subprojects/msnumpress/src/main/R/RMSNumpress/src/include/MSNumpress.hpp b/subprojects/msnumpress/src/main/R/RMSNumpress/src/include/MSNumpress.hpp new file mode 100644 index 0000000..525104d --- /dev/null +++ b/subprojects/msnumpress/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/subprojects/msnumpress/src/main/R/RMSNumpress/tests/testthat.R b/subprojects/msnumpress/src/main/R/RMSNumpress/tests/testthat.R new file mode 100644 index 0000000..895d5b2 --- /dev/null +++ b/subprojects/msnumpress/src/main/R/RMSNumpress/tests/testthat.R @@ -0,0 +1,4 @@ +library(testthat) +library(RMSNumpress) + +test_check("RMSNumpress") diff --git a/subprojects/msnumpress/src/main/R/RMSNumpress/tests/testthat/test_RMSNumpress.R b/subprojects/msnumpress/src/main/R/RMSNumpress/tests/testthat/test_RMSNumpress.R new file mode 100644 index 0000000..9658a4a --- /dev/null +++ b/subprojects/msnumpress/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/subprojects/msnumpress/src/main/cpp/MSNumpress.cpp b/subprojects/msnumpress/src/main/cpp/MSNumpress.cpp new file mode 100644 index 0000000..480c743 --- /dev/null +++ b/subprojects/msnumpress/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/subprojects/msnumpress/src/main/cpp/MSNumpress.hpp b/subprojects/msnumpress/src/main/cpp/MSNumpress.hpp new file mode 100644 index 0000000..985820f --- /dev/null +++ b/subprojects/msnumpress/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/subprojects/msnumpress/src/main/cpp/MSNumpressTest.cpp b/subprojects/msnumpress/src/main/cpp/MSNumpressTest.cpp new file mode 100644 index 0000000..a56543a --- /dev/null +++ b/subprojects/msnumpress/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/subprojects/msnumpress/src/main/csharp/MSNumpress.cs b/subprojects/msnumpress/src/main/csharp/MSNumpress.cs new file mode 100644 index 0000000..4b511c5 --- /dev/null +++ b/subprojects/msnumpress/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/subprojects/msnumpress/src/main/csharp/MSNumpressTest.cs b/subprojects/msnumpress/src/main/csharp/MSNumpressTest.cs new file mode 100644 index 0000000..9bdc71d --- /dev/null +++ b/subprojects/msnumpress/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/subprojects/msnumpress/src/main/java/ms/numpress/IntDecoder.java b/subprojects/msnumpress/src/main/java/ms/numpress/IntDecoder.java new file mode 100644 index 0000000..1d3e857 --- /dev/null +++ b/subprojects/msnumpress/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/subprojects/msnumpress/src/main/java/ms/numpress/MSNumpress.java b/subprojects/msnumpress/src/main/java/ms/numpress/MSNumpress.java new file mode 100644 index 0000000..97ffdf8 --- /dev/null +++ b/subprojects/msnumpress/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/subprojects/msnumpress/src/main/python/setup.py b/subprojects/msnumpress/src/main/python/setup.py new file mode 100644 index 0000000..a8de752 --- /dev/null +++ b/subprojects/msnumpress/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/subprojects/msnumpress/src/main/python/test_pymsnumpress.py b/subprojects/msnumpress/src/main/python/test_pymsnumpress.py new file mode 100644 index 0000000..89bab64 --- /dev/null +++ b/subprojects/msnumpress/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/subprojects/msnumpress/src/test/java/ms/numpress/MSNumpressTest.java b/subprojects/msnumpress/src/test/java/ms/numpress/MSNumpressTest.java new file mode 100644 index 0000000..45459fc --- /dev/null +++ b/subprojects/msnumpress/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 + +#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); +} 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()); 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"); }