diff --git a/CHANGELOG.rst b/CHANGELOG.rst index 56fcc917..716b3432 100644 --- a/CHANGELOG.rst +++ b/CHANGELOG.rst @@ -1,3 +1,14 @@ +Unreleased +---------- + +- Compression filters and libraries: + + * Added **Terse** filter: `TERSE/PROLIX (TRPX) v2.0 `_ (commit 93cf8e4) + +- Build: + + * Added support of C++20 compile flag handling for filter extensions (previously only for static libraries) + 7.1.0: 01/09/02026 ------------------ diff --git a/README.rst b/README.rst index 1167fae3..46775c73 100644 --- a/README.rst +++ b/README.rst @@ -4,7 +4,7 @@ hdf5plugin .. image:: https://zenodo.org/badge/DOI/10.5281/zenodo.7257761.svg :target: https://doi.org/10.5281/zenodo.7257761 -*hdf5plugin* provides `HDF5 compression filters `_ (namely: Blosc, Blosc2, BitShuffle, BZip2, Htj2k, FciDecomp, LZ4, Sperr, SZ, SZ3, Zfp, ZStd) and makes them usable from `h5py `_. +*hdf5plugin* provides `HDF5 compression filters `_ (namely: Blosc, Blosc2, BitShuffle, BZip2, Htj2k, FciDecomp, LZ4, Sperr, SZ, SZ3, TERSE/PROLIX, Zfp, ZStd) and makes them usable from `h5py `_. See `documentation `_. diff --git a/doc/contribute.rst b/doc/contribute.rst index a3c4687f..07ca3d71 100644 --- a/doc/contribute.rst +++ b/doc/contribute.rst @@ -282,6 +282,20 @@ compression_opts: - **psnr high** (big endian float64) - **psnr low** +terse +..... + +No options can be provided by the user, but the filter stores the following ``cd_values``: + +- **type_code**: Code of the dataset's integer data type: + + * 0: int16 + * 1: uint16 + * 2: int32 + * 3: uint32 + * 4: int8 + * 5: uint8 + zfp ... diff --git a/doc/index.rst b/doc/index.rst index aab2a825..93460012 100644 --- a/doc/index.rst +++ b/doc/index.rst @@ -1,7 +1,7 @@ hdf5plugin ========== -*hdf5plugin* provides `HDF5 compression filters `_ (namely: Blosc, Blosc2, BitShuffle, BZip2, FciDecomp, Htj2k, LZ4, Sperr, SZ, SZ3, Zfp, ZStd) and makes them usable from `h5py `_. +*hdf5plugin* provides `HDF5 compression filters `_ (namely: Blosc, Blosc2, BitShuffle, BZip2, FciDecomp, Htj2k, LZ4, Sperr, SZ, SZ3, TERSE/PROLIX, Zfp, ZStd) and makes them usable from `h5py `_. * Supported operating systems: Linux, Windows, macOS. diff --git a/doc/information.rst b/doc/information.rst index 8d66a70e..5f640444 100644 --- a/doc/information.rst +++ b/doc/information.rst @@ -83,6 +83,8 @@ HDF5 compression filters and compression libraries sources were obtained from: * `SZ3 plugin `_ (commit `4bbe9df7e4bcb `_) using `SZ3 `_ and ZStd. +* `TERSE/PROLIX (TRPX v2.0) plugin `_ + (commit `93cf8e4 `_). * `HDF5-ZFP plugin `_ (v1.1.1) using ZFP. * `ZStd plugin `_ (commit `5573db8 `_) using ZStd. @@ -120,6 +122,7 @@ Please read the different licenses: * SPERR: `lib/H5Z-SPERR/LICENSE `_ and `lib/SPERR/LICENSE `_ * SZ: `lib/SZ/copyright-and-BSD-license.txt `_ * SZ3: `lib/SZ3/copyright-and-BSD-license.txt `_ + * TERSE/PROLIX: `lib/trpx2.0/LICENSE `_ * zfp: `lib/H5Z-ZFP/LICENSE `_ * zstd: `lib/HDF5Plugin-Zstandard/LICENSE `_ diff --git a/doc/usage.rst b/doc/usage.rst index 8ca04494..84505df6 100644 --- a/doc/usage.rst +++ b/doc/usage.rst @@ -30,7 +30,7 @@ As for reading compressed datasets, ``import hdf5plugin`` is required to enable To create a compressed dataset use `h5py.Group.create_dataset`_ and set the ``compression`` and ``compression_opts`` arguments. -``hdf5plugin`` provides helpers to prepare those compression options: `Bitshuffle`_, `Blosc`_, `Blosc2`_, `BZip2`_, `FciDecomp`_, `Htj2k`_, `LZ4`_, `Sperr`_, `SZ`_, `SZ3`_, `Zfp`_, `Zstd`_. +``hdf5plugin`` provides helpers to prepare those compression options: `Bitshuffle`_, `Blosc`_, `Blosc2`_, `BZip2`_, `FciDecomp`_, `Htj2k`_, `LZ4`_, `Sperr`_, `SZ`_, `SZ3`_, `Terse`_, `Zfp`_, `Zstd`_. Sample code: @@ -133,6 +133,14 @@ SZ3 :inherited-members: FilterRefBase, Mapping :undoc-members: +Terse +===== + +.. autoclass:: Terse + :members: + :inherited-members: FilterRefBase, Mapping + :undoc-members: + Zfp === diff --git a/lib/trpx2.0/.gitignore b/lib/trpx2.0/.gitignore new file mode 100644 index 00000000..e41d1f8d --- /dev/null +++ b/lib/trpx2.0/.gitignore @@ -0,0 +1,5 @@ +.DS_Store +.vscode/ +build/ + + diff --git a/lib/trpx2.0/CMakeLists.txt b/lib/trpx2.0/CMakeLists.txt new file mode 100755 index 00000000..e430005e --- /dev/null +++ b/lib/trpx2.0/CMakeLists.txt @@ -0,0 +1,22 @@ +cmake_minimum_required(VERSION 3.15) +project(pyterse_and_tersecodec LANGUAGES CXX) + + +option(BUILD_PYTERSE "Build the pyterse Python extension" ON) +option(BUILD_TERSECODEC "Build the HDF5 filter library (tersecodec) plugin" ON) + + +set(CMAKE_CXX_STANDARD 20) +set(CMAKE_CXX_STANDARD_REQUIRED ON) +if(NOT CMAKE_BUILD_TYPE) + set(CMAKE_BUILD_TYPE Release) +endif() + + +if(BUILD_PYTERSE) + add_subdirectory(pyterse) +endif() + +if(BUILD_TERSECODEC) + add_subdirectory(tersecodec) +endif() diff --git a/lib/trpx2.0/DIALS_integration/FormatTRPX.py b/lib/trpx2.0/DIALS_integration/FormatTRPX.py new file mode 100644 index 00000000..dcc11e3f --- /dev/null +++ b/lib/trpx2.0/DIALS_integration/FormatTRPX.py @@ -0,0 +1,171 @@ +from __future__ import annotations + +import os +import re +import xml.etree.ElementTree as ET +import numpy as np +import sys + + +sys.path.append(os.path.join(os.getcwd(), 'build')) + +from pyterse import Terse + +from scitbx.array_family import flex +from dxtbx.format.Format import Format +from dxtbx.model import ScanFactory +from dxtbx.model.beam import Probe + + +class FormatTRPX(Format): + @staticmethod + def understand(image_file): + try: + with FormatTRPX.open_file(image_file, "rb") as fh: + header = fh.read(24) + print(f"Header: {header}") + + except OSError: + return False + + # Check if header starts with '' + raise RuntimeError("Failed to find the end of the header.") + header_data += chunk + if b'/>' in header_data: + end_index = header_data.find(b'/>') + 2 + header_data = header_data[:end_index] + break + + try: + header_str = header_data.decode('utf-8') + except UnicodeDecodeError as e: + raise RuntimeError(f"Failed to decode header data as UTF-8: {e}") + + try: + root = ET.fromstring(header_str) + except ET.ParseError as e: + raise RuntimeError(f"Failed to parse XML header: {e}") + + hd['prolix_bits'] = int(root.attrib.get('prolix_bits', '12')) + hd['signed'] = int(root.attrib.get('signed', '0')) + hd['block'] = int(root.attrib.get('block', '12')) + hd['memory_size'] = int(root.attrib.get('memory_size', '0')) + hd['number_of_values'] = int(root.attrib.get('number_of_values', '0')) + hd['number_of_frames'] = int(root.attrib.get('number_of_frames', '1')) + + hd['pixel_size'] = (0.055, 0.055) + hd['trusted_range'] = (0, 65535) + hd['distance'] = 478.0 + hd['start_tilt'] = 0.0 + hd['delta_tilt'] = 0.0 + hd['exposure'] = 1.0 + + if 'dimensions' in root.attrib: + dimensions_str = root.attrib['dimensions'] + dimensions = list(map(int, dimensions_str.strip().split())) + hd['image_size'] = dimensions + else: + num_values = hd['number_of_values'] + image_dimension = int(np.sqrt(num_values)) + hd['image_size'] = (image_dimension, image_dimension) + + hd['distance'] = float(root.attrib.get('distance', '478.0')) + pixel_size_str = root.attrib.get('pixel_size', '0.055 0.055') + hd['pixel_size'] = tuple(map(float, pixel_size_str.strip().split())) + trusted_range_str = root.attrib.get('trusted_range', '0 65535') + hd['trusted_range'] = tuple(map(float, trusted_range_str.strip().split())) + + return hd + + + def _start(self): + """Open the image file and read useful metadata into an internal dictionary""" + self._header_dictionary = self._read_metadata(self._image_file) + + def _goniometer(self): + """Dummy goniometer, 'vertical' as the images are viewed.""" + goniometer = self._goniometer_factory.known_axis((0, -1, 0)) + return goniometer + + def _detector(self): + beam_centre = [ + (p * i) / 2 + for p, i in zip( + self._header_dictionary["pixel_size"], + self._header_dictionary["image_size"], + ) + ] + d = self._detector_factory.simple( + "PAD", + self._header_dictionary["distance"], + beam_centre, + "+x", + "-y", + self._header_dictionary["pixel_size"], + self._header_dictionary["image_size"], + self._header_dictionary["trusted_range"], + ) + return d + + + def _beam(self): + """Unpolarized beam, default energy 200 keV""" + beam = self._beam_factory.make_polarized_beam( + sample_to_source=(0.0, 0.0, 1.0), + wavelength=0.02508, + polarization=(0, 1, 0), + polarization_fraction=0.5, + probe=Probe.electron, + ) + return beam + + def _scan(self): + """Scan model for this image, filling out any unavailable items with dummy values""" + alpha = self._header_dictionary.get("start_tilt", 0.0) + dalpha = self._header_dictionary.get("delta_tilt", 1.0) + exposure = self._header_dictionary.get("exposure", 0.0) + oscillation = (alpha, dalpha) + fname = os.path.split(self._image_file)[-1] + s = fname.split("_")[-1].split(".")[0] + try: + index = int(re.match(".*?([0-9]+)$", s).group(1)) + except AttributeError: + index = 1 + scan = ScanFactory.make_scan((index, index), exposure, oscillation, {index: 0}) + return scan + + def get_raw_data(self): + terse = Terse.load(self._image_file) + decompressed_data = terse.prolix() + raw_data_flex = flex.double(decompressed_data) + return raw_data_flex + +if __name__ == "__main__": + import sys + if len(sys.argv) < 2: + print("Usage: python script.py image_file.trpx") + sys.exit(1) + image_file = sys.argv[1] + try: + if FormatTRPX.understand(image_file): + format_instance = FormatTRPX(image_file) + format_instance._start() + raw_data = format_instance.get_raw_data() + print("Raw data:", raw_data) + else: + print("File format not understood.") + except Exception as e: + print(f"An error occurred: {e}") diff --git a/lib/trpx2.0/Doxygen/Doxyfile b/lib/trpx2.0/Doxygen/Doxyfile new file mode 100644 index 00000000..25e2db8d --- /dev/null +++ b/lib/trpx2.0/Doxygen/Doxyfile @@ -0,0 +1,180 @@ +PROJECT_NAME = "TRPX compression algorithm" +PROJECT_BRIEF = "TRPX is a fast and efficient algorithm for compressing diffraction data" +PROJECT_LOGO = +OUTPUT_DIRECTORY = docs +CREATE_SUBDIRS = YES +ALLOW_UNICODE_NAMES = YES +OUTPUT_LANGUAGE = English +BRIEF_MEMBER_DESC = YES +REPEAT_BRIEF = YES +ALWAYS_DETAILED_SEC = NO +INLINE_INHERITED_MEMB = NO +FULL_PATH_NAMES = YES +SHORT_NAMES = NO +JAVADOC_AUTOBRIEF = YES +QT_AUTOBRIEF = YES +MULTILINE_CPP_IS_BRIEF = NO +INHERIT_DOCS = YES +SEPARATE_MEMBER_PAGES = NO +TAB_SIZE = 4 +OPTIMIZE_OUTPUT_FOR_C = NO +OPTIMIZE_OUTPUT_JAVA = NO +OPTIMIZE_FOR_FORTRAN = NO +OPTIMIZE_OUTPUT_VHDL = NO +MARKDOWN_SUPPORT = YES +TOC_INCLUDE_HEADINGS = 5 +AUTOLINK_SUPPORT = YES +BUILTIN_STL_SUPPORT = YES +CPP_CLI_SUPPORT = YES +SIP_SUPPORT = YES +IDL_PROPERTY_SUPPORT = YES +DISTRIBUTE_GROUP_DOC = NO +GROUP_NESTED_COMPOUNDS = NO +SUBGROUPING = YES +INLINE_GROUPED_CLASSES = NO +INLINE_SIMPLE_STRUCTS = NO +TYPEDEF_HIDES_STRUCT = NO +LOOKUP_CACHE_SIZE = 0 +EXTRACT_ALL = YES +EXTRACT_PRIVATE = YES +EXTRACT_PRIV_VIRTUAL = NO +EXTRACT_PACKAGE = NO +EXTRACT_STATIC = YES +EXTRACT_LOCAL_CLASSES = YES +EXTRACT_LOCAL_METHODS = YES +EXTRACT_ANON_NSPACES = YES +HIDE_UNDOC_MEMBERS = NO +HIDE_UNDOC_CLASSES = NO +HIDE_FRIEND_COMPOUNDS = NO +HIDE_IN_BODY_DOCS = NO +INTERNAL_DOCS = NO +CASE_SENSE_NAMES = YES +HIDE_SCOPE_NAMES = NO +HIDE_COMPOUND_REFERENCE= NO +SHOW_INCLUDE_FILES = YES +SHOW_GROUPED_MEMB_INC = NO +FORCE_LOCAL_INCLUDES = NO +INLINE_INFO = YES +SORT_MEMBER_DOCS = YES +SORT_BRIEF_DOCS = NO +SORT_MEMBERS_CTORS_1ST = NO +SORT_GROUP_NAMES = NO +SORT_BY_SCOPE_NAME = NO +STRICT_PROTO_MATCHING = NO +GENERATE_TODOLIST = YES +GENERATE_TESTLIST = YES +GENERATE_BUGLIST = YES +GENERATE_DEPRECATEDLIST= YES +MAX_INITIALIZER_LINES = 30 +SHOW_USED_FILES = YES +SHOW_FILES = YES +SHOW_NAMESPACES = YES +QUIET = NO +WARNINGS = YES +WARN_IF_UNDOCUMENTED = YES +WARN_IF_DOC_ERROR = YES +WARN_NO_PARAMDOC = NO +WARN_FORMAT = "$file:$line: $text" +INPUT = ./include ./src ./pyterse ./tersecodec +INPUT_ENCODING = UTF-8 +FILE_PATTERNS = *.c *.cc *.cxx *.cpp *.c++ *.h *.hh *.hxx *.hpp *.h++ +RECURSIVE = YES +EXCLUDE = build +EXCLUDE_SYMLINKS = NO +EXCLUDE_PATTERNS = */build/* */test/* +EXCLUDE_SYMBOLS = +EXAMPLE_PATH = +EXAMPLE_PATTERNS = * +EXAMPLE_RECURSIVE = NO +IMAGE_PATH = +INPUT_FILTER = +FILTER_PATTERNS = +FILTER_SOURCE_FILES = NO +FILTER_SOURCE_PATTERNS = +USE_MDFILE_AS_MAINPAGE = README.md +SOURCE_BROWSER = YES +INLINE_SOURCES = YES +STRIP_CODE_COMMENTS = NO +REFERENCED_BY_RELATION = YES +REFERENCES_RELATION = YES +REFERENCES_LINK_SOURCE = YES +SOURCE_TOOLTIPS = YES +USE_HTAGS = NO +VERBATIM_HEADERS = YES +CLANG_ASSISTED_PARSING = YES +ALPHABETICAL_INDEX = YES +GENERATE_HTML = YES +HTML_OUTPUT = html +HTML_FILE_EXTENSION = .html +HTML_HEADER = +HTML_FOOTER = +HTML_STYLESHEET = +HTML_EXTRA_STYLESHEET = doxygen-awesome.css +HTML_EXTRA_FILES = +HTML_COLORSTYLE_HUE = 220 +HTML_COLORSTYLE_SAT = 100 +HTML_COLORSTYLE_GAMMA = 80 +HTML_TIMESTAMP = NO +HTML_DYNAMIC_MENUS = YES +HTML_DYNAMIC_SECTIONS = NO +HTML_INDEX_NUM_ENTRIES = 100 +GENERATE_DOCSET = NO +GENERATE_HTMLHELP = NO +GENERATE_QHP = NO +GENERATE_ECLIPSEHELP = NO +GENERATE_TREEVIEW = YES +ENUM_VALUES_PER_LINE = 4 +TREEVIEW_WIDTH = 250 +EXT_LINKS_IN_WINDOW = NO +FORMULA_FONTSIZE = 10 +FORMULA_TRANSPARENT = YES +USE_MATHJAX = NO +SEARCHENGINE = YES +SERVER_BASED_SEARCH = NO +EXTERNAL_SEARCH = NO +SEARCHDATA_FILE = searchdata.xml +GENERATE_LATEX = NO +GENERATE_RTF = NO +GENERATE_MAN = NO +GENERATE_XML = NO +GENERATE_DOCBOOK = NO +GENERATE_AUTOGEN_DEF = NO +GENERATE_PERLMOD = NO +ENABLE_PREPROCESSING = YES +MACRO_EXPANSION = YES +EXPAND_ONLY_PREDEF = NO +SEARCH_INCLUDES = YES +PREDEFINED = +SKIP_FUNCTION_MACROS = YES +TAGFILES = +ALLEXTERNALS = NO +EXTERNAL_GROUPS = YES +EXTERNAL_PAGES = YES +CLASS_DIAGRAMS = YES +DIA_PATH = +HIDE_UNDOC_RELATIONS = YES +HAVE_DOT = YES +DOT_NUM_THREADS = 0 +DOT_FONTNAME = Helvetica +DOT_FONTSIZE = 10 +DOT_FONTPATH = +CLASS_GRAPH = YES +COLLABORATION_GRAPH = YES +GROUP_GRAPHS = YES +UML_LOOK = NO +UML_LIMIT_NUM_FIELDS = 10 +TEMPLATE_RELATIONS = YES +INCLUDE_GRAPH = YES +INCLUDED_BY_GRAPH = YES +CALL_GRAPH = YES +CALLER_GRAPH = YES +GRAPHICAL_HIERARCHY = YES +DIRECTORY_GRAPH = YES +DOT_IMAGE_FORMAT = svg +INTERACTIVE_SVG = YES +DOT_GRAPH_MAX_NODES = 50 +MAX_DOT_GRAPH_DEPTH = 0 +DOT_TRANSPARENT = NO +DOT_MULTI_TARGETS = NO +GENERATE_LEGEND = YES +DOT_CLEANUP = YES \ No newline at end of file diff --git a/lib/trpx2.0/Doxygen/doxygen-awesome.css b/lib/trpx2.0/Doxygen/doxygen-awesome.css new file mode 100644 index 00000000..8d27c8c9 --- /dev/null +++ b/lib/trpx2.0/Doxygen/doxygen-awesome.css @@ -0,0 +1,195 @@ +:root { + --primary-color: #1779c4; + --primary-dark-color: #0e5789; + --primary-light-color: #4da3e5; + --secondary-color: #37474f; + --accent-color: #03a9f4; + --background-color: #ffffff; + --alt-background-color: #f7f7f7; + --text-color: #37474f; + --link-color: var(--primary-color); + --link-hover-color: var(--primary-dark-color); + --code-background: #f5f5f5; + --code-color: #1f1f1f; + --border-color: #e0e0e0; + --font-family: -apple-system, BlinkMacSystemFont, "Segoe UI", Roboto, Oxygen, Ubuntu, Cantarell, "Fira Sans", "Droid Sans", "Helvetica Neue", Arial, sans-serif; + --monospace-font-family: "SFMono-Regular", Consolas, "Liberation Mono", Menlo, monospace; + } + + html { + font-size: 16px; + font-family: var(--font-family); + color: var(--text-color); + background-color: var(--background-color); + } + + body, table, div, p, dl { + font-family: var(--font-family); + font-size: 1rem; + line-height: 1.6; + } + + /* Header styling */ + #titlearea { + background-color: var(--primary-color); + color: white; + border-bottom: none; + padding: 15px 0; + margin-bottom: 20px; + } + + #projectname { + font-size: 2.2rem; + font-weight: bold; + } + + #projectbrief { + font-size: 1.2rem; + margin-top: 5px; + } + + /* Navigation tree */ + #nav-tree { + background-color: var(--alt-background-color); + border-right: 1px solid var(--border-color); + } + + #nav-tree .selected { + background-image: none; + background-color: var(--primary-light-color); + color: white; + text-shadow: none; + } + + /* Links */ + a { + color: var(--link-color); + text-decoration: none; + } + + a:hover { + color: var(--link-hover-color); + text-decoration: underline; + } + + /* Code blocks */ + div.fragment, pre.fragment { + background-color: var(--code-background); + border: 1px solid var(--border-color); + border-radius: 4px; + padding: 12px; + margin: 16px 0; + overflow: auto; + font-family: var(--monospace-font-family); + font-size: 0.9rem; + } + + /* Syntax highlighting */ + span.keyword { + color: #0066cc; + font-weight: bold; + } + + span.keywordtype { + color: #007020; + font-weight: bold; + } + + span.comment { + color: #60a0b0; + font-style: italic; + } + + span.preprocessor { + color: #007020; + } + + span.stringliteral { + color: #d32f2f; + } + + /* Tables */ + table.memberdecls, table.fieldtable { + border-spacing: 0; + width: 100%; + margin-bottom: 20px; + } + + table.memberdecls td, table.fieldtable td { + padding: 8px; + border-bottom: 1px solid var(--border-color); + } + + table.memberdecls th, table.fieldtable th { + background-color: var(--primary-color); + color: white; + padding: 8px; + text-align: left; + } + + /* Memitem blocks */ + .memitem { + margin-bottom: 24px; + border: 1px solid var(--border-color); + border-radius: 4px; + overflow: hidden; + } + + .memproto { + background-color: var(--alt-background-color); + padding: 12px; + border-bottom: 1px solid var(--border-color); + font-weight: bold; + } + + .memdoc { + padding: 12px; + } + + /* Navigation */ + .tabs, .tabs2, .tabs3 { + background-image: none; + background-color: var(--primary-color); + color: white; + } + + .tablist li { + background-image: none !important; + } + + .tablist a { + background-image: none !important; + color: white !important; + text-shadow: none !important; + } + + .tablist a:hover { + background-color: var(--primary-dark-color) !important; + } + + .tablist li.current a { + background-color: var(--primary-dark-color) !important; + } + + /* Mobile responsiveness */ + @media (max-width: 768px) { + #top { + display: flex; + flex-direction: column; + } + + #side-nav { + width: 100% !important; + height: auto !important; + margin-bottom: 20px; + } + + #nav-tree { + width: 100% !important; + } + + #doc-content { + margin-left: 0 !important; + height: auto !important; + } + } \ No newline at end of file diff --git a/lib/trpx2.0/ImageJ/TRPX_Reader.java b/lib/trpx2.0/ImageJ/TRPX_Reader.java new file mode 100644 index 00000000..b6bfe14c --- /dev/null +++ b/lib/trpx2.0/ImageJ/TRPX_Reader.java @@ -0,0 +1,547 @@ +import java.io.BufferedReader; +import java.io.FileReader; +import java.io.FileInputStream; +import java.io.FileNotFoundException; +import java.io.IOException; +import java.math.BigInteger; +import java.nio.file.Paths; +import java.nio.file.Path; +import java.util.Arrays; +import java.util.Scanner; +import java.util.regex.Pattern; +import ij.IJ; +import ij.ImagePlus; +import ij.process.ImageProcessor; +import ij.gui.NewImage; +import ij.ImageStack; +import ij.plugin.PlugIn; +import ij.io.OpenDialog; + +public class TRPX_Reader implements PlugIn { + + public void run(String arg) { + + + OpenDialog dialog = new OpenDialog("Select a Terse/Prolix data file", null); + String filePath = dialog.getPath(); + if (filePath == null) return; + if (!filePath.toLowerCase().endsWith(".trpx")) return; + + // Parameters required for unpacking the file + long prolixBits = 0; + long signed=1; + long block=0; + long terseDataSize = 0; + long imageSize = 0; + long dataStartIndex; + int dim0 = 0; + int dim1 = 0; + int nFrames = 1; + long[] frameSizes = null; + long[] metadataSizes = null; + + + try (BufferedReader br = new BufferedReader(new FileReader(filePath))) { + String line; + dataStartIndex = 0; + while ((line = br.readLine()) != null) { + int indexTerse = line.indexOf("", indexTerse) + 2; + if (indexEndTag != -1) { + dataStartIndex += indexEndTag; + try (Scanner scanner = new Scanner(line.substring(indexTerse, indexEndTag))) { + scanner.findWithinHorizon("<.*?prolix_bits=\"(\\d+)\".*?/>", 0); + prolixBits = Long.parseLong(scanner.match().group(1)); + } + try (Scanner scanner = new Scanner(line.substring(indexTerse, indexEndTag))) { + scanner.findWithinHorizon("<.*?signed=\"(\\d+)\".*?/>", 0); + signed = Long.parseLong(scanner.match().group(1)); + } + try (Scanner scanner = new Scanner(line.substring(indexTerse, indexEndTag))) { + scanner.findWithinHorizon("<.*?block=\"(\\d+)\".*?/>", 0); + block = Long.parseLong(scanner.match().group(1)); + } + try (Scanner scanner = new Scanner(line.substring(indexTerse, indexEndTag))) { + scanner.findWithinHorizon("<.*?memory_size=\"(\\d+)\".*?/>", 0); + terseDataSize = Long.parseLong(scanner.match().group(1)); + } + try (Scanner scanner = new Scanner(line.substring(indexTerse, indexEndTag))) { + scanner.findWithinHorizon("<.*?number_of_values=\"(\\d+)\".*?/>", 0); + imageSize = Long.parseLong(scanner.match().group(1)); + } + try (Scanner scanner = new Scanner(line.substring(indexTerse, indexEndTag))) { + if (scanner.findWithinHorizon("<.*?number_of_frames=\"(\\d+)\".*?/>", 0) != null) + nFrames = Integer.parseInt(scanner.match().group(1)); + } + try (Scanner scanner = new Scanner(line.substring(indexTerse, indexEndTag))) { + if (scanner.findWithinHorizon("<.*?dimensions=\"(\\d+)(?:\\s+(\\d+))?(?:\\s+(\\d+))?\".*?/>", 0) != null) { + dim0 = Integer.parseInt(scanner.match().group(1)); + dim1 = scanner.match().group(2) != null ? Integer.parseInt(scanner.match().group(2)) : dim1; + } else { + dim0 = (int) Math.sqrt(imageSize); + dim1 = dim0; + } + } + // Parse memory_sizes_of_frames if present + try (Scanner scanner = new Scanner(line.substring(indexTerse, indexEndTag))) { + if (scanner.findWithinHorizon("<.*?memory_sizes_of_frames=\"([^\"]+)\".*?/>", 0) != null) { + String[] sizes = scanner.match().group(1).trim().split("\\s+"); + frameSizes = new long[sizes.length]; + for (int i = 0; i < sizes.length; i++) { + frameSizes[i] = Long.parseLong(sizes[i]); + } + } + } + // Parse metadata_string_sizes if present + try (Scanner scanner = new Scanner(line.substring(indexTerse, indexEndTag))) { + if (scanner.findWithinHorizon("<.*?metadata_string_sizes=\"([^\"]+)\".*?/>", 0) != null) { + String[] sizes = scanner.match().group(1).trim().split("\\s+"); + metadataSizes = new long[sizes.length]; + for (int i = 0; i < sizes.length; i++) { + metadataSizes[i] = Long.parseLong(sizes[i]); + } + } + } + } + break; + } + } + } + catch (FileNotFoundException e) {return;} + catch (IOException e) {return;} + + + if (signed != 0 || prolixBits > 16) { + System.err.println("Error: Invalid .trpx images for this plugin: images must be unsigned 16 bit."); + System.err.println("Data of " + filePath + " are " + ((signed == 0) ? "un" : "") + "signed " + prolixBits + " bit"); + return; + } + + if (terseDataSize <= 0 || imageSize <= 0 || dim0 <= 0 || dim1 <= 0) { + System.err.println("Error: Invalid parameters in .trpx file - data size: " + terseDataSize + + ", image size: " + imageSize + ", dimensions: " + dim0 + "x" + dim1); + return; + } + + long totalMetadataSize = 0; + if (metadataSizes != null) { + for (long size : metadataSizes) { + totalMetadataSize += size; + } + } + + + System.out.println("TRPX File Parameters:"); + System.out.println(" Prolix bits: " + prolixBits); + System.out.println(" Signed: " + signed); + System.out.println(" Block size: " + block); + System.out.println(" Image size: " + imageSize); + System.out.println(" Dimensions: " + dim0 + "x" + dim1); + System.out.println(" Frames: " + nFrames); + System.out.println(" Terse data size: " + terseDataSize); + System.out.println(" Data start index: " + dataStartIndex); + System.out.println(" Total metadata size: " + totalMetadataSize); + + + ImagePlus imageStack = NewImage.createShortImage(filePath, dim0, dim1, nFrames, NewImage.FILL_BLACK); + + + dTerseData = new byte[(int) (dataStartIndex + totalMetadataSize + terseDataSize + 2)]; + Path path = Paths.get(filePath); + try (FileInputStream fileInputStream = new FileInputStream(path.toFile())) { + fileInputStream.read(dTerseData); + } + catch (IOException e) { return; } + + + dBitStart = (dataStartIndex + totalMetadataSize) * 8; + + + long currentByteOffset = dataStartIndex + totalMetadataSize; + + for (int frameNumber = 1; frameNumber <= nFrames; ++frameNumber) { + System.out.println("Processing frame " + frameNumber + " at bit position: " + dBitStart); + + // Detect compression mode + int compressionMode = detectCompressionMode(); + System.out.println("Frame " + frameNumber + " compression mode: " + + (compressionMode == MODE_SIGNED ? "Signed" : + compressionMode == MODE_UNSIGNED ? "Unsigned" : "Small_unsigned")); + + imageStack.setSlice(frameNumber); + ImageProcessor ip = imageStack.getProcessor(); + short[] pixels = (short[]) ip.getPixels(); + + if (compressionMode == MODE_UNSIGNED) { + decompressUnsignedFrame(pixels, (int)imageSize, (int)block, (int)prolixBits); + } else if (compressionMode == MODE_SMALL_UNSIGNED) { + decompressSmallUnsignedFrame(pixels, (int)imageSize, (int)block, (int)prolixBits); + } else { + + decompressSignedFrame(pixels, (int)imageSize, (int)block); + } + + // Set start of next frame using frame sizes from XML if available + if (frameSizes != null && frameNumber < nFrames) { + currentByteOffset += frameSizes[frameNumber - 1]; + dBitStart = currentByteOffset * 8; + System.out.println("Next frame starts at byte " + currentByteOffset + " (bit " + dBitStart + ")"); + } + + + } + + + imageStack.show(); + IJ.run(imageStack, "Enhance Contrast", "saturated=0.35"); + imageStack.updateAndDraw(); + } + + private long dBitStart; + private byte[] dTerseData; + + + private static final int MODE_SIGNED = 0; + private static final int MODE_UNSIGNED = 1; + private static final int MODE_SMALL_UNSIGNED = 2; + + + private int detectCompressionMode() { + // Peek at the next 18 bits to check if it's a mode prefix + long savedBitStart = dBitStart; + int prefix = ToUnsignedInt(18); + + System.out.println("Checking 18-bit value at bit " + savedBitStart + ": 0b" + + String.format("%18s", Integer.toBinaryString(prefix)).replace(' ', '0') + + " (decimal: " + prefix + ")"); + + if (prefix == 0b111111111111111100) { + // Small_unsigned mode - prefix consumed, don't rewind + System.out.println("Detected: Small_unsigned mode (prefix consumed)"); + return MODE_SMALL_UNSIGNED; + } else if (prefix == 0b111111111111111000) { + // Unsigned mode - prefix consumed, don't rewind + System.out.println("Detected: Unsigned mode (prefix consumed)"); + return MODE_UNSIGNED; + } else { + // Signed mode - NO prefix exists, so rewind to start of frame data + System.out.println("Detected: Signed mode (no prefix, rewinding to bit " + savedBitStart + ")"); + dBitStart = savedBitStart; + return MODE_SIGNED; + } + } + + + private void decompressSignedFrame(short[] pixels, int imageSize, int block) { + short significant_bits = 0; + for (int from = 0; from < imageSize; from += block) { + if (0 == ToShort(1)) + if (7 == (significant_bits = ToShort(3))) + if (10 == (significant_bits += ToShort(2))) + significant_bits += ToShort(6); + int to = (int) Math.min(imageSize, from + block); + if (significant_bits == 0) + Arrays.fill(pixels, from, to, (short)0); + else + for (int j = from; j < to; ++j) + pixels[j] = ToShort(significant_bits); + } + } + + + private void decompressUnsignedFrame(short[] pixels, int imageSize, int block, int prolixBits) { + short significant_bits = 0; + short masked_bits = 0; + + for (int from = 0; from < imageSize; from += block) { + if (0 == ToShort(1)) + if (7 == (significant_bits = ToShort(3))) + if (10 == (significant_bits += ToShort(2))) + significant_bits += ToShort(6); + + int to = (int) Math.min(imageSize, from + block); + + if (significant_bits != prolixBits) { + if (significant_bits == 0) + Arrays.fill(pixels, from, to, (short)0); + else + for (int j = from; j < to; ++j) + pixels[j] = ToShort(significant_bits); + } else { + if (0 == ToShort(1)) + if (7 == (masked_bits = ToShort(3))) + if (10 == (masked_bits += ToShort(2))) + masked_bits += ToShort(6); + + + for (int j = from; j < to; ++j) { + int value = ToShort(masked_bits) - 1; + pixels[j] = (short) Math.max(0, value); + } + } + } + } + + private short ToShort(int s) { + int indx = (int) dBitStart >> 3; + + if (indx + 2 >= dTerseData.length) { + System.err.println("Warning: Trying to read beyond data bounds. Index: " + indx + + ", Array length: " + dTerseData.length + ", Bit position: " + dBitStart); + dBitStart += s; + return 0; + } + + int trpx = ((dTerseData[indx] & 0xFF) + + ((dTerseData[1 + indx] & 0xFF) << 8) + + ((dTerseData[2 + indx] & 0xFF) << 16)); + short result = (short)((trpx >>> (dBitStart & 7)) & ((1 << s) - 1)); + dBitStart += s; + return result; + } + + + private int ToUnsignedInt(long s) { + int indx = (int) dBitStart >> 3; + + if (indx + 4 >= dTerseData.length) { + System.err.println("Warning: Trying to read beyond data bounds for unsigned int. Index: " + indx + + ", Array length: " + dTerseData.length + ", Bit position: " + dBitStart); + dBitStart += s; + return 0; + } + + long trpx = ((dTerseData[indx] & 0xFFL) + + ((dTerseData[1 + indx] & 0xFFL) << 8) + + ((dTerseData[2 + indx] & 0xFFL) << 16) + + ((dTerseData[3 + indx] & 0xFFL) << 24)); + int result = (int)((trpx >>> (dBitStart & 7)) & ((1L << s) - 1)); + dBitStart += s; + return result; + } + + private BigInteger ToBigInteger(long numBits) { + if (numBits <= 0) return BigInteger.ZERO; + if (numBits <= 32) return BigInteger.valueOf(ToUnsignedInt(numBits)); + + BigInteger result = BigInteger.ZERO; + long bitsRemaining = numBits; + int shift = 0; + + while (bitsRemaining > 0) { + long bitsToRead = Math.min(32, bitsRemaining); + long chunk = ToUnsignedInt(bitsToRead); + result = result.or(BigInteger.valueOf(chunk).shiftLeft(shift)); + shift += (int)bitsToRead; + bitsRemaining -= bitsToRead; + } + + return result; + } + + // Small_unsigned decompression + private void decompressSmallUnsignedFrame(short[] pixels, int imageSize, int block, int prolixBits) { + + + int effectiveBlock = Math.min(block, 24); + int bits = 0; + int max = 0; + + for (int from = 0; from < imageSize; from += effectiveBlock) { + int to = Math.min(imageSize, from + effectiveBlock); + int blockSize = to - from; + + + int[] result = getMaxValueSmallUnsigned(max, bits); + max = result[0]; + bits = result[1]; + + + + switch (max) { + case 0: + // All zeros + Arrays.fill(pixels, from, to, (short)0); + break; + case 1: + // 1-bit per value + for (int j = from; j < to; j++) { + pixels[j] = (short)ToShort(1); + } + break; + case 2: + // Base-3 encoding + decompressBase3BlockSmall(pixels, from, to, effectiveBlock); + break; + case 3: + // 2-bit per value + for (int j = from; j < to; j++) { + pixels[j] = (short)ToShort(2); + } + break; + case 7: + // 3-bit per value + for (int j = from; j < to; j++) { + pixels[j] = (short)ToShort(3); + } + break; + default: + if (max < 7) { + // Base-(max+1) encoding for small values + decompressBaseNBlockSmall(pixels, from, to, max, effectiveBlock); + } else if (bits == prolixBits) { + + decompressSmallUnsignedMaskedData(pixels, from, imageSize, effectiveBlock, prolixBits, max, bits); + return; + } else { + // Regular bit encoding for larger values + for (int j = from; j < to; j++) { + pixels[j] = (short)ToShort(bits); + } + max = 0xFFFF / 2; + } + break; + } + } + } + + private int[] getMaxValueSmallUnsigned(int max, int bits) { + int flag = ToShort(1); + if (flag == 1 && max == 0) { + return new int[]{max, bits}; + } + + flag = (flag << 1) + ToShort(1); + if (flag == 0b11) { + // Same as previous block + return new int[]{max, bits}; + } + if (flag == 0b10) { + + bits--; + max--; + } else if (flag == 0b01) { + // Increment from previous + bits++; + max = (max == 6) ? max - 2 : max + 1; + } else { + // flag == 0b00: read new value + max = bits = ToShort(3); + if (bits == 7) { + max = 0xFFFF / 2; + bits = 3 + ToShort(3); + if (bits == 10) { + bits += ToShort(3); + if (bits == 17) { + bits += ToShort(6); + } + } + } + } + return new int[]{max, bits}; + } + private void decompressBase3BlockSmall(short[] pixels, int from, int to, int effectiveBlock) { + int blockSize = to - from; + long numBits = (long)Math.ceil(effectiveBlock * Math.log(3) / Math.log(2)); + BigInteger val = ToBigInteger(numBits); + BigInteger three = BigInteger.valueOf(3); + for (int i = 0; i < blockSize; i++) { + pixels[from + i] = (short)(val.mod(three).intValue()); + val = val.divide(three); + } + } + + + private void decompressBaseNBlockSmall(short[] pixels, int from, int to, int max, int effectiveBlock) { + int blockSize = to - from; + int base = max + 1; + long numBits = (long)Math.ceil(effectiveBlock * Math.log(base) / Math.log(2)); + BigInteger val = ToBigInteger(numBits); + BigInteger bigBase = BigInteger.valueOf(base); + for (int i = 0; i < blockSize; i++) { + pixels[from + i] = (short)(val.mod(bigBase).intValue()); + val = val.divide(bigBase); + } + } + + + + private void decompressSmallUnsignedMaskedData(short[] pixels, int from, int imageSize, int effectiveBlock, int prolixBits, int max, int bits) { + int to = Math.min(imageSize, from + effectiveBlock); + + for (int j = from; j < to; j++) { + int value = ToShort(bits) - 1; + pixels[j] = (short) Math.max(0, value); + } + + from += effectiveBlock; + + max = 0xFFFF / 2; + bits = prolixBits + 2; + + + while (from < imageSize) { + to = Math.min(imageSize, from + effectiveBlock); + + int[] result = getMaxValueSmallUnsigned(max, bits); + max = result[0]; + bits = result[1]; + + if (max >= 7) { + // Read values with 'bits' bits each, then subtract 1 + for (int j = from; j < to; j++) { + int value = ToShort(bits) - 1; + pixels[j] = (short) Math.max(0, value); + } + max = 0xFFFF / 2; + } else { + + int blockSize = to - from; + int base = max + 1; + long numBits = (long)Math.ceil(effectiveBlock * Math.log(base) / Math.log(2)); + BigInteger val = ToBigInteger(numBits); + BigInteger bigBase = BigInteger.valueOf(base); + for (int i = 0; i < blockSize; i++) { + int decodedVal = val.mod(bigBase).intValue() - 1; + pixels[from + i] = (short) Math.max(0, decodedVal); + val = val.divide(bigBase); + } + } + + from += effectiveBlock; + + if (from < imageSize) { + int masked = ToShort(1); + if (masked == 0) { + max = 0; + bits = 0; + break; + } + } + } + } + + + private long integerPower(int base, int exp) { + long result = 1; + long LongBase = base; + while (exp > 0) { + if ((exp & 1) == 1) result *= LongBase; + LongBase *= LongBase; + exp >>= 1; + } + return result; + } + + private int mostSignificantBit(long val) { + int r = 0; + while (val > 0) { + val >>= 1; + r++; + } + return r; + } + } diff --git a/lib/trpx2.0/LICENSE b/lib/trpx2.0/LICENSE new file mode 100644 index 00000000..1ca1f1ac --- /dev/null +++ b/lib/trpx2.0/LICENSE @@ -0,0 +1,21 @@ +MIT License + +Copyright (c) 2026 Senik Matinyan + +Permission is hereby granted, free of charge, to any person obtaining a copy +of this software and associated documentation files (the "Software"), to deal +in the Software without restriction, including without limitation the rights +to use, copy, modify, merge, publish, distribute, sublicense, and/or sell +copies of the Software, and to permit persons to whom the Software is +furnished to do so, subject to the following conditions: + +The above copyright notice and this permission notice shall be included in all +copies or substantial portions of the Software. + +THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR +IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, +FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE +AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER +LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, +OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE +SOFTWARE. diff --git a/lib/trpx2.0/README.md b/lib/trpx2.0/README.md new file mode 100644 index 00000000..03154835 --- /dev/null +++ b/lib/trpx2.0/README.md @@ -0,0 +1,245 @@ +# TRPX V2.0 data compression algorithm + +This repository provides two related components: + +- HDF5 **terse** filter – an HDF5 plugin that compresses diffraction data using the next-generation TERSE/PROLIX (TRPX) algorithm. +- **Pyterse** – a Python package that provides direct Python bindings (via pybind11) to the TERSE/PROLIX library. + +Both are powered by the same TERSE/PROLIX algorithm. +You can build only the HDF5 filter, only the Pyterse Python package, or both. + +# Pyterse python package + + +The pyterse python package provides Python bindings for the C++ TERSE/PROLIX(TRPX) compression algorithm scheme ([https://github.com/senikm/trpx](https://github.com/senikm/trpx)). This version of TRPX is parallel, has better memory management, and compression modes. + +## Prerequisites + +Before using pyterse, ensure your data meets these requirements: +- Signed or unsigned integral type data +- Grayscale data +- Preferably has high dynamic range + +## Installation + +Create a virtual environment: +```bash +conda create -n pyterse python pip numpy pillow +conda activate pyterse +``` + +Install the package: +```python +pip install pyterse +``` + +## Usage guide + +### Basic operations + +#### Creating a Terse object + +There are multiple ways to create a Terse object: + +```python +import pyterse + +# Empty constructor +terse = pyterse.Terse() + +# From NumPy array +terse = pyterse.Terse(data) # data can be nD NumPy array or slice + +# With custom compression mode +terse = pyterse.Terse(data, pyterse.TerseMode.SIGNED) # Available modes: SIGNED, UNSIGNED, SMALL_UNSIGNED, DEFAULT +``` + +#### Inserting and managing data + +Add data to an existing Terse object: +```python +# Append data at the end +terse.push_back(data) # Data must match existing shape (terse.dim()) + +# Insert at specific position +terse.insert(pos, data) # pos is the frame index +``` + +#### File operations + +Save and load compressed data: +```python +# Save to file +terse.save('filename.trpx') + +# Load from file +loaded_terse = pyterse.Terse.load('filename.trpx') +``` + +#### Data Decompression + +Decompress data: +```python +# Decompress all data +decompressed_data = terse.prolix() + +# Decompress specific frame +frame = terse.at(0) +decompressed_frame = frame.prolix() +``` + +#### Metadata Management + +```python +# Set metadata for a frame +terse.set_metadata(frame, "metadata string") + +# Get metadata from a frame +metadata = terse.metadata(frame) +``` + +#### Utility Methods + +```python +# Data information +terse.dim() # Dimensions of one frame +terse.size # Number of elements per frame +terse.number_of_frames # Number of frames +terse.number_of_bytes # Size in bytes of compressed data +terse.bits_per_val # Bits used per value +terse.is_signed # Whether data is signed + +# Data management +terse.erase(pos) # Remove frame at position +terse.shrink_to_fit() # Optimize memory usage + +# Compression settings +terse.set_block_size(size) # Set compression block size (before adding frames) +terse.set_fast(bool) # Toggle fast compression mode +terse.set_small(bool) # Toggle small data optimization +terse.dop() # Degree of parallelism +terse.set_dop(value) # Set degree of parallelism (0.0 to 1.0) + +# Float data compression +terse.set_fractional_precision (value) +terse.fractional_precision() +``` + +## Building from source + +```bash +conda create -n pyterse python pip numpy pillow +conda activate pyterse +git clone https://github.com/senikm/trpx2.0.git +cd trpx2.0 +mkdir build +cd build +cmake .. -DCMAKE_PREFIX_PATH=$CONDA_PREFIX +make +python3 py_tests/pyterse_test.py # Will run unittests +``` + +# HDF5 tersecodec + +## Prerequisites + +Before proceeding, ensure you have the following installed on your system: + +1. **Create a new conda environment** + +```bash +conda create -n hdf5_terse_env python=3.10 +conda activate hdf5_terse_env +``` +2. **Install CMake, HDF5 and other dependencies** + +```bash +conda install cmake hdf5 numpy h5py hdf5plugin pillow +cmake --version +``` + +3. **Clone the repository and build the filter** + +``` +git clone https://github.com/senikm/trpx2.0.git +cd trpx2.0 +mkdir build +cd build +cmake .. -DCMAKE_PREFIX_PATH=$CONDA_PREFIX +make +``` +Note: The CMakeLists.txt provided uses the CONDA_PREFIX environment variable to locate HDF5. Ensure the Conda environment is active when running CMake. + +4. **Running performance tests** + +> Run the test.py in py_tests folder. +``` +python3 py_tests/filter_test.py +``` +Note: The conda environment is activated and by default the script loads the filter from `build/src/terse.so`. + +> Sample output + +```sql +HDF5_PLUGIN_PATH: None +h5py version: 3.12.1 +HDF5 version: 1.14.3 +TERSE filter registered successfully with ID 32029 +TERSE filter is available. +Original data shape: (30, 512, 512) +Original data size: 30720.00 KB +Compressed HDF5 file size: 1044.42 KB +Data successfully compressed and decompressed with TERSE filter. +``` + +5. **You can make the filter path permanent** + +``` +echo 'export HDF5_PLUGIN_PATH=/path/to/hdf5_filter/build/src' >> ~/.bashrc +source ~/.bashrc +``` + +# Fiji/ImageJ plugin for .trpx format files + + +> For compilation, use the Java version that came with Fiji, to ensure Java compatibility. Also make sure the ij-1.??.jar package is included in the compilation: +> For example, compile with: + +```bash + Applications/Fiji.app/java/macosx/zulu8.60.0.21-ca-fx-jdk8.0.322-macosx_x64/jre/Contents/Home/bin/javac -cp /Applications/Fiji.app/jars/ij-1.53t.jar TRPX_Reader.java +``` +> Then create the .jar files with: +```bash + Applications/Fiji.app/java/macosx/zulu8.60.0.21-ca-fx-jdk8.0.322-macosx_x64/jre/Contents/Home/bin/jar -cvf Terse_Reader.jar TRPX_Reader*.class +``` +> Then copy Terse_Reader.jar to the "plugins" directory of Fiji: +```bash + cp TRPX_Reader.jar /Applications/Fiji.app/plugins/. +``` +> Then restart Fiji, and Terse Reader is in the plugins menu. + +# Reading TRPX files with DIALS + +> Please, install Pyterse package with the following command before reading single TRPX files in DIALS: + +```bash +dials.python -m pip install pyterse +``` + +## Citation + +**If you use our software in your research, please cite our paper using the following BibTeX entry**: + +```bibtex +@article{MATINYAN2026114298, + title = {TRPXv2.0: superfast, parallel compression of diffraction patterns and images, with native Python and HDF5 support}, + journal = {Ultramicroscopy}, + volume = {281}, + pages = {114298}, + year = {2026}, + issn = {0304-3991}, + doi = {https://doi.org/10.1016/j.ultramic.2025.114298}, + url = {https://www.sciencedirect.com/science/article/pii/S0304399125001962}, + author = {S. Matinyan and P. Filipcik and D.G. Waterman and C.D. Owen and J.P. Abrahams} +} +``` \ No newline at end of file diff --git a/lib/trpx2.0/build.sh b/lib/trpx2.0/build.sh new file mode 100755 index 00000000..9038b42c --- /dev/null +++ b/lib/trpx2.0/build.sh @@ -0,0 +1,49 @@ +#!/bin/bash + + +set -e + + +set -x + + +echo "Creating/cleaning build directory..." +rm -rf build +mkdir -p build + +echo "Entering build directory..." +cd build + + +if [ -n "$CONDA_PREFIX" ]; then + echo "Conda environment detected at: $CONDA_PREFIX" + PYTHON_PATH="$CONDA_PREFIX/bin/python" +else + echo "No Conda environment detected, using system Python" + PYTHON_PATH=$(which python3) +fi + + +echo "Configuring project with CMake..." +cmake .. \ + -DCMAKE_BUILD_TYPE=Release \ + -DCMAKE_CXX_FLAGS="-O3" \ + -DPython_EXECUTABLE="$PYTHON_PATH" \ + -DPYTHON_EXECUTABLE="$PYTHON_PATH" \ + -DPython_ROOT_DIR="$CONDA_PREFIX" \ + -DCMAKE_PREFIX_PATH="$CONDA_PREFIX" + + +echo "Building project..." +make VERBOSE=1 + + +cd .. + + +echo "Running tests..." +python3 py_tests/pyterse_test.py -v +python3 py_tests/filter_test.py -v + + +echo "Build and test process completed!" \ No newline at end of file diff --git a/lib/trpx2.0/include/Bitqueue.hpp b/lib/trpx2.0/include/Bitqueue.hpp new file mode 100644 index 00000000..162fa1ae --- /dev/null +++ b/lib/trpx2.0/include/Bitqueue.hpp @@ -0,0 +1,408 @@ +// +// Bitqueue.h +// Bitqueue +// +// Created by Jan Pieter Abrahams on 14.12.2024. +// + +#ifndef Bitqueue_h +#define Bitqueue_h + +#include "Shared_array.hpp" +#include +#include +#include +#include +#include + +namespace jpa { + +/** + * @brief Concept to define a Bitqueue-like container. + * + * A container that: + * - Provides sequential access to its data. + * - Has a `data()` member returning a pointer to the underlying memory (`std::uint8_t*`). + * - Offers `begin()` and `end()` methods for iterating over the memory sequentially. + * - Ensures compatibility between the raw pointer returned by `data()` and the iterator returned by `begin()`. + * + * This concept can be used to constrain generic code to only accept containers that satisfy the + * required properties of a Bitqueue. + * + * **Requirements:** + * - `data()` must return a `std::uint8_t*`. + * - `begin()` and `end()` must return iterators, and `*begin()` must yield a `std::uint8_t&`. + * - The iterator returned by `begin()` must satisfy the `std::contiguous_iterator` concept. + * - `data()` and `std::to_address(begin())` must return the same address. + * + * This ensures that the container stores its values sequentially in memory and allows efficient access. + */ +template +concept Bitqueue = requires(T obj) { + { obj.data() } -> std::same_as; + { obj.begin() } -> std::forward_iterator; + { obj.end() } -> std::forward_iterator; + { *obj.begin() } -> std::same_as; + requires std::contiguous_iterator; + requires std::is_same_v; +}; + +/** + * @brief A utility class for inserting integal values with bit widths up to 64 bits into a Bitqueue-like container. + * + * This class writes data to a container that satisfies the `Bitqueue` concept, treating the container + * as a contiguous stream of bits. It supports pushing single or multiple values of arbitrary bit widths. + * + * **Usage:** + * - Instantiate the class with a `Bitqueue` container. + * - Use the `push_back` methods to insert values of specific bit widths into the buffer. + * + * **Features:** + * - Supports fixed-width pushes at compile-time (`push_back`). + * - Supports dynamic-width pushes at runtime (`push_back(std::size_t B)`). + * - Handles both signed and unsigned integral types. + * - Ensures that remaining buffered bits are flushed upon destruction. + * + * Example: + * @code + * Bitqueue_push_back pusher(bitqueue); + * pusher.push_back<5>(value); // Push 5 bits of value into the buffer + * @endcode + * + */ +class Bitqueue_push_back { +public: + + /** + * @brief Constructor that initializes the Bitqueue_push_back instance with a writable data buffer. + * @param container Bitqueue container that stores the data buffer to write to. + */ + template + Bitqueue_push_back(CONTAINER& container) noexcept : d_data(container.data()) { + assert(&*container.begin() == d_data); + assert((container.size() * sizeof(*d_data)) % sizeof(std::size_t) == 0 && "Container size must be a multiple of d_buffer size"); + assert(reinterpret_cast(d_data) % alignof(std::size_t) == 0 && "Container data must be aligned to the size of d_buffer"); + } + + /** + * @brief Destructor that flushes the remaining bits to the data buffer. + */ + ~Bitqueue_push_back() noexcept { flush_buffer(); } + + /** + * @brief Relocates the writable data buffer to a new position. + * @param data Pointer to the new writable data buffer. + */ + constexpr void relocate(std::uint8_t* const data) noexcept { d_data = data; } + + /** + * @brief flushes any buffered bits to memory + */ + void flush_buffer() noexcept { std::memcpy(d_data, &d_buffer, sizeof(decltype(d_buffer))); } + + /** + * @brief Retrieves the current writable data pointer. + * @return Pointer to the current writable data buffer. + */ + constexpr std::uint8_t* data() const noexcept { return d_data; } + + /** + * @brief Pushes a single value of a compile-time fixed bit width into the buffer. + * @tparam B Number of bits to push. + * @tparam T Integral type of the value to push. + * @param val The value to push into the buffer. + */ + template + constexpr void push_back(T const val) noexcept { + if constexpr (B != 0) + push_back(std::span(&val, 1)); + } + + /** + * @brief Pushes a single value of a specified bit width into the buffer. + * @tparam T Integral type of the value to push. + * @param B Number of bits to push. + * @param val The value to push into the buffer. + */ + template + constexpr void push_back(std::uint8_t const B, T const val) { + if (B != 0) + push_back(B, std::span(&val, 1)); + } + + /** + * @brief Pushes multiple values of a compile-time fixed bit width into the buffer. + * @tparam B Number of bits to push for each value. + * @tparam T Integral type of the values to push. + * @tparam N Size of the span containing the values. + * @param values Span containing the values to push into the buffer. + */ + template + constexpr void push_back(std::span const values) noexcept { + static_assert(sizeof(T) <= sizeof(std::size_t), "Type too large"); + if constexpr (B != 0) { + decltype(d_buffer) const mask = ((1ul << B) - 1); + for (auto const& val : values) { + if constexpr (std::is_signed_v) + d_buffer |= static_cast(static_cast>(val) & mask) << d_buffered_bits; + else + d_buffer |= static_cast(val) << d_buffered_bits; + d_buffered_bits += B; + if (d_buffered_bits >= sizeof(decltype(d_buffer)) * 8) { + std::memcpy(d_data, &d_buffer, sizeof(decltype(d_buffer))); + d_data += sizeof(decltype(d_buffer)); + d_buffered_bits -= sizeof(decltype(d_buffer)) * 8; + if constexpr (std::is_signed_v) + d_buffer = static_cast(static_cast>(val) & mask) >> (B - d_buffered_bits); + else + d_buffer = static_cast(val) >> (B - d_buffered_bits); + } + } + } + } + + /** + * @brief Pushes multiple values of a specified bit width into the buffer. + * @tparam T Integral type of the values to push. + * @tparam N Size of the span containing the values. + * @param B Number of bits to push for each value. + * @param values Span containing the values to push into the buffer. + */ + template + constexpr void push_back(std::uint8_t const B, std::span const values) noexcept { + static_assert(sizeof(T) <= sizeof(std::size_t), "Type too large"); + assert(B <= sizeof(T) * 8); + if (B != 0) { + decltype(d_buffer) const mask = ((1ul << B) - 1); + for (auto const& val : values) { + if constexpr (std::is_signed_v) + d_buffer |= static_cast(static_cast>(val) & mask) << d_buffered_bits; + else + d_buffer |= static_cast(val) << d_buffered_bits; + d_buffered_bits += B; + if (d_buffered_bits >= sizeof(decltype(d_buffer)) * 8) { + std::memcpy(d_data, &d_buffer, sizeof(decltype(d_buffer))); + d_data += sizeof(decltype(d_buffer)); + d_buffered_bits -= sizeof(decltype(d_buffer)) * 8; + if constexpr (std::is_signed_v) + d_buffer = static_cast(static_cast>(val) & mask) >> (B - d_buffered_bits); + else + d_buffer = static_cast(val) >> (B - d_buffered_bits); + } + } + } + } + + /** + * @brief Computes the distance in bytes between the current buffer position and the given pointer. + * @param p Pointer to compare the current position against. + * @return Distance in bytes between the current position and the pointer. + */ + constexpr std::ptrdiff_t operator-(std::uint8_t const* p) const noexcept { return d_data + (d_buffered_bits + 7) / 8 - p; } + +private: + std::uint8_t *d_data; + std::uint8_t d_buffered_bits = 0; + std::uint64_t d_buffer = 0; +}; + +/** + * @brief A utility class for extracting integal values with bit widths up to 64 bits from a Bitqueue-like container. + * + * This class reads data from a container that satisfies the `Bitqueue` concept, interpreting the data as + * a contiguous stream of bits. It supports popping single or multiple values of arbitrary bit widths. + * + * **Usage:** + * - Instantiate the class with a `Bitqueue` container. + * - Use the `pop` methods to extract values of specific bit widths from the buffer. + * - Use the `skip` methods to skip over bits without extracting them. + * + * **Features:** + * - Supports fixed-width pops at compile-time (`pop`). + * - Supports dynamic-width pops at runtime (`pop(std::size_t B)`). + * - Handles both signed and unsigned integral types. + * + * Example: + * @code + * Bitqueue_pop popper(bitqueue); + * auto value = popper.pop<5, std::uint8_t>(); // Extract 5 bits into a uint8_t + * @endcode + */ +class Bitqueue_pop { +public: + /** + * @brief Constructor that initializes the Bitqueue_pop instance with data to process. + * @param container Bitqueue container that stores the data buffer to read from. + */ + template + Bitqueue_pop(CONTAINER const& container) noexcept : d_data(container.data()) { + assert(&*container.begin() == d_data); + assert((container.size() * sizeof(*d_data)) % sizeof(std::size_t) == 0 && "Container size must be a multiple of d_buffer size"); + assert(reinterpret_cast(d_data) % alignof(std::size_t) == 0 && "Container data must be aligned to the size of d_buffer"); + std::memcpy(&d_buffer, d_data, std::min(container.size() * sizeof(*d_data), sizeof(std::size_t))); + } + /** + * @brief Retrieves the current readable data pointer. + * @return Pointer to the current readable data buffer. + */ + constexpr std::uint8_t const* data() const noexcept { return d_data; } + + /** + * @brief Pops a single value of specified bit width from the buffer. + * @tparam B Number of bits to pop. + * @tparam T Integral type to store the popped value. + * @return The value popped from the buffer. + */ + template + constexpr T pop() noexcept { + static_assert(sizeof(T) <= sizeof(decltype(d_buffer)), "Type too large"); + if (B == 0) + return 0; + else { + T val; + pop(std::span(&val, 1)); + return val; + } + } + + /** + * @brief Pops a single value with a specified bit width from the buffer. + * @tparam T Integral type to store the popped value. + * @param B Number of bits to pop. + * @return The value popped from the buffer. + */ + template + constexpr T pop(std::uint8_t B) noexcept { + static_assert(sizeof(T) <= sizeof(std::size_t), "Type too large"); + if (B == 0) + return 0; + else { + T val; + pop(B, std::span(&val, 1)); + return val; + } + } + + /** + * @brief Pops multiple values of a specified bit width from the buffer. + * @tparam T Integral type to store the popped values. + * @tparam N Size of the span for storing the values. + * @param B Number of bits to pop for each value. + * @param values Span to store the popped values. + */ + template + constexpr void pop(std::uint8_t const B, std::span values) noexcept { + assert(B <= 8 * sizeof(T)); // More bits requested than depth of T: B too large + if (B==0) + std::fill(values.begin(), values.end(), 0); + else { + T const mask = static_cast((1ul << B) - 1); + T const is_negative = std::is_unsigned_v ? T(0) : T(T(1) << (B - 1)); + for (auto& p : values) { + p = static_cast(d_buffer); + if (B <= d_buffered_bits){ + if (B == 8 * sizeof(decltype(d_buffer))) + d_buffer = d_buffered_bits = 0; + else { + d_buffer >>= B; + d_buffered_bits -= B; + } + } + else { + d_data += sizeof(decltype(d_buffer)); + std::memcpy(&d_buffer, d_data, sizeof(decltype(d_buffer))); + p |= static_cast(d_buffer << d_buffered_bits); + d_buffer >>= B - d_buffered_bits; + d_buffered_bits = sizeof(decltype(d_buffer)) * 8 + d_buffered_bits - B; + } + if (B != 8 * sizeof(T)) { + if (p & is_negative) p |= ~mask; + else p &= mask; + } + } + } + } + + /** + * @brief Pops multiple values of a compile-time fixed bit width from the buffer. + * @tparam B Number of bits to pop. + * @tparam T Integral type to store the popped values. + * @tparam N Size of the span for storing the values. + * @param values Span to store the popped values. + */ + template requires (!std::is_same_v) + constexpr void pop(std::span values) noexcept { + assert(B <= 8 * sizeof(T)); // More bits requested than depth of T: B too large + if constexpr (B==0) + std::fill(values.begin(), values.end(), 0); + else { + for (auto& p : values) { + p = static_cast(d_buffer); + if (B <= d_buffered_bits){ + if constexpr (B == 8 * sizeof(decltype(d_buffer))) + d_buffered_bits = d_buffer = 0; + else { + d_buffer >>= B; + d_buffered_bits -= B; + } + } + else { + d_data += sizeof(decltype(d_buffer)); + std::memcpy(&d_buffer, d_data, sizeof(decltype(d_buffer))); + p |= static_cast(d_buffer << d_buffered_bits); + d_buffer >>= B - d_buffered_bits; + d_buffered_bits = sizeof(decltype(d_buffer)) * 8 + d_buffered_bits - B; + } + if constexpr (B != 8 * sizeof(T)) { + constexpr T is_negative = std::is_unsigned_v ? T(0) : T(T(1) << (B - 1)); + constexpr auto mask = static_cast>((1ul << B) - 1); + if (p & is_negative) p |= static_cast(~mask); + else p &= static_cast(mask); + } + } + } + } + + /** + * @brief Skips a fixed number of bits in the buffer. + * @tparam B Number of bits to skip. + */ + template + constexpr void skip() { + if (B <= d_buffered_bits) + d_buffered_bits -= B; + else { + d_data += sizeof(std::size_t) + B / (sizeof(std::size_t) * 8) ; + d_buffered_bits = (sizeof(std::size_t) * 8) - (B % (sizeof(std::size_t) * 8)); + } + } + + /** + * @brief Skips a specified number of bits in the buffer. + * @param B Number of bits to skip. + */ + constexpr void skip(std::size_t const B) { + if (B <= d_buffered_bits) + d_buffered_bits -= B; + else { + d_data += sizeof(std::size_t) + B / (sizeof(std::size_t) * 8) ; + d_buffered_bits = (sizeof(std::size_t) * 8) - (B % (sizeof(std::size_t) * 8)); + } + } + + /** + * @brief Computes the distance in bytes between the current buffer position and the given pointer. + * @param p Pointer to compare the current position against. + * @return Distance in bytes between the current position and the pointer. + */ + constexpr std::ptrdiff_t operator-(std::uint8_t const* p) const noexcept { return d_data + (d_buffered_bits + 7) / 8 - p; } + +private: + std::uint8_t const* d_data; + std::size_t d_buffer; + std::uint8_t d_buffered_bits = sizeof(decltype(d_buffer)) * 8; +}; + +}// end namespace jpa + +#endif /* Bitqueue_h */ diff --git a/lib/trpx2.0/include/Concurrent.hpp b/lib/trpx2.0/include/Concurrent.hpp new file mode 100644 index 00000000..ea6277cc --- /dev/null +++ b/lib/trpx2.0/include/Concurrent.hpp @@ -0,0 +1,412 @@ +#ifndef Concurrent_h +#define Concurrent_h + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +namespace jpa { +/** + * @brief Class to control and represent the degree of parallelism. + * + * The `Deg_of_parallelism` class allows you to control the degree of parallelism + * for multi-threaded applications. It dynamically determines how many + * cores should be available for execution and ensures that parallelism is within + * a defined range. + * + * This class is used both in concurrent programming (managing different tasks simultaneously) + * and in parallel processing (dividing a task into subtasks processed concurrently). + */ +class Deg_of_parallelism { +public: + /** + * @brief Default constructor initializes the degree of parallelism to 0.0 (sequential execution). + */ + constexpr Deg_of_parallelism() : d_dop(0) {} + + /** + * @brief Constructor that initializes the degree of parallelism to a given value. + * + * @param val A value between 0 and 1, where 0 implies sequential execution and + * 1 implies full use of available cores. Values outside this range will be clamped. + */ + constexpr Deg_of_parallelism(double const val) : d_dop(std::clamp(val, 0.0, 1.0)) {} + + /** + * @brief Cast operator to double, allowing implicit conversion to a degree of parallelism value. + * + * @return The degree of parallelism as a double between 0 and 1. + */ + constexpr operator double() const { return d_dop; } + + /** + * @brief Set the degree of parallelism to a new value. + * + * @param val A value between 0 and 1, where 0 implies sequential execution and + * 1 implies full use of available cores. Values outside this range will be clamped. + * + * @return A reference to the updated `Deg_of_parallelism` object. + */ + constexpr Deg_of_parallelism& operator=(double const val) { + d_dop = std::clamp(val, 0.0, 1.0); + return *this; + } + + /** + * @brief Set the degree of parallelism based on a given number of cores. + * + * @param val The desired number of cores. It will be clamped to the range of 1 to the total number of available cores. + * + * @return A reference to the updated `Deg_of_parallelism` object. + */ + Deg_of_parallelism& cores(unsigned const val) { + d_dop = std::clamp(val, 1u, d_cores) / double(d_cores); + return *this; + } + + /** + * @brief Get the number of cores currently assigned based on the degree of parallelism. + * + * @return The number of cores available for concurrent processing, based on the current degree of parallelism. + */ + unsigned cores() const { + return std::clamp(unsigned(d_dop * d_cores), 1u, d_cores); + } + +private: + double d_dop; ///< Degree of parallelism, a value between 0 (sequential) and 1 (full parallelism). + static inline unsigned const d_cores = std::max(1u, std::thread::hardware_concurrency()); ///< Total number of available cores on the system. +}; + +/** + * @brief Default degree of parallelism is sequential execution. + */ +static Deg_of_parallelism degree_of_parallelism_default = Deg_of_parallelism(); + +/** + * @brief A concurrent task manager for dispatching and executing different tasks using a thread_local thread pool. + * + * The `Concurrent` class is designed for managing the concurrent execution of multiple, potentially different tasks. + * It allows you to run tasks asynchronously, facilitating concurrent programming where tasks can run simultaneously. + * This is distinct from parallel processing of a single task divided into subtasks (handled by separate algorithms). + * + * The class manages a pool of threads shared by all instances of `Concurrent` to optimize resource use. + * You can run background tasks using `std::future` for standard use cases. + * It also accepts **functors** for advanced control, such as monitoring, pausing, or stopping tasks while they run. + * + * ### Example Usage: + * + * **Standard usage:** + * @code + * jpa::Concurrent tasks; + * auto result1 = tasks.background(my_function1, arg1, arg2); + * auto result2 = tasks.background(my_function2, arg3, arg4); + * // Do other work here... + * int final_result = result1.get() + result2.get(); // Blocks execution until the tasks are finished + * @endcode + * This approach allows for **fire-and-forget** or **fire-and-expect** behavior, where you dispatch tasks and retrieve + * their results later. + * + * ### Functor-based Task Control: + * + * You can pass a **functor** to the `background()` function, where the functor has state variables + * like `progress`, `kill`, etc. This allows you to monitor or control the task as it runs. + * + * **Example:** + * @code + * // Define a functor that tracks progress and can be killed + * struct TaskControl { + * int progress = 0; + * bool kill = false; + * + * // Define the task (operator()) to be run in the background + * void operator()(int duration) { + * for (int i = 0; i <= 100; ++i) { + * if (kill) { + * std::cout << "Task killed at " << progress << "%!" << std::endl; + * return; + * } + * progress = i; + * std::this_thread::sleep_for(std::chrono::milliseconds(duration)); + * } + * std::cout << "Task completed!" << std::endl; + * } + * }; + * + * // Instantiate the functor and start the task in the background + * jpa::Concurrent tasks; + * TaskControl task; + * auto result = tasks.background(task, 50); + * + * // Monitor the task's progress and kill it when necessary + * while (task.progress < 100) { + * std::cout << "Progress: " << task.progress << "%" << std::endl; + * std::this_thread::sleep_for(std::chrono::milliseconds(500)); + * + * if (task.progress > 50) { + * task.kill = true; // Stop the task + * } + * } + * + * result.get(); // Blocks until the task completes or is killed + * @endcode + * + * ### Note on Concurrency vs. Parallelism: + * + * The `Concurrent` class focuses on concurrent programming—running different tasks simultaneously. + * For parallel processing of a single task (like processing elements of a container in parallel), + * refer to the algorithms provided in `Parallel.hpp`. + * + * ### Note on Nested Concurrent Calls: + * + * The Concurrent class is designed to handle nested calls to ensure that concurrency occurs at the highest + * abstraction level while avoiding deadlocks. If multiple instances of Concurrent are created (for example, + * within different functions), only the first instance created will manage all its tasks asynchronously via the thread + * pool. All subsequent Concurrent instances that are created during the lifetime of the first instance will + * only run asynchronously if there are free threads in the thread pool. If there are no free threads, they will + * execute their tasks sequentially in the current thread. This ensures that tasks within tasks are executed, even + * if there are top-level tasks in the queue. Otherwise a deadlock could result, when top level tasks are queued + * with a higher priority than sub-tasks, yet requiring the sub-tasks to finish first. + * + * This design allows for efficient task management, where lower-level tasks run sequentially when dispatched + * by higher-level concurrent operations, thus preventing potential deadlocks caused by excessive nesting and + * ensuring that concurrency is handled where it’s most needed. + */ +class Concurrent { +public: + /** + * @brief Constructs the task manager, specifying the degree of parallelism (dop). + * + * @param dop The degree of parallelism, representing the proportion of available cores to use. + * It is a number between 0 and 1, where 0 implies sequential execution and + * 1 implies full use of available cores. + */ + Concurrent(Deg_of_parallelism dop) : d_dop(dop) { + if (d_dop > 0) + d_unique_id = d_thread_pool.f_get_unique_id(); + } + + /** + * @brief Backgrounds a task for asynchronous execution. + * + * This function creates a packaged task from the function and arguments, pushes it into the task queue, + * and manages the number of threads allocated for this task manager. + * + * @tparam Func The function type. + * @tparam Args The argument types for the function. + * @param func The function to be executed in the background. + * @param args The arguments to be passed to the function. + * @return A future object representing the result of the task. + */ + template + auto background(Func&& func, Args&&... args) -> std::future> { + using return_type = typename std::invoke_result_t; + auto task = std::make_shared> + ([captured_func = std::forward(func), ... captured_args = std::forward(args)]() mutable { + if constexpr (std::is_void_v) + captured_func(std::forward(captured_args)...); + else + return captured_func(std::forward(captured_args)...); + }); + std::future result = task->get_future(); + if (d_dop == 0.0 || !d_thread_pool.f_add_task(*this, [task]() { (*task)(); })) + (*task)(); // Execute the task immediately + return result; + } + + /** + * @brief Waits until all backgrounded tasks have been completed. + * + * This method blocks the caller until the task queue is empty and no active threads are running. + */ + void finish() { + while (running_tasks() > 0) + std::this_thread::yield(); // Yield control to other threads + } + + /** + * @brief Returns the number of tasks submitted by this instance of `Concurrent` that are either in the queue or being processed. + * + * @return The number of running tasks. + */ + unsigned running_tasks() const { + unsigned task_count = 0; + { + std::lock_guard lock(d_thread_pool.d_mutex); + for (const auto& task : d_thread_pool.d_task_queue) + if (std::get<2>(task) == d_unique_id) + ++task_count; + if (d_thread_pool.d_running_tasks.count(d_unique_id)) + task_count += d_thread_pool.d_running_tasks.at(d_unique_id); + } + return task_count; + } + + /** + * @brief Returns the degree of partiality of this instance of `Concurrent`. + * + * @return The Deg_of_parallelism: 0 for sequential processing, 1 for full paralllel processing. + */ + constexpr Deg_of_parallelism dop() const noexcept { return d_dop; } + + +private: + Deg_of_parallelism d_dop; + std::size_t d_unique_id = 0; + + inline static class c_Global_thread_pool { + friend class Concurrent; + + unsigned const d_max_threads = std::max(1u, std::thread::hardware_concurrency() - 1); + std::atomic d_unique_id_generator = 0; + std::mutex d_mutex; + std::list, Deg_of_parallelism, std::size_t>> d_task_queue; + std::unordered_map d_running_tasks; + std::vector d_workers; + std::condition_variable d_condition; + std::atomic d_stop {false}; + + c_Global_thread_pool() { + d_workers.reserve(d_max_threads); + for (unsigned i = 0; i < d_max_threads; ++i) + d_workers.emplace_back(&c_Global_thread_pool::f_worker_thread, this); + } + + ~c_Global_thread_pool() { + d_stop = true; + d_condition.notify_all(); + for (auto& worker : d_workers) + if (worker.joinable()) { + try { worker.join(); } + catch (const std::system_error& e) { std::cerr << "Error joining thread: " << e.what() << std::endl; } + catch (...) { std::cerr << "Unknown error occurred while joining thread." << std::endl; } + } + } + + bool f_add_task(Concurrent const& concurrent_instance, std::function const& task) { + if (d_max_threads == 0) + return false; + auto top_level_id = f_get_top_level_id(); + if (top_level_id.has_value() && + concurrent_instance.d_unique_id != *top_level_id && + d_running_tasks.size() == d_max_threads) + return false; + { + std::lock_guard lock(d_mutex); + d_task_queue.emplace_back(std::move(task), concurrent_instance.d_dop, concurrent_instance.d_unique_id); + } + d_condition.notify_one(); + return true; + } + + std::size_t f_get_unique_id() { return d_unique_id_generator++; } + + std::optional f_get_top_level_id() { + std::lock_guard lock(d_mutex); + std::optional top_level_id; + for (const auto& [id, count] : d_running_tasks) + if (!top_level_id || id < *top_level_id) + top_level_id = id; + for (const auto& [task, dop, id] : d_task_queue) + if (!top_level_id || id < *top_level_id) + top_level_id = id; + return top_level_id; + } + + void f_worker_thread() { + while (true) { + std::function task; + std::size_t concurrent_id; + { + std::unique_lock lock(d_mutex); + d_condition.wait(lock, [this]() { return d_stop || !d_task_queue.empty(); }); + if (d_task_queue.empty() && d_stop) + return; + auto it = d_task_queue.begin(); + while (it != d_task_queue.end()) { + Deg_of_parallelism const& dop = get<1>(*it); + concurrent_id = get<2>(*it); + if (d_running_tasks[concurrent_id] < dop.cores()) { + task = std::move(get<0>(*it)); + d_task_queue.erase(it); + ++d_running_tasks[concurrent_id]; + break; + } + ++it; + } + if (!task) + continue; + } + try { task(); } + catch (const std::exception& e) { std::cerr << "Exception in worker thread: " << e.what() << std::endl; } + catch (...) { std::cerr << "Unknown exception in worker thread." << std::endl; } + { + std::lock_guard lock(d_mutex); + if (d_running_tasks[concurrent_id] == 1) + d_running_tasks.erase(concurrent_id); + else + --d_running_tasks[concurrent_id]; + d_condition.notify_all(); + } + } + } + } d_thread_pool; +}; + +/** + * @brief Backgrounds a task for asynchronous execution, returning a future representing the result. + * + * This function provides a simple and intuitive interface for executing tasks in the background. + * It works similarly to `std::async()`, allowing you to run tasks asynchronously in a **global thread pool**, + * ensuring efficient resource management without excessive thread creation. + * + * **Example:** + * @code + * auto result = jpa::background(my_function, arg1, arg2); + * // Do other stuff + * int final_result = result.get(); // Blocks until my_function completes + * @endcode + * + * ### When to Use `background()` vs `std::async`: + * + * - **Use `background()`** for tasks that are **CPU-bound** and need efficient concurrent execution, especially when + * you need to execute multiple different tasks concurrently using a centralized thread pool. It prevents the creation of + * excessive threads and manages resources effectively. + * + * - **Use `std::async`** for tasks that are expected to spend most of their time waiting, such as: + * - Network requests (waiting for responses). + * - File I/O (waiting for disk read/write operations). + * - Database queries (waiting for results). + * + * Also, use `std::async` if you need explicit control over task execution policies, such as choosing between launching + * a new thread or using deferred execution. `std::async` is appropriate for lightweight or occasional tasks where + * thread management is not a primary concern. + * + * ### Note on Concurrency vs. Parallelism: + * + * The `background()` function and the `Concurrent` class are designed for concurrent programming—running multiple + * different tasks simultaneously. For parallel processing of a single task (e.g., operating on elements of a container concurrently), + * refer to the parallel algorithms provided in a separate header (e.g., `Parallel.hpp`). + * + * @tparam Func The function type. + * @tparam Args The argument types for the function. + * @param func The function to be executed in the background. + * @param args The arguments to be passed to the function. + * @return A future object representing the result of the background task. + */ +template +auto background(Func&& func, Args&&... args) -> std::future> { + Concurrent task(1); + return task.background(std::forward(func), std::forward(args)...); +} +} // namespace jpa + +#endif /* Concurrent_h */ diff --git a/lib/trpx2.0/include/Grey_tif.hpp b/lib/trpx2.0/include/Grey_tif.hpp new file mode 100644 index 00000000..39d9a9df --- /dev/null +++ b/lib/trpx2.0/include/Grey_tif.hpp @@ -0,0 +1,882 @@ +// +// Grey_tif.h +// Grey_tif +// +// Created by Jan Pieter Abrahams on 23.05.23. +// + +#ifndef Grey_tif_h +#define Grey_tif_h + + +#include +#include +#include +#include +#include +#include +#include +#include + + +// DISCLAIMER: This is not a general-purpose TIFF library! + +namespace jpa { +/** + * @brief Grey_tif: A specialized class for two-dimensional scientific greyscale TIFF files, images & stacks + * + * Key Features: + * - Reads TIFF files that contain a single image or a stack of images. + * - Seamless appending and extraction of image data from/to TIFF files, STL containers and other Grey_tif objects. + * - Compact API that is powerful, intuitive, flexible and easy-to-learn/use. + * - Image data of any image on the stack can be accessed through std::span, an STL compatible range view. + * - Full read and write access to all pixel values. + * - Multi-dimensional subscript operator (C++23 feature) for convenient indexing. + * - Compact code, optimized for high-speed performance and minimal memory foorprint. + * - Accommodating images with varying types and sizes. + * - Support for signed and unsigned integral values, and float and double values. + * - Flexible bit depths: 8, 16, 32, or (for double precision values) 64 bits. + * - Accomodates raw TIFF data with runtinme pixel type identification. + * - Allows regularizing pixel types to a compile-time determined type. + * - The first (or only) image managed by a Grey_tif object is directly available as a readonly std::span + * + * Limitations: Grey_tif does not support: + * - TIFF files with fragmented images. + * - TIFF files with single-bit (black & white) or color images. + * - TIFF files that contain compressed images. + * + * The Grey_tif API has 7 member functions: + * - A default constructor() and constructors using a TIFF stream or an STL container. + * - write(): write a TIFF file to a stream. + * - clear(): remove all images. + * - swap(): swap image contants with another Grey_tif object, reacsting the pixel types of both Grey_tif containers if required. + * - image_stack_size(): get the number of images managed by the Grey_tif object. + * - push_back(): append one or more images from an STL container or other Grey_tif object. + * - const_image(): returns a const Grey_tif_image from the tif stack for readonly access to images form a TIFF stack. + * - image(): returns a Grey_tif_image object for read/write access to images form a TIFF stack. + * + * The Grey_tif_image API is publicly derived from std::span, and has two additional member functions: + * - type(): returns a POD_type_traits object providing information about the size, signedness, and integral nature of the image's pixel type. + * - operator[int, int]: the multi-dimensional subscript operator (C++23 feature) is available if supported by the compiler. + * + * Example that tests all features of Grey_tif.hpp: + * + * @code{.cpp} + * #include + * #include + * #include "Grey_tif.hpp" + * + * int main(int argc, const char * argv[]) { + * using namespace jpa; + * + * // Generate a Grey_tif object from a container + * std::vector img0({42,1,2,3,4,5,6,7,8,9,10,11,12,13,14,15}); + * Grey_tif tif(img0, {4, 4}); + * + * // The pixel type is defaults to the type of the initializing container + * Grey_tif tif_uint16(img0, {4, 4}); + * assert (tif.type().is() && tif_uint16.type().is()); + * + * // Append a Grey_tif object using a TIFF file; pixel format of all images is int, regardless of + * // the encoding in the TIFF file + * std::ifstream istream("/Users/abrahams/Dropbox/develop/jpa/Grey_tif/Grey_tif/infile.tif", std::ios::binary); + * tif.push_back(Grey_tif(istream)); + * + * // Readonly access of pixels of the first (or only) image of a Grey_tif as an STL range view (a std::span) + * assert(tif[0] == 42); + * + * // Read/write access is achieved through the image() member function, for constant images, use const_image() + * tif.image(0)[0] = tif.image(1)[0] = 43; + * // tif.const_image(1) = 44; <--- does not compile + * assert (tif.const_image(1)[0] == 43 && tif.image(1)[0] == 43); + * + * + * // A Grey_tif object can be written to disk using its 'write(ostream)' member function. + * std::filesystem::path output = "/Users/abrahams/Dropbox/develop/jpa/Grey_tif/Grey_tif/outfile.tif"; + * std::ofstream ostream(output, std::ios::binary); + * tif.write(ostream); + * + * // Raw Grey_tif objects can storore images with different pixel types + * Grey_tif tiff_raw(img0, {4, 4}); + * tiff_raw.push_back({10,12}); // You can also pushback a zero-initialized image by only providing its dimensions + * + * // Access to pixel data of raw Grey_tif only possible with the correct template parameter + * // assert(tif[0] == 42 && tif.image(0)[0] == 42); <--- does not compile + * assert(tiff_raw.image(0)[0] == 42); + * assert((static_cast>(tiff_raw))[0] == 42); // alternatively, a typecast is possible. + * + * // Grey_tif objects can be swapped. If a raw Grey_tif is swapped with a typed Grey_tiff, type conversion + * // of the pixel values occurs + * Grey_tif tif_float; // construct an empty Grey_tif object with float pixel values + * tif_float.swap(tiff_raw); // swapping converts the pixel values of all raw images from tif_raw to float + * tif_float.swap(tiff_raw); //swapping back preserves the new pixel values as float + * assert(tiff_raw.image(0).size() == 16); + * + * // tif objects can also be cleared: + * assert(tif.image_stack_size() == 2); + * assert(tif.clear().image_stack_size() == 0); + * } + * @endcode + * + * @tparam T The default / assumed pixel data type (default: std::uint16_t). + */ +template class Grey_tif; + +/** + * @brief POD_type_traits represents traits information about Plain Old Data (POD) types. + * + * This struct provides information about the size, signedness, and integral nature of a type. + */ +struct POD_type_traits { + /** + * @brief The size in bytes of the type. + */ + std::size_t size; + + /** + * @brief Flag indicating whether the type is signed. + */ + bool is_signed; + + /** + * @brief Flag indicating whether the type is integral. + */ + bool is_integral; + + /** + * @brief Check if type T matches the specified criteria. + * + * This function checks if a given type T matches the size, signedness, and integral nature + * criteria specified by this POD_type_traits object. + * + * @tparam T The type to check. + * @return True if the type matches the criteria, false otherwise. + */ + template + bool is() const { return sizeof(T) == size && std::is_integral_v == is_integral && std::is_signed_v == is_signed; } + + /** + * @brief Compare this POD_type_traits object with another for equality. + * + * @param other The POD_type_traits object to compare with. + * @return True if the objects are equal, false otherwise. + */ + bool operator ==(POD_type_traits const& other) const noexcept { + return size == other.size && is_signed == other.is_signed && is_integral == other.is_integral; + } +}; + +/** + * @brief API to an image managed by a regularized Grey_tif object with a compile time defined pixel type. + * + * This class provides methods for manipulating and accessing pixel data of an image in a Grey_tif object. It is an STL compatible range view, and inherits and provides all the features of std::span. + * It has no public constructors and can only be created by a Grey_tif object to allow access to its TIFF image data. If the In case of raw TIFF data (stored in a Grey_tif object), the size() member function returns zero if the template parameter does not correspond to the type encoded in the raw TIFF data. + * + * @tparam T The type of pixel data in the image. + * @see std::span + * @see Grey_tif + */ +template +class Grey_tif_image : public std::span { + + template friend class Grey_tif; + template friend class Grey_tif_image; + +public: + +#ifdef __cpp_multidimensional_subscript + /** + * @brief Accessor for image data using multidimensional subscripting. + * + * Allows read-only access to individual pixel values within the TIFF images using multidimensional subscripting without bounds checking. + * + * @param i0 The first dimension index. + * @param i1 The second dimension index. + * @return A reference to the pixel value at the specified indices. + */ + constexpr auto const& operator [] (std::size_t i0, std::size_t i1) const noexcept { return (*this)[d_dim[0] * i0 + i1]; } + + /** + * @brief Accessor for image data using multidimensional subscripting. + * + * Allows read/write access to individual pixel values within the TIFF images using multidimensional subscripting without bounds checking. + * + * @param i0 The first dimension index. + * @param i1 The second dimension index. + * @return A reference to the pixel value at the specified indices. + */ + constexpr auto & operator [] (std::size_t i0, std::size_t i1) noexcept { return (*this)[d_dim[0] * i0 + i1]; } + + /** + * @brief Accessor for image data using one-dimensional subscripting. + * + * Allows read-only access to individual pixel values within the TIFF images using single dimenion subscripting without bounds checking. + * + * @param i0 The first dimension index. + * @return A reference to the pixel value at the specified indices. + */ + constexpr auto const& operator [] (std::size_t i0) const noexcept {return std::span::operator[](i0);} + + /** + * @brief Accessor for image data using one-dimensional subscripting. + * + * Allows read/write access to individual pixel values within the TIFF images using single dimenion subscripting without bounds checking. + * + * @param i0 The first dimension index. + * @return A reference to the pixel value at the specified indices. + */ + constexpr auto & operator [] (std::size_t i0) noexcept {return std::span::operator[](i0);} +#endif + + /** + * @brief Get the dimensions of the TIFF image. + * + * @return An array containing the width and height of the image. + */ + std::array const& dim() const noexcept { return d_dim; } + + /** + * @brief Get the pixel data type of the TIFF image. + * + * @return The pixel data type of the image as a POD_type_traits. + * @see POD_type_traits + */ + POD_type_traits const& type() const noexcept {return d_image_type;} + +private: + Grey_tif_image(POD_type_traits const& image_type, std::array const& dim, std::span const& data) : + std::span(data), + d_image_type(image_type), + d_dim(dim) {} + + POD_type_traits d_image_type; + std::array d_dim; +}; + + +/** + * @brief API to image parameters of a raw Grey_tif object with a runtime defined pixel type. + * + * This class allows determining the runtime pixel type and dimensions of an image managed by a raw Grey_tif object. + * But it does not allow access to the pixel values. In order to address its pixel values, it must first be cast to + * a Grey_tif_image of the correct pixel type T, which is the case when raw_tif.type().is() is true. + * It has no public constructors and can only be created by a Grey_tif object. The size() member function returns zero. + */ +template requires std::is_same_v, std::byte> +class Grey_tif_image { + + template friend class Grey_tif; + +public: + /** + * @brief Always returns zero to indicate that the image pixel type must be determined at runtime. + * + * @return zero. + */ + std::size_t size() const noexcept { return 0; } + + /** + * @brief Get the dimensions of the raw TIFF image. + * + * @return An array containing the width and height of the image. + */ + std::array const& dim() const noexcept { return d_dim; } + + /** + * @brief Get the pixel data type of the TIFF image. + * + * @return The pixel data type of the image as a POD_type_traits. + * @see POD_type_traits + */ + POD_type_traits const& type() const noexcept {return d_image_type;} + + /** + * @brief Typecasting a Grey_tif_image to its runtime POD type + * + * Grey_tif_image and Grey_tif_image can only be typecast to a type that corresponds for which type().is() is true. + * + * @tparam T0 POD type that is encoded by the POD type as returned by type(). + * @return Grey_tif_image + */ + template + operator Grey_tif_image() const { + static_assert ((std::is_integral_v && sizeof(T0) <= 4 && !std::is_same_v) || + (std::is_floating_point_v && sizeof(T0) <= 8), + "Only signed & unsigned 8-, 16 and 32-bit integers, floats and doubles allowed in TIFF image data"); + if constexpr (std::is_const_v) + static_assert(std::is_const_v, "Can only cast a Grey_tif_image to a Grey_tif_image, when T is const"); + assert(type().template is()); + return Grey_tif_image(type(), dim(), std::span(reinterpret_cast(d_data), dim()[0] * dim()[1])); + } + +private: + Grey_tif_image(POD_type_traits const& image_type, std::array const& dim, std::span const& data) : + d_data(data.size() > 0 ? reinterpret_cast(&data.begin()[0]) : nullptr), + d_image_type(image_type), + d_dim(dim) { } + + T* d_data; + POD_type_traits d_image_type; + std::array d_dim; +}; + +template +class Grey_tif : public Grey_tif_image { + static_assert ((std::is_integral_v && sizeof(T) <= 4 && !std::is_same_v) || + std::is_same_v || + (std::is_floating_point_v && sizeof(T) <= 8), + "Only signed & unsigned 8-, 16 and 32-bit integers, floats and doubles allowed in TIFF image data"); + + template friend class Grey_tif; + +public: + using value_type = T; + + /** + * @brief Default constructor for Grey_tif with no images. + */ + Grey_tif() noexcept : + Grey_tif_image(POD_type_traits{.size = sizeof(T), .is_signed = std::is_signed_v, .is_integral = std::is_integral_v}, + {0, 0}, std::span()) { + d_tif = std::vector(8, std::byte(0x00)); + d_tif[0] = d_tif[1] = std::byte((std::endian::native == std::endian::little) ? 'I' : 'M'); + reinterpret_cast(d_tif[2]) = 42; + d_last_ifd_offset = 4; + } + + /** + * @brief Constructor for Grey_tif of a single image with specified dimensions, populating the data with the container + * + * If the container has a member function dim(), the second parameter may be omitted, and the dimensions of the container will be used. + * + * @tparam C The container type. + * @param container The container containing pixel data. + * @param dim The dimensions of the image in the container, may be omitted if the container has a member function dim(). + */ + template requires (!std::derived_from) + Grey_tif(C const& container, std::array const& dim = {-1, -1}) noexcept : Grey_tif() { + push_back(container, dim); + } + + /** + * @brief Constructor for Grey_tif that populates it from a TIFF input stream. + * + * @param is The input stream containing TIFF data. + */ + Grey_tif(std::istream& is) : Grey_tif() { + is.seekg(0, std::ios::end); + std::streampos fileSize = is.tellg(); + is.seekg(0, std::ios::beg); + d_tif.resize(static_cast(fileSize)); + if (!is.read(reinterpret_cast(d_tif.data()), fileSize) || + (d_tif[0] != std::byte('I') && d_tif[0] != std::byte('M')) || + d_tif[0] != d_tif[1] || std::byte(42) != (d_tif[0] == std::byte('I') ? d_tif[2] : d_tif[3]) ) + is.setstate(std::ios::failbit); + else + f_scan_images(); + } + + /** + * @brief swapping the contents of two TIFF stacks + * + * @tparam Tother The type of the Grey_tif object to be swapped. + * @param other The Grey_tif object to be swapped. + */ + template + void swap(Grey_tif&& other) { + d_last_ifd_offset = other.d_last_ifd_offset = 4; + std::swap(d_tif, other.d_tif); + f_scan_images(); + other.f_scan_images(); + if constexpr(!std::is_same_v) + f_regularize(); + if constexpr(!std::is_same_v) + other.f_regularize(); + } + + /** + * @brief swapping the contents of two TIFF stacks + * + * Swaps the contents of this Grey_tif with that of another. If required, the data are regularized, so that they conform to the pixel type of the Grey_tif objects. + * + * @tparam Tother The type of the Grey_tif object to be swapped. + * @param other The Grey_tif object to be swapped. + */ + template + void swap(Grey_tif& other) { swap(std::move(other)); } + + /** + * @brief Read/write accessor to images from a single or multi-image Grey_tif object. + * + * @param i The index of the image in the stack. + * @return A reference to the requested image. + */ + auto const& image(std::size_t const i = 0) const { return d_img.at(i); } + + /** + * @brief Read accessor to const images from a single or multi-image Grey_tif object. + * + * @param i The index of the image in the stack. + * @return A reference to the requested image. + */ + auto const& const_image(int const i = 0) const { return d_img_const.at(i); } + + /** + * @brief Read/write accessor to images from a single or multi-image raw Grey_tif object. + * + * Read/write access by index. If the pixel type T0 is incorrect, the image has a zero size. + * + * @tparam T0 The desired pixel data type. + * @param i The index of the image in the stack. + * @return Image with the specified pixel type. + */ + template requires std::is_same_v + auto image(int const i) const { return static_cast>(d_img.at(i)); } + + /** + * @brief Accessor for retrieving a const image from a raw Grey_tif object, of which the pixels cannot be edited. + * + * Read access by index. If the pixel type T0 is incorrect, the image has a zero size. + * + * @tparam T0 The desired pixel data type. + * @param i The index of the image in the stack. + * @return Image with the specified pixel type. + */ + template requires std::is_same_v + auto const_image(int const i) const { return Grey_tif_image(d_img_const.at(i)); } + + /** + * @brief Get the number of images in the TIFF stack. + * + * @return The number of images in the TIFF stack. + */ + std::size_t image_stack_size() const noexcept {return d_img.size();} + + /** + * @brief Get the number of bytes of the TIFF data. + * + * @return The number of images bytes of the TIFF data. + */ + std::size_t raw_data_size() const noexcept {return d_tif.size();} + + /** + * @brief Write the TIFF data to an output stream. + * + * @param os The output stream where TIFF data will be written. + */ + void write(std::ostream& os) const { os.write(reinterpret_cast(d_tif.data()), d_tif.size());} + + /** + * @brief Appends a new image from an STL container or a Grey_tif object to the end of the image stack. + * + * For a Grey_tif object (encodes raw TIFF data) that is appended, the pixel values in the input container must be either + * signed or unsigned 8-, 16-, 32-bit integrals, float, or double. + * + * @tparam C The type of the container. + * @param container The image container to push into the Grey_tif object. + * @param dim The dimension of the image as a 2D array of long, can be omitted if the container has a member function dim(). + */ + template requires (requires (C const& c) {std::begin(c), std::end(c), std::size(c);} || std::is_same_v>) + void push_back(C const& container, std::array dim = {-1,-1}) { + if constexpr (is_Grey_tif::value) { + for (int i = 0; i < container.image_stack_size(); ++i) { + if constexpr (!std::is_same_v>) + push_back(container.d_img[i]); + else { + POD_type_traits const& img_type = container.image(i).type(); + if (img_type.is()) push_back(Grey_tif_image(container.d_img[i])); + else if (img_type.is()) push_back(Grey_tif_image(container.d_img[i])); + else if (img_type.is()) push_back(Grey_tif_image(container.d_img[i])); + else if (img_type.is()) push_back(Grey_tif_image(container.d_img[i])); + else if (img_type.is()) push_back(Grey_tif_image(container.d_img[i])); + else if (img_type.is()) push_back(Grey_tif_image(container.d_img[i])); + else if (img_type.is()) push_back(Grey_tif_image(container.d_img[i])); + else if (img_type.is()) push_back(Grey_tif_image(container.d_img[i])); + } + } + } + else { + if constexpr (requires(C const& c) { c.dim().size();}) + if (dim[0] == -1 && dim[1] == -1) { + assert(container.dim().size() == 2); // Grey_tif.pushback(container) only accepts 2D containers + dim[0] = container.dim()[0]; + dim[1] = container.dim()[1]; + } + using CT = typename C::value_type; + static_assert ((std::is_arithmetic_v && sizeof(CT) <= 8 && !std::is_same_v && !(sizeof(CT) == 8 && std::is_integral_v)) || + std::is_same_v || std::is_same_v); + assert(container.size() == dim[0] * dim[1]); + assert((d_tif.size() + dim[0] * dim[1] + 6 * 12 + 4) < std::numeric_limits::max()); + std::uint32_t index = static_cast(d_tif.size()); + std::uint32_t data_start = index; + if constexpr (std::is_same_v) { + d_tif.resize(d_tif.size() + dim[0] * dim[1] * sizeof(CT) + 7 * 12 + 6); + for (auto const val : container) { + reinterpret_cast(d_tif[index]) = static_cast(val); + index += sizeof(CT); + } + } + else { + d_tif.resize(d_tif.size() + dim[0] * dim[1] * sizeof(T) + 7 * 12 + 6); + for (auto const val : container) { + reinterpret_cast(d_tif[index]) = static_cast(val); + index += sizeof(T); + } + } + if ((d_tif.size() & 1) == 1) { + d_tif.push_back(std::byte(0)); + ++index; + } + reinterpret_cast(d_tif[d_last_ifd_offset]) = index; + reinterpret_cast(d_tif[index]) = 7; + index += 2; + f_set_ifd(index, 0x0100, 3, static_cast(dim[1])); + f_set_ifd(index, 0x0101, 3, static_cast(dim[0])); + f_set_ifd(index, 0x0102, 3, 8 * (std::is_same_v ? sizeof(CT) : sizeof(T))); + f_set_ifd(index, 0x0103, 3, 1); + f_set_ifd(index, 0x0106, 3, 1); + f_set_ifd(index, 0x0111, 4, data_start); + if constexpr (std::is_same_v) { + if (std::is_unsigned_v ) + f_set_ifd(index, 0x0153, 3, 1); + else if (std::is_integral_v) + f_set_ifd(index, 0x0153, 3, 2); + else if (std::is_floating_point_v) + f_set_ifd(index, 0x0153, 3, 3); + } + else { + if (std::is_unsigned_v ) + f_set_ifd(index, 0x0153, 3, 1); + else if (std::is_integral_v) + f_set_ifd(index, 0x0153, 3, 2); + else if (std::is_floating_point_v) + f_set_ifd(index, 0x0153, 3, 3); + } + d_tif[index] = d_tif[index + 1] = d_tif[index + 2] = d_tif[index + 3] = std::byte(0); + d_last_ifd_offset = static_cast(index); + f_scan_images(); + } + } + + /** + * @brief Appends an empty image to the end of the image stack. + * + * @tparam C The pixel type of the empty image. + * @param dim The dimension of the image as a 2D array of long. + */ + template requires (std::is_same_v || (std::is_same_v && + ((std::is_integral_v && sizeof(T) <= 4 && !std::is_same_v) || + (std::is_floating_point_v && sizeof(C) <= 8)))) + void push_back(std::array const dim) { + f_push_back(dim, POD_type_traits{.size = sizeof(C), .is_signed = std::is_signed_v, .is_integral = std::is_integral_v}); + } + + /** + * @brief Wipes all images, resetting to an empty Grey_tif object + */ + Grey_tif& clear() { + Grey_tif tmp; + this->swap(tmp); + this->d_image_type = {.size = sizeof(T), .is_signed = std::is_signed_v, .is_integral = std::is_integral_v}; + this->d_dim = {0, 0}; + std::span::operator=({}); + return *this; + } + +private: + std::uint32_t d_last_ifd_offset = 0; // location of the last tif tag + std::vector d_tif; // raw tiff data as std::byte + std::vector> d_img; // vector of images that retrieve their pixel values from d_tif + std::vector> d_img_const; // vector of images that retrieve their pixel values from d_tif + + template struct is_Grey_tif : std::false_type {}; + template struct is_Grey_tif> : std::true_type {}; + + bool f_all_same_type(POD_type_traits const& check_type) const noexcept { + for (auto it = d_img.begin(); it != d_img.end(); ++it) + if (it->type() != check_type) + return false; + return true; + } + + // Push_back an all-zero image of the specified type + void f_push_back(std::array const& dim, POD_type_traits const& new_type) { + assert((d_tif.size() + static_cast(dim[0] * dim[1]) + 6 * 12 + 4) < std::numeric_limits::max()); + assert((this->type().template is() || std::is_same_v)); + std::uint32_t index = static_cast(d_tif.size()); + std::uint32_t data_start = index; + d_tif.resize(d_tif.size() + static_cast(dim[0] * dim[1]) * new_type.size + 7 * 12 + 6); + index += static_cast(dim[0] * dim[1]) * new_type.size; + if ((d_tif.size() & 1) == 1) { + d_tif.push_back(std::byte(0)); + ++index; + } + reinterpret_cast(d_tif[d_last_ifd_offset]) = index; + reinterpret_cast(d_tif[index]) = 7; + index += 2; + f_set_ifd(index, 0x0100, 3, static_cast(dim[0])); + f_set_ifd(index, 0x0101, 3, static_cast(dim[1])); + f_set_ifd(index, 0x0102, 3, 8 * static_cast(new_type.size)); + f_set_ifd(index, 0x0103, 3, 1); + f_set_ifd(index, 0x0106, 3, 1); + f_set_ifd(index, 0x0111, 4, data_start); + f_set_ifd(index, 0x0153, 3, new_type.is_integral ? (!new_type.is_signed ? 1 : 2) : 3); + d_tif[index] = d_tif[index + 1] = d_tif[index + 2] = d_tif[index + 3] = std::byte(0); + d_last_ifd_offset = static_cast(index); + f_scan_images(); + } + + Grey_tif& f_regularize() noexcept { + static_assert(!std::is_same_v); + POD_type_traits new_type = {.size = sizeof(T), .is_signed = std::is_signed_v, .is_integral = std::is_integral_v }; + if (f_all_same_type(new_type)) + return *this; + else { + bool same_size = this->type().size == new_type.size; + for (auto it = d_img.begin() + 1; it != d_img.end() && same_size == true; ++it) + same_size = it->type().size == new_type.size; + if (same_size) { + for (auto& img : d_img) { + if (!new_type.is_integral && img.type().is_integral && img.type().is_signed) + for (auto& pixel : img) + reinterpret_cast(pixel) = static_cast(reinterpret_cast(pixel)); + else if (!new_type.is_integral && img.type().is_integral && !img.type().is_signed) + for (auto& pixel : img) + reinterpret_cast(pixel) = static_cast(reinterpret_cast(pixel)); + else if (new_type.is_integral && new_type.is_signed && !img.type().is_integral) + for (auto& pixel : img) + reinterpret_cast(pixel) = static_cast(reinterpret_cast(pixel)); + else if (new_type.is_integral && !new_type.is_signed && !img.type().is_integral) + for (auto& pixel : img) + reinterpret_cast(pixel) = static_cast(reinterpret_cast(pixel)); + img = Grey_tif_image(new_type, img.dim(), img); + } + } + else { + Grey_tif new_tif; + for (std::size_t i=0; i != image_stack_size(); ++i) { + new_tif.f_push_back(image(i).dim(), new_type); + POD_type_traits const& old_type = image(i).type(); + std::size_t size = static_cast(image(i).dim()[0] * image(i).dim()[1]); + auto new_begin = new_tif.image(i).begin(); + if (old_type.is()) std::copy_n(reinterpret_cast(&d_img[i][0]) , size, new_begin); + else if (old_type.is()) std::copy_n(reinterpret_cast(&d_img[i][0]) , size, new_begin); + else if (old_type.is()) std::copy_n(reinterpret_cast(&d_img[i][0]) , size, new_begin); + else if (old_type.is()) std::copy_n(reinterpret_cast(&d_img[i][0]) , size, new_begin); + else if (old_type.is()) std::copy_n(reinterpret_cast(&d_img[i][0]) , size, new_begin); + else if (old_type.is()) std::copy_n(reinterpret_cast(&d_img[i][0]) , size, new_begin); + else if (old_type.is()) std::copy_n(reinterpret_cast(&d_img[i][0]) , size, new_begin); + else if (old_type.is()) std::copy_n(reinterpret_cast(&d_img[i][0]) , size, new_begin); + } + std::swap(new_tif, *this); + } + return *this; + } + } + + void f_scan_images() { + d_img.clear(); + d_img_const.clear(); + this->d_dim = {0,0}; + bool native = d_tif[0] == std::byte('I') && std::endian::native == std::endian::little; + if (!native) + d_tif[0] = d_tif[1] = (d_tif[0] == std::byte('M')) ? std::byte('I'): std::byte('M'); + reinterpret_cast(d_tif[2]) = 42; + std::uint32_t index = reinterpret_cast(d_tif[4]); + if (!native) { + index = (index << 24) | ((index & 0xff00) << 8) | ((index >> 8) & 0xff00) | (index >> 24); + reinterpret_cast(d_tif[4]) = index; + } + if (native) + while (index != 0) { + f_make_Image(index); + index = reinterpret_cast(d_tif[d_last_ifd_offset]); + } + else + while (index != 0) { + f_make_Image(index); + index = reinterpret_cast(d_tif[d_last_ifd_offset]); + } + if (image_stack_size() != 0) { + if constexpr(!std::is_same_v) + std::span::operator=(d_img_const[0]); + else + this->d_data = d_img_const[0].d_data; + this->d_dim = d_img[0].dim(); + this->d_image_type = d_img[0].type(); + } + if constexpr (!std::is_same_v) + f_regularize(); + } + + template + void f_make_Image(std::uint32_t& index) { + bool compatible_tif = true; + std::array dim = {0,0}; + std::size_t bits_per_pixel = 0; + std::byte* cursor = d_tif.data() + index; + int const tag_count = f_int16(cursor); +// std::int32_t rows_per_strip; + std::vector strip_offsets(1); + std::vector strip_byte_counts(1); + bool signed_pixels = false; + bool int_pixels = true; + for (int i = 0; i != tag_count; ++i) { + uint16_t const tag = f_int16(cursor); + uint16_t const type = f_int16(cursor); + uint32_t const count = f_int32(cursor); + uint32_t val(0); + double double_val(0); + if ((type == 1) || (type == 2) || (type == 6) || (type == 7)) { + val = f_int8(cursor); + cursor += 3; + } + else if ((type == 3) || (type == 8)) { + val = f_int16(cursor); + cursor += 2; + } + else if ((type == 4) || (type == 9)) + val = f_int32(cursor); + else if (type == 5 || type == 10) { + std::byte* p = d_tif.data() + f_int32(cursor); + double_val = f_int32(p); + double_val /= f_int32(p); + } + else if (type == 11) { + val = f_int32(cursor); + double_val = reinterpret_cast(val); + val = 0; + } + else if (type == 12) { + std::byte* p = d_tif.data() + f_int32(cursor); + std::uint64_t tmp = f_int64(p); + double_val = reinterpret_cast(tmp); + } + if (tag == 0x0100) dim[0] = val; + else if (tag == 0x0101) dim[1] = val; + else if (tag == 0x0102 ) { + if ((val == 8 || val == 16 || val == 32 || val == 64)) + bits_per_pixel = val; + else { + std::cerr << "Warning: Grey_tif can only read greyscale tiff files with 8-, 16-, 32-, or 64-bit pixels" << std::endl; + compatible_tif = false; + } + } + else if (tag == 0x0103 && val != 1) { + std::cerr << "Warning: Grey_tif cannot read compressed tiff files"<< std::endl; + compatible_tif = false; + } + else if (tag == 0x0106 && val > 1) { + std::cerr << "Warning: Grey_tif cannot read colour tiff files"<< std::endl; + compatible_tif = false; + } + else if ((tag == 0x0107 || tag == 0x0108 || tag == 0x0109 || tag == 0x010A) && val != 1) { + std::cerr << "Warning: Grey_tif cannot read black & white tiff files"<< std::endl; + compatible_tif = false; + } + else if (tag == 0x0111) { + if (count == 1) strip_offsets[0] = static_cast(val); + else { + strip_offsets.resize(count); + std::byte* p = d_tif.data() + val; + for (std::size_t j = 0; j != static_cast(count); ++j) + strip_offsets[j] = static_cast(f_int32(p)); + } + } + else if (tag == 0x0115 && val != 1) { + std::cerr << "Warning: Grey_tif cannot read RGB colour tiff files"<< std::endl; + compatible_tif = false; + } +// else if (tag == 0x0116) rows_per_strip = static_cast(val); + else if (tag == 0x0117) { + if (count == 1) strip_byte_counts[0] = static_cast(val); + else { + strip_byte_counts.resize(count); + std::byte* p = d_tif.data() + val; + for (std::size_t j=0; j != count; ++j) + strip_byte_counts[j] = static_cast(f_int32(p)); + } + } + else if (tag == 0x0153) { + if (val != 1) signed_pixels = true; + if (val == 3) int_pixels = false; + } + } + for (std::size_t i=0; i != strip_offsets.size() - 1; ++i) + if (strip_byte_counts[i] != strip_offsets[i + 1] - strip_offsets[i]) { + std::cerr << "Warning: Grey_tif cannot read tiff files with non-consecutive strips" << std::endl; + std::cerr << " most likely the tiff file is corrupted" << std::endl; + compatible_tif = false; + break; + } + if (!compatible_tif) + throw std::runtime_error("Incompatible TIFF file\n"); + d_last_ifd_offset = static_cast(cursor - d_tif.data()); + index = f_int32(cursor); + if constexpr (!NATIVE) { + std::byte* begin = d_tif.data() + strip_offsets[0]; + std::byte* end = begin + static_cast(dim[0] * dim[1]) * bits_per_pixel / 8; + if (bits_per_pixel == 16) + while (begin != end) f_int16(begin); + else if (bits_per_pixel == 32) + while (begin != end) f_int32(begin); + else if (bits_per_pixel == 64) + while (begin != end) f_int64(begin); + } + POD_type_traits image_type = { .size = bits_per_pixel / 8, .is_signed = signed_pixels, .is_integral = int_pixels }; + d_img.push_back(Grey_tif_image(image_type, dim, + std::span(reinterpret_cast(&d_tif[static_cast(strip_offsets[0])]), + static_cast(dim[0] * dim[1])))); + d_img_const.push_back(Grey_tif_image(image_type, dim, + std::span(reinterpret_cast(&d_tif[static_cast(strip_offsets[0])]), + static_cast(dim[0] * dim[1])))); + }; + + void f_set_ifd(std::uint32_t& index, uint16_t tag, uint16_t type, uint32_t val) { + reinterpret_cast(d_tif[index]) = tag; + reinterpret_cast(d_tif[index + 2]) = type; + reinterpret_cast(d_tif[index + 4]) = 1; + if (type == 1) + reinterpret_cast(d_tif[index + 8]) = static_cast(val); + else if (type == 3) + reinterpret_cast(d_tif[index + 8]) = static_cast(val); + else if (type == 4) + reinterpret_cast(d_tif[index + 8]) = val; + index += 12; + }; + + std::uint8_t f_int8(std::byte* &cursor) { + return reinterpret_cast(*cursor++); + } + + template + constexpr std::uint16_t f_int16(std::byte* &cursor) { + std::uint16_t r; + std::memcpy(&r, cursor, sizeof(std::uint16_t)); + cursor += 2; + return NATIVE ? r : r = static_cast(r << 8) | static_cast(r >> 8); + } + + template + constexpr std::uint32_t f_int32(std::byte* &cursor) const { + std::uint32_t r; + std::memcpy(&r, cursor, sizeof(std::uint32_t)); + cursor += 4; + return NATIVE ? r : r = (r << 24) | ((r & 0x00FF0000) >> 8) | ((r & 0x0000FF00) << 8) | (r >> 24); + } + + template + constexpr std::uint64_t f_int64(std::byte* &cursor) { + std::uint64_t r; + std::memcpy(&r, cursor, sizeof(std::uint64_t)); + cursor += 8; + return NATIVE ? r : r = ((r << 56) | ((r & 0xff00) << 40) | ((r & 0xff0000) << 24) | ((r & 0xff000000) << 8) | + ((r >> 8) & 0xff000000) | ((r >> 24) & 0xff0000) | ((r >> 40) & 0xff00) | (r >> 56)); + } +}; + +// Constructor deduction guide for Grey_tif +template Grey_tif(C const& container, std::array const& dim) -> Grey_tif; + +} // namespace jpa +#endif /* Grey_tif_h */ diff --git a/lib/trpx2.0/include/Operators.hpp b/lib/trpx2.0/include/Operators.hpp new file mode 100755 index 00000000..9d1892c6 --- /dev/null +++ b/lib/trpx2.0/include/Operators.hpp @@ -0,0 +1,194 @@ +// +// Operators.hpp +// Operators +// +// Created by Jan Pieter Abrahams on 13/04/2019. +// Copyright © 2019 Jan Pieter Abrahams. All rights reserved. +// + +#ifndef Operators_h +#define Operators_h +#include + +// Some useful operators and constants in a namespace Operator:: that are not defined in the C++ standard: +// +// static const bool big_endian +// Runtime function returning true if running on a big-endian machine, otherwise returns false. Will +// be superceded in C++20 by the class std::endian. +// swap_bytes(T val) +// Swaps byte order of 'val'. Note that 'val' can also be of type std::complex +// euclidian_remainder(T v, Ti p) +// Returns the Euclidian remainer (positive modulus). Faster than operator%() or fmod() when 'v' is +// within the interval [-p; 2p] and 'p' is positive. When 'v' outside this interval, or when (p<0), +// use ((v %= p) < 0) ? v+abs(p) : v for integral types, or ((v = fmod(v,p) < 0) ? v+abs(p) : v for +// real-valued types. +// euclidian_remainder(T v) +// Returns the Euclidian remainer for N is a positive power of 2 and 'v' is an integral type. +// euclidian_division(T v, Ti p) +// Returns an anonymous struct with members (anonymous)::quot and (anonymous)::rem, which +// are the integral divisor and the positive remainder of v / p. The type of the quotient and +// remainder is decltype(T*Ti). +// euclidian_division(T v) +// Fast euclidian division for N is a positive power of 2 and 'v' is an integral type. If N is not +// a positive power of 2, euclidian_division(v, N) is called instead. +// int highest_set_bit(T val) +// Returns the number of the most significant bit that is set in an integral type T. If no bits are +// set, 0 is returned. If T is signed, the sign bit is added, unless (val == 0). So +// (highest_set_bit(uint8_t(0b00100000)) == 6), whilst for a signed argument: +// (highest_set_bit( int8_t(0b00100000)) == 7). +// int lowest_set_bit(T val) +// Returns the number of the least significant bit that is set in an integral type T. If no bits +// are set, returns 0. So (lowest_set_bit(uint8_t(0b00100100)) == 3) and (lowest_set_bit(0) == 0). +// bool is_bounded(std::initializer_list const& list) +// Returns true if the values of 'list' between the first and last item on the list, are at least +// equal to the first item of the list and smaller than the last item of the list. +// std::string filename_extension(std::string const& filename) +// Returns the filename extension: string starting at the last occurrence of the character '.' +// For instance: filename_extension("foo.bar") == ".bar". +// std::string strip_filename_extension(std::string const& filename) +// Returns the filename without extension: string starting up to the last occurrence of the character '.' +// For instance: filename_extension("foo.bar") == "foo". + +namespace Operator { + +static const bool big_endian = []() { + int i=1; + return ! *(reinterpret_cast(&i)); +}(); + +template +constexpr inline T swap_bytes(T val) noexcept { + uint8_t *s = reinterpret_cast(&val); + uint8_t *e = s + sizeof(T) - 1; + while (s < e) std::swap(*s++, *e--); + return val; +} + +template +constexpr inline std::complex swap_bytes(std::complex const& val) noexcept { + T r = val.real(); + T i = val.imag(); + return std::complex(swap_bytes(r), swap_bytes(i)); +} + +template +constexpr inline decltype(T0()*T1()) euclidian_remainder(T0 v, T1 const p) noexcept { + if constexpr (std::is_signed_v) { + if (((v < 0) && ((v += p) < 0)) || + ((v >= p) && ((v -= p) >= p))) { + if constexpr (std::is_integral_v) { + if constexpr (std::is_signed_v) + return ((v -= (v / p) * p) < 0) ? v+abs(p) : v; // not use % as this may have a similar optimisation + else if (v < 0) + return p - ((-v) - ((-v) / p) * p); + else + return v - (v / p) * p; + } + else { + if constexpr (std::is_signed_v) + return ((v = fmod(v,p)) < 0) ? v+abs(p) : v; + else + return ((v = fmod(v,p)) < 0) ? v+p : v; + } + } + } + else if constexpr (std::is_signed_v) { + if ((p > 0) && (v >= p) && ((v -= p) >= p)) { + if constexpr (std::is_integral_v) + return v - (v / p) * p; + else + return v = fmod(v,p); + } + else if ((v += p) >= (-p)) { + if constexpr (std::is_integral_v) + return v - (v / (-p)) * (-p); + else + return v = fmod(v,p); + } + } + else if ((v >= p) && ((v -= p) >= p)) + return v - (v / p) * p; + return v; +} + +template +constexpr inline T0 euclidian_remainder(T0 const v) noexcept { + if constexpr (std::is_integral_v && ((N & (N - 1)) == 0)) + return v & N-1; + else + return euclidian_remainder(v, N); +} + +template +constexpr inline auto euclidian_division(T0 const v, T1 const p) noexcept { + if constexpr (std::is_integral_v) { + auto r = std::div( v, p); + if (r.rem < 0) { + r.rem += p; + --r.quot; + } + return r; + } + else { + struct {decltype(T0()*T1()) quot; decltype(T0()*T1()) rem;} r {}; + r.rem = euclidian_remainder(v, p); + r.quot = (v - r.rem) / p; + return r; + } +} + +template +constexpr inline auto euclidian_division(T0 const v) noexcept { + if constexpr (std::is_integral_v && ((N & (N - 1)) == 0)) { + struct {T0 quot; T0 rem;} r {}; + r.rem = v & T0(N-1); + r.quot = (v - r.rem) / T0(N); + return r; + } + else + return euclidian_division(v, N); +} + +template +constexpr inline int highest_set_bit(T val) noexcept { + static_assert(std::is_integral_v , "most_significant_bit() requires an integral argument."); + if constexpr (std::is_signed_v) + return (val == 0) ? 0 : 1 + highest_set_bit(std::make_unsigned_t (abs(val))); + else { + int r=0; + for ( ; val; val>>=1, ++r); + return r; + } +} + +template +constexpr inline int lowest_set_bit(T val) noexcept { + static_assert(std::is_integral_v , "least_significant_bit() requires an integral argument."); + if (val) { + int r=1; + val = (val ^ (val - 1)) >> 1; + for ( ; val; val>>=1, ++r); + return r; + } + else return 0; +} + +template +constexpr inline bool is_bounded(std::initializer_list const& list) noexcept { + auto const last = list.end() - 1; + for (auto p = list.begin() + 1; p != last; ++p) + if ((*p < *list.begin()) || (*p >= *last)) + return false; + return true; +} + +inline std::string filename_extension(std::string const& filename) noexcept { + return filename.substr(filename.find_last_of('.')); +} + + inline std::string strip_filename_extension(std::string const& filename) noexcept { + return filename.substr(0, filename.find_last_of('.')); + } + +} +#endif /* Operators_h */ diff --git a/lib/trpx2.0/include/Shared_array.hpp b/lib/trpx2.0/include/Shared_array.hpp new file mode 100644 index 00000000..3945850b --- /dev/null +++ b/lib/trpx2.0/include/Shared_array.hpp @@ -0,0 +1,177 @@ +// +// Shared_array.hpp +// Shared_array +// +// Created by Jan Pieter Abrahams on 23.01.2024. +// +// The `Shared_array` class is a light-weight C++ container designed for managing arrays with shared ownership +// semantics, combining aspects of `std::shared_ptr` and `std::span`. It acts as a smart pointer for arrays, +// offering automatic memory management with shared ownership, alongside convenient array-like access. This class +// is particularly useful in scenarios where multiple parts of a program need to manage and access the same array +// safely and efficiently. +// +// ### Key Features: +// 1. **Shared Ownership**: Unlike `Unique_array` or 'std::vector', `Shared_array` allows multiple instances to +// share ownership of the underlying array. The array is automatically deallocated when the last `Shared_array` +// owning it goes out of scope, thus preventing memory leaks. +// +// 2. **Slicing and Viewing Data Structure**: This feature complements the shared ownership model by allowing +// users to create slices or views of the array. These views can access specific segments of the data without +// duplication, all while maintaining the guarantee that the underlying memory won't be released as long +// as any view or slice remains active. This ensures safety and consistency when accessing subsets of the array. +// +// 3. **Array-like Access**: Inheriting from `std::span`, `Shared_array` provides direct, bounds-checked access +// to its elements, similar to standard C++ arrays and containers like `std::vector`. +// +// 4. **Support for Dynamic and Fixed-Size Arrays**: The class template can be used for both dynamic arrays +// (where the size is determined at runtime) and fixed-size arrays (with compile-time size). +// +// 5. **Copy and Move Semantics**: `Shared_array` supports both copying and moving, making it flexible for +// use in a variety of contexts where shared data management is needed. +// +// 6. **Reinterpreting Data**: This feature complements the shared ownership model by allowing users to reinterpret +// data, for instance, to access a Shared_array of bytes as a Shared_array of unsigned longs. +// +// ### Use Cases: +// - **Collaborative Data Management**: Ideal for applications where multiple components or threads need to +// access and potentially modify the same array data. +// +// - **Complex Data Structures in Multi-Threaded Environments**: Useful for managing shared resources in +// multi-threaded applications, facilitating safe and efficient access to shared data. +// +// - **Slicing and Viewing Data Structure**: Useful for extracting slices or views that allow access to common +// data, with the guarantee that associated memory is not released until all views or slices have gone out of +// scope. +// +// - **Graphs and Trees with Shared Nodes**: In data structures like graphs or trees, `Shared_array` can +// efficiently manage arrays of shared nodes or edges. +// +// - **Cache and Buffer Sharing**: For scenarios where different parts of an application need to share access +// to a common cache or buffer. +// +// ### Example: +// jpa::Shared_array sharedData1(1000); // Shared array of 1000 doubles +// auto sharedData2 = sharedData1; // sharedData2 now shares ownership with sharedData1 +// +// // Both sharedData1 and sharedData2 can access and modify the array +// for (int i = 0; i < 1000; ++i) { +// sharedData1[i] = computeSomeValue(i); +// } +// +// In this example, `Shared_array` is used to create a shared array of doubles. Modifications through +// either reference will be reflected in the other. The array's lifetime is managed through shared ownership, +// ensuring that it exists as long as any `Shared_array` instance owns it. +// +// In summary, `Shared_array` is a valuable tool for managing arrays in C++, particularly useful in shared +// ownership scenarios. Its design offers a balance between safety, efficiency, and shared data management, +// suitable for a variety of applications, especially those involving collaborative data access or complex +// data structures. + +#ifndef SHARED_ARRAY_HPP +#define SHARED_ARRAY_HPP + +#include +#include +#include +#include + +namespace jpa { + +/// \brief A template class for managing a shared array. +/// \details Shared_array combines shared ownership (like std::shared_ptr) with array-like access (like std::span). +/// It can manage either dynamic arrays (with std::dynamic_extent) or fixed-size arrays. +/// \tparam T The type of the elements in the array. +/// \tparam N The size of the array, defaults to std::dynamic_extent for a dynamic array. +template +class Shared_array : +private std::shared_ptr, +public std::span { + template friend class Shared_array; + static_assert(std::is_default_constructible_v, "Type T must be default constructible for Shared_array"); +public: + /// \brief Constructor for a dynamic array. + /// \param size The size of the array to create. + /// \note This constructor is only available for dynamic arrays (N == std::dynamic_extent). + explicit Shared_array(std::size_t size) noexcept requires (N == std::dynamic_extent) : + std::shared_ptr(std::make_shared_for_overwrite(size * sizeof(T))), + std::span(reinterpret_cast(this->get()), size) {} + + /// \brief Default constructor for a fixed-size array. + /// \note This constructor is only available for fixed-size arrays (N != std::dynamic_extent). + constexpr Shared_array() noexcept requires (N != std::dynamic_extent) : + std::shared_ptr(std::make_shared_for_overwrite(N * sizeof(T))), + std::span(reinterpret_cast&>(*this->get())) {} + + /// \brief Constructor to create a dynamic Shared_array from another fixed Shared_array of a different type. + /// \tparam U The type parameter of the other Shared_array. + /// \param other The other Shared_array. + template + constexpr Shared_array(Shared_array const& other) noexcept requires (N == std::dynamic_extent) : + std::shared_ptr(other), + std::span(reinterpret_cast(this->get()), other.size() * sizeof(U) / sizeof(T)) { + static_assert(sizeof(U) % sizeof(T) == 0 || sizeof(T) % sizeof(U) == 0, + "Incompatible type sizes for Shared_array reinterpretation."); + } + + /// \brief Constructor to create a fixed-size Shared_array from another fixed Shared_array of a different type. + /// \tparam U The type parameter of the other Shared_array. + /// \param other The other Shared_array. + template + constexpr Shared_array(Shared_array const& other) noexcept requires (N != std::dynamic_extent) : + std::shared_ptr(other), + std::span(reinterpret_cast&>(*this->get())) { + static_assert(N * sizeof(T) == M * sizeof(U), + "Incompatible type sizes for Shared_array reinterpretation."); + } + + /// Shared_array is copyable and movable. + constexpr Shared_array(const Shared_array&) = default; + constexpr Shared_array(Shared_array&&) = default; + constexpr Shared_array& operator=(const Shared_array&) = default; + constexpr Shared_array& operator=(Shared_array&&) = default; + + using std::span::operator[]; + + /// \brief Checks if the current instance is the only owner of the array. + /// \return True if this is the only owner, false otherwise. + constexpr bool is_unique() const noexcept { + return use_count() == 1; + } + + /// \brief Creates a slice of the Shared_array + /// \param offset The starting index of the slice + /// \param length The number of elements in the slice + /// \return A new Shared_array instance that shares ownership of the array but only accesses a subset of it. + constexpr Shared_array slice(std::size_t offset, std::size_t length) const { + if (offset + length > this->size()) + throw std::out_of_range("Slice range is out of bounds"); + return Shared_array(reinterpret_cast&>(*this), this->data() + offset, length); + } + + /// \brief Creates a fixed-sized slice of the Shared_array + /// \param offset The starting index of the slice + /// \tparam M The number of elements in the slice + /// \return A new fixed-size Shared_array instance that shares ownership of the array but only accesses a subset of it. + template + constexpr Shared_array slice(std::size_t offset) const { + if (offset + M > this->size()) + throw std::out_of_range("Slice range is out of bounds"); + return Shared_array(reinterpret_cast&>(*this), this->data() + offset); + } + +private: + // Private constructor for creating dynamic slices + constexpr Shared_array(const std::shared_ptr& sharedArray, T* slicedPtr, std::size_t length) requires (N == std::dynamic_extent) : + std::shared_ptr(sharedArray, reinterpret_cast(slicedPtr)), + std::span(slicedPtr, length) {} + + // Private constructor for creating dynamic slices + constexpr Shared_array(const std::shared_ptr& sharedArray, T* slicedPtr) requires (N != std::dynamic_extent) : + std::shared_ptr(sharedArray, reinterpret_cast(slicedPtr)), + std::span(reinterpret_cast&>(*this->get())) {} +}; + +} // namespace jpa + + +#endif /* SHARED_ARRAY_HPP */ diff --git a/lib/trpx2.0/include/Terse.hpp b/lib/trpx2.0/include/Terse.hpp new file mode 100644 index 00000000..66c8b778 --- /dev/null +++ b/lib/trpx2.0/include/Terse.hpp @@ -0,0 +1,1603 @@ +// +// Terse.hpp +// Terse +// +// Created by Jan Pieter Abrahams on 30/04/2019. +// Copyright © 2019 Jan Pieter Abrahams. All rights reserved. +// + +#ifndef Terse_h +#define Terse_h + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include "Bitqueue.hpp" +#include "Unique_array.hpp" +#include "XML_element.hpp" + +// Terse allows efficient and fast compression of integral diffraction data and other integral greyscale +// data into a Terse object that can be decoded by the member function prolix(). The +// prolix(iterator) member function decompresses the data starting at the location defined by 'iterator' +// (which can also be a pointer). A Terse object is constructed by supplying it with uncompressed data or a +// stream that contains compressed Terse data. +// +// A Terse object may contain compressed data of multiple frames, and data of a particular frame can be +// extracted by indexing. All frames must have the same size and dimensions, and must all be signed or unsigned. +// +// The template parameter C is either void or Concurrent. If C is Concurrent, decompressing the multiple frames +// into a std::vector is multi-threaded, and compressing is multithreaded if the data to be compressed are +// passed as an rvalue container (e.g. Terse terse(std::move(data_vector));). This empties the +// rvalue container, thus optimizing memory use. +// +// A Terse object can be unpacked into any arithmetic type T, including float and double. Unpacking into +// values of a type with fewer bits than the original data is not allowed. Compressing as unsigned yields a tighter +// compression, because the sign bit does not need to be encoded. +// +// A Terse object can be written or appended to any stream. The resulting file is independent of the endian-nes +// of the machine: both big- and small-endian machines produce identical files, making data transfer optimally +// transparent. +// +// Terse data in a file are immediately preceded by a small header, which is encoded in standard XML as follows: +// +// - "n" is the number of bits required for the most extreme value in the Terse data, representing the bit depth of +// the original, uncompressed data. +// - "s" is "0" for unsigned data, "1" for signed data. +// - "b" is the block size of the stretches of data values that are encoded (by default 12 values). +// - "v" is the number of elements of a single frame of a stack (it is the product of the dimensions, if these +// are provided). +// - "f" is the number of frames encoded in the file. +// - "m" is the total memory size of all frames, excluding the header and metadata. +// - "d [...]" is optional and encodes the dimensions of a single frame. Frames can have any number of dimensions. +// - "ms" is optional and encodes the sizes of metadata strings associated with each frame. +// - "mf" is optional and encodes the memory sizes of individual frames in the stack. +// Here is an example: +// +// +// The default algorithm is a run-length encoding type. Each data block (by default 12 integral values in the +// constructor, but this can be changed) is preceded by one or more data block header bits. The values in the +// data block are stripped of their most significant bits, provided they are all zero (for unsigned values), +// or either all zero or all one (for signed values). In the latter case, the sign bit is maintained. So, for a +// block size of 3 with values 3, 4, 2, the encoded bits would be: 011 (denoting 3) 100 (denoting 4) 010 +// (denoting 2). Hence 011100010 would be pushed into the Terse object. In case of signed values -3, 4, 2, the +// encoded bits would be 1011 (denoting -3) 0100 (denoting +4) 0010 (denoting +2), resulting in a data block +// 101101000010. So if the values that need to be encoded are all positive or zero, they should be encoded as +// unsigned for optimal compression: it saves 1 bit per encoded value. +// +// The header bits define how the values are encoded. They have the following following structure: +// bit 1: If the first bit of the block header is set, then there are no more bits in the block header, +// and the parameters of the previous block header are used. +// bit 2-4: The first header bit is 0. The three bits 2 to 4 define how many bits are used per value of the +// encoded block. If bits 2 to 4 are all set, 7 or more bits per value are required and the +// header is expanded by a further 2 bits. +// bit 5-6: The first 4 header bits are 0111. The number encoded bits 5 and 6 is added to 7 and this defines +// how many bits are used to encode the block. So if bits 5 & 6 are 00 then 7 bits are used, 01 means +// 8 bits, 10 means 9 bits and 11 means at least 10 bits. If bits 5 & 6 are both set, the header is +// expanded by another 6 bits. +// bit 7-12: The first 6 header bits are 011111. The number encoded by bits 7 to 12 is added to decimal 10 and +// this defines how many bits are used to encode the block. So if bits 7 to 12 are 000000, then the +// number of bits per value in the data block is decimal 10. If bits 7 to 12 are 110110, then the +// number of bits per value in the data block is 10 + 54 = 64. +// The default compression mode is compatible with ealier versions of Terse and available through the flag Terse_mode::Signed. +// +// For unsigned data that contains overloads (all bits set to 1), a separate encoding scheme is available through +// the flag Terse_mode::Unsigned. In this scheme, the algorithm increments all values in the block by 1 prior to encoding. +// This optimization saves space when many overloads are present. However, earlier versions of Terse will not be able +// to read this data. +// +// An additional algorithm improves compression for very weak data available through the flag Terse_mode::Small_unsigned. +// Blocks with a maximum value of 2, 4, 5, or 6 are encoded using an alternative approach, where values are represented +// in a number base system with radices 3, 5, 6, and 7, respectively. This method reduces the number of bits required +// for encoding these blocks. Furthermore, the algorithm leverages patterns between consecutive blocks. If all values in +// the current block are smaller than 7, it checks whether the block differs by 1 from the previous block. For blocks with +// values greater than 6, it checks if the significant number of bits differs by 1 from the previous block. In these cases, +// the new block size or precision can be inferred from the previous block, further reducing the overhead. Additionally, +// blocks with all zero values are encoded using a single bit, maximizing efficiency in cases with sparse data. Earlier +// versions of Terse will not be able to read this data. +// +// Constructors: +// Terse(std::ifstream& istream) +// Reads in a Terse object that has been written to a file by Terse::write(...). +// Terse(container_type&& data, Terse_mode const mode = Terse_mode::Signed) +// Creates a Terse object from data (which can be a std::vector, Field, etc.). Only containers of +// integral types are allowed. If the container has a member function dim(), that will set the dimensions +// of the Terse object. Otherwise the dimensions can be set once using the dim(vector const&) member function. +// If the the 'data' parameter is an rvalue and the Terse template parameter C is Concurrent, compression is +// branched to a different thread and proceeds concurrently. In this case, the 'data' container is emptied. +// Terse(iterator begin, std::size_t size, Terse_mode const mode = Terse_mode::Signed) +// Creates a Terse object given a starting iterator or pointer and the number of elements that need to be +// encoded. +// +// Member functions: +// void insert(std::size_t const pos, Iterator const data, size_t const size, Terse_mode const mode = Terse_mode::Signed) noexcept +// Adds another frame to the Terse object at the position defined by 'pos'. The new frame is defined by +// its begin iterator and size. +// void insert(std::size_t const pos, C&& data, Terse_mode const mode = Terse_mode::Signed) noexcept +// Inserts a frame into the Terse object at the position defined by 'pos'. The new frame is defined by a reference +// to a container. The size and dimensions of the input container must be the same as that of the first frame that +// was used for creating the Terse object, unless the Terse object is empty, in which case the inserted container +// determines the size and dimension of subsequent frames that are added. +// If the input container is passed as an r-value, it will be consumed, and will be empty after the call. +// If the input container is passed as an r-value, and the include file Concurrent.hpp is included, compression +// will be performed concurrently in the background. +// void push_back(Iterator const begin, size_t const size, Terse_mode const mode = Terse_mode::Signed) +// Adds another frame to the Terse object, given a starting iterator or pointer and the number of elements +// that need to be encoded. The 'size' must be identical to that of existing frames. +// void push_back(container_type&& container, Terse_mode const mode = Terse_mode::Signed) +// Adds another frame to the Terse object. The new frame is defined by its begin iterator and size, or by a +// reference to a container. The size must be the same as that of the frame used to create the Terse object. +// If the container has a member function dim(), that must return the same dimension as provided for the first +// frame. +// If the the 'data' parameter is an rvalue and the Terse template parameter C is Concurrent, compression is +// branched to a different thread and proceeds concurrently. In this case, the 'data' container is emptied. +// void erase(std::size_t pos) noexcept +// Removes the frame with index 'pos' from the Terse object. Also waits until concurrent compression has finished, +// and releases unused storage to the heap. +// Terse at(std::size_t pos) noexcept +// Returns the frame with index 'pos' as a Terse object. +// void prolix(iterator begin, std::size_t const pos = 0) +// Unpacks the Terse frame with index 'pos', storing it from the location defined by 'begin'. +// Terse integral signed data cannot be unpacked into integral unsigned data. Terse data cannot be decompressed +// into elements that have fewer bits than bits_per_val(), but can be decompressed into larger values. Terse data +// can always be unpacked into signed integral, double and float data and will have the correct sign (with one +// exception: an unsigned overflowed - all 1's - value will be unpacked as -1 signed value. As all other values are +// positive in this case). +// void prolix(container_type& container, std::size_t const pos = 0) +// Unpacks the Terse frame with index 'pos' and stores it in the provided container. Also checks the container is +// large enough. +// std::size_t size() const noexcept +// Returns the number of encoded elements. +// std::size_t const number_of_frames() const +// Returns the number of frames stored in the Terse object +// std::vector const& dim() const noexcept +// Returns the dimensions of the frames (all frames must have the same dimensions). +// void dim(std::vector const& dim) noexcept +// Sets the dimensions of the Terse frames. Since all frames must have the same dimensions, they can be set only once. +// bool is_signed() const noexcept +// Returns true if the encoded data are signed, false if unsigned. Signed data cannot be decompressed into +// unsigned data. +// unsigned bits_per_val() const noexcept +// Returns the maximum number of bits per element that can be expected. So for uncompressed uint_16 type +// data, bits_per_val() returns 16. Terse data cannot be decompressed into a container type with elements +// that are smaller in bits than bits_per_val(). +// std::size_t terse_size() noexcept +// Returns the number of bytes used for encoding the Terse data. +// void write(std::ostream& ostream) +// Writes Terse data to 'ostream'. The Terse data are preceded by an XML element containing the parameters +// that are required for constructing a Terse object from the stream. Data are written as a byte stream +// and are therefore independent of endian-ness. A small-endian memory lay-out produces the a Terse file +// that is identical to a big-endian machine. +// void shrink_to_fit() noexcept +// Releases unused buffer storage to heap memory. This increases available heap memory when Terse objects were constructed +// from uncompressed data sources held in memory. If compression is performed concurrently, also waits for all compression +// processes to finish. It has no effect when Terse object are read from a stream. +// std::string const& metadata(std::size_t frame = 0) const noexcept +// Returns the metadata that are associated with the specified frame. +// void metadata(std::string data, std::size_t frame = 0) noexcept +// Sets / overwrites any optional metadata that are associated with the specified frame. Metadata are not compressed. +// Metadata will be written to, and read from a stream as part of the Terse object. +// +// Example: +// +// std::vector numbers(1000); // Uncompressed data location +// std::iota(numbers.begin(), numbers.end(), -500); // Fill with numbers -500, -499, ..., 499 +// Terse trpx(numbers); // Compress the data to less than 30% of memory +// std::cout << "compression rate " << float(compressed.terse_size()) / (numbers.size() * sizeof(unsigned)) << std::endl; +// std::ofstream outfile("junk.trpx"); +// compressed.write(outfile); // Write Terse data to disk +// std::ifstream infile("junk.trpx"); +// Terse from_file(infile); // Read it back in again +// std::vector uncompressed(1000); +// from_file.prolix(uncompressed.begin()); // Decompress the data... +// for (int i=0; i != 5; ++i) +// std::cout << uncompressed[i] << std::endl; +// for (int i=995; i != 1000; ++i) +// std::cout << uncompressed[i] << std::endl; +// +// Produces as output: +// +// compression rate 0.29 +// -500 +// -499 +// -498 +// -497 +// -496 +// 495 +// 496 +// 497 +// 498 +// 499 + + +namespace jpa { +template +concept Container = requires(T& c) { + { c.begin() } -> std::input_iterator; + { c.end() } -> std::sentinel_for; + { c.size() } -> std::convertible_to; +// { c.data() } -> std::convertible_to::type*>; +}; + +/** + * @class Terse + * @brief Terse allows efficient and fast compression of integral diffraction data and other integral greyscale data. + * + * The Terse class is used for compression and decompression of integral diffraction data and other integral greyscale data. + * It supports efficient compression and decompression, as well as data extraction and output to streams. + * + * A Terse object may contain compressed data of multiple frames, and data of a particular frame can be + * extracted by indexing. All frames in a Terse object must have the same size and dimensions, and must all be signed or unsigned. + * + * A Terse object can be unpacked into any arithmetic type, including float and double. Unpacking into + * values of type T with fewer bits than the original data results in truncation of overflowed data to + * std::numeric_limits::max() and to std::numeric_limits::min() for underflowed signed types. + * Unpacking signed data into unsigned values is not allowed. Compressing as unsigned yields a tighter compression. + * + * A Terse object can be written or appended to any stream. The Terse data in the file is independent of the endian-ness + * of the machine: both big-endian and small-endian machines produce identical files, making data transfer optimally + * transparent. + * + * Terse data in a file are immediately preceded by a small header, which is encoded in standard XML as follows: + * + * - "n" is the number of bits required for the most extreme value in the Terse data, representing the bit depth of + * the original, uncompressed data. + * - "s" is "0" for unsigned data, "1" for signed data. + * - "b" is the block size of the stretches of data values that are encoded (by default 12 values). + * - "v" is the number of elements of a single frame of a stack (it is the product of the dimensions, if these are provided). + * - "f" is the number of frames encoded in the file. + * - "m" is the total memory size of all frames, excluding the header and metadata. + * - "d [...]" is optional and encodes the dimensions of a single frame. Frames can have any number of dimensions. + * - "ms" is optional and encodes the sizes of metadata strings associated with each frame. + * - "mf" is optional and encodes the memory sizes of individual frames in the stack. + * + * Here is an example of a Terse file with two frames of 512x512 pixels: + * 
+ *       
+ *
+ * + * Example of usage: + * \code{.cpp} + * std::vector numbers(1000); // Uncompressed data location + * std::iota(numbers.begin(), numbers.end(), -500); // Fill with numbers -500, -499, ..., 499 + * Terse trpx(numbers); // Compress the data to less than 30% of memory + * std::cout << "compression rate " << float(compressed.terse_size()) / (numbers.size() * sizeof(unsigned)) << std::endl; + * std::ofstream outfile("junk.trpx"); + * compressed.write(outfile); // Write Terse data to disk + * std::ifstream infile("junk.trpx"); + * Terse from_file(infile); // Read it back in again + * std::vector uncompressed(1000); + * from_file.prolix(uncompressed.begin()); // Decompress the data... + * for (int i = 0; i != 5; ++i) + * std::cout << uncompressed[i] << std::endl; + * for (int i = 995; i != 1000; ++i) + * std::cout << uncompressed[i] << std::endl; + * \endcode + * + * Produces the following output: + * 
+ *    compression rate 0.29
+ *    -500
+ *    -499
+ *    -498
+ *    -497
+ *    -496
+ *    495
+ *    496
+ *    497
+ *    498
+ *    499
+ *
+ */ + +class Concurrent; + +/** + * @enum Terse_mode + * @brief Defines the modes used for encoding data in a Terse object. + * + * The `Terse_mode` enum specifies whether specific optimizations are applied for (small) unsigned values. + * It also allows backwards compatibility with earlier versions of Terse by selecting Signed. + * + * Enumerators: + * - **Signed**: The fastest and most general compression mode: the sign bit is maintained for signed data, + * enabling proper decompression of signed values. Unsigned data is compressed without sign bit, allowing + * tighter compression. Data stored as Signed can be read by earlier versions of Terse. + * - **Unsigned**: Data is treated as unsigned for compression; overloads (values with all bits set) are encoded + * efficiently. Unsigned compression is just as fast as the Signed compression. Earlier versions of Terse cannot + * read Unsigned data. + * - **Small_unsigned**: Optimized encoding mode for small unsigned values, providing tighter compression + * for values that require fewer bits. However, it is a bit slower than Unsigned or Signed compression. Earlier + * versions of Terse cannot read Small_unsigned data. + * + * Usage example: + * \code{.cpp} + * Terse_mode mode = Terse_mode::Signed; // Specify signed mode for data compression. + * \endcode + */ +enum class Terse_mode { + Unsigned, + Small_unsigned, + Signed, + Default +}; + +template +class Terse { + static_assert(std::is_same_v || std::is_same_v, "Terse is either Concurrent or not"); +#ifndef Concurrent_h + static_assert(std::is_same_v, "Concurrent.hpp must be included before #include \"Terse.hpp\""); +#endif + + template friend class Terse; + +public: + /** + * @brief Initializes an empty Terse object. + * + * Data can be appended to an empty Terse object. The first dataset to be pushed in + * determines size, signedness and dimensions of the remaining datasets that can be pushed in. + */ + Terse() noexcept {}; + + /** + * @brief Creates a Terse object from data (which can be a std::vector, Field, etc.). + * Only containers of integral types are allowed. If the container has a member function dim(), + * that will set the dimensions of the Terse object. Otherwise, the dimensions can be set once + * using the dim(vector const&) member function. + * If the input container is passed as an r-value, it will be consumed, and will be empty after the call. + * If the input container is passed as an r-value, and the include file Concurrent.hpp was included before + * the Terse.hpp file, compression will be performed concurrently in the background. + * + * @tparam C The type of the container containing integral data. + * @param data The container containing integral data. + */ + template + Terse(C&& data, Terse_mode const mode = Terse_mode::Default) noexcept { + insert(0, std::forward(data), mode); + }; + + /** + * @brief Creates a Terse object given a starting iterator or pointer and the number of elements that need to be encoded. + * + * @tparam Iterator The type of the iterator. + * @param data The starting iterator or pointer to the data. + * @param size The number of elements to be encoded. + * @param block The block size for compression (default is 12). + */ + template + Terse(Iterator const data, size_t const size, unsigned int const block=12, Terse_mode mode = Terse_mode::Default) noexcept : + d_signed(std::is_signed_v::value_type>), + d_block(block), + d_size(size) { + push_back(data, size, mode); + } + + /** + * @brief Reads in a Terse object that has been written to a stream . + * + * Scans the stream for the Terse XML header, then reads the binary Terse data, leaving the stream position exactly one byte beyond the binary Terse data. + * + * @tparam STREAM The type of the input stream where the Terse data are stored (e.g. a file stream or a string stream). + * @param istream The input stream containing Terse data. + */ + template requires std::derived_from + Terse(STREAM& istream) : Terse(istream, XML_element(istream, "Terse")) {}; + + /** + * @brief Adds another frame to the Terse object. The new frame is defined by its begin iterator and size. + * + * The size must be the same as that of the first frame that was used for creating the Terse object. + * + * @tparam Iterator The type of the iterator. + * @param pos The location where the data need to be inserted. + * @param data The starting iterator or pointer to the data. + * @param size The number of elements to be encoded. + */ + template + void insert(std::size_t const pos, Iterator const data, size_t const size, Terse_mode mode = Terse_mode::Default) noexcept { + insert(pos, std::span(data, data + static_cast(size)), mode); + } + + /** + * @brief Inserts a frame into the Terse object. The new frame is defined by a reference to a container. + * The size and dimensions of the input container must be the same as that of the first frame that was used + * for creating the Terse object, unless the Terse object is empty, in which case the inserted container determines + * the size and dimension of subsequent frames that are added. + * If the input container is passed as an r-value, it will be consumed, and will be empty after the call. + * If the input container is passed as an r-value, and the include file Concurrent.hpp was included before + * the Terse.hpp file, compression will be performed concurrently in the background. + * + * @tparam C The type of the container containing integral data. + * @param pos The location where the data need to be inserted. + * @param data The container containing integral data. + * @throws std::invalid_argument If the dimensions of this Terse object and the provided frame differ. + */ + template + void insert(std::size_t const pos, C&& data, Terse_mode mode = Terse_mode::Default) { + bool dim_ok = true; + if constexpr (requires(C& c) { c.dim(); }) { + for (std::size_t i = 0; i != data.dim().size(); ++i) + if (number_of_frames() == 0) + d_dim.push_back(static_cast(data.dim()[i])); + else + dim_ok = dim_ok && d_dim[i] == static_cast(data.dim()[i]); + } + else if constexpr(requires (C &c) {c.request().shape;}) { + auto shape = data.request().shape; + if (number_of_frames() == 0) + d_dim = std::vector(shape.begin(), shape.end()); + else + dim_ok = dim_ok && d_dim == std::vector(shape.begin(), shape.end()); + } + if (!dim_ok) + throw(std::invalid_argument("The provided container and the requested frame have different dimensions")); + using T = std::remove_cv_t>; + if (std::is_signed_v) + mode = Terse_mode::Signed; + d_prolix_bits = std::max(d_prolix_bits, 8 * static_cast(sizeof(T))); + if (number_of_frames() == 0) { + d_size = data.size(); + d_signed = std::is_signed_v; + } + d_metadata.insert(d_metadata.begin() + static_cast(pos), ""); + auto at = d_terse_frames.begin() + static_cast(pos); + if constexpr (std::is_lvalue_reference_v) + d_terse_frames.insert(at, f_compress(mode, data.data())); + else if constexpr (std::is_same_v) + d_terse_frames.insert(at, d_concurrent->background([this, d = std::move(data), mode] { return f_compress(mode, d.data()); })); + else { + d_terse_frames.insert(at, f_compress(mode, data.data())); + auto local_data = std::move(data); + } + } + + /** + * @brief Insert a (potentially multiframe) Terse object at the specified position. + * The size and dimensions of both Terse objects must be the same as that of the first frame that was used + * for creating the Terse object, unless the Terse object is empty, in which case the inserted container determines + * the size and dimension of subsequent frames that are added. + * + * @tparam T Either void or Concurrent. + * @param pos The location where the data need to be inserted. + * @param trs The Terse object to be appended. + * @throws std::out_of_range If the provided Terse object index is greater than or equal to the number of frames. + * @throws std::invalid_argument If the provided Terse object differs in signedness, size, block size of dimensionality from this Terse object. + */ + template + void insert(std::ptrdiff_t const pos, Terse& trs) { + if (number_of_frames() == 0) { + d_signed = trs.d_signed; + d_block = trs.d_block; + d_size = trs.d_size; + d_dim = trs.d_dim; + } + if (static_cast(pos) > number_of_frames()) throw std::out_of_range("Frame index is out of range."); + if (trs.d_dim != d_dim) + throw(std::invalid_argument("Dimension mismatch of the provided Terse object")); + if (trs.d_signed != d_signed) + throw(std::invalid_argument("Sign mismatch of the provided Terse object")); + if (trs.d_block != d_block) + throw(std::invalid_argument("Blocksize mismatch of the provided Terse object")); + if (trs.d_size != d_size) + throw(std::invalid_argument("Size mismatch of the provided Terse object")); + d_metadata.insert(d_metadata.begin() + pos, trs.d_metadata.begin(), trs.d_metadata.end()); + shrink_to_fit(); + trs.shrink_to_fit(); + if constexpr (std::is_same_v) + d_terse_frames.insert(d_terse_frames.begin() + pos, trs.d_terse_frames.begin(), trs.d_terse_frames.end()); + else { + std::size_t i_end = trs.d_terse_frames.size(); + for (size_t i = 0; i < i_end; ++i) { + auto& frame = trs.d_terse_frames[i]; + if (std::holds_alternative>(frame)) + d_terse_frames.insert(d_terse_frames.begin() + pos + i, std::get>(frame)); + else + d_terse_frames.insert(d_terse_frames.begin() + pos + i, std::get>>(frame).get()); + } + } + } + + /** + * @brief Appends a frame to the Terse object. + * + * @tparam Iterator The type of the iterator pointing to the data. + * @param it The iterator pointing to the data. + * @param size The number of elements in the frame (must be equal to the size of previously inserted frames). + */ + template + void push_back(Iterator const it, size_t const size, Terse_mode mode = Terse_mode::Default) noexcept { + insert(d_terse_frames.size(), it, size, mode); + } + + /** + * @brief Appends a frame to the Terse object. The new frame is defined by a reference to a container. + * The size and dimensions of the input container must be the same as that of the first frame that was used + * for creating the Terse object, unless the Terse object is empty, in which case the inserted container determines + * the size and dimension of subsequent frames that are added. + * If the input container is passed as an r-value, it will be consumed, and will be empty after the call. + * If the input container is passed as an r-value, and the include file Concurrent.hpp was included before + * the Terse.hpp file, compression will be performed concurrently in the background. + * + * @tparam C The type of the container containing integral data. + * @param data The container containing integral data. + */ + template + void push_back(C&& data, Terse_mode mode = Terse_mode::Default) noexcept { + insert(number_of_frames(), std::forward(data), mode); + } + + /** + * @brief Appends a (potentially multiframe) Terse object. + * The size and dimensions of bothe Terse objects must be the same as that of the first frame that was used + * for creating the Terse object, unless the Terse object is empty, in which case the inserted container determines + * the size and dimension of subsequent frames that are added. + * + * @tparam T Either void or Concurrent. + * @param trs The Terse object to be appended. + */ + template + void push_back(Terse& trs) noexcept { insert(static_cast(number_of_frames()), trs); } + + /** + * @brief Removes one of the frames from the Terse object. Also waits until concurrent compression has finished, and release unused storage to the heap. + * + * @param i The index of the frame to be removed. + */ + void erase(std::ptrdiff_t i) noexcept { + shrink_to_fit(); + d_metadata.erase(d_metadata.begin() + i); + d_terse_frames.erase(d_terse_frames.begin() + i); + } + + /** + * @brief Returns a selected frame as a Terse object. + * + * @param pos The index of the selected frame. + * @throws std::out_of_range If the provided frame index is greater than or equal to the number of frames. + */ + Terse at(std::size_t pos) { + if (pos >= number_of_frames()) throw std::out_of_range("Frame index is out of range."); + Terse result; + result.d_signed = d_signed; + result.d_block = d_block; + result.d_size = d_size; + result.d_prolix_bits = d_prolix_bits; + result.d_dim = d_dim; + result.d_metadata.push_back(d_metadata[pos]); + result.d_terse_frames.push_back(f_get_frame(pos)); + return result; + } + + /** + * @brief Unpacks the Terse data of the requested frame and stores it in the provided container. + * + * @tparam C The type of the container. + * @param container The container where the data will be stored. + * @param frame The index of the frame to unpack (default is 0). + * @throws std::invalid_argument If the provided container does not have enough space to unpack all frames, + * or the number of dimensions of the requested frame and the provided container differ. + */ + template + C&& prolix(C&& container, std::size_t frame) { + if (this->size() != container.size()) + throw(std::invalid_argument("The provided container not have enough space for the requested frame of the Terse object")); + bool dim_ok = true; + if constexpr(requires (C &c) {c.dim();}) + for (int i = 0; i != d_dim.size(); ++i) + dim_ok = dim_ok && d_dim[i] == container.dim()[i]; + else if constexpr(requires (C &c) {c.request().shape;}) + dim_ok = dim_ok && d_dim == container.request().shape; + if (!dim_ok) + throw(std::invalid_argument("The provided container and the requested frame have different dimensions")); + prolix(container.data(), frame); + return std::forward(container); + } + + /** + * @brief Unpacks all the Terse data and stores these consecutively in the provided container of + * numerical values. + * + * @tparam C The type of the container. + * @param container The container where the data will be stored. + * @throws std::invalid_argument If the provided container does not have enough space to unpack all frames. + */ + template requires std::is_arithmetic_v::value_type> + C&& prolix(C&& container) { + if (size() == container.size() && number_of_frames() != 0) + return prolix(std::forward(container), 0); + if (size() * number_of_frames() != container.size()) + throw(std::invalid_argument("The provided container not have enough space for unpacking all frames of the Terse object")); + auto* data_ptr = [&]() { + if constexpr (requires (C &c) { c.mutable_data(); }) return container.mutable_data(); + else return container.data(); + }(); + if constexpr (std::is_same_v) + for (std::size_t i = 0; i != number_of_frames(); ++i) + prolix(data_ptr + i * size(), i); + else { + std::vector> futures; + for (std::size_t i = 0; i != number_of_frames(); ++i) + futures.push_back(d_concurrent->background([this, i, &data_ptr] { + prolix(data_ptr + i * size(), i); + })); + for (auto& future : futures) future.get(); + } + return std::forward(container); + } + + /** + * @brief Unpacks all the Terse data and stores it in the provided container of containers. + * Each frame is stored separately in each sub-container. + * + * @tparam C The type of the container. + * @param container The container where the data will be stored. + * @throws std::invalid_argument If the provided container does not have enough sub-containers to unpack all frames. + */ + template requires Container::value_type> + C&& prolix(C&& container) { + if (number_of_frames() == 1) + prolix(container[0], 0); + else { + if (container.size() < number_of_frames()) + throw(std::invalid_argument("The provided container does not have enough subcontainers for unpacking all frames of the Terse object")); + if constexpr (std::is_same_v) + for (std::size_t i = 0; i != number_of_frames(); ++i) + prolix(container[i], i); + else { + std::vector> futures; + for (std::size_t i = 0; i != number_of_frames(); ++i) + futures.push_back(d_concurrent->background([this, &container, i] { + prolix(container[i], i); + })); + for (auto& future : futures) future.get(); + } + } + return std::forward(container); + } + + /** + * @brief Unpacks the Terse data, storing the unpacked data from the location defined by 'begin'. + * + * Unpacks with bounds checking. + * + * @tparam Iterator The type of the iterator. + * @param begin The starting iterator or pointer where the data will be stored. + * @param frame The index of the frame to unpack (default is 0). + * @throws std::out_of_range If the provided frame index is greater than or equal to the number of frames. + * @throws std::invalid_argument If the iterator refers unsigned values, when the Terse object contains signed values. + */ + template requires requires (Iterator& i) {*i;} + void prolix(Iterator begin, std::size_t frame = 0) { + if (frame >= number_of_frames()) throw std::out_of_range("Frame index is out of range."); + if (is_signed() && std::unsigned_integral::value_type>) + throw std::invalid_argument("Cannot decompress signed data into an unsigned container."); + switch (d_block) { + case(8) : return f_prolix<8> (begin, frame); + case(9) : return f_prolix<9> (begin, frame); + case(10) : return f_prolix<10> (begin, frame); + case(11) : return f_prolix<11> (begin, frame); + case(12) : return f_prolix<12> (begin, frame); + case(13) : return f_prolix<13> (begin, frame); + case(14) : return f_prolix<14> (begin, frame); + case(15) : return f_prolix<15> (begin, frame); + case(16) : return f_prolix<16> (begin, frame); + case(20) : return f_prolix<20> (begin, frame); + case(24) : return f_prolix<24> (begin, frame); + case(32) : return f_prolix<32> (begin, frame); + default : return f_prolix<0> (begin, frame); + } + } + + /** + * @brief Returns the number of encoded elements of a single frame (all frames in a Terse object must have the + * same number of elements). + * + * @return The number of encoded elements. + */ + std::size_t size() const noexcept { return d_size; } + + /** + * @brief Returns the number of frames stored in the Terse object. + * + * @return The number of frames stored in the Terse object. + */ + std::size_t number_of_frames() const noexcept { return d_terse_frames.size(); } + + /** + * @brief Returns the dimensions of each of the Terse frames (all frames in a Terse object must have the same dimensions). + * + * @return The dimensions of each of the Terse frames. + */ + std::vector const& dim() const noexcept { return d_dim; } + + /** + * @brief Sets the dimensions of the Terse frames. All frames must have the same dimensions. If the Terse object is empty, + * the dimensions can be set to any value. Otherwise, the the total number of elements defined per frame is fixed. + * + * @param dim The dimensions to set for the Terse frames. + * @throws std::invalid_argument If the provided dimensions differ from dimensions that were set earlier. + */ + void dim(std::vector const& dim) { + if (d_dim.size() == 0 || size() == std::accumulate(dim.begin(), dim.end(), 1ul, std::multiplies<>())) + d_dim = dim; + if (dim != d_dim) + throw std::invalid_argument("Dimensions already set."); + } + + /** + * @brief Returns the block size that is used for compression. + * + * @return The block size that is used for compression. + */ + std::size_t const& block_size() const noexcept { return d_block; } + + /** + * @brief Sets the the block size used for compression, provided no frames have been compressed yet. + * + * @param block The block size that is used for compression. If the Terse object is not empty, the block size + * cannot be changed, and the original block size is returned. + */ + void block_size(std::size_t const& block) noexcept { + d_block = number_of_frames() == 0 ? block : d_block; + } + + /** + * @brief Returns the degree of partiality of a Terse object, or zero for Terse.. + * + * @return The degree of partiality (0 for sequential compression, 1 for using all cores). + */ + double dop() const noexcept { + if constexpr (std::is_same_v) + return d_concurrent->dop(); + else + return 0; + } + + /** + * @brief Sets the degree of partiality for a Terse object. + * + * @param new_dop The new degree of partiality (0 for sequential compression, 1 for using all cores). + */ + void dop(double new_dop) noexcept requires std::is_same_v { + for (std::size_t i = 0; i!= d_terse_frames.size(); ++i) + f_get_frame(i); + d_concurrent = Concurrent(new_dop); + } + + /** + * @brief Returns the fractional precision for lossy floating point compression. + * + * @return The fractional precision for lossy floating point compression; a number between 1 (loss of all information) + * and 0 (lossless floating point compression). For fractional precsion of 0.01 indicates a precision of 1%. + */ + double fractional_precision() const noexcept { return static_cast(1) / (static_cast(1) << (d_binary_precision - 1)); } + + /** + * @brief Sets the fractional precision for lossy floating point compression. + * + * @param frac The fractional precision for lossy floating point compression; a number between 1 (loss of all information) + * and 0 (lossless floating point compression). For fractional precision of 0.01 indicates a precision of 1%. + */ + void fractional_precision(double frac) { + auto precision = 1 + std::ceil(std::log2(1 / frac)); + d_binary_precision = static_cast(std::min(precision, 54.0)); + } + + /** + * @brief Checks if at least one of the frames in the Terse object has floating point values. + * + * @return True any of the frames has floating point values. + */ + bool is_float() { + for (std::size_t i = 0; i != d_terse_frames.size(); ++i) + if (is_float(i)) return true; + return false; + } + + /** + * @brief Checks if the frame with the specified index has floatinging point values. + * + * @return True if the frame with the specified index has floating point values + */ + bool is_float(std::size_t frame_num) { return (Bitqueue_pop(f_get_frame(frame_num)).pop<18, std::size_t>() == 0b111111111111111010); } + + /** + * @brief Returns if the default mode of compression for unsigned data is set to Terse_mode::Small_unsigned. + * + * @return true if the default mode of compression for unsigned data is set to Terse_mode::Small_unsigned. + */ + bool small() const noexcept { return d_small; } + + /** + * @brief Sets/resets the default mode of compression for unsigned to Terse_mode::Small_unsigned. + * + * @param val true: sets default mode for unsigned compression to Terse_mode::Small_unsigned; false sets it to Terse_mode::Unsigned + */ + void small(bool val) { d_small = val; } + + /** + * @brief Returns if the default mode of compression for unsigned data is set to Terse_mode::Unsigned. + * + * @return true if the default mode of compression for unsigned data is set to Terse_mode::Unsigned. + */ + bool fast() const noexcept { return !small(); } + + /** + * @brief Sets/resets the default mode of compression for unsigned to Terse_mode::Unsigned. + * + * @param val true: sets default mode for unsigned compression to Terse_mode::Unsigned; false sets it to Terse_mode::Small_unsigned + */ + void fast(bool val) { small(!val); } + + /** + * @brief Returns true if the encoded data are signed, false if unsigned. Signed data cannot be decompressed into unsigned data. + * + * @return True if the encoded data are signed, false otherwise. + */ + bool is_signed() const noexcept {return d_signed;} + + /** + * @brief Returns the bit depth of the data before compression. + * + * @return The maximum number of bits per element. + */ + unsigned bits_per_val() const noexcept {return d_prolix_bits;} + + /** + * @brief Returns the number of bytes used for encoding the Terse data. + * + * Mainly useful for reporting compression rates. + * + * @return The number of bytes used for encoding the Terse data. + */ + std::size_t terse_size() noexcept { + std::size_t memsize = 0; + for (std::size_t i = 0; i!= d_terse_frames.size(); ++i) + memsize += f_get_frame(i).size(); + return memsize; + } + + /** + * @brief Returns the number of bytes that a Terse file of the object would contain. + * + * Mainly useful for codecs that use Terse. + * + * @return The number of bytes used for storing the Terse file. + */ + std::size_t file_size() noexcept { + if (number_of_frames() == 0) + return 0; + std::ostringstream oss; + f_write_metadata(oss); + return oss.str().size() + terse_size(); + } + + /** + * @brief Write the Terse object to the specified output stream. + * + * This member function first writes a small XML header with data required to unpack the Terse object, then writes + * any provided Terse metadata, then writes the compressed data to the provided output stream. + * + * @param ostream The output stream to which the Terse data will be written. + */ + void write(std::ostream& ostream) { + if (number_of_frames() == 0) return; + f_write_metadata(ostream); + for (std::size_t i = 0; i!= d_terse_frames.size(); ++i) + ostream.write(reinterpret_cast(f_get_frame(i).data()), static_cast(f_get_frame(i).size())); + ostream.flush(); + } + + /** + * @brief Releases unused buffer storage to heap memory. This increases available heap memory when Terse objects were constructed + * from uncompressed data sources held in memory. If compression is performed concurrently, also waits for all compression + * processes to finish. It has no effect when Terse object are read from a stream. + */ + void shrink_to_fit() noexcept { + for (std::size_t i = 0; i!= d_terse_frames.size(); ++i) + f_get_frame(i).shrink_to_fit(); + } + + /** + * @brief Sets / overwrites any optional metadata that are associated with the specified frame. Metadata are not compressed. + * Metadata will be written to, and read from a stream as part of the Terse object. + * + * @param data The metadata as a string object. + * @param pos The number of the frame pertaining to the metadata. + * @throws std::out_of_range If the provided frame index is greater than or equal to the number of frames. + */ + void metadata(std::size_t const pos, std::string const data) { + if (pos >= number_of_frames()) throw std::out_of_range("Frame index is out of range."); + d_metadata[pos] = data; + } + + void metadata(std::string data) { + if (number_of_frames() == 0) throw std::invalid_argument("Cannot add metadata, as this Terse object has no frames."); + d_metadata[0] = data; + } + + /** + * @brief Returns the metadata that are associated with the specified frame. + * + * @param frame The number of the frame pertaining to the metadata. + * @return A constant reference to the metadata as a std::string. If no metadata are available, an empty string is returned. + */ + std::string const& metadata(std::size_t frame = 0) const noexcept { + return d_metadata[frame]; + } + +private: + using FrameStorage = std::conditional_t< + std::is_same_v, + std::vector>, std::vector>>, + std::vector>>; + + FrameStorage d_terse_frames; + bool d_signed; + bool d_small = true; + std::size_t d_block = 12; + std::size_t d_size = 0; + unsigned d_prolix_bits = 0; + std::uint8_t d_binary_precision = 54; + std::vector d_dim; + std::vector d_metadata; + std::optional d_concurrent = []() -> std::optional { + if constexpr (std::is_same_v) return Concurrent(1); + else return std::nullopt; + }(); + + template requires std::derived_from + Terse(STREAM& istream, XML_element const& xmle) { + if (!istream.eof()) { + d_prolix_bits = unsigned(std::stoul(xmle.attribute("prolix_bits"))); + d_signed = std::stoul(xmle.attribute("signed")); + d_block = std::stoul(xmle.attribute("block")); + d_size = std::stoull(xmle.attribute("number_of_values")); + std::istringstream dim_str(xmle.attribute("dimensions")); + unsigned int d; + while (dim_str >> d) + d_dim.push_back(d); + if (xmle.attribute("metadata_string_sizes") != "") { + std::istringstream meta_str(xmle.attribute("metadata_string_sizes")); + unsigned int val; + while (meta_str >> val) { + d_metadata.push_back(std::string(val, ' ')); + istream.read(reinterpret_cast(&d_metadata.back()), val); + } + } + d_terse_frames.resize(std::stoull(xmle.attribute("number_of_frames"))); + for (auto& frame : d_terse_frames) + frame = std::vector(); // Initialize each element to an empty vector + if (xmle.attribute("memory_sizes_of_frames") == "") + f_fill_terse_frames(istream, std::stoul(xmle.attribute("memory_size"))); + else { + std::istringstream frame_sizes_str(xmle.attribute("memory_sizes_of_frames")); + unsigned int val; + for (std::size_t i = 0; i != d_terse_frames.size(); ++i) { + frame_sizes_str >> val; + f_get_frame(i).resize(val); + istream.read(reinterpret_cast(&(f_get_frame(i)[0])), static_cast(f_get_frame(i).size())); + } + } + } + } + + void f_write_metadata(std::ostream& ostream) { + std::size_t memory_size = terse_size(); + XML_element xml(""); + xml.add_attribute("prolix_bits", d_prolix_bits); + xml.add_attribute("signed", d_signed); + xml.add_attribute("block", d_block); + xml.add_attribute("number_of_values", size()); + if (!d_dim.empty()) xml.add_attribute("dimensions", d_dim); + xml.add_attribute("number_of_frames", d_terse_frames.size()); + std::vector frame_sizes; + for (size_t i = 0; i != d_terse_frames.size(); ++i) + frame_sizes.push_back(f_get_frame(i).size()); + xml.add_attribute("memory_sizes_of_frames", frame_sizes); + xml.add_attribute("memory_size", memory_size); + if (!d_metadata.empty()) { + std::vector metadata_sizes; + for (size_t i = 0; i != d_metadata.size(); ++i) + metadata_sizes.push_back(d_metadata[i].size()); + xml.add_attribute("metadata_string_sizes", metadata_sizes); + } + ostream << xml.XML(); + for (auto& str : d_metadata) ostream << str; + ostream.flush(); + } + + std::vector& f_get_frame(std::size_t index) noexcept { + if constexpr (std::is_same_v) + return d_terse_frames[index]; + else { + if (std::holds_alternative>>(d_terse_frames[index])) + d_terse_frames[index] = std::get>>(d_terse_frames[index]).get(); + return std::get>(d_terse_frames[index]); + } + } + + template + auto f_compress(Terse_mode mode, Iterator const data_begin) { + if constexpr (std::is_signed_v>) + return f_compress(data_begin); + else switch (mode) { + case Terse_mode::Unsigned: return f_compress(data_begin); + case Terse_mode::Small_unsigned: return f_compress(data_begin); + case Terse_mode::Signed: return f_compress(data_begin); + default: return d_small ? f_compress(data_begin) : f_compress(data_begin); + } + } + +#define PUSH_BITS(prevbits, bits) \ + if (prevbits == bits) bitqueue.push_back<1>(0b1); \ + else if (bits < 7) bitqueue.push_back<4>(bits << 1); \ + else if (bits < 10) bitqueue.push_back<6>(0b1110 + ((bits - 7) << 4)); \ + else bitqueue.push_back<12>(0b111110 + ((bits - 10) << 6)); \ + prevbits = bits; + + template requires ((MODE == Terse_mode::Signed || MODE == Terse_mode::Unsigned) && + std::integral::value_type>) + std::vector const f_compress(Iterator data) noexcept { + using T = std::iterator_traits::value_type; + std::size_t terse_frame_size = static_cast(d_size * sizeof(decltype(*data)) * 0.01) + 2 * (d_block * sizeof(decltype(*data)) + 4); + terse_frame_size = (terse_frame_size + 7) & ~std::size_t(7); + std::vector terse_frame(terse_frame_size); + Unique_array> buffer(d_block); + Bitqueue_push_back bitqueue(terse_frame); + if constexpr (MODE != Terse_mode::Signed) + bitqueue.push_back<18>(0b111111111111111000); + std::size_t prevbits = 0; + std::size_t prevmasked_bits = 0; + for (std::size_t from = 0; from < d_size; from += d_block) { + std::size_t index = static_cast(bitqueue.data() - terse_frame.data()); + if (index > terse_frame.size() - d_block * sizeof(decltype(*data)) - 2) { + terse_frame_size = static_cast(1.1 * terse_frame.size() * d_size / from); + terse_frame_size = (terse_frame_size + 7) & ~std::size_t(7); + terse_frame.resize(terse_frame_size); + bitqueue.relocate(terse_frame.data() + index); + } + std::span data_block(data + static_cast(from), std::min(d_size - from, d_block)); + std::uint8_t significant_bits = f_most_significant_bit(data_block); + PUSH_BITS(prevbits, significant_bits); + if constexpr (MODE == Terse_mode::Signed) + bitqueue.push_back(significant_bits, data_block); + else { + if (significant_bits != sizeof(T) * 8) + bitqueue.push_back(significant_bits, data_block); + else { // masked data + std::transform(data_block.begin(), data_block.end(), buffer.begin(), [](T val) {return val + 1;}); + std::uint8_t masked_bits = f_most_significant_bit(std::span(buffer.begin(), buffer.end())); + PUSH_BITS(prevmasked_bits, masked_bits); + bitqueue.push_back(masked_bits, buffer); + } + } + } + terse_frame_size = static_cast(bitqueue - terse_frame.data()); + terse_frame_size = (terse_frame_size + 7) & ~std::size_t(7); + terse_frame.resize(terse_frame_size); + return terse_frame; + } + + template requires (std::floating_point::value_type>) + std::vector const f_compress(Iterator data) noexcept { + using T = std::iterator_traits::value_type; + std::uint8_t binary_precision = std::min(d_binary_precision, static_cast(std::numeric_limits::digits +1)); + std::size_t terse_frame_size = static_cast(d_size * sizeof(T) * 0.01) + 2 * (d_block * sizeof(decltype(*data)) + 4); + terse_frame_size = (terse_frame_size + 7) & ~std::size_t(7); + std::vector terse_frame(terse_frame_size); + Unique_array> buffer(d_block); + Bitqueue_push_back bitqueue(terse_frame); + bitqueue.push_back<18>(0b111111111111111010); + bitqueue.push_back<6>(binary_precision); + std::size_t prevbits_exponents = 0; + Unique_array mantissas(d_block); + Unique_array exponents(d_block); + std::size_t mantissa_bits = (static_cast(1) << (binary_precision - 1)); + for (std::size_t from = 0; from < d_size; from += d_block) { + std::size_t index = static_cast(bitqueue.data() - terse_frame.data()); + if (index > terse_frame.size() - d_block * sizeof(decltype(*data)) - 2) { + terse_frame_size = static_cast(1.1 * terse_frame.size() * d_size / from); + terse_frame_size = (terse_frame_size + 7) & ~std::size_t(7); + terse_frame.resize(terse_frame_size); + bitqueue.relocate(terse_frame.data() + index); + } + std::span data_block(data + static_cast(from), std::min(d_size - from, d_block)); + bool unsigned_block = true; + for (std::size_t i=0; i != data_block.size(); ++i) { + T mantissa_float = std::frexp(data_block[i], &exponents[i]); + unsigned_block = unsigned_block && (mantissa_float >= 0); + if (!std::isnan(mantissa_float)) + mantissas[i] = static_cast(std::round(mantissa_float * mantissa_bits)); + else { + mantissas[i] = 0; + exponents[i] = 1 << 12; + } + } + bitqueue.push_back<1>(unsigned_block); + std::uint8_t significant_bits_exponents = f_most_significant_bit(std::span(exponents.begin(), data_block.size())); + PUSH_BITS(prevbits_exponents, significant_bits_exponents); + bitqueue.push_back(binary_precision - unsigned_block, std::span(mantissas.begin(), data_block.size())); + bitqueue.push_back(significant_bits_exponents, std::span(exponents.begin(), data_block.size())); + } + terse_frame_size = static_cast(bitqueue - terse_frame.data()); + terse_frame_size = (terse_frame_size + 7) & ~std::size_t(7); + terse_frame.resize(terse_frame_size); + return terse_frame; + } + +#undef PUSH_BITS + + template requires std::is_unsigned_v + constexpr void f_compress_weak_block(std::span const data_block, + Bitqueue_push_back& bitqueue, + T const maxval, + T& previous_maxval) { + if (previous_maxval == 0 && maxval == 0) + bitqueue.push_back<1>(0b1); + else if (previous_maxval == maxval) + bitqueue.push_back<2>(0b11); + else if (previous_maxval + 1 == maxval) + bitqueue.push_back<2>(0b10); + else if (maxval != 6 && previous_maxval - 1 == maxval) + bitqueue.push_back<2>(0b01); + else if (previous_maxval == 6 && maxval == 4) + bitqueue.push_back<2>(0b10); + else + bitqueue.push_back<5>(maxval << 2); + switch (maxval) { + case 0: break; + case 1: bitqueue.push_back<1>(data_block); break; + case 3: bitqueue.push_back<2>(data_block); break; + default: + std::size_t mult = 1; + std::size_t compact = 0; + for (std::size_t i = 0; i != data_block.size(); ++i) { + compact += mult * data_block[i]; + mult *= (maxval + 1); + } + bitqueue.push_back(f_most_significant_bit(mult - 1), compact); + } + previous_maxval = maxval; + } + + template + constexpr void f_compress_strong_block(std::span const data_block, + Bitqueue_push_back& bitqueue, + std::uint8_t const significant_bits, + std::size_t& previous_significant_bits) { + if (previous_significant_bits == significant_bits) + bitqueue.push_back<2>(0b11); + else if (previous_significant_bits + 1 == significant_bits) + bitqueue.push_back<2>(0b10); + else if (previous_significant_bits - 1 == significant_bits) + bitqueue.push_back<2>(0b01); + else if (significant_bits < 10) + bitqueue.push_back<8>(0b11100 + ((significant_bits - 3) << 5)); + else if (significant_bits < 17) + bitqueue.push_back<11>(0b11111100 + ((significant_bits - 10) << 8)); + else + bitqueue.push_back<17>(0b11111111100 + ((significant_bits - 17) << 11)); + bitqueue.push_back(significant_bits, data_block); + previous_significant_bits = significant_bits; + } + + template requires (MODE == Terse_mode::Small_unsigned) + std::vector const f_compress(Iterator data) noexcept { + using T = std::iterator_traits::value_type; + static_assert(std::is_unsigned_v, "Cannot compress signed data with Terse_mode::small"); + //std::size_t const block = std::min(d_block, 24ul); + std::size_t block = std::min(d_block, std::size_t(24)); + std::size_t terse_frame_size = static_cast(d_size * sizeof(decltype(*data)) * 0.01) + 2 * (d_block * sizeof(decltype(*data)) + 4); + terse_frame_size = (terse_frame_size + 7) & ~std::size_t(7); + std::vector terse_frame(terse_frame_size); + Bitqueue_push_back bitqueue(terse_frame); + bitqueue.push_back<18>(0b111111111111111100); + T prevmax = 0; + std::size_t prevbits = 0; + for (std::size_t from = 0; from < d_size; from += block) { + std::size_t index = static_cast(bitqueue.data() - terse_frame.data()); + if (index > terse_frame.size() - block * sizeof(decltype(*data)) - 2) { + terse_frame_size = static_cast(1.1 * terse_frame.size() * d_size / from); + terse_frame_size = (terse_frame_size + 7) & ~std::size_t(7); + terse_frame.resize(terse_frame_size); + bitqueue.relocate(terse_frame.data() + index); + } + std::span data_block(data + static_cast(from), std::min(d_size - from, d_block)); + T max = std::ranges::max(data_block); + if (max < 7) { + f_compress_weak_block(data_block, bitqueue, max, prevmax); + prevbits = 65; + } + else { + std::uint8_t const significant_bits = f_most_significant_bit(max); + if (significant_bits == sizeof(T) * 8) + f_compress_masked(data, terse_frame, from, bitqueue, max, prevmax, prevbits); + else { + f_compress_strong_block(data_block, bitqueue, significant_bits, prevbits); + prevmax = std::numeric_limits::max() / 2; + } + } + } + terse_frame_size = static_cast(bitqueue - terse_frame.data()); + terse_frame_size = (terse_frame_size + 7) & ~std::size_t(7); + terse_frame.resize(terse_frame_size); + return terse_frame; + } + + template + void f_compress_masked(Iterator data, + std::vector& terse_frame, + std::size_t& from, + Bitqueue_push_back& bitqueue, T& max, T& prevmax, std::size_t& prevbits) noexcept { + //std::size_t const block = std::min(d_block, 24ul); + std::size_t block = std::min(d_block, std::size_t(24)); + Unique_array buffer(block); + if constexpr (sizeof(T) * 8 < 10) + bitqueue.push_back<8>(0b11100 + ((sizeof(T) * 8 - 3) << 5)); + else if constexpr (sizeof(T) * 8 < 17) + bitqueue.push_back<11>(0b11111100 + ((sizeof(T) * 8 - 10) << 8)); + else + bitqueue.push_back<17>(0b11111111100 + ((sizeof(T) * 8 - 17) << 11)); + prevmax = std::numeric_limits::max(); + prevbits = sizeof(T) * 8 + 2; + for ( ; from < d_size; from += block) { + std::size_t index = static_cast(bitqueue.data() - terse_frame.data()); + if (index > terse_frame.size() - block * sizeof(decltype(*data)) - 2) { + terse_frame.resize(static_cast(1.1 * terse_frame.size() * d_size / from)); + bitqueue.relocate(terse_frame.data() + index); + } + std::size_t to = std::min(d_size, from + block); + std::transform(data + static_cast(from), data + static_cast(to), buffer.begin(), [](T val) {return val + 1;}); + std::span unmasked(buffer.begin(), buffer.begin() + static_cast(to - from)); + max = std::ranges::max(unmasked); + if (max < 7) + f_compress_weak_block(unmasked, bitqueue, max, prevmax); + else + f_compress_strong_block(unmasked, bitqueue, f_most_significant_bit(max), prevbits); + if (to != d_size) { + std::span test_for_masked(data + static_cast(to), std::min(d_size - to, block)); + if (std::ranges::max(test_for_masked) != std::numeric_limits::max()) { + bitqueue.push_back<1>(0); // End masking + break; + } + bitqueue.push_back<1>(1); + } + } + } + +#define POP(T, bits, begin, num) \ + do { if constexpr (N == 0) \ + bitqueue.pop(bits, std::span(begin, begin + num)); \ + else \ + switch (bits + (num != static_cast(d_block) ? 9 : 0)) { \ + case 0: bitqueue.template pop<0>(std::span(begin, N)); break; \ + case 1: bitqueue.template pop<1>(std::span(begin, N)); break; \ + case 2: bitqueue.template pop<2>(std::span(begin, N)); break; \ + case 3: bitqueue.template pop<3>(std::span(begin, N)); break; \ + case 4: bitqueue.template pop<4>(std::span(begin, N)); break; \ + case 5: bitqueue.template pop<5>(std::span(begin, N)); break; \ + case 6: bitqueue.template pop<6>(std::span(begin, N)); break; \ + case 7: bitqueue.template pop<7>(std::span(begin, N)); break; \ + case 8: bitqueue.template pop<8>(std::span(begin, N)); break; \ + default: bitqueue.pop(bits, std::span(begin, begin + num)); \ +} } while (0) + +#define BITS(bits) \ +if (bitqueue.template pop<1,uint8_t>() == 0) { \ + bits = bitqueue.template pop<3,uint8_t>(); \ + if (bits == 7) { \ + bits += bitqueue.template pop<2,uint8_t>(); \ + if (bits == 10) \ + bits += bitqueue.template pop<6,uint8_t>(); \ + } \ +} + + template + void f_prolix(Iterator const begin, std::size_t const frame) noexcept { + Bitqueue_pop bitqueue(f_get_frame(frame)); + std::size_t flag = bitqueue.pop<18, std::size_t>(); + switch (flag) { + case 0b111111111111111100: return f_prolix_small_unsigned(begin, frame); + case 0b111111111111111000: return f_prolix_unsigned(begin, frame); + case 0b111111111111111010: return f_prolix_float(begin, frame); + default: return f_prolix_signed(begin, frame); + } + } + + template + void f_prolix_signed(Iterator begin, std::size_t frame) noexcept { + using T = std::iterator_traits::value_type; + Bitqueue_pop bitqueue(f_get_frame(frame)); + uint8_t significant_bits = 0; + for (std::size_t from = 0; from < d_size; from += d_block) { + auto const to = std::min(d_size, from + d_block); + BITS(significant_bits); + if constexpr (std::is_integral_v) + POP(T, significant_bits, begin + static_cast(from), static_cast(to - from)); + else + f_pop_into(bitqueue, significant_bits, begin + from, static_cast(to - from)); + } + } + + template + void f_prolix_float(Iterator begin, std::size_t frame) noexcept { + using T = std::iterator_traits::value_type; + Bitqueue_pop bitqueue(f_get_frame(frame)); + bitqueue.pop<18, std::size_t>(); + std::uint8_t const binary_precision = bitqueue.pop<6, std::uint8_t>(); + std::uint8_t significant_bits_exponents = 0; + Unique_array mantissas(d_block); + Unique_array exponents(d_block); + std::size_t mantissa_bits = (static_cast(1) << (binary_precision - 1)); + for (size_t from = 0; from < d_size; from += d_block) { + auto const to = std::min(d_size, from + d_block); + std::uint8_t unsigned_block = bitqueue.pop<1, std::uint8_t>(); + BITS(significant_bits_exponents); + if (unsigned_block) + POP(std::uint64_t, binary_precision - 1, reinterpret_cast(mantissas.data()), static_cast(to - from)); + else + POP(std::int64_t, binary_precision, mantissas.begin(), static_cast(to - from)); + POP(std::int16_t, significant_bits_exponents, exponents.begin(), static_cast(to - from)); + for (std::size_t i = 0; i != static_cast(to - from); ++i) { + if (exponents[i] == 1 << 12) { + if constexpr (std::integral) *begin++ = std::numeric_limits::max(); + else *begin++ = std::numeric_limits::quiet_NaN(); + } + else { + double const val = std::ldexp(static_cast(mantissas[i]) / mantissa_bits, exponents[i]); + if constexpr (std::floating_point) + *begin++ = static_cast(val); + else + *begin++ = static_cast(std::round(std::max(std::min(static_cast(std::numeric_limits::max()), val), val))); + } + } + } + } + + template + void f_prolix_unsigned(Iterator begin, std::size_t frame) noexcept { + using T = std::iterator_traits::value_type; + Bitqueue_pop bitqueue(f_get_frame(frame)); + bitqueue.pop<18, std::size_t>(); + uint8_t significant_bits = 0; + uint8_t masked_bits = 0; + for (std::size_t from = 0; from < d_size; from += d_block) { + auto const to = std::min(d_size, from + d_block); + BITS(significant_bits); + if (significant_bits != d_prolix_bits) { + if constexpr (std::is_integral_v) + POP(T, significant_bits, begin + static_cast(from), static_cast(to - from)); + else + f_pop_into(bitqueue, significant_bits, begin + static_cast(from), static_cast(to - from)); + } + else { + BITS(masked_bits); + if constexpr (std::is_integral_v) + POP(T, masked_bits, begin + static_cast(from), static_cast(to - from)); + else + f_pop_into(bitqueue, masked_bits, begin + static_cast(from), static_cast(to - from)); + for (std::ptrdiff_t i = static_cast(from); i != static_cast(to); ++i) + --begin[i]; + } + } + } + + template + requires std::floating_point::value_type> + void f_prolix_small_unsigned(Iterator begin, std::size_t frame) noexcept { + using F = typename std::iterator_traits::value_type; + auto prolix_float = [&](auto type_tag) { + using T = decltype(type_tag); + std::vector buffer(d_size); + prolix(buffer.begin(), frame); + std::transform(buffer.begin(), buffer.end(), begin, [](T val) { return static_cast(val); }); + }; + if (d_prolix_bits <= 8) prolix_float(uint8_t{}); + else if (d_prolix_bits <= 16) prolix_float(uint16_t{}); + else if (d_prolix_bits <= 32) prolix_float(uint32_t{}); + else prolix_float(uint64_t{}); + } + + constexpr std::size_t f_integer_power(std::size_t base, std::size_t exp) { + std::size_t result = 1; + while (exp) { + if (exp & 1) result *= base; // If the exponent is odd, multiply by base + base *= base; // Square the base + exp >>= 1; // Divide exponent by 2 + } + return result; + } + + template requires std::integral::value_type> + void f_prolix_small_unsigned(Iterator begin, std::size_t frame) noexcept { + using T = std::iterator_traits::value_type; + //std::size_t block = std::min(d_block, 24ul); + std::size_t block = std::min(d_block, std::size_t(24)); + Bitqueue_pop bitqueue(f_get_frame(frame)); + bitqueue.pop<18, std::size_t>(); + uint8_t bits = 0; + T max = 0; + for (size_t from = 0; from < d_size; from += block) { + auto const to = std::min(d_size, from + block); + f_get_max(bitqueue, max, bits); + switch (max) { + case 0: POP(T, 0, begin + static_cast(from), static_cast(to - from)); break; + case 1: POP(T, 1, begin + static_cast(from), static_cast(to - from)); break; + case 2: { + std::span data_block(begin + static_cast(from), begin + static_cast(to)); + std::size_t val = bitqueue.pop(f_most_significant_bit(f_integer_power(3, block) - 1)); + for (std::size_t i = 0; i != data_block.size(); ++i) { + data_block[i] = val % (3); + val /= 3; + } + break; + } + case 3: POP(T, 2, begin + static_cast(from), static_cast(to - from)); break; + case 7: POP(T, 3, begin + static_cast(from), static_cast(to - from)); break; + default: { + if (max < 7) { + std::span data_block(begin + static_cast(from), begin + static_cast(to)); + std::size_t const mult = f_integer_power(static_cast(max) + 1, block) - 1; + std::size_t val = bitqueue.pop(f_most_significant_bit(mult)); + for (std::size_t i = 0; i != data_block.size(); ++i) { + auto [quot, rem] = std::div(static_cast(val), static_cast(max + 1)); + data_block[i] = static_cast(rem); + val = static_cast(quot); + } + } + else if (bits == d_prolix_bits) + f_prolix_small_unsigned_masked(bitqueue, begin, from, max, bits); + else + POP(T, bits, begin + static_cast(from), static_cast(to - from)); + } + } + } + } + + template + void f_prolix_small_unsigned_masked(Bitqueue_pop& bitqueue, Iterator begin, std::size_t& from, auto& max, auto& bits) noexcept { + using T = std::iterator_traits::value_type; + //std::size_t block = std::min(d_block, 24ul); + std::size_t block = std::min(d_block, std::size_t(24)); + Unique_array buffer(block); + bits = 0; + for (uint8_t masked = true; from < d_size; from += block) { + auto const to = std::min(d_size, from + block); + f_get_max(bitqueue, max, bits); + if (max >= 7) { + POP(T, bits, begin + static_cast(from), static_cast(to - from)); + max = std::ranges::max(std::span(begin + static_cast(from), to - from)); + for (std::ptrdiff_t i = static_cast(from); i != static_cast(to); ++i) + --begin[i]; + } + else { + std::span data_block(begin + static_cast(from), begin + static_cast(to)); + std::size_t const mult = f_integer_power(static_cast(max) + 1, block) - 1; + std::size_t val = bitqueue.pop(f_most_significant_bit(mult)); + for (std::size_t i = 0; i != data_block.size(); ++i) { + auto [quot, rem] = std::div(static_cast(val), static_cast(max + 1)); + data_block[i] = static_cast(rem) - 1; + val = static_cast(quot); + } + } + masked = bitqueue.pop<1, uint8_t>(); + if (!masked) + break; + } + } + + template + void f_pop_into(Bitqueue_pop& bitqueue, uint8_t significant_bits, Iterator begin, std::ptrdiff_t const num) noexcept { + std::vector tmp(static_cast(num)); + if (is_signed()) + POP(std::int64_t, significant_bits, tmp.data(), num); + else + POP(std::uint64_t, significant_bits, reinterpret_cast(tmp.data()), num); + std::copy(tmp.begin(), tmp.end(), begin); + } + +#undef POP +#undef BITS + + template + inline void f_get_max(Bitqueue_pop& bitqueue, T& max, std::uint8_t& bits) const noexcept { + uint8_t flag = bitqueue.template pop<1,uint8_t>(); + if (flag == 1 && max == 0) return; + flag = static_cast(flag << 1) + bitqueue.template pop<1,uint8_t>(); + if (flag == 0b11) return; + if (flag == 0b10) { + --bits; + --max; + } + else if (flag == 0b01) { + ++bits; + max = (max == 6) ? max - 2 : max + 1; + } + else { + max = static_cast(bits = bitqueue.template pop<3,std::uint8_t>()); + if (bits == 7) { + max = std::numeric_limits::max() / 2; + bits = 3 + bitqueue.template pop<3,std::uint8_t>(); + if (bits == 10) { + bits += bitqueue.template pop<3,std::uint8_t>(); + if (bits == 17) + bits += bitqueue.template pop<6,std::uint8_t>(); + } + } + } + } + + template requires std::integral + constexpr inline std::uint8_t f_most_significant_bit(std::span const values) const noexcept { + std::size_t setbits = 0; + if constexpr (std::is_unsigned_v) + for (auto const& val : values) + setbits |= val; + else + for (auto const& val : values) + setbits |= (val == -1) ? std::size_t(1) : static_cast(std::abs(val)) << 1; + std::size_t r=0; + for ( ; setbits; setbits>>=1, ++r); + return static_cast(std::min(r, sizeof(T0) * 8)); + } + + template + constexpr inline std::uint8_t f_most_significant_bit(T0 val) const noexcept { + if constexpr (std::is_signed_v) + val = (val == -1) ? 1 : std::abs(val) << 1; + std::size_t r=0; + for ( ; val; val>>=1, ++r); + return static_cast(std::min(r, sizeof(T0) * 8)); + } + + template requires std::derived_from + void f_fill_terse_frames(STREAM& istream, std::size_t const number_of_bytes) { + f_get_frame(0).resize(number_of_bytes); + istream.read(reinterpret_cast(&(f_get_frame(0)[0])), static_cast(number_of_bytes)); + std::vector terse_sizes(d_terse_frames.size()); + std::uint8_t const* terse_begin = f_get_frame(0).data(); + for (std::size_t current_frame = 0; current_frame != d_terse_frames.size(); ++current_frame) { + Bitqueue_pop bitqueue(f_get_frame(0)); + uint8_t significant_bits = 0; + for (std::size_t from = 0; from < size(); from += d_block) { + if (bitqueue.template pop<1,uint8_t>() == 0) { + significant_bits = bitqueue.template pop<3,uint8_t>(); + if (significant_bits == 7) { + significant_bits += bitqueue.template pop<2,uint8_t>(); + if (significant_bits == 10) + significant_bits += bitqueue.template pop<6,uint8_t>(); + } + } + bitqueue.skip(significant_bits * std::min(static_cast(d_block), size() - from)); + } + terse_sizes[current_frame] = static_cast(bitqueue - terse_begin); + terse_begin += terse_sizes[current_frame]; + } + std::size_t data_start = terse_sizes[0]; + for (std::size_t frame_num = 1; frame_num != d_terse_frames.size(); ++frame_num) { + terse_sizes[frame_num] = ((terse_sizes[frame_num] + sizeof(size_t) - 1) / sizeof(size_t)) * sizeof(size_t); + f_get_frame(frame_num).resize(terse_sizes[frame_num]); + std::memcpy(f_get_frame(frame_num).data(), f_get_frame(0).data() + data_start, terse_sizes[frame_num]); + data_start += terse_sizes[frame_num]; + } + f_get_frame(0).resize(terse_sizes[0]); + f_get_frame(0).shrink_to_fit(); + } +}; +} // end namespace jpa + +#endif /* Terse_h */ diff --git a/lib/trpx2.0/include/Unique_array.hpp b/lib/trpx2.0/include/Unique_array.hpp new file mode 100644 index 00000000..d6d7046c --- /dev/null +++ b/lib/trpx2.0/include/Unique_array.hpp @@ -0,0 +1,184 @@ +// +// Unique_array.hpp +// Unique_array +// +// Created by Jan Pieter Abrahams on 23.01.2024. +// +// The `Unique_array` class is a light-weight C++ container designed for managing arrays with unique ownership +// semantics, combining aspects of `std::unique_ptr` and `std::span`. It is essentially a smart pointer for +// arrays, offering automatic memory management alongside convenient array-like access. This class is +// particularly useful in scenarios where you need the dynamic allocation of an array with the safety and +// functionality of a standard container. +// +// ### Key Features: +// 1. **Unique Ownership**: Like `std::unique_ptr`, `Unique_array` ensures that there is only one owner of +// the underlying array at any time. This guarantees that the array is automatically deallocated when the +// `Unique_array` object goes out of scope, preventing memory leaks. +// +// 2. **Array-like Access**: Thanks to inheriting from `std::span`, `Unique_array` provides direct, +// bounds-checked access to its elements, similar to standard C++ arrays and containers like `std::vector`. +// +// 3. **Support for Dynamic and Fixed-Size Arrays**: The class template can be instantiated for either +// dynamic arrays (where the size is determined at runtime) or fixed-size arrays (with compile-time size). +// +// 4. **Move Semantics**: `Unique_array` supports move semantics, allowing efficient transfer of ownership +// of the array between `Unique_array` instances. +// +// 5. **No Copy Semantics**: In line with unique ownership, copying of `Unique_array` instances is not +// allowed, preventing unintentional duplication of array ownership. +// +// ### Use Cases: +// - **Resource Management in High-Level Code**: When dealing with raw dynamic arrays in C++, there is +// always a risk of memory leaks or dangling pointers. `Unique_array` encapsulates this complexity, ensuring +// that memory is automatically managed, making it a safer choice for managing dynamically allocated arrays. +// +// - **Performance-Critical Applications**: In scenarios where performance is critical, such as real-time +// systems or high-performance computing, the deterministic destruction and lack of overhead in copying that +// `Unique_array` provides can be very beneficial. +// +// - **Safe Array Manipulation**: For applications that require safe access to array elements, `Unique_array` +// offers bounds-checked access, reducing the risk of out-of-bounds errors. +// +// - **Complex Data Structures**: In applications involving complex data structures, such as trees or graphs +// with dynamically allocated nodes, `Unique_array` can manage arrays of pointers or objects efficiently, +// ensuring proper cleanup. +// +// - **Scientific Computing and Image Processing**: In fields like scientific computing or image processing, +// where large arrays of data are common, `Unique_array` provides an efficient way to handle these arrays with +// automatic memory management and easy access. +// +// ### Example: +// jpa::Unique_array data(1000); // Dynamically allocated array of 1000 doubles +// for (int i = 0; i < 1000; ++i) { +// data[i] = computeSomeValue(i); // Safe access to elements +// } +// +// In this example, `Unique_array` is used to create a dynamically allocated array of doubles. The array is +// automatically deallocated when it goes out of scope, preventing memory leaks. +// +// Advantages of Using Unique_array over Direct new T[] and new std::array: +// +// ### Safety and Automatic Memory Management: +// - **Automatic Deallocation**: `Unique_array` automatically deallocates the array when it goes out of scope, +// preventing memory leaks, which is a common issue with manual memory management using `new` and `delete`. +// - **Exception Safety**: `Unique_array` provides strong exception safety. If an exception is thrown after +// the array is allocated, `Unique_array` ensures proper cleanup, a challenge with raw pointers. +// +// ### Modern C++ Best Practices: +// - **Resource Management**: Following the RAII (Resource Acquisition Is Initialization) idiom, +// `Unique_array` ensures that resources are tied to the lifespan of objects, leading to more predictable +// and maintainable code. +// - **Avoiding Raw Pointers**: Modern C++ guidelines advise against using raw pointers for ownership and +// memory management. `Unique_array` encapsulates these best practices. +// +// ### Ease of Use and Convenience: +// - **Array-like Access**: Provides a familiar array-like interface, including element access via the `[]` +// operator and bounds checking (when using `std::span`). +// - **Compatibility with Standard Library**: Easier integration with C++ standard library algorithms and +// functions compared to raw pointers. +// - **Size Retrieval**: Unlike raw dynamic arrays, `Unique_array` allows for easy retrieval of the array +// size, which is crucial for safe array operations and iteration. +// +// ### Flexibility: +// - **Move Semantics**: Supports move semantics, enabling efficient transfer and ownership management, +// not directly available with raw array pointers. +// - **Support for Dynamic and Fixed-Size Arrays**: Useful for both dynamically allocated and fixed-size +// arrays, providing a unified interface. +// +// In summary, `Unique_array` is a versatile and safe container for managing arrays in C++, particularly +// useful when both unique ownership and array-like access are desired. Its design simplifies memory +// management and enhances code safety and clarity, especially in complex or performance-critical +// applications. `Unique_array` is recommended over direct use of `new T[]` and `new std::array`, +// as it leads to safer, more readable, and maintainable code. + +#ifndef UNIQUE_ARRAY_HPP +#define UNIQUE_ARRAY_HPP + +#include +#include +#include + +namespace jpa { + +/// \brief A template class for managing a unique array. +/// Can be either dynamic (with std::dynamic_extent) or fixed-size. +/// \tparam T Type of the elements in the array. +/// \tparam N Size of the array. Defaults to std::dynamic_extent for a dynamic array. +template +class Unique_array : +private std::unique_ptr>>, +public std::span { +public: + /// Deleted copy constructor. + Unique_array(Unique_array const&) = delete; + + /// Deleted copy assignment operator. + Unique_array& operator=(const Unique_array&) = delete; + + /// Constructor for dynamic array. + /// \param size The size of the array to create. + Unique_array(std::size_t size) noexcept requires (N == std::dynamic_extent) : + std::unique_ptr(std::make_unique(size)), + std::span(this->get(), size) + {} + + /// Default constructor for fixed-size array. + Unique_array() noexcept requires (N != std::dynamic_extent) : + std::unique_ptr>(new std::array), + std::span(*this->get()) + {} + + /// Move constructor for dynamic array. + /// \param other Another Unique_array to move from. + Unique_array(Unique_array&& other) noexcept requires (N == std::dynamic_extent) : + std::unique_ptr(std::move((other))), + std::span(this->get(), other.size()) + {} + + /// Move constructor for fixed-size array. + /// \param other Another Unique_array to move from. + Unique_array(Unique_array&& other) noexcept requires (N != std::dynamic_extent) : + std::unique_ptr>(std::move(other)), + std::span(*this->get()) + {} + + /// Move assignment operator. + /// \param other Another Unique_array to assign from. + /// \return Reference to this Unique_array. + Unique_array& operator=(Unique_array&& other) noexcept { + if (this != &other) { + std::unique_ptr>>::operator=(std::move(other)); + if constexpr (N == std::dynamic_extent) + this->std::span::operator=(std::span(this->get(), other.size())); + else + this->std::span::operator=(std::span(*this->get())); + } + return *this; + } + + using std::span::operator[]; +}; + +/// Compare two Unique_array objects for equality. +/// \param lhs Left-hand side Unique_array to compare. +/// \param rhs Right-hand side Unique_array to compare. +/// \return True if equal, false otherwise. +template +bool operator==(Unique_array const& lhs, Unique_array const& rhs) { + return std::equal(lhs.begin(), lhs.end(), rhs.begin(), rhs.end()); +} + +/// Compare two Unique_array objects using three-way comparison. +/// \param lhs Left-hand side Unique_array to compare. +/// \param rhs Right-hand side Unique_array to compare. +/// \return std::strong_ordering result of the comparison. +template +std::strong_ordering operator<=>(const Unique_array& lhs, const Unique_array& rhs) { + for (std::size_t i = 0; i < lhs.size() && i < rhs.size(); ++i) + if (auto cmp = lhs[i] <=> rhs[i]; cmp != 0) return cmp; + return lhs.size() <=> rhs.size(); +} + +} // namespace jpa + +#endif /* UNIQUE_ARRAY_HPP */ diff --git a/lib/trpx2.0/include/XML_element.hpp b/lib/trpx2.0/include/XML_element.hpp new file mode 100644 index 00000000..85aa1cee --- /dev/null +++ b/lib/trpx2.0/include/XML_element.hpp @@ -0,0 +1,665 @@ +// +// XML_element.hpp +// XML_element +// +// Created by Jan Pieter Abrahams on 17/05/2019. +// Copyright © 2019 Jan Pieter Abrahams. All rights reserved. +// + +#ifndef XML_element_h +#define XML_element_h + +#include +#include +#include +#include + +// The XML_element class is derived from std::string. It scans an istream for an XML element with a +// specified tag. The element can contain attibutes and data. The tag itself is stripped and its +// attributes and data can be accessed as a string or parsed further by constructing another XML_element. +// +// After reading the specified Element from istream, the istream is left in the condition that exists +// after having read the Element, but not further. This allows mixing XML straightforwardly with +// binary data. For instance, after having read a XML metadata header, the file is positioned where the +// binary data starts. +// +// Constructor: +// template +// XML_element(std::istream& is, std::string const& tag) +// XML_element(std::istream&& is, std::string const& tag) +// istr is skipped until the specified XML tag ('<' + tag + '>') is encountered. Then all the XML data +// defined by that tag, are parsed into the XML_element, so until '' is encountered, +// and scanning the file stops after having read the end tag. Any attributes of the tag are parsed into +// XML_element. Tags can be empty (when they end with “/>”), in which case only the attributes of tag +// are parsed into the XML_element. +// XML_element(std::istream& is) +// XML_element(std::istream&& is) +// istr is skipped until the first occurence of a valid XML tag is encountered. Then this tag is used for +// creating the XML_element. +// XML_element(std::string const& xml_string) +// XML_element(std::string const& xml_string, std::string const& tag) +// Delegates to the constructor requiring a std::istream, after creating a std::istringstream(xml_string). +// +// Member functions: +// std::string const attribute(std::string const& att_name) const +// Returns the attribute of the tag of XML_element that has the name 'att_name'. An XML_element line does +// not have to have attributes. If none of the attributes is named 'att_name', an empty string is returned. +// std::string const attribute(std::size_t n) const +// Returns the n-th attribute of an XML_element. If n is larger than the number of attributes, an empty +// string is returned. +// std::string const tag() +// Returns the tag name of the XML_element that was used for extracting it from the XML data. If no tag name +// was specified when the XML_element was created, an empty string is returned. +// std::string const tag(std::string const& tag) +// Returns the XML data encoded by the first occurrence of requested tag as an XML_element. +// std::vector const multitag(std::string const& tag) +// Returns all requested tags as a vector, in the order in which they occur in original XML data. +// std::vector const value(std::string const& tag) +// Returns a vector of the requested type, containing the numerical values of the requested tag of the +// XML_element. +// std::vector> const multivalue(std::string const& tag) +// Returns a vector of all the requested values. +// +// An example: +// +// int main(int argc, const char * argv[]) { +// std::string xml("\n \n\n"); +// std::cout << "xml stream:\n" << xml << std::endl <<"Parsing produces:\n"; +// std::cout << "Element 1: " << XML_element(xml, "Element1") << std::endl; +// std::cout << "Attribute att0 of element 1: " << XML_element(xml, "Element1").attribute("att0") << std::endl; +// std::cout << "Attribute att1 of element 1 does not exist, so empty string is returned: " << '\"' << XML_element(xml, "Element1").attribute("att1") << '\"' << std::endl; +// std::cout << "Attribute att1 of element 2: " << XML_element(xml, "Element2").attribute("att1") << std::endl; +// std::cout << "Attribute att2 of element 2: " << XML_element(xml, "Element2").attribute("att2") << std::endl; +// std::istringstream xml_stream(xml); +// XML_element(xml_stream, "Element2"); +// std::string rest; +// xml_stream >> rest; +// std::cout << "After reading Element2, the xml_stream is now located ath the end of Element2. \nSubsequent output from theis stream produces:\n"<< rest << std::endl; +// +// return 0; +// } +// +// Produces the following output: +// +// xml stream: +// +// +// 1 2 3 4 +// +// +// +// Parsing produces: +// Element 1: +// +// 1 2 3 4 +// +// Attribute att0 of element 1: plop +// Attribute att1 of element 1 does not exist, so an empty string is returned: "" +// Attribute att1 of element 2: 1 +// Attribute att2 of element 2: 22 +// After reading Element2, the xml_stream is now located at the end of Element2. +// Subsequent output from this stream produces: +// + +namespace jpa { + +/** + * @class XML_element + * @brief Represents an XML element with attributes and data. + * + * The XML_element class is derived from std::string and is used to parse and store XML data. + * It scans an input stream for an XML element with a specified tag, extracts its attributes, + * and stores the its XML data defined inside the element. The tag itself is stripped, and attributes + * and data can be accessed as strings. Or data can be parsed if they contain XML_elements. + * + * After reading the specified element from the input stream, the stream is left in a state that + * allows mixing XML with binary data. The stream is positioned at the first location after the closing '/>' tag + * of the requested XML_element. This can be the location where the binary data starts. For instance, + * after having read a XML metadata header, the file is positioned where the binary data starts. + * + * An example: + * @code{.cpp} + * int main(int argc, const char * argv[]) { + * std::string xml("\n \n\n"); + * std::cout << "xml stream:\n" << xml << std::endl <<"Parsing produces:\n"; + * std::cout << "Element 1: " << XML_element(xml, "Element1") << std::endl; + * std::cout << "Attribute att0 of element 1: " << XML_element(xml, "Element1").attribute("att0") << std::endl; + * std::cout << "Attribute att1 of element 1 does not exist, so an empty string is returned: " << '\"' << XML_element(xml, "Element1").attribute("att1") << '\"' << std::endl; + * std::cout << "Attribute att1 of element 2: " << XML_element(xml, "Element2").attribute("att1") << std::endl; + * std::cout << "Attribute att2 of element 2: " << XML_element(xml, "Element2").attribute("att2") << std::endl; + * std::istringstream xml_stream(xml); + * XML_element(xml_stream, "Element2"); + * std::string rest; + * xml_stream >> rest; + * std::cout << "After reading Element2, the xml_stream is now located at the end of Element2. \nSubsequent output from this stream produces:\n"<< rest << std::endl; + * return 0; + * } + * @endcode + * + * Produces the following output: + * 
+ *    xml stream:
+ *    
+ *      
+ *       1 2 3 4 
+ *      
+ *    
+ *
+ *    Parsing produces:
+ *    Element 1:
+ *      
+ *       1 2 3 4 
+ *      
+ *    Attribute att0 of element 1: plop
+ *    Attribute att1 of element 1 does not exist, so an empty string is returned: ""
+ *    Attribute att1 of element 2: 1
+ *    Attribute att2 of element 2: 22
+ *    After reading Element2, the xml_stream is now located at the end of Element2.
+ *    Subsequent output from this stream produces:
+ *    
+ *
+ */ + +class XML_element : public std::string { + + class s_param { + friend XML_element; + s_param(std::string const& str, std::size_t frm, std::size_t siz) : s(str), from(frm), size(siz) {} + friend void swap(s_param& left, s_param &right) {std::swap(left.from, right.from); std::swap(left.size, right.size);} + std::string const& s; + std::size_t from; + std::size_t size; + + std::size_t to() const noexcept {return from + size;} + + public: + s_param(s_param const& other) : s(other.s), from(other.from), size(other.size) {} + constexpr auto& operator=(s_param const& other) noexcept {from = other.from; size = other.size; return *this;} + + template + operator const T() const { + T r; + std::string const s1(s, from, size); + std::istringstream(s1) >> r; + return r; + } + + operator const std::string() const { return std::string(s, from, size); } + + std::string const tag() const { + if (s[from] != '<') return ""; + std::size_t tagsize = s.find_first_of(" >", from) - from - 1; + if (s[from + tagsize] == '/') --tagsize; + return std::string(s, from + 1, tagsize); + } + + operator XML_element() const { return XML_element(std::string(s, from, size), tag()); } + }; + +public: + + /** + * @brief Constructor for XML_element with specified tag. + * + * Constructs an XML_element by parsing XML data from the input stream 'istr' with the specified XML tag. + * When a character sequence '<' + tag + '>' is encountered, it checks if the element is empty (which + * is the case when the last two characters of the tag are '/>'), in which case it only stores the attributes. If + * the the tag is not empty, it scans for the end tag: '', and stores the data between the + * beginning and end tag. The input stream is left at the position immediately after the closing '>' character + * of the specified XML element. If no element is found, the file pointer is 'eof', and no data are stored. + * + * @param istr The input stream containing the XML data. + * @param tag The tag to search for in the XML data. + */ + XML_element(std::istream&& istr, std::string const& tag) : + std::string(), + d_tag(tag), + d_attributes(f_get_attributes(istr)) { + if (istr.good() && !f_prev_is(istr, "/>") && d_tag_found) { // only do if this is not an empty element + static_cast(*this) = f_read_element(istr, tag); + d_param = f_find_all_params(std::string(*this)); + } + } + + /** + * @brief Constructor for XML_element with specified tag. + * + * @see XML_element(std::istream&& istr, std::string const& tag) + * @param istr The input stream containing the XML data. + * @param tag The tag to search for in the XML data. + */ + XML_element(std::istream& istr, std::string const& tag) : XML_element(std::move(istr), tag) {} + + /** + * @brief Constructor for XML_element with specified tag. + * + * Constructs an XML_element by parsing XML data from the input string 'xml_string' with the specified + * XML tag. + * + * @see XML_element(std::istream&& istr, std::string const& tag) + * @param xml_string The input string containing the XML data. + * @param tag The tag to search for in the XML data. + */ + XML_element(std::string const& xml_string, std::string const& tag) : XML_element(std::istringstream(xml_string), tag) {} + + + /** + * @brief Constructor for XML_element without specifying the tag. + * + * Constructs an XML_element by parsing XML data from the input stream 'istr' without specifying a tag. + * The tag is determined by searching for the first valid XML tag in the input stream. + * + * @see XML_element(std::istream&& istr, std::string const& tag) + * @param istr The input stream containing the XML data. + */ + XML_element(std::istream&& istr) : + XML_element(std::move(istr), std::string([&] { + f_read_upto(istr, '<'); + auto tag_start = istr.tellg() - std::basic_istream::pos_type(1); + std::string tag = f_read_upto_any_of(istr, "/> "); + istr.seekg(tag_start); + return tag; + }())) {} + + /** + * @brief Constructor for XML_element without specifying the tag. + * + * Constructs an XML_element by parsing XML data from the input stream 'istr' without specifying a tag. + * The tag is determined by searching for the first valid XML tag in the input stream. + * + * @see XML_element(std::istream&& istr, std::string const& tag) + * @param istr The input stream containing the XML data. + */ + XML_element(std::istream& istr) : XML_element(std::move(istr)) {} + + /** + * @brief Constructor for XML_element without specifying the tag. + * + * Constructs an XML_element by parsing XML data from the input string 'istr' without specifying a tag. + * The tag is determined by searching for the first valid XML tag in the input stream. + * + * @see XML_element(std::istream&& istr, std::string const& tag) + * @param xml_string The input string containing the XML data. + */ + XML_element(std::string const& xml_string) : XML_element(std::istringstream(xml_string)) {} + + /** + * @brief Retrieves the attribute value by name. + * + * Returns the attribute value with the specified 'name' from the XML element's attributes. + * + * @param name The name of the attribute to retrieve. + * @return The attribute value or an empty string if the attribute is not found. + */ + std::string const attribute(std::string const& name) const noexcept { + for (std::size_t i = 0; i <= (d_attributes.size() - name.size() - 3); ++i) { + if ((d_attributes[i + name.size()] == '=') && (name == d_attributes.substr(i, name.size()))) { + auto strt = 1 + (i += name.size() + 1); + for (char q = d_attributes[i]; d_attributes[++i] != q;); + return d_attributes.substr(strt, i - strt); + } + else if (d_attributes[i] == '=') + for (char q = d_attributes[++i]; d_attributes[++i] != q;); + } + return ""; + } + + /** + * @brief Retrieves an attribute value by index. + * + * Returns the attribute value at the specified index 'n' from the XML element's attributes. + * + * @param n The index of the attribute to retrieve. + * @return The attribute value or an empty string if the index is out of bounds. + */ + std::string const attribute(std::size_t const& n) const noexcept { + std::string r; + std::size_t i=0; + std::size_t found = 0; + for ( ; (i <= d_attributes.size()) && (found != (n+1)); ++i, ++found) + i = d_attributes.find('='); + if (found == (n+1)) { + i = d_attributes.find_first_of("'\"", i); + for (std::size_t j = i+1 ; d_attributes[j] != d_attributes[i]; ++j) + r.push_back(d_attributes[j]); + } + return r; + } + + /** + * @brief Appends an attribute to the tag. + * + * @param attribute The attribute name. + * @param val The optional attribute value or an empty string. + * @return A reference to the updated XML_element instance. + */ + XML_element& add_attribute(std::string const& attribute, std::string const& val = "") { + if (d_attributes.size() != 0) + d_attributes.append(" "); + d_attributes.append(attribute); + if (val.size() != 0) + d_attributes.append("=\"" + val + "\""); + return *this; + } + + /** + * @brief Appends an attribute with a numeric value to a tag. + * + * @param attribute The attribute name. + * @param val The arithmetic value. + * @return A reference to the updated XML_element instance. + */ + template requires std::is_arithmetic_v + XML_element& add_attribute(std::string const& attribute, T const val) { + return add_attribute(attribute, std::to_string(val)); + } + + /** + * @brief Appends an attribute with a range of numeric values to a tag. + * + * @param attribute The attribute name. + * @param container The range of arithmetic values. + * @return A reference to the updated XML_element instance. + */ + template requires (requires (T c) { std::to_string(*c.begin()); }) + XML_element& add_attribute(std::string const& attribute, T const container) { + std::string values; + for (auto it = container.begin(); it != container.end(); ++it) { + values.append(std::to_string(*it)); + if (it != container.end() + 1) + values.append(" "); + } + return add_attribute(attribute, values); + } + + /** + * @brief Appends an element. + * + * @param elem The element to be added. + * @return A reference to the updated XML_element instance. + */ + XML_element& add_element(std::string const& elem) { + static_cast(*this) = *this + elem; + d_param = f_find_all_params(std::string(*this)); + return *this; + } + + /** + * @brief Appends a range of arithmetic elements. + * + * @param container The range of elements to be added. + * @return A reference to the updated XML_element instance. + */ + template requires (requires (T c) { std::to_string(*c.begin()); }) + XML_element& add_element(std::string const& tag, T const container) { + std::string s = *this + '<' + tag + '>'; + for (auto it = container.begin(); it != container.end(); ++it) + s = s + ' ' + std::to_string(*it); + static_cast(*this) = s + " ';; + return *this; + } + + /** + * @brief Returns the encompassing XML string. + * @return An XML string, tagged with the tag (and tag attributes) of the XML_element object. + */ + std::string XML() const noexcept { + std::string xml = '<' + d_tag; + if (d_attributes.size() != 0) + xml = xml + ' ' + d_attributes; + if (d_param.size() == 0) + xml = xml + "/>"; + else + xml = xml + '>' + *this + "'; + return xml; + } + + /** + * @brief Retrieves the tag name of the XML_element. + * + * Returns the tag name of the XML_element that was used to extract it from the XML data. + * + * @return The tag name of the XML element or an empty string if no tag was specified. + */ + std::string const& tag() const noexcept {return d_tag;} + + /** + * @brief Retrieves an XML element by tag name. + * + * Returns the first XML element with the specified 'tag' as an XML_element. + * + * @param tag The tag name to search for. + * @return An XML_element containing the first matching element. + */ + XML_element const tag(std::string const& tag) const noexcept { return multitag(tag, 1)[0]; } + + /** + * @brief Retrieves multiple XML elements with the same tag. + * + * Returns a vector containing all XML elements with the specified 'tag' found in the parsed data, + * up to a maximum of 'max_tags', if this is provided. If max_tags is not provided, all + * XML_elements in the stream are returned. + * + * @param tag The tag name to search for. + * @param max_tags The maximum number of elements to retrieve. Default is -1 (no limit). + * @return A vector of XML elements with the specified tag. + */ + std::vector const multitag(std::string const& tag, std::size_t max_tags = static_cast(-1)) const noexcept { + std::vector r; + for (std::size_t i=0; (max_tags != 0) && (i != d_param.size()); ++i) + if (d_param[i].tag() == tag) { + r.push_back(d_param[i]); + --max_tags; + } + return r; + } + + /** + * @brief Retrieves values from an XML element that are defined by the first occurence of a specified tag. + * + * Returns a vector containing the values of type 'T' parsed from the XML element with the specified 'tag'. + * If no matching elements are found, an empty vector is returned. For example: + * @code{.cpp} + * std::string xml = " 512 512 "; + * std::vector val = XML_element(xml).value("size"); + * assert(val[0] == 512 && val[1] == 512); + * @endcode + * + * @tparam T The type of values to retrieve. + * @param tag The tag name to search for. + * @return A vector of values of type 'T'. + */ + template + std::vector const value(std::string const& tag) const noexcept { + auto tmp = multitag(tag, 1); + if (tmp.size() == 0) + return std::vector(); + else + return tmp[0]; + } + + /** + * @brief Retrieves values from an XML element that are defined by some or all instances of a specified tag. + * + * Returns a vector of vectors containing the values of type 'T' parsed from the XML element with the specified 'tag'. + * If no matching elements are found, an empty vector is returned. For example: + * + * @code{.cpp} + * std::string xml = " 2 50 3 49 " + * std::vector> vals = XML_element(xml).multivalue("dead_pix"); + * assert(vals[0][0] == 2 && vals[1][1] == 49); + * @endcode + * + * @tparam T The type of values to retrieve. + * @param tag The tag name to search for. + * @param max_tags The maximum instances of 'tag' to consider; if -1 or not specidfies, all occurences are considered. + * @return A vector of vectors of values of type 'T'. + */ + template + std::vector> const multivalue(std::string const& tag, std::size_t max_tags = static_cast(-1)) const noexcept { + std::vector> r; + auto tmp = multitag(tag, max_tags); + for (auto const& t: tmp) + r.push_back(t); + return r; + } + + +private: + std::vector d_param; + std::string d_tag; + bool d_tag_found = false; + std::string d_attributes; + + std::string f_get_attributes(std::istream& istr) { + if (d_tag == "") + return ""; + f_find_tag(istr, d_tag); + if (d_tag_found) { + std::string r = std::string(f_read_upto(istr, '>'), d_tag.size() + 1); + if (r.size() >= 1) + r.pop_back(); + if ((r.size() > 1) && (r.back() == '/')) + r.pop_back(); + return r; + } + return ""; + } + + void f_find_tag(std::istream& istr, std::string const& tag) { + static auto const all = std::numeric_limits::max(); + for (istr.ignore(all, '<'); istr.good() && !f_next_is(istr, tag); istr.ignore(all, '<')) { + if (f_next_is(istr, "![CDATA[")) + for (f_read_upto(istr, ']'); istr.good() && !f_next_is(istr, "]>"); f_read_upto(istr, ']')); + else if (f_next_is(istr, "!--")) + for (f_read_upto(istr, '-'); istr.good() && !f_next_is(istr, "->"); f_read_upto(istr, '-')); + } + if (istr.good()) + d_tag_found = true; + } + + std::vector f_find_comments_and_data(std::string const& s) { + std::vector r; + for (s_param p = f_find_interval(s, ""); p.from != std::string::npos; p = f_find_interval(s, "", p.to())) + r.push_back(p); + for (s_param p = f_find_interval(s, ""); p.from != std::string::npos; p = f_find_interval(s, "", p.to())) + r.push_back(p); + std::sort(r.begin(), r.end(), [](s_param const& left, s_param const& right) {return left.from < right.from;}); + for (std::size_t i = 0; i!=r.size(); ++i) + for (std::size_t j = i+1; (j!=r.size()) && (r[j].from < (r[i].from + r[i].size)); ++j) + r.erase(r.begin() + static_cast(j--)); + return r; + } + + std::vector f_find_all_params(std::string const& s) { + std::vector skip = f_find_comments_and_data(s); + std::vector r; + for (s_param param(*this, f_ignore(s, 0, skip), 1); param.from < s.size(); param.from = f_ignore(s, param.to(), skip)) { + if (s[param.from] != '<') { + for (param.size = 1; !f_is_white(s[param.to()]) && (s[param.to()] != '<'); ++param.size); + r.push_back(param); + } + else { + param.size = s.find('>', param.from) + 1 - param.from; + if (s[param.to() - 2] == '/') + r.push_back(param); + else { + std::size_t end_tag_size = s.find_first_of(" >", param.from) - param.from - 1; + std::string const end_tag = "'; + param.size = s.find(end_tag, param.to()) - param.from + end_tag.size(); + for (auto& ignore: skip) + if ((param.to() > ignore.from) && (param.to() < ignore.to())) + param.size = s.find(end_tag, param.to()) - param.from; + r.push_back(param); + } + } + } + r.insert(r.begin(), skip.begin(), skip.end()); + std::sort(r.begin(), r.end(), [](s_param const& left, s_param const& right) {return left.from < right.from;}); + return r; + } + + s_param f_find_interval(std::string const& s, std::string const& from, std::string const& to, std::size_t pos = 0) { + std::size_t param_start = s.find(from, pos); + if (param_start == std::string::npos) + return s_param(*this, param_start, 0); + else + return s_param(*this, param_start, s.find(to, param_start + from.size()) - param_start + to.size()); + } + + std::size_t f_ignore(std::string const& s, std::size_t pos, std::vector skip = std::vector()) { + std::size_t start = pos; + std::size_t i = 0; + do { + pos = start; + while (f_is_white(s[start]) && (start < s.size())) + ++start; + for ( ; (i != skip.size()) && (start >= skip[i].from); ++i) + if ((start >= skip[i].from) && (start < skip[i].to())) + start = skip[i].to(); + } while ((pos != start) && (start < s.size())); + return start; + } + + std::string f_read_upto(std::istream& istr, char const c) { + std::string r; + int ch; + while (istr.good() && (ch = istr.get()) != std::char_traits::eof()) { + r.push_back(static_cast(ch)); + if (ch == c) break; + } + return r; + } + + std::string f_read_upto_any_of(std::istream& istr, std::string const& stop) { + std::string r; + int ch; + while (istr.good() && (ch = istr.get()) != std::char_traits::eof()) { + if (stop.find(static_cast(ch)) != std::string::npos) break; // ✅ Safe conversion + r.push_back(static_cast(ch)); + } + return r; + } + + std::string f_read_element(std::istream& istr, std::string const& tag) { + std::string element = f_read_upto(istr, '<'); + for (; !f_next_is(istr, '/' + tag + '>'); element += f_read_upto(istr, '<')) { + if (f_next_is(istr, "![CDATA[")) + for (element += f_read_upto(istr, '>'); !element.compare(element.size() - 3, 3, "]]>"); element += f_read_upto(istr, '>')); + else if (f_next_is(istr, "!--")) + for (element += f_read_upto(istr, '>'); !element.compare(element.size() - 3, 3,"-->"); element += f_read_upto(istr, '>')); + } + element.pop_back(); + return element; + } + + bool f_is_white(char const c) const noexcept {return (c == 0x20) || (c==0x09) || (c==0x0D) || (c==0x0A);} + + bool f_next_is(std::istream& istr, std::string const& s) { + std::string tst(s.size(), 0); + if (istr.read(&tst[0], static_cast(s.size()))) { + istr.seekg(-std::ssize(s), std::ios_base::cur); + return d_tag_found = tst == s; + } + return d_tag_found = false; // Unable to read from the stream, so it's not a match + } + + bool f_prev_is(std::istream& istr, std::string const& s) { + std::string tst(s.size(), 0); + istr.seekg(-std::ssize(s), std::ios_base::cur); + istr.read(&tst[0], static_cast(s.size())); + return tst == s; + } + + template + operator std::vector() const { + std::vector r; + std::istringstream str(std::string(*this)); + for (T tmp; (str >> tmp).good(); ) r.push_back(tmp); + return r; + } +}; + +} + +#endif /* XML_element_h */ diff --git a/lib/trpx2.0/notebook/pyterse.ipynb b/lib/trpx2.0/notebook/pyterse.ipynb new file mode 100644 index 00000000..7a32e8df --- /dev/null +++ b/lib/trpx2.0/notebook/pyterse.ipynb @@ -0,0 +1,486 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Terse Python Binder\n", + "\n", + "## Constructors\n", + "\n", + "- **Default Constructor:** `Terse()`\n", + "- **From NumPy Array with Mode:** `Terse(py::array data, Terse_mode mode)`\n", + "- **From Python Stream:** `Terse(py::object stream)`\n", + "\n", + "## Methods\n", + "\n", + "- **`insert(pos, data, mode)`** \n", + " *Insert data at a specific frame position.*\n", + " \n", + "- **`push_back(data, mode)`** \n", + " *Append data to the end.*\n", + " \n", + "- **`prolix(data)`** \n", + " *Decompress into a provided NumPy array.*\n", + " \n", + "- **`prolix()`** \n", + " *Decompress and return a new NumPy array.*\n", + " \n", + "- **`at(pos)`** \n", + " *Retrieve a specific frame as a new Terse object.*\n", + " \n", + "- **`erase(pos)`** \n", + " *Remove a specific frame.*\n", + " \n", + "- **`write(stream)`** \n", + " *Write compressed data to a Python file-like object.*\n", + "\n", + "## Properties\n", + "\n", + "### Read-only\n", + "\n", + "- **`size`** \n", + " *Number of elements in each frame.*\n", + " \n", + "- **`number_of_frames`** \n", + " *Number of frames stored in the Terse object.*\n", + " \n", + "- **`is_signed`** \n", + " *Indicates if the data is signed (`True`) or unsigned (`False`).*\n", + " \n", + "- **`bits_per_val`** \n", + " *Number of bits per value (e.g., 8, 16, 32, 64).*\n", + " \n", + "- **`terse_size`** \n", + " *Size of the compressed data in bytes.*\n", + " \n", + "- **`file_size`** \n", + " *Size of the Terse file in bytes.*\n", + "\n", + "### Metadata Management\n", + "\n", + "- **`metadata(frame)`** \n", + " *Retrieve metadata for the specified frame.*\n", + " \n", + "- **`set_metadata(frame, metadata)`** \n", + " *Set metadata for the specified frame.*\n", + "\n", + "### Dimension Management\n", + "\n", + "- **`dim()`** \n", + " *Get the dimensions of each Terse frame.*\n", + " \n", + "- **`set_dim(dims)`** \n", + " *Set the dimensions of each Terse frame.*\n", + "\n", + "### Configuration\n", + "\n", + "- **`block_size()`** \n", + " *Get the block size used for Terse compression.*\n", + " \n", + "- **`set_block_size(block_size)`** \n", + " *Set the block size for Terse compression.*\n", + " \n", + "- **`fast()`** \n", + " *Check if the default compression mode is `UNSIGNED`.*\n", + " \n", + "- **`set_fast(value)`** \n", + " *Enable (`True`) or disable (`False`) the default compression mode to `UNSIGNED`.*\n", + " \n", + "- **`small()`** \n", + " *Check if the default compression mode is `SMALL_UNSIGNED`.*\n", + " \n", + "- **`set_small(value)`** \n", + " *Enable (`True`) or disable (`False`) the default compression mode to `SMALL_UNSIGNED`.*\n", + " \n", + "- **`dop()`** \n", + " *Get the degree of parallelism for compression or decompression (0 for sequential, 1 for all cores, or a fractional value).*\n", + " \n", + "- **`set_dop(value)`** \n", + " *Set the degree of parallelism for compression or decompression.*\n", + "\n", + "### Memory Management\n", + "\n", + "- **`shrink_to_fit()`** \n", + " *Release unused memory and shrink the storage to fit the data.*\n", + "\n" + ] + }, + { + "cell_type": "code", + "execution_count": 1, + "metadata": {}, + "outputs": [], + "source": [ + "import sys\n", + "import os\n", + "import numpy as np\n", + "from PIL import Image\n", + "import matplotlib.pyplot as plt\n", + "import tifffile as tiff\n", + "import os" + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "build_dir \n" + ] + } + ], + "source": [ + "build_dir = os.path.join('') #Provide the path to the build directory\n", + "if build_dir not in sys.path:\n", + " sys.path.append(build_dir)\n", + "print('build_dir', build_dir)\n", + "import pyterse\n", + "from pyterse import Terse, TerseMode" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Basic usage" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "\n", + "\n", + "tiff_stack = tiff.imread('') #nD TiFF stack\n", + "\n", + "print(f\"TIFF stack shape: {tiff_stack.shape}\")\n", + "print(f\"TIFF stack dtype: {tiff_stack.dtype}\")" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Block size" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "terse = Terse() #empty Terse object" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "terse.block_size()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "terse.set_block_size(6)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "terse = Terse(tiff_stack, TerseMode.SIGNED) #Create Terse object from nD numpy array with signed compression mode" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Degree of parallelism" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "terse.dop() #1.0 means use all available system resources" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "terse.set_dop(0.5)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Terse compression mode::unsigned" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "terse.set_fast(True) # This mode will use fast compression" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "terse.push_back(tiff_stack) #push_back() method to add numpy arrays to the Terse object" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Properties" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "terse.dim()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "terse.size" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "terse.number_of_frames" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "terse.number_of_bytes" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "terse.bits_per_val" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "terse.is_signed" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Terse compression mode::small_unsigned" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "terse = Terse()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "terse.set_small(True)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "terse.push_back(tiff_stack)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "terse.number_of_bytes" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Save and load" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "terse.save('filename.trpx')" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "terse = Terse.load('filename.trpx')" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "original_data = terse.prolix()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "specific_frame = terse.at(0).prolix()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "specific_frame.shape" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Erase and insert" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "terse.erase(0)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "terse.insert(0, tiff_stack)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Metadata" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "terse.set_metadata(0, \"Some metadata\")\n", + "md = terse.metadata(0)\n", + "md" + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "hdf5-env", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.9.19" + } + }, + "nbformat": 4, + "nbformat_minor": 2 +} diff --git a/lib/trpx2.0/py_tests/filter_test.py b/lib/trpx2.0/py_tests/filter_test.py new file mode 100644 index 00000000..9bc0805e --- /dev/null +++ b/lib/trpx2.0/py_tests/filter_test.py @@ -0,0 +1,86 @@ +import os +import h5py +import numpy as np +import h5py.h5pl +import h5py.h5z as h5z +import sys +import hdf5plugin + + +plugin_dir = os.path.abspath(os.path.join(os.path.dirname(__file__), '../build', 'tersecodec/')) +#plugin_dir = os.path.join(sys.prefix, "lib", "hdf5", "plugins") +os.environ["HDF5_PLUGIN_PATH"] = plugin_dir + +# Enable HDF5 plugin debugging +os.environ["HDF5_DEBUG"] = "all" + +print("HDF5_PLUGIN_PATH:", os.environ.get("HDF5_PLUGIN_PATH")) +print("h5py version:", h5py.__version__) +print("HDF5 version:", h5py.version.hdf5_version) + +TERSE_FILTER_ID = 32029 + +if not h5z.filter_avail(TERSE_FILTER_ID): + print("Automatic loading failed, trying manual registration...") + import ctypes + lib_path = os.path.join(plugin_dir, "libh5terse.so") + try: + compression_lib = ctypes.CDLL(lib_path) + register_filter = compression_lib.register_terse_filter + register_filter.restype = ctypes.c_int + if register_filter() < 0: + raise RuntimeError("Failed to register Terse filter") + print("Manual registration successful") + except Exception as e: + raise RuntimeError(f"Failed to load and register filter: {str(e)}") + +if not h5z.filter_avail(TERSE_FILTER_ID): + raise RuntimeError("TERSE filter not available after registration attempts") + +print("TERSE filter is available.") +file_name = "test_terse.h5" +dataset_name = "compressed_data" + + +data = np.random.choice([0, 1], size=(30, 512, 512), p=[0.5, 0.5]).astype(np.uint32) + + +print(f"Original data shape: {data.shape}") +print(f"Original data size: {data.nbytes / 1024:.2f} KB") + +with h5py.File(file_name, "w") as f: + space = h5py.h5s.create_simple(data.shape) + dcpl = h5py.h5p.create(h5py.h5p.DATASET_CREATE) + dcpl.set_chunk(data.shape) + + + numpy_to_type_code = { + np.int16: 0, # TYPE_CODE_INT16 + np.uint16: 1, # TYPE_CODE_UINT16 + np.int32: 2, # TYPE_CODE_INT32 + np.uint32: 3, # TYPE_CODE_UINT32 + } + + data_type_code = numpy_to_type_code.get(data.dtype.type) + if data_type_code is None: + raise ValueError(f"Data type {data.dtype} not supported by TERSE filter") + + + dcpl.set_filter(TERSE_FILTER_ID, h5py.h5z.FLAG_OPTIONAL, (data_type_code,)) + + dtype = h5py.h5t.py_create(data.dtype) + dset_id = h5py.h5d.create(f.id, dataset_name.encode(), dtype, space, dcpl=dcpl) + dset = h5py.Dataset(dset_id) + dset[...] = data + + +file_size = os.path.getsize(file_name) +print(f"Compressed HDF5 file size: {file_size / 1024:.2f} KB") + +with h5py.File(file_name, "r") as f: + dset = f[dataset_name] + read_data = dset[:] + +assert np.array_equal(data, read_data), "Data verification failed" + +print("Data successfully compressed and decompressed with TERSE filter.") \ No newline at end of file diff --git a/lib/trpx2.0/py_tests/pyterse_test.py b/lib/trpx2.0/py_tests/pyterse_test.py new file mode 100644 index 00000000..48273c2f --- /dev/null +++ b/lib/trpx2.0/py_tests/pyterse_test.py @@ -0,0 +1,160 @@ +import unittest +import numpy as np +import io +import shutil +import sys +import os + +for root, dirs, files in os.walk('.'): + if '__pycache__' in dirs: + shutil.rmtree(os.path.join(root, '__pycache__')) + +sys.path.append(os.path.join(os.getcwd(), 'build', 'pyterse/')) + +from pyterse import Terse, TerseMode + +class TestTerseLibrary(unittest.TestCase): + def setUp(self): + # Common test data + self.test_data_1d = np.array([1, 2, 3, 4, 5], dtype=np.int32) + self.test_data_2d = np.array([[1, 2, 3], [4, 5, 6]], dtype=np.int32) + self.test_data_3d = np.array([[[1, 2], [3, 4]], [[5, 6], [7, 8]]], dtype=np.int32) + + def test_basic_construction(self): + """Test basic Terse object construction""" + terse = Terse() + self.assertEqual(terse.number_of_frames, 0) + self.assertEqual(terse.size, 0) + + def test_construction_with_data(self): + """Test Terse construction with initial data""" + terse = Terse(self.test_data_1d, TerseMode.SIGNED) + self.assertEqual(terse.size, len(self.test_data_1d)) + self.assertTrue(terse.is_signed) + + def test_data_types(self): + """Test different numpy data types""" + for dtype in [np.int8, np.uint8, np.int16, np.uint16, np.int32, np.uint32, np.int64, np.uint64]: + data = np.array([1, 2, 3], dtype=dtype) + terse = Terse(data, TerseMode.SIGNED) + result = terse.prolix() + np.testing.assert_array_equal(result, data) + self.assertEqual(result.dtype, data.dtype) + + + def test_dimensions(self): + """Test handling of different dimensional arrays""" + # 1D + terse_1d = Terse(self.test_data_1d) + np.testing.assert_array_equal(terse_1d.prolix(), self.test_data_1d) + + # 2D + terse_2d = Terse(self.test_data_2d) + np.testing.assert_array_equal(terse_2d.prolix(), self.test_data_2d) + + # 3D + terse_3d = Terse(self.test_data_3d) + np.testing.assert_array_equal(terse_3d.prolix(), self.test_data_3d) + + def test_compression_modes(self): + """Test different compression modes""" + data = np.array([1, 2, 3], dtype=np.int32) + for mode in [TerseMode.SIGNED, TerseMode.UNSIGNED, TerseMode.SMALL_UNSIGNED, TerseMode.DEFAULT]: + terse = Terse(data, mode) + np.testing.assert_array_equal(terse.prolix(), data) + + def test_frame_operations(self): + """Test frame-based operations""" + # Insert frames + terse = Terse() + terse.push_back(self.test_data_1d) + terse.insert(0, self.test_data_1d) + self.assertEqual(terse.number_of_frames, 2) + + # Access frame + frame = terse.at(0) + np.testing.assert_array_equal(frame.prolix(), self.test_data_1d) + + # Erase frame + terse.erase(0) + self.assertEqual(terse.number_of_frames, 1) + + def test_metadata(self): + """Test metadata operations""" + terse = Terse(self.test_data_1d) + test_metadata = "test" + terse.set_metadata(0, test_metadata) + self.assertEqual(terse.metadata(0), test_metadata) + + def test_dimension_operations(self): + """Test dimension-related operations""" + terse = Terse(self.test_data_2d) + original_dims = terse.dim() + + # Test dimension setting + new_dims = [6, 1] # Same total size as 2x3 + terse.set_dim(new_dims) + self.assertEqual(terse.dim(), new_dims) + + # Test invalid dimension setting + with self.assertRaises(ValueError): + terse.set_dim([2, 2]) # Invalid total size + + def test_stream_operations(self): + """Test stream reading and writing""" + # Write to stream + stream = io.BytesIO() + terse_write = Terse(self.test_data_1d) + terse_write.write(stream) + + # Read from stream + stream.seek(0) + terse_read = Terse(stream) + np.testing.assert_array_equal(terse_read.prolix(), self.test_data_1d) + + def test_compression_settings(self): + """Test compression settings""" + terse = Terse() + + # Test block size + terse.set_block_size(1024) + self.assertEqual(terse.block_size(), 1024) + + # Test fast mode + terse.set_fast(True) + self.assertTrue(terse.fast()) + + # Test small mode + terse.set_small(True) + self.assertTrue(terse.small()) + + # Test degree of parallelism + terse.set_dop(0.5) + self.assertEqual(terse.dop(), 0.5) + + def test_error_handling(self): + """Test error handling""" + terse = Terse(self.test_data_1d) + + # Test invalid frame access + with self.assertRaises(IndexError): + terse.at(999) + + # Test invalid frame insertion + with self.assertRaises(ValueError): + terse.insert(999, self.test_data_1d) + + # Test dimension mismatch + different_shape = np.array([[1, 2], [3, 4]], dtype=np.int32) + with self.assertRaises(ValueError): + terse.push_back(different_shape) + + def test_memory_management(self): + """Test memory management features""" + terse = Terse(self.test_data_1d) + original_size = terse.terse_size + terse.shrink_to_fit() + self.assertLessEqual(terse.terse_size, original_size) + +if __name__ == '__main__': + unittest.main(verbosity=2) \ No newline at end of file diff --git a/lib/trpx2.0/pyterse/CMakeLists.txt b/lib/trpx2.0/pyterse/CMakeLists.txt new file mode 100755 index 00000000..6a2509ec --- /dev/null +++ b/lib/trpx2.0/pyterse/CMakeLists.txt @@ -0,0 +1,47 @@ +cmake_minimum_required(VERSION 3.15) +project(pyterse LANGUAGES CXX) + +set(CMAKE_CXX_STANDARD 20) +set(CMAKE_CXX_STANDARD_REQUIRED ON) +set(CMAKE_VERBOSE_MAKEFILE ON) + +if(NOT CMAKE_BUILD_TYPE) + set(CMAKE_BUILD_TYPE Release) +endif() + +include_directories(${CMAKE_CURRENT_SOURCE_DIR}/../include) + +find_package(pybind11 CONFIG REQUIRED) + +message(STATUS "pybind11_INCLUDE_DIRS: ${pybind11_INCLUDE_DIRS}") +message(STATUS "pybind11_VERSION: ${pybind11_VERSION}") + +pybind11_add_module(pyterse src/pyterse.cpp) + +target_include_directories(pyterse PRIVATE + ${CMAKE_CURRENT_SOURCE_DIR}/../include +) + +set_target_properties(pyterse PROPERTIES OUTPUT_NAME "pyterse") +set_target_properties(pyterse PROPERTIES PREFIX "") + + +if(MSVC) + target_compile_options(pyterse PRIVATE /O2) +else() + target_compile_options(pyterse PRIVATE -O3 -fpermissive) +endif() + + +if(APPLE) + set_target_properties(pyterse PROPERTIES LINK_FLAGS "-undefined dynamic_lookup") +endif() + +target_link_libraries(pyterse PRIVATE pybind11::pybind11) + +# Install the module +install(TARGETS pyterse + LIBRARY DESTINATION . + ARCHIVE DESTINATION . + RUNTIME DESTINATION . +) diff --git a/lib/trpx2.0/pyterse/src/pyterse.cpp b/lib/trpx2.0/pyterse/src/pyterse.cpp new file mode 100644 index 00000000..2bb3a9c6 --- /dev/null +++ b/lib/trpx2.0/pyterse/src/pyterse.cpp @@ -0,0 +1,422 @@ + + /** + * @file pyterse.cpp + * @author Senik Matinyan + * @date 2024 + * @brief Python bindings for Terse compression library using pybind11 + * + * This file provides Python bindings for the Terse compression library, + * allowing seamless integration with NumPy arrays and Python I/O streams. + */ + + + #include + #include + #include "Concurrent.hpp" + #include "Terse.hpp" + #include + #include + #include + #include + #include + #include + + + PYBIND11_MODULE(pyterse, m) + { + m.doc() = "Python bindings for Terse"; + + using namespace jpa; + namespace py = pybind11; + + /** + * @class python_streambuf + * @brief Custom stream buffer implementation that wraps Python file-like objects + * + * This class allows C++ streams to interact with Python file-like objects by + * implementing the necessary streambuf interface methods. + */ + class python_streambuf : public std::streambuf { + public: + /** + * @brief Construct a new python_streambuf with a Python stream + * @param py_stream Python file-like object with read/write methods + */ + explicit python_streambuf(py::object py_stream) : py_stream(std::move(py_stream)), buffer(1024) {} + + protected: + /** + * @brief Write characters to the underlying Python stream + * @param s Pointer to the character buffer + * @param n Number of characters to write + * @return Number of characters successfully written + * @throws Various exceptions if writing fails + */ + std::streamsize xsputn(char const* s, std::streamsize const n) override { + try { + py_stream.attr("write")(py::bytes(s, static_cast(n))); + return n; + } + catch (const py::error_already_set& e) { throw; } + catch (const std::exception& e) { throw py::value_error(std::string("[xsputn] Error: ") + e.what()); } + catch (...) { throw std::runtime_error("[xsputn] Unknown error occurred"); } + } + + /** + * @brief Read characters from the underlying Python stream + * @param s Pointer to the destination buffer + * @param n Maximum number of characters to read + * @return Number of characters successfully read + * @throws Various exceptions if reading fails + */ + std::streamsize xsgetn(char* s, std::streamsize const n) override { + if (n == 0) return 0; + try { + py::bytes data = py_stream.attr("read")(n); + std::string str_data = data; + if (str_data.empty()) throw py::value_error("[xsgetn] No data read; stream may be at EOF."); + std::streamsize size = std::min(n, str_data.size()); + std::copy(str_data.begin(), str_data.begin() + size, s); + return size; + } + catch (const py::error_already_set& e) { throw; } + catch (const std::exception& e) { throw py::value_error(std::string("[xsgetn] Error: ") + e.what()); } + catch (...) { throw std::runtime_error("[xsgetn] Unknown error occurred"); } + } + + /** + * @brief Seek to a position within the Python stream + * @param off Offset to seek to + * @param way Direction of seeking (beginning, current, end) + * @param which I/O mode (not used) + * @return New absolute position after seeking + * @throws Various exceptions if seeking fails + */ + std::streampos seekoff(std::streamoff const off, std::ios_base::seekdir const way, std::ios_base::openmode const which) override { + try { + int whence = (way == std::ios_base::beg) ? 0 : + (way == std::ios_base::cur) ? 1 : + (way == std::ios_base::end) ? 2 : -1; + if (whence == -1) throw py::value_error("Invalid seek direction"); + py_stream.attr("seek")(off, whence); + auto new_pos = py_stream.attr("tell")().cast(); + return std::streampos(new_pos); + } + catch (const py::error_already_set& e) { throw; } + catch (const std::exception& e) { throw py::value_error(std::string("[seekoff] Error: ") + e.what()); } + catch (...) { throw std::runtime_error("[seekoff] Unknown error occurred"); } + } + + /** + * @brief Read a character when the input buffer is empty + * @return The next character as an int, or EOF if no more characters + */ + int underflow() override { + char c; + if (xsgetn(&c, 1) == 0) return traits_type::eof(); + setg(&c, &c, &c + 1); + return traits_type::to_int_type(c); + } + + private: + py::object py_stream; ///< Python file-like object + std::vector buffer; ///< Internal buffer for streambuf operations + }; + + /** + * @class python_ostream + * @brief C++ output stream that wraps a Python file-like object + */ + class python_ostream : public std::ostream { + public: + /** + * @brief Construct a new python_ostream with a Python stream + * @param py_stream Python file-like object with write method + */ + explicit python_ostream(py::object py_stream) : std::ostream(nullptr), buffer(std::move(py_stream)) { this->init(&buffer); } + private: + python_streambuf buffer; ///< Stream buffer that interfaces with Python + }; + + /** + * @class python_istream + * @brief C++ input stream that wraps a Python file-like object + */ + class python_istream : public std::istream { + public: + /** + * @brief Construct a new python_istream with a Python stream + * @param py_stream Python file-like object with read method + */ + explicit python_istream(py::object py_stream) : std::istream(nullptr), buffer(std::move(py_stream)) { this->init(&buffer); } + private: + python_streambuf buffer; ///< Stream buffer that interfaces with Python + }; + + /** + * @brief Select appropriate Terse function based on NumPy array dtype + * @param data NumPy array to analyze + * @param func Function to call with the determined C++ type + * @return Result of the function call + * @throws py::type_error if the data type is not supported + */ + auto select_terse_func = [](py::array& data, auto&& func) { + // Handle float types + if (data.dtype().kind() == 'f') { + switch (data.dtype().itemsize()) { + case 4: return func(float{}); + case 8: return func(double{}); + default: + throw py::type_error("Terse supports float32 and float64 types"); + } + } + else if (data.dtype().kind() == 'i' || data.dtype().kind() == 'u') { + switch (data.dtype().itemsize()) { + case 1: return data.dtype().kind() == 'i' ? func(std::int8_t {}) : func(std::uint8_t {}); + case 2: return data.dtype().kind() == 'i' ? func(std::int16_t{}) : func(std::uint16_t{}); + case 4: return data.dtype().kind() == 'i' ? func(std::int32_t{}) : func(std::uint32_t{}); + case 8: return data.dtype().kind() == 'i' ? func(std::int64_t{}) : func(std::uint64_t{}); + default: + throw py::type_error("Terse only supports compression of integral values with 1, 2, 4 or 8 bytes"); + } + } + else { + throw py::type_error("Terse only supports integral and floating point types"); + } + }; + + /** + * @brief Determine appropriate NumPy dtype based on Terse object properties + * @param terse Terse object to analyze + * @return Corresponding NumPy dtype + */ + auto pydtype_of_terse = [&] (Terse& terse) { + if (terse.is_float()) { + return py::dtype::of(); + } + switch (terse.bits_per_val()) { + case 8 : return terse.is_signed() ? py::dtype::of() : py::dtype::of(); + case 16: return terse.is_signed() ? py::dtype::of() : py::dtype::of(); + case 32: return terse.is_signed() ? py::dtype::of() : py::dtype::of(); + default: return terse.is_signed() ? py::dtype::of() : py::dtype::of(); + } + }; + + /** + * @brief Calculate the shape of a NumPy array from a Terse object + * @param terse Terse object to analyze + * @return Vector of dimensions for the NumPy array + */ + auto pyshape_of_terse = [] (Terse& terse) { + if (terse.dim().empty()) + return std::vector({terse.size()}); + if (terse.number_of_frames() == 1) + return terse.dim(); + std::vector shape = terse.dim(); + shape.insert(shape.begin(), terse.number_of_frames()); + return shape; + }; + + /** + * @brief Insert data from a NumPy array into a Terse object + * @tparam T C++ data type of the array elements + * @param terse Target Terse object + * @param pos Position to insert at + * @param data Source NumPy array + * @param mode Terse compression mode + * @param type Type tag used for template deduction + * @throws py::value_error if dimensions do not match + */ + auto insert = [] (Terse& terse, std::size_t pos, py::array& data, Terse_mode mode, T type) { + auto const buf = data.request(); + auto const shape = std::vector(buf.shape.begin(), buf.shape.end()); + auto dim = std::vector(shape.begin() + (shape.size() <= 2 ? 0 : 1), shape.end()); + if (terse.number_of_frames() == 0) terse.dim(dim); + else if (dim != terse.dim()) + throw py::value_error("Dimension mismatch: Terse cannot insert data because of shape mismatch."); + size_t const frame_size = std::max(static_cast(shape.back()), + static_cast(std::accumulate(dim.begin(), dim.end(), 1ul, std::multiplies<>()))); + size_t const num_frames = shape.size() <= 2 ? 1 : shape[0]; + T* base_ptr = static_cast(buf.ptr); + for (std::size_t i = 0; i < num_frames; ++i) + terse.insert(pos + i, std::span(base_ptr + i * frame_size, frame_size), mode); + terse.shrink_to_fit(); + }; + + /** + * @brief Decompress Terse data into a NumPy array + * @tparam T C++ data type of the array elements + * @param terse Source Terse object + * @param data Target NumPy array + * @param type Type tag used for template deduction + * @throws py::value_error if dimensions do not match + */ + auto prolix = [&] (Terse& terse, py::array& data, T type) { + auto const buf = data.request(); + auto const shape = std::vector(buf.shape.begin(), buf.shape.end()); + if (terse.dim().empty()) { + if (shape.size() != 1 || shape[0] != terse.size()) + throw py::value_error("Dimension mismatch: Terse cannot expand data into array because of shape mismatch."); + } + else if (terse.number_of_frames() == 1) { + if (terse.dim() != shape) + throw py::value_error("Dimension mismatch: Terse cannot expand data into array because of shape mismatch."); + } + else if (terse.dim() != std::vector(shape.begin() + 1, shape.end()) && shape[0] != terse.number_of_frames()) + throw py::value_error("Dimension mismatch: Terse cannot expand data into array because of shape mismatch."); + T* base_ptr = static_cast(buf.ptr); + std::vector> frames; + for (std::size_t i = 0; i != terse.number_of_frames(); ++i) + frames.push_back(std::span(base_ptr + i * terse.size(), terse.size())); + terse.prolix(frames); + }; + + + /** + * @brief Python enum for Terse compression modes + */ + py::enum_(m, "TerseMode") + .value("SIGNED", Terse_mode::Signed, "Signed integer mode") + .value("UNSIGNED", Terse_mode::Unsigned, "Unsigned integer mode") + .value("SMALL_UNSIGNED", Terse_mode::Small_unsigned, "Small unsigned integer mode") + .value("DEFAULT", Terse_mode::Default, "Default mode based on data type") + .export_values(); + + /** + * @brief Python bindings for the Terse class + */ + py::class_, std::shared_ptr>>(m, "Terse") + .def(py::init<>(), "Create an empty Terse object") + + .def(py::init([&](py::array data, Terse_mode mode) -> std::shared_ptr> { + auto terse_ptr = std::make_shared>(); + select_terse_func(data, [&](auto Type) { return insert(*terse_ptr, 0, data, mode, Type); }); + return terse_ptr; + }), py::arg("data"), py::arg("mode") = Terse_mode::Default, + "Initialize a Terse object with a NumPy array and a Terse_mode") + + .def(py::init([](py::object py_stream) -> std::shared_ptr> { + python_istream stream(py_stream); + return std::make_shared>(stream); + }), py::arg("stream"), + "Create a Terse object from a binary input stream. The stream must contain Terse data preceded by the required XML header.") + + .def("insert", [&](Terse& terse, std::size_t pos, py::array data, Terse_mode mode) { + if (pos > terse.number_of_frames()) + throw std::invalid_argument("Cannot insert at this position: the terse object has fewer frames that the position requested"); + select_terse_func(data, [&](auto Type) { return insert(terse, pos, data, mode, Type); }); + }, py::arg("pos"), py::arg("data"), py::arg("mode") = Terse_mode::Default, + "Insert data into the Terse object at the specified position.") + + .def("push_back", [&](Terse& terse, py::array data, Terse_mode mode) { + select_terse_func(data, [&](auto Type) { return insert(terse, terse.number_of_frames(), data, mode, Type); }); + }, py::arg("data"), py::arg("mode") = Terse_mode::Default, + "Append data at the end of the Terse object.") + + .def("prolix", [&](Terse& terse, py::array& data) { + select_terse_func(data, [&](auto Type) { return prolix(terse, data, Type); }); + return data; + }, py::arg("data"), + "Decompress data into the provided data array.") + + .def("prolix", [&](Terse& terse) -> py::array { + py::array data(pydtype_of_terse(terse), pyshape_of_terse(terse)); + select_terse_func(data, [&](auto Type) { return prolix(terse, data, Type); }); + return data; + }, + "Decompress data into a new data array.") + + .def("at", [](Terse& self, std::size_t pos) -> std::shared_ptr> { + if (pos >= self.number_of_frames()) + throw py::index_error("Requested frame not present: index too high."); + return std::make_shared>(self.at(pos)); + }, py::arg("pos"), + "Return a Terse object representing a single frame at the specified position.") + + .def("erase", &Terse::erase, py::arg("pos"), + "Erase the frame at the specified position.") + + .def("write", [](Terse& self, py::object stream) { + python_ostream out(stream); + self.write(out); + }, py::arg("stream"), + "Write Terse data to a binary output stream.") + + .def("save", [](Terse& self, const std::string& filename) { + std::ofstream outfile(filename, std::ios::binary); + if (!outfile.is_open()) { + throw std::runtime_error("Failed to open file for writing"); + } + self.write(outfile); + outfile.close(); + }, + py::arg("filename"), + "Save Terse data to a file.") + + .def_static("load", [](const std::string& filename) -> std::shared_ptr> { + std::ifstream infile(filename, std::ios::binary); + if (!infile.is_open()) { + throw std::runtime_error("Failed to open file for reading"); + } + auto obj = std::make_shared>(infile); + return obj; + }, py::arg("filename"), + "Load Terse data from a file.") + + .def_property_readonly("size", &Terse::size, + "Get the number of values in each frame.") + .def_property_readonly("number_of_frames", &Terse::number_of_frames, + "Get the number of frames in the Terse object.") + .def_property_readonly("is_signed", &Terse::is_signed, + "Check if the data is stored as signed values.") + .def_property_readonly("bits_per_val", &Terse::bits_per_val, + "Get the number of bits used per value in the compressed data.") + .def_property_readonly("terse_size", &Terse::terse_size, + "Get the size of the compressed data in bytes.") + .def_property_readonly("file_size", &Terse::file_size, + "Get the total file size including headers in bytes.") + .def_property_readonly("number_of_bytes", &Terse::terse_size, + "Alias for terse_size.") + + .def("metadata", py::overload_cast(&Terse::metadata, py::const_), py::arg("frame") = 0, + "Get the metadata for the specified frame.") + .def("set_metadata", py::overload_cast(&Terse::metadata), + py::arg("frame"), py::arg("metadata"), + "Set the metadata for the specified frame.") + .def("dim", py::overload_cast<>(&Terse::dim, py::const_), + "Get the dimensions of the frames.") + .def("set_dim", [](Terse& self, const std::vector& dims) { + if (!self.dim().empty() && std::accumulate(dims.begin(), dims.end(), 1ul, std::multiplies<>()) != self.size()) + throw py::value_error("The total number of elements must remain unchanged when setting dimensions, unless the Terse object is empty."); + self.dim(dims); + }, py::arg("dimensions"), + "Set the dimensions of the frames.") + .def("fractional_precision", py::overload_cast<>(&Terse::fractional_precision, py::const_), + "Get the fractional precision for floating-point data.") + .def("set_fractional_precision", [](Terse& self, double precision) { + self.fractional_precision(precision); + }, py::arg("precision"), + "Set the fractional precision for floating-point data.") + .def("block_size", py::overload_cast<>(&Terse::block_size, py::const_), + "Get the block size used for compression.") + .def("set_block_size", [](Terse& self, std::size_t block_size) { + if (self.number_of_frames() > 0) + throw py::value_error("Cannot set the block size after frames have been added to the Terse object."); + self.block_size(block_size); + }, py::arg("block_size"), + "Set the block size used for compression.") + .def("fast", py::overload_cast<>(&Terse::fast, py::const_), + "Check if fast mode is enabled.") + .def("set_fast", py::overload_cast(&Terse::fast), py::arg("value"), + "Enable or disable fast mode.") + .def("small", py::overload_cast<>(&Terse::small, py::const_), + "Check if small mode is enabled.") + .def("set_small", py::overload_cast(&Terse::small), py::arg("value"), + "Enable or disable small mode.") + .def("dop", py::overload_cast<>(&Terse::dop, py::const_), + "Get the degree of parallelism.") + .def("set_dop", py::overload_cast(&Terse::dop), py::arg("value"), + "Set the degree of parallelism.") + .def("shrink_to_fit", &Terse::shrink_to_fit, + "Reduce memory usage by freeing unused capacity."); + } \ No newline at end of file diff --git a/lib/trpx2.0/tersecodec/CMakeLists.txt b/lib/trpx2.0/tersecodec/CMakeLists.txt new file mode 100755 index 00000000..6189be25 --- /dev/null +++ b/lib/trpx2.0/tersecodec/CMakeLists.txt @@ -0,0 +1,56 @@ +cmake_minimum_required(VERSION 3.12) +project(tersecodec) + +set(CMAKE_CXX_STANDARD 20) +set(CMAKE_BUILD_TYPE Release) +set(CMAKE_CXX_FLAGS_RELEASE "${CMAKE_CXX_FLAGS_RELEASE} -O3 ") + + +cmake_policy(SET CMP0074 NEW) + +set(HDF5_ROOT "$ENV{HDF5_ROOT}") +if(NOT HDF5_ROOT) + set(HDF5_ROOT "$ENV{CONDA_PREFIX}") +endif() + +set(TERSE_INCLUDE_DIR "${CMAKE_CURRENT_SOURCE_DIR}/../include") +include_directories(${TERSE_INCLUDE_DIR}) + + +find_package(HDF5 REQUIRED COMPONENTS C CXX HL) +if(HDF5_FOUND) + message(STATUS "HDF5 include dirs: ${HDF5_INCLUDE_DIRS}") + message(STATUS "HDF5 library dirs: ${HDF5_LIBRARY_DIRS}") +else() + message(FATAL_ERROR "HDF5 not found in ${HDF5_ROOT}") +endif() + +include_directories(${HDF5_INCLUDE_DIRS}) + + +set(PLUGIN_NAME "h5terse") +add_library(${PLUGIN_NAME} SHARED src/filter.cpp) + + +target_include_directories(${PLUGIN_NAME} PUBLIC ${TERSE_INCLUDE_DIR}) +target_compile_options(${PLUGIN_NAME} PRIVATE -fPIC) +target_link_libraries(${PLUGIN_NAME} ${HDF5_LIBRARIES}) + + +set_target_properties(${PLUGIN_NAME} PROPERTIES + OUTPUT_NAME "${PLUGIN_NAME}" + PREFIX "lib" + SUFFIX ".so" + BUILD_WITH_INSTALL_RPATH TRUE + INSTALL_RPATH "${HDF5_LIBRARY_DIRS}" + INSTALL_RPATH_USE_LINK_PATH TRUE +) + +# Install to standard HDF5 plugin location +if(CMAKE_INSTALL_PREFIX_INITIALIZED_TO_DEFAULT) + set(CMAKE_INSTALL_PREFIX "${HDF5_ROOT}" CACHE PATH "Default install path" FORCE) +endif() + +install(TARGETS ${PLUGIN_NAME} + LIBRARY DESTINATION lib/hdf5/plugins +) \ No newline at end of file diff --git a/lib/trpx2.0/tersecodec/src/filter.cpp b/lib/trpx2.0/tersecodec/src/filter.cpp new file mode 100644 index 00000000..bfc7a68c --- /dev/null +++ b/lib/trpx2.0/tersecodec/src/filter.cpp @@ -0,0 +1,366 @@ +// +// Written by: Senik Matinyan, 2024 +// + +#include "hdf5.h" +#include "H5PLextern.h" + +#include +#include +#include +#include +#include +#include +#include +#include + +#include "Concurrent.hpp" +#include "Terse.hpp" + +#define TERSE_FILTER_ID 32029 +#define TERSE_DEFAULT_CHUNK_SIZE (1 << 18) + +#define TYPE_CODE_INT16 0 +#define TYPE_CODE_UINT16 1 +#define TYPE_CODE_INT32 2 +#define TYPE_CODE_UINT32 3 +#define TYPE_CODE_INT8 4 +#define TYPE_CODE_UINT8 5 + +template +std::size_t hdf5_buffer_to_terse( + void** data, + std::size_t* size, + std::size_t chunk_size) { + struct Fixed_ostreambuf : public std::streambuf { + Fixed_ostreambuf(char* buffer, std::size_t size) { + setp(buffer, buffer + size); + } + }; + + std::span buffer(static_cast(*data), *size / sizeof(T)); + + jpa::Terse terse_chunks; + jpa::Terse terse_rest; + + if (chunk_size == 0) { + terse_chunks.push_back(buffer); + } else { + std::size_t pos = 0; + + for (; pos + chunk_size < buffer.size(); pos += chunk_size) { + terse_chunks.push_back( + std::span(buffer.data() + pos, chunk_size)); + } + + terse_rest.push_back( + std::span(buffer.data() + pos, buffer.size() - pos)); + } + + std::size_t terse_size = + terse_chunks.file_size() + terse_rest.file_size(); + + void* new_buf = malloc(terse_size); + if (!new_buf) { + std::cerr + << "Allocation of memory failed while compressing data in " + "hdf5_buffer_to_terse" + << std::endl; + return 0; + } + + Fixed_ostreambuf oss_buffer( + static_cast(new_buf), + terse_size); + std::ostream oss(&oss_buffer); + + terse_chunks.write(oss); + terse_rest.write(oss); + + free(buffer.data()); + *data = new_buf; + + return *size = terse_size; +} + +template +std::size_t terse_to_hdf5_buffer(void** data, std::size_t* size) { + struct Fixed_istreambuf : public std::streambuf { + char* buffer; + std::size_t buffer_size; + + Fixed_istreambuf(char* buffer, std::size_t size) + : buffer(buffer), buffer_size(size) { + setg(buffer, buffer, buffer + size); + } + + std::streampos seekoff( + std::streamoff off, + std::ios_base::seekdir way, + std::ios_base::openmode which = + std::ios_base::in) override { + switch (way) { + case std::ios_base::beg: + setg(buffer, buffer + off, buffer + buffer_size); + break; + case std::ios_base::end: + setg( + buffer, + buffer + buffer_size + off, + buffer + buffer_size); + break; + default: + gbump(static_cast(off)); + break; + } + + return gptr() - eback(); + } + + std::streampos seekpos( + std::streampos pos, + std::ios_base::openmode which = + std::ios_base::in) override { + return seekoff(pos, std::ios_base::beg, which); + } + }; + + std::span buffer(static_cast(*data), *size); + + Fixed_istreambuf iss_buffer( + static_cast(buffer.data()), + buffer.size()); + std::istream iss(&iss_buffer); + + jpa::Terse terse_chunks(iss); + std::size_t chunk_size = terse_chunks.size(); + + jpa::Terse terse_rest; + if (!iss.eof()) { + terse_rest = jpa::Terse(iss); + } + + std::size_t prolix_size = + terse_chunks.size() * terse_chunks.number_of_frames() + + terse_rest.size(); + + void* new_buf = malloc(prolix_size * sizeof(T)); + if (!new_buf) { + std::cerr + << "Allocation of memory failed while expanding data in " + "terse_to_hdf5_buffer" + << std::endl; + return 0; + } + + free(*data); + *data = new_buf; + + std::span prolix( + static_cast(new_buf), + prolix_size); + + terse_chunks.prolix( + std::span( + prolix.data(), + chunk_size * terse_chunks.number_of_frames())); + + terse_rest.prolix( + std::span( + prolix.data() + + chunk_size * terse_chunks.number_of_frames(), + terse_rest.size())); + + return *size = prolix_size * sizeof(T); +} + +size_t Terse_filter_cpp( + unsigned int flags, + size_t cd_nelmts, + const unsigned int cd_values[], + size_t* bufsize, + void** buf) { + if (!buf + || !bufsize + || cd_nelmts == 0 + || cd_values == nullptr + || cd_values[0] > 5) { + std::cerr << "Invalid arguments provided to Terse_filter" + << std::endl; + return 0; + } + + unsigned int data_type_code = cd_values[0]; + + if (flags & H5Z_FLAG_REVERSE) { + switch (data_type_code) { + case TYPE_CODE_INT16: + return terse_to_hdf5_buffer(buf, bufsize); + case TYPE_CODE_UINT16: + return terse_to_hdf5_buffer(buf, bufsize); + case TYPE_CODE_INT32: + return terse_to_hdf5_buffer(buf, bufsize); + case TYPE_CODE_UINT32: + return terse_to_hdf5_buffer(buf, bufsize); + case TYPE_CODE_INT8: + return terse_to_hdf5_buffer(buf, bufsize); + case TYPE_CODE_UINT8: + return terse_to_hdf5_buffer(buf, bufsize); + } + } else { + switch (data_type_code) { + case TYPE_CODE_INT16: + return hdf5_buffer_to_terse( + buf, + bufsize, + TERSE_DEFAULT_CHUNK_SIZE); + case TYPE_CODE_UINT16: + return hdf5_buffer_to_terse( + buf, + bufsize, + TERSE_DEFAULT_CHUNK_SIZE); + case TYPE_CODE_INT32: + return hdf5_buffer_to_terse( + buf, + bufsize, + TERSE_DEFAULT_CHUNK_SIZE); + case TYPE_CODE_UINT32: + return hdf5_buffer_to_terse( + buf, + bufsize, + TERSE_DEFAULT_CHUNK_SIZE); + case TYPE_CODE_INT8: + return hdf5_buffer_to_terse( + buf, + bufsize, + TERSE_DEFAULT_CHUNK_SIZE); + case TYPE_CODE_UINT8: + return hdf5_buffer_to_terse( + buf, + bufsize, + TERSE_DEFAULT_CHUNK_SIZE); + } + } + + return 0; +} + +extern "C" size_t Terse_filter( + unsigned int flags, + size_t cd_nelmts, + const unsigned int cd_values[], + size_t nbytes, + size_t* bufsize, + void** buf) { + return Terse_filter_cpp( + flags, + cd_nelmts, + cd_values, + bufsize, + buf); +} + +static herr_t Terse_set_local( + hid_t dcpl_id, + hid_t type_id, + hid_t space_id) { + if (H5Tget_class(type_id) != H5T_INTEGER) { + std::cerr + << "Terse filter only supports integer datasets" + << std::endl; + return -1; + } + + H5T_sign_t sign = H5Tget_sign(type_id); + if (sign == H5T_SGN_ERROR) { + return -1; + } + + bool is_signed = sign != H5T_SGN_NONE; + unsigned int data_type_code; + + switch (H5Tget_size(type_id)) { + case 1: + data_type_code = + is_signed ? TYPE_CODE_INT8 : TYPE_CODE_UINT8; + break; + case 2: + data_type_code = + is_signed ? TYPE_CODE_INT16 : TYPE_CODE_UINT16; + break; + case 4: + data_type_code = + is_signed ? TYPE_CODE_INT32 : TYPE_CODE_UINT32; + break; + default: + std::cerr + << "Terse filter only supports 8-, 16- and 32-bit " + "integer datasets" + << std::endl; + return -1; + } + + unsigned int flags; + size_t cd_nelmts = 1; + unsigned int cd_values[1] = {0}; + + if (H5Pget_filter_by_id2( + dcpl_id, + TERSE_FILTER_ID, + &flags, + &cd_nelmts, + cd_values, + 0, + nullptr, + nullptr) < 0) { + return -1; + } + + cd_values[0] = data_type_code; + + if (H5Pmodify_filter( + dcpl_id, + TERSE_FILTER_ID, + flags, + 1, + cd_values) < 0) { + return -1; + } + + return 1; +} + +extern "C" herr_t register_terse_filter() { + H5Z_class2_t filter_class = { + H5Z_CLASS_T_VERS, + TERSE_FILTER_ID, + 1, + 1, + "TERSE", + nullptr, + Terse_set_local, + Terse_filter + }; + + return H5Zregister(&filter_class); +} + +extern "C" const void* H5PLget_plugin_info(void) { + static H5Z_class2_t filter_class = { + H5Z_CLASS_T_VERS, + TERSE_FILTER_ID, + H5Z_FILTER_CONFIG_ENCODE_ENABLED + | H5Z_FILTER_CONFIG_DECODE_ENABLED, + 1, + "TERSE", + nullptr, + Terse_set_local, + Terse_filter + }; + + return static_cast(&filter_class); +} + +extern "C" H5PL_type_t H5PLget_plugin_type(void) { + return H5PL_TYPE_FILTER; +} diff --git a/setup.py b/setup.py index 24a8f1dd..b5e4bf16 100644 --- a/setup.py +++ b/setup.py @@ -705,6 +705,17 @@ def build_extensions(self): arg for arg in e.extra_link_args if arg.startswith(prefix) ] + # Patch C++20 flag for older gcc/clang support + if "-std=c++20" in e.extra_compile_args: + if config.cpp20_compile_arg is None: + raise RuntimeError( + "Cannot compile C++20 code with current compiler" + ) + e.extra_compile_args = [ + config.cpp20_compile_arg if arg == "-std=c++20" else arg + for arg in e.extra_compile_args + ] + build_ext.build_extensions(self) @@ -1646,6 +1657,34 @@ def _get_sperr_plugin(): PLUGIN_LIB_DEPENDENCIES["sperr"] = ("sperr", "sperr_filter") +def _get_terse_plugin(): + """TERSE/PROLIX (TRPX v2.0) plugin build config + + Source from https://github.com/senikm/trpx2.0 + """ + trpx_dir = "lib/trpx2.0" + + extra_compile_args = ["-O3", "-std=c++20", "-pthread"] + extra_compile_args += ["/Ox", "/std:c++20"] + # Concurrent.hpp uses std::clamp without including , + # which fails with some libstdc++ versions: force include it. + # Use joined forms since compile args must be single tokens starting with - or / + extra_compile_args += ["-includealgorithm", "/FIalgorithm"] + extra_link_args = ["-pthread"] + + return HDF5PluginExtension( + "hdf5plugin.plugins.libh5terse", + sources=[f"{trpx_dir}/tersecodec/src/filter.cpp"], + include_dirs=[f"{trpx_dir}/include"], + extra_compile_args=extra_compile_args, + extra_link_args=extra_link_args, + cpp20_required=True, + ) + + +PLUGIN_LIB_DEPENDENCIES["terse"] = () + + _EMBEDDED_PLUGIN_EXTENSIONS = { "blosc": _get_blosc_plugin, "blosc2": _get_blosc2_plugin, @@ -1657,6 +1696,7 @@ def _get_sperr_plugin(): "sperr": _get_sperr_plugin, "sz": _get_sz_plugin, "sz3": _get_sz3_plugin, + "terse": _get_terse_plugin, "zfp": _get_h5zfp_plugin, "zstd": _get_zstd_plugin, } diff --git a/src/hdf5plugin/__init__.py b/src/hdf5plugin/__init__.py index 1255ddaa..e652a8e1 100644 --- a/src/hdf5plugin/__init__.py +++ b/src/hdf5plugin/__init__.py @@ -42,6 +42,7 @@ SZ3, SZ3_ID, SZ_ID, + TERSE_ID, ZFP_ID, ZSTD_ID, Bitshuffle, @@ -51,6 +52,7 @@ FciDecomp, Htj2k, Sperr, + Terse, Zfp, Zstd, ) @@ -75,6 +77,7 @@ "Sperr", "SZ", "SZ3", + "Terse", "Zfp", "Zstd", # Utilities diff --git a/src/hdf5plugin/_filters.py b/src/hdf5plugin/_filters.py index ea8ace21..6617837a 100644 --- a/src/hdf5plugin/_filters.py +++ b/src/hdf5plugin/_filters.py @@ -73,6 +73,9 @@ SPERR_ID = 32028 """SPERR filter ID""" +TERSE_ID = 32029 +"""TERSE/PROLIX (TRPX) filter ID""" + class FilterBase(h5py.filters.FilterRefBase): filter_id: int @@ -1072,6 +1075,47 @@ def __unpack_options(cls, meta: int, ret: int) -> tuple[int, float, bool, int]: return mode, quality, swap, missing_value_mode +class Terse(FilterBase): + """``h5py.Group.create_dataset``'s compression arguments for using the TERSE/PROLIX (TRPX) filter. + + .. code-block:: python + + f = h5py.File('test.h5', 'w') + f.create_dataset( + 'terse', + data=numpy.arange(100, dtype=numpy.uint16), + compression=hdf5plugin.Terse()) + f.close() + + This filter compresses integer data losslessly using the + `TERSE/PROLIX v2.0 `_ algorithm, + which is tailored for diffraction data and other greyscale images with + a high dynamic range. + + It only supports datasets of signed or unsigned 8-, 16- or 32-bit integers + (``int8``, ``uint8``, ``int16``, ``uint16``, ``int32``, ``uint32``). + """ + + filter_name = "terse" + filter_id = TERSE_ID + + def __init__(self) -> None: + if not build_config.cpp20: + logger.error( + "The Terse filter is not available as hdf5plugin was not built with C++20.\n" + "You may need to reinstall hdf5plugin with a recent version of pip, or rebuild it with a newer compiler." + ) + super().__init__(filter_options=(), config={}) + + @classmethod + def _from_filter_options(cls, filter_options: tuple[int, ...]) -> Terse: + """Returns compression arguments from HDF5 compression filters options + + :raises ValueError: Unsupported filter_options + """ + return cls() + + def _sz_pack_float64(value: float) -> tuple[int, int]: # Pack as big-endian IEEE 754 double packed = struct.pack(">d", value) @@ -1386,6 +1430,7 @@ def _from_filter_options(cls, filter_options: tuple[int, ...]) -> Zstd: Sperr, SZ, SZ3, + Terse, Zfp, Zstd, ) diff --git a/src/hdf5plugin/test.py b/src/hdf5plugin/test.py index 69620f91..f59cbdef 100644 --- a/src/hdf5plugin/test.py +++ b/src/hdf5plugin/test.py @@ -68,6 +68,7 @@ def should_test(filter_name: str) -> bool: "sperr": hdf5plugin.Sperr, "sz": hdf5plugin.SZ, "sz3": hdf5plugin.SZ3, + "terse": hdf5plugin.Terse, "zfp": hdf5plugin.Zfp, "htj2k": hdf5plugin.Htj2k, "zstd": hdf5plugin.Zstd, @@ -381,6 +382,29 @@ def testSZ3(self): with self.subTest(options=options, dtype=dtype): self._test("sz3", dtype=dtype, lossless=False, options=options) + @unittest.skipUnless(should_test("terse"), "Terse filter not available") + def testTerse(self): + """Write/read test with terse filter plugin""" + # Test with supported datatypes and their type code + for dtype, type_code in ( + (numpy.int16, 0), + (numpy.uint16, 1), + (numpy.int32, 2), + (numpy.uint32, 3), + (numpy.int8, 4), + (numpy.uint8, 5), + ): + with self.subTest(dtype=dtype): + filter_ = self._test("terse", dtype=dtype) + # The filter stores the type code of the dataset's dtype + self.assertEqual(filter_[2], (type_code,)) + + # Test that unsupported datatypes are rejected at dataset creation + for unsupported_dtype in (numpy.float32, numpy.int64): + with self.subTest(dtype=unsupported_dtype): + with self.assertRaises(ValueError): + self._test("terse", dtype=unsupported_dtype) + @unittest.skipUnless(should_test("zfp"), "ZFP filter not available") def testZfp(self): """Write/read test with zfp filter plugin""" @@ -694,6 +718,10 @@ def testSZ3(self): compression_filter.filter_options, expected_filter.filter_options ) + def testTerse(self): + compression_filter = hdf5plugin.Terse._from_filter_options((1,)) + self.assertEqual(compression_filter.filter_options, ()) + def testZfp(self): for filter_options, expected_filter in ( ( @@ -821,6 +849,11 @@ def testSZ3(self): data = numpy.arange(256**2, dtype=numpy.float32).reshape(256, 256) self._test(hdf5plugin.SZ3(), data) + @unittest.skipUnless(should_test("terse"), "Terse filter not available") + def testTerse(self): + data = numpy.arange(256**2, dtype=numpy.uint16).reshape(256, 256) + self._test(hdf5plugin.Terse(), data) + @unittest.skipUnless(should_test("zfp"), "Zfp filter not available") def testZfp(self): data = numpy.arange(256**2, dtype=numpy.float32).reshape(256, 256) diff --git a/test/terse.h5 b/test/terse.h5 new file mode 100644 index 00000000..07322c8b Binary files /dev/null and b/test/terse.h5 differ diff --git a/test/test.py b/test/test.py index c6ddd582..5c34c67a 100644 --- a/test/test.py +++ b/test/test.py @@ -202,6 +202,47 @@ def testSperrV0_2MissingValue(self): self.assertTrue(numpy.isnan(compressed[0, 10])) self.assertTrue(numpy.isnan(compressed[1, 36])) + @unittest.skipUnless( + h5py.h5z.filter_avail(hdf5plugin.TERSE_ID), "Terse filter not available" + ) + def testTerse(self): + """Test reading TERSE/PROLIX compressed data (lossless) + + The file was written with h5py's low-level API and the raw filter + ``cd_values`` (i.e., the data type code) as in the upstream + ``py_tests/filter_test.py``, without using ``hdf5plugin.Terse``. + """ + dirname = os.path.abspath(os.path.dirname(__file__)) + fname = os.path.join(dirname, "terse.h5") + self.assertTrue(os.path.exists(fname), "Cannot find %s file" % fname) + expected = { + "uint16": numpy.arange(60 * 30, dtype=numpy.uint16).reshape(60, 30), + "int32": numpy.arange(-900, 900, dtype=numpy.int32).reshape(60, 30), + "frames_uint16": numpy.arange(4 * 32 * 32, dtype=numpy.uint16).reshape( + 4, 32, 32 + ) + % 1000, + } + with h5py.File(fname, "r") as h5: + for name, original in expected.items(): + with self.subTest(dataset=name): + dataset = h5[name] + compressed = dataset[()] + self.assertEqual(compressed.shape, original.shape) + self.assertEqual(compressed.dtype, original.dtype) + self.assertTrue( + numpy.array_equal(compressed, original), "Incorrect values" + ) + # Check filter configuration + plist = dataset.id.get_create_plist() + self.assertEqual(plist.get_nfilters(), 1) + filter_id, _, filter_options, _ = plist.get_filter(0) + self.assertEqual(filter_id, hdf5plugin.TERSE_ID) + self.assertEqual( + hdf5plugin.from_filter_options(filter_id, filter_options), + hdf5plugin.Terse(), + ) + @unittest.skipUnless( h5py.h5z.filter_avail(hdf5plugin.SZ_ID), "SZ filter not available" )