From 2882416b599514f5a7e60b07c6a20b53a32f9027 Mon Sep 17 00:00:00 2001 From: Daniel Spangberg Date: Wed, 15 May 2013 14:31:28 +0200 Subject: Added tng_compress trajectory compression algorithms diff --git a/include/CMakeLists.txt b/include/CMakeLists.txt new file mode 100644 index 0000000..74d02ac --- /dev/null +++ b/include/CMakeLists.txt @@ -0,0 +1 @@ +add_subdirectory(compression) diff --git a/include/compression/bwlzh.h b/include/compression/bwlzh.h new file mode 100644 index 0000000..3e58d2d --- /dev/null +++ b/include/compression/bwlzh.h @@ -0,0 +1,74 @@ +/* This code is part of the tng compression routines. + * + * Written by Daniel Spangberg + * Copyright (c) 2010, 2013, The GROMACS development team. + * + * + * This program is free software; you can redistribute it and/or + * modify it under the terms of the GNU Lesser General Public License + * as published by the Free Software Foundation; either version 2.1 + * of the License, or (at your option) any later version. + */ + + +#ifndef BWLZH_H +#define BWLZH_H + +/* Compress the integers (positive, small integers are preferable) + using bwlzh compression. The unsigned char *output should be + allocated to be able to hold worst case. You can obtain this length + conveniently by calling comp_get_buflen() +*/ +void DECLSPECDLLEXPORT bwlzh_compress(unsigned int *vals, int nvals, + unsigned char *output, int *output_len); + +void DECLSPECDLLEXPORT bwlzh_compress_no_lz77(unsigned int *vals, int nvals, + unsigned char *output, int *output_len); + +int DECLSPECDLLEXPORT bwlzh_get_buflen(int nvals); + +void DECLSPECDLLEXPORT bwlzh_decompress(unsigned char *input, int nvals, + unsigned int *vals); + + +/* The routines below are mostly useful for testing, and for internal + use by the library. */ + +void DECLSPECDLLEXPORT bwlzh_compress_verbose(unsigned int *vals, int nvals, + unsigned char *output, int *output_len); + +void DECLSPECDLLEXPORT bwlzh_compress_no_lz77_verbose(unsigned int *vals, int nvals, + unsigned char *output, int *output_len); + +void DECLSPECDLLEXPORT bwlzh_decompress_verbose(unsigned char *input, int nvals, + unsigned int *vals); + +/* Compress the integers (positive, small integers are preferable) + using huffman coding, with automatic selection of how to handle the + huffman dictionary. The unsigned char *huffman should be allocated + to be able to hold worst case. You can obtain this length + conveniently by calling comp_huff_buflen() +*/ +void Ptngc_comp_huff_compress(unsigned int *vals, int nvals, + unsigned char *huffman, int *huffman_len); + +int Ptngc_comp_huff_buflen(int nvals); + +void Ptngc_comp_huff_decompress(unsigned char *huffman, int huffman_len, + unsigned int *vals); + + +/* the value pointed to by chosen_algo should be + sent as -1 for autodetect. */ +void Ptngc_comp_huff_compress_verbose(unsigned int *vals, int nvals, + unsigned char *huffman, int *huffman_len, + int *huffdatalen, + int *huffman_lengths,int *chosen_algo, + int isvals16); + +#define N_HUFFMAN_ALGO 3 +char *Ptngc_comp_get_huff_algo_name(int algo); +char *Ptngc_comp_get_algo_name(int algo); + + +#endif diff --git a/include/compression/bwt.h b/include/compression/bwt.h new file mode 100644 index 0000000..e95f7ee --- /dev/null +++ b/include/compression/bwt.h @@ -0,0 +1,28 @@ +/* This code is part of the tng compression routines. + * + * Written by Daniel Spangberg + * Copyright (c) 2010, 2013, The GROMACS development team. + * + * + * This program is free software; you can redistribute it and/or + * modify it under the terms of the GNU Lesser General Public License + * as published by the Free Software Foundation; either version 2.1 + * of the License, or (at your option) any later version. + */ + + +#ifndef BWT_H +#define BWT_H + +void Ptngc_comp_to_bwt(unsigned int *vals, int nvals, + unsigned int *output, int *index); + +void Ptngc_comp_from_bwt(unsigned int *input, int nvals, int index, + unsigned int *vals); + +void Ptngc_bwt_merge_sort_inner(int *indices, int nvals,unsigned int *vals, + int start, int end, + unsigned int *nrepeat, + int *workarray); + +#endif diff --git a/include/compression/coder.h b/include/compression/coder.h new file mode 100644 index 0000000..5cef38a --- /dev/null +++ b/include/compression/coder.h @@ -0,0 +1,42 @@ +/* This code is part of the tng compression routines. + * + * Written by Daniel Spangberg + * Copyright (c) 2010, 2013, The GROMACS development team. + * + * + * This program is free software; you can redistribute it and/or + * modify it under the terms of the GNU Lesser General Public License + * as published by the Free Software Foundation; either version 2.1 + * of the License, or (at your option) any later version. + */ + +#ifndef CODER_H +#define CODER_H + +struct coder +{ + unsigned int pack_temporary; + int pack_temporary_bits; + int stat_overflow; + int stat_numval; +}; + +struct coder *Ptngc_coder_init(void); +void Ptngc_coder_deinit(struct coder *coder); +unsigned char *Ptngc_pack_array(struct coder *coder,int *input, int *length, int coding, int coding_parameter, int natoms, int speed); +int Ptngc_unpack_array(struct coder *coder,unsigned char *packed,int *output, int length, int coding, int coding_parameter, int natoms); +unsigned char *Ptngc_pack_array_xtc2(struct coder *coder,int *input, int *length); +int Ptngc_unpack_array_xtc2(struct coder *coder,unsigned char *packed,int *output, int length); +unsigned char *Ptngc_pack_array_xtc3(int *input, int *length, int natoms, int speed); +int Ptngc_unpack_array_xtc3(unsigned char *packed,int *output, int length, int natoms); + +void Ptngc_out8bits(struct coder *coder, unsigned char **output); +void Ptngc_pack_flush(struct coder *coder,unsigned char **output); +void Ptngc_write_pattern(struct coder *coder,unsigned int pattern, int nbits, unsigned char **output); + +void Ptngc_writebits(struct coder *coder,unsigned int value,int nbits, unsigned char **output_ptr); +void Ptngc_write32bits(struct coder *coder,unsigned int value,int nbits, unsigned char **output_ptr); +void Ptngc_writemanybits(struct coder *coder,unsigned char *value,int nbits, unsigned char **output_ptr); + + +#endif diff --git a/include/compression/dict.h b/include/compression/dict.h new file mode 100644 index 0000000..e29b17f --- /dev/null +++ b/include/compression/dict.h @@ -0,0 +1,23 @@ +/* This code is part of the tng compression routines. + * + * Written by Daniel Spangberg + * Copyright (c) 2010, 2013, The GROMACS development team. + * + * + * This program is free software; you can redistribute it and/or + * modify it under the terms of the GNU Lesser General Public License + * as published by the Free Software Foundation; either version 2.1 + * of the License, or (at your option) any later version. + */ + + +#ifndef DICT_H +#define DICT_H + +void Ptngc_comp_canonical_dict(unsigned int *dict, int *ndict); + +void Ptngc_comp_make_dict_hist(unsigned int *vals, int nvals, + unsigned int *dict, int *ndict, + unsigned int *hist); + +#endif diff --git a/include/compression/fixpoint.h b/include/compression/fixpoint.h new file mode 100644 index 0000000..d20d5e1 --- /dev/null +++ b/include/compression/fixpoint.h @@ -0,0 +1,45 @@ +/* This code is part of the tng compression routines. + * + * Written by Daniel Spangberg + * Copyright (c) 2010, 2013, The GROMACS development team. + * + * + * This program is free software; you can redistribute it and/or + * modify it under the terms of the GNU Lesser General Public License + * as published by the Free Software Foundation; either version 2.1 + * of the License, or (at your option) any later version. + */ + +#ifndef FIXPOINT_H +#define FIXPOINT_H + +#if HAVE_CONFIG_H +#include +#endif + +#include "my64bit.h" + +/* There are at least 32 bits available in a long. */ +typedef unsigned long fix_t; + +/* Positive double to 32 bit fixed point value */ +fix_t Ptngc_ud_to_fix_t(double d,double max); + +/* double to signed 32 bit fixed point value */ +fix_t Ptngc_d_to_fix_t(double d,double max); + +/* 32 bit fixed point value to positive double */ +double Ptngc_fix_t_to_ud(fix_t f, double max); + +/* signed 32 bit fixed point value to double */ +double Ptngc_fix_t_to_d(fix_t f, double max); + +/* Convert a floating point variable to two 32 bit integers with range + -2.1e9 to 2.1e9 and precision to somewhere around 1e-9. */ +void Ptngc_d_to_i32x2(double d, fix_t *hi, fix_t *lo); + +/* Convert two 32 bit integers to a floating point variable + -2.1e9 to 2.1e9 and precision to somewhere around 1e-9. */ +double Ptngc_i32x2_to_d(fix_t hi, fix_t lo); + +#endif diff --git a/include/compression/huffman.h b/include/compression/huffman.h new file mode 100644 index 0000000..c36e94d --- /dev/null +++ b/include/compression/huffman.h @@ -0,0 +1,35 @@ +/* This code is part of the tng compression routines. + * + * Written by Daniel Spangberg + * Copyright (c) 2010, 2013, The GROMACS development team. + * + * + * This program is free software; you can redistribute it and/or + * modify it under the terms of the GNU Lesser General Public License + * as published by the Free Software Foundation; either version 2.1 + * of the License, or (at your option) any later version. + */ + + +#ifndef HUFFMAN_H +#define HUFFMAN_H + +void Ptngc_comp_conv_to_huffman(unsigned int *vals, int nvals, + unsigned int *dict, int ndict, + unsigned int *prob, + unsigned char *huffman, + int *huffman_len, + unsigned char *huffman_dict, + int *huffman_dictlen, + unsigned int *huffman_dict_unpacked, + int *huffman_dict_unpackedlen); + +void Ptngc_comp_conv_from_huffman(unsigned char *huffman, + unsigned int *vals, int nvals, + int ndict, + unsigned char *huffman_dict, + int huffman_dictlen, + unsigned int *huffman_dict_unpacked, + int huffman_dict_unpackedlen); + +#endif diff --git a/include/compression/lz77.h b/include/compression/lz77.h new file mode 100644 index 0000000..e811256 --- /dev/null +++ b/include/compression/lz77.h @@ -0,0 +1,27 @@ +/* This code is part of the tng compression routines. + * + * Written by Daniel Spangberg + * Copyright (c) 2010, 2013, The GROMACS development team. + * + * + * This program is free software; you can redistribute it and/or + * modify it under the terms of the GNU Lesser General Public License + * as published by the Free Software Foundation; either version 2.1 + * of the License, or (at your option) any later version. + */ + + +#ifndef LZ77_H +#define LZ77_H + +void Ptngc_comp_to_lz77(unsigned int *vals, int nvals, + unsigned int *data, int *ndata, + unsigned int *len, int *nlens, + unsigned int *offsets, int *noffsets); + +void Ptngc_comp_from_lz77(unsigned int *data, int ndata, + unsigned int *len, int nlens, + unsigned int *offsets, int noffsets, + unsigned int *vals, int nvals); + +#endif diff --git a/include/compression/merge_sort.h b/include/compression/merge_sort.h new file mode 100644 index 0000000..970d5ee --- /dev/null +++ b/include/compression/merge_sort.h @@ -0,0 +1,22 @@ +/* This code is part of the tng compression routines. + * + * Written by Daniel Spangberg + * Copyright (c) 2010, 2013, The GROMACS development team. + * + * + * This program is free software; you can redistribute it and/or + * modify it under the terms of the GNU Lesser General Public License + * as published by the Free Software Foundation; either version 2.1 + * of the License, or (at your option) any later version. + */ + + +#ifndef MERGE_SORT_H +#define MERGE_SORT_H + +void Ptngc_merge_sort(void *base, size_t nmemb, size_t size, + int (*compar)(const void *v1,const void *v2,const void *private), + void *private); + + +#endif diff --git a/include/compression/mtf.h b/include/compression/mtf.h new file mode 100644 index 0000000..9c5e175 --- /dev/null +++ b/include/compression/mtf.h @@ -0,0 +1,37 @@ +/* This code is part of the tng compression routines. + * + * Written by Daniel Spangberg + * Copyright (c) 2010, 2013, The GROMACS development team. + * + * + * This program is free software; you can redistribute it and/or + * modify it under the terms of the GNU Lesser General Public License + * as published by the Free Software Foundation; either version 2.1 + * of the License, or (at your option) any later version. + */ + + +#ifndef MTF_H +#define MTF_H + +void Ptngc_comp_conv_to_mtf(unsigned int *vals, int nvals, + unsigned int *dict, int ndict, + unsigned int *valsmtf); + +void Ptngc_comp_conv_from_mtf(unsigned int *valsmtf, int nvals, + unsigned int *dict, int ndict, + unsigned int *vals); + +void Ptngc_comp_conv_to_mtf_partial(unsigned int *vals, int nvals, + unsigned int *valsmtf); + +void Ptngc_comp_conv_from_mtf_partial(unsigned int *valsmtf, int nvals, + unsigned int *vals); + +void Ptngc_comp_conv_to_mtf_partial3(unsigned int *vals, int nvals, + unsigned char *valsmtf); + +void Ptngc_comp_conv_from_mtf_partial3(unsigned char *valsmtf, int nvals, + unsigned int *vals); + +#endif diff --git a/include/compression/my64bit.h b/include/compression/my64bit.h new file mode 100644 index 0000000..fb2d5d0 --- /dev/null +++ b/include/compression/my64bit.h @@ -0,0 +1,67 @@ +/* This code is part of the tng compression routines. + * + * Written by Daniel Spangberg + * Copyright (c) 2010, 2013, The GROMACS development team. + * + * + * This program is free software; you can redistribute it and/or + * modify it under the terms of the GNU Lesser General Public License + * as published by the Free Software Foundation; either version 2.1 + * of the License, or (at your option) any later version. + */ + + +/* Here presence of 64 bit integers are checked for. If on windows + define USE_WINDOWS if it does not get defined automatically. + If on UNIX, define the sizes of SIZEOF_INT, SIZEOF_LONG, SIZEOF_LONG_LONG + If none of these symbols are defined, stdint.h is included (if on C99 + or USE_STDINT is defined) and 64 bit integers are assumed to be + present. If none of this is fulfilled, 64 bit integers are not + available. */ + +#ifndef MY64BIT_H +#define MY64BIT_H + +#if HAVE_CONFIG_H +#include +#endif + +#if HAVE_STDINT_H +#include +#endif + +/* The USE_WINDOWS symbol should be automatically defined in tng_compress.h */ +#include "tng_compress.h" + +#ifdef USE_WINDOWS +typedef __int64 my_int64_t; +typedef unsigned __int64 my_uint64_t; +#define HAVE64BIT +#else /* USE_WINDOWS */ +/* UNIX. Use config.h */ +#if SIZEOF_INT >= 8 +typedef int my_int64_t; +typedef unsigned int my_uint64_t; +#define HAVE64BIT +#else /* SIZEOF_INT */ +#if SIZEOF_LONG >= 8 +typedef long my_int64_t; +typedef unsigned long my_uint64_t; +#define HAVE64BIT +#else /* SIZEOF_LONG */ +#if SIZEOF_LONG_LONG >= 8 +typedef long long my_int64_t; +typedef unsigned long long my_uint64_t; +#define HAVE64BIT +#else /* SIZEOF_LONG_LONG */ +#if HAVE_STDINT_H +typedef int64_t my_int64_t; +typedef uint64_t my_uint64_t; +#define HAVE64BIT +#endif /* STDINT_H */ +#endif /* SIZEOF_LONG_LONG */ +#endif /* SIZEOF_LONG */ +#endif /* SIZEOF_INT */ +#endif /* USE_WINDOWS */ + +#endif diff --git a/include/compression/rle.h b/include/compression/rle.h new file mode 100644 index 0000000..3adf8dc --- /dev/null +++ b/include/compression/rle.h @@ -0,0 +1,24 @@ +/* This code is part of the tng compression routines. + * + * Written by Daniel Spangberg + * Copyright (c) 2010, 2013, The GROMACS development team. + * + * + * This program is free software; you can redistribute it and/or + * modify it under the terms of the GNU Lesser General Public License + * as published by the Free Software Foundation; either version 2.1 + * of the License, or (at your option) any later version. + */ + + +#ifndef RLE_H +#define RLE_H + +void Ptngc_comp_conv_to_rle(unsigned int *vals, int nvals, + unsigned int *rle, int *nrle, + int min_rle); + +void Ptngc_comp_conv_from_rle(unsigned int *rle, + unsigned int *vals, int nvals); + +#endif diff --git a/include/compression/tng_compress.h b/include/compression/tng_compress.h new file mode 100644 index 0000000..eefcfaa --- /dev/null +++ b/include/compression/tng_compress.h @@ -0,0 +1,165 @@ +/* This code is part of the tng compression routines. + * + * Written by Daniel Spangberg + * Copyright (c) 2010, 2013, The GROMACS development team. + * + * + * This program is free software; you can redistribute it and/or + * modify it under the terms of the GNU Lesser General Public License + * as published by the Free Software Foundation; either version 2.1 + * of the License, or (at your option) any later version. + */ + +#ifndef TNG_COMPRESS_H +#define TNG_COMPRESS_H + +#ifndef USE_WINDOWS +#if defined(WIN32) || defined(_WIN32) || defined(WIN64) || defined(_WIN64) +#define USE_WINDOWS +#endif /* win32... */ +#endif /* not defined USE_WINDOWS */ + +#ifdef USE_WINDOWS +#define DECLSPECDLLEXPORT __declspec(dllexport) +#else +#define DECLSPECDLLEXPORT +#endif + +#ifdef __cplusplus + extern "C" { +#endif + +/* tng_compress_pos expects positions to have the order: + first xyz, then sorted in atom order + then all the frames repeated, i.e.: + nframes * [ + natoms* [ + x, y, z + ] + ] + desired_precision what to round the numbers to, i.e. integers will be created as: + round(pos[]/desired_precision). + + algo should first be determined by calling + tng_compress_pos_find_algo + + The compressed data is returned in a malloced pointer (so free can + be called to free the memory), the number of chars in the compressed + data is put into *nitems. +*/ + +char DECLSPECDLLEXPORT *tng_compress_pos(double *pos, int natoms, int nframes, + double desired_precision, + int speed, int *algo, + int *nitems); + + +/* The tng_compress_pos_find_algo works the same as tng_compress_pos, but + it performs benchmarking to find the algorithms with the best + compression ratio. + The search is controlled by giving speed: + speed=1: Fast algorithms only. This excludes all BWLZH algorithms and + the XTC3 algorithm. + speed=2: Same as 1 and also includes the XTC3 algorithm using base compression + only. + speed=3: Same as 2 and also includes the XTC3 algorithm which will use BWLZH + compression when it seems likely to give better + compression. Also includes the interframe BWLZH algorithm for + coordinates and velocities. + speed=4: Enable the inter frame BWLZH algorithm for the coordinates. + The one-to-one BWLZH algorithm is enabled for velocities. + speed=5: Enable the LZ77 part of the BWLZH algorithm. + speed=6: Enable the intra frame BWLZH algorithm for the coordinates. Always try + the BWLZH compression in the XTC3 algorithm. + + Set speed=0 to allow tng_compression to set the default speed (which is currently 2). + For very good compression it makes sense to choose speed=4 or speed=5 + + The number of items required in the algorithm array can be found + by calling tng_compress_nalgo +*/ + +char DECLSPECDLLEXPORT *tng_compress_pos_find_algo(double *pos, int natoms, int nframes, + double desired_precision, + int speed, + int *algo, + int *nitems); + +/* This returns the number of integers required for the storage of the algorithm + with the best compression ratio. */ +int DECLSPECDLLEXPORT tng_compress_nalgo(void); + +/* The following two routines does the same as the compression of the + positions, but compresses velocities instead. The algorithm + selection for velocities is different, so the position and + velocities routines should not be mixed. */ + +char DECLSPECDLLEXPORT *tng_compress_vel(double *vel, int natoms, int nframes, + double desired_precision, + int speed, int *algo, + int *nitems); + +char DECLSPECDLLEXPORT *tng_compress_vel_find_algo(double *vel, int natoms, int nframes, + double desired_precision, + int speed, + int *algo, + int *nitems); + +/* From a compressed block, obtain information about + whether it is a position or velocity block: + *vel=1 means velocity block, *vel=0 means position block. + It also gives info about the number of atoms, + frames, and the precision used to compress the block, and the algorithms used to + compress the block. The return value=0 if the block looks like a tng compressed block, + and 1 otherwise. If the return value is 1 no information is returned. */ +int DECLSPECDLLEXPORT tng_compress_inquire(char *data,int *vel, int *natoms, + int *nframes, double *precision, + int *algo); + +/* Uncompresses any tng compress block, positions or velocities. It determines whether it is positions or velocities from the data buffer. The return value is 0 if ok, and 1 if not. +*/ +int DECLSPECDLLEXPORT tng_compress_uncompress(char *data,double *posvel); + + + /* Compression algorithms (matching the original trajng + assignments) The compression backends require that some of the + algorithms must have the same value. */ + +#define TNG_COMPRESS_ALGO_STOPBIT 1 +#define TNG_COMPRESS_ALGO_TRIPLET 2 +#define TNG_COMPRESS_ALGO_BWLZH1 8 +#define TNG_COMPRESS_ALGO_BWLZH2 9 + +#define TNG_COMPRESS_ALGO_POS_STOPBIT_INTER TNG_COMPRESS_ALGO_STOPBIT +#define TNG_COMPRESS_ALGO_POS_TRIPLET_INTER TNG_COMPRESS_ALGO_TRIPLET +#define TNG_COMPRESS_ALGO_POS_TRIPLET_INTRA 3 +#define TNG_COMPRESS_ALGO_POS_XTC2 5 +#define TNG_COMPRESS_ALGO_POS_TRIPLET_ONETOONE 7 +#define TNG_COMPRESS_ALGO_POS_BWLZH_INTER TNG_COMPRESS_ALGO_BWLZH1 +#define TNG_COMPRESS_ALGO_POS_BWLZH_INTRA TNG_COMPRESS_ALGO_BWLZH2 +#define TNG_COMPRESS_ALGO_POS_XTC3 10 +#define TNG_COMPRESS_ALGO_VEL_STOPBIT_ONETOONE TNG_COMPRESS_ALGO_STOPBIT +#define TNG_COMPRESS_ALGO_VEL_TRIPLET_INTER TNG_COMPRESS_ALGO_TRIPLET +#define TNG_COMPRESS_ALGO_VEL_TRIPLET_ONETOONE 3 +#define TNG_COMPRESS_ALGO_VEL_STOPBIT_INTER 6 +#define TNG_COMPRESS_ALGO_VEL_BWLZH_INTER TNG_COMPRESS_ALGO_BWLZH1 +#define TNG_COMPRESS_ALGO_VEL_BWLZH_ONETOONE TNG_COMPRESS_ALGO_BWLZH2 +#define TNG_COMPRESS_ALGO_MAX 11 + + + +/* Obtain strings describing the actual algorithms. These point to static memory, so should + not be freed. */ +char DECLSPECDLLEXPORT *tng_compress_initial_pos_algo(int *algo); +char DECLSPECDLLEXPORT *tng_compress_pos_algo(int *algo); +char DECLSPECDLLEXPORT *tng_compress_initial_vel_algo(int *algo); +char DECLSPECDLLEXPORT *tng_compress_vel_algo(int *algo); + + + +#ifdef __cplusplus + } +#endif + + +#endif diff --git a/include/compression/vals16.h b/include/compression/vals16.h new file mode 100644 index 0000000..a8acdf4 --- /dev/null +++ b/include/compression/vals16.h @@ -0,0 +1,23 @@ +/* This code is part of the tng compression routines. + * + * Written by Daniel Spangberg + * Copyright (c) 2010, 2013, The GROMACS development team. + * + * + * This program is free software; you can redistribute it and/or + * modify it under the terms of the GNU Lesser General Public License + * as published by the Free Software Foundation; either version 2.1 + * of the License, or (at your option) any later version. + */ + + +#ifndef VALS16_H +#define VALS16_H + +void Ptngc_comp_conv_to_vals16(unsigned int *vals,int nvals, + unsigned int *vals16, int *nvals16); + +void Ptngc_comp_conv_from_vals16(unsigned int *vals16,int nvals16, + unsigned int *vals, int *nvals); + +#endif diff --git a/include/compression/warnmalloc.h b/include/compression/warnmalloc.h new file mode 100644 index 0000000..c2eb28d --- /dev/null +++ b/include/compression/warnmalloc.h @@ -0,0 +1,32 @@ +/* This code is part of the tng compression routines. + * + * Written by Daniel Spangberg + * Copyright (c) 2010, 2013, The GROMACS development team. + * + * + * This program is free software; you can redistribute it and/or + * modify it under the terms of the GNU Lesser General Public License + * as published by the Free Software Foundation; either version 2.1 + * of the License, or (at your option) any later version. + */ + + +#ifndef WARNMALLOC_H +#define WARNMALLOC_H + +#if HAVE_CONFIG_H +#include +#endif + +#include "tng_compress.h" + +void DECLSPECDLLEXPORT *Ptngc_warnmalloc_x(size_t size, char *file, int line); + +#define warnmalloc(size) Ptngc_warnmalloc_x(size,__FILE__,__LINE__) + +void DECLSPECDLLEXPORT *Ptngc_warnrealloc_x(void *old, size_t size, char *file, int line); + +#define warnrealloc(old,size) Ptngc_warnrealloc_x(old,size,__FILE__,__LINE__) + + +#endif diff --git a/include/compression/widemuldiv.h b/include/compression/widemuldiv.h new file mode 100644 index 0000000..b7574fa --- /dev/null +++ b/include/compression/widemuldiv.h @@ -0,0 +1,33 @@ +/* This code is part of the tng compression routines. + * + * Written by Daniel Spangberg + * Copyright (c) 2010, 2013, The GROMACS development team. + * + * + * This program is free software; you can redistribute it and/or + * modify it under the terms of the GNU Lesser General Public License + * as published by the Free Software Foundation; either version 2.1 + * of the License, or (at your option) any later version. + */ + + +#ifndef WIDEMULDIV_H +#define WIDEMULDIV_H + +/* Multiply two 32 bit unsigned integers returning a 64 bit unsigned value (in two integers) */ +void Ptngc_widemul(unsigned int i1, unsigned int i2, unsigned int *ohi, unsigned int *olo); + +/* Divide a 64 bit unsigned value in hi:lo with the 32 bit value i and + return the result in result and the remainder in remainder */ +void Ptngc_widediv(unsigned int hi, unsigned int lo, unsigned int i, unsigned int *result, unsigned int *remainder); + +/* Add a unsigned int to a largeint. */ +void Ptngc_largeint_add(unsigned int v1, unsigned int *largeint, int n); + +/* Multiply v1 with largeint_in and return result in largeint_out */ +void Ptngc_largeint_mul(unsigned int v1, unsigned int *largeint_in, unsigned int *largeint_out, int n); + +/* Return the remainder from dividing largeint_in with v1. Result of the division is returned in largeint_out */ +unsigned int Ptngc_largeint_div(unsigned int v1, unsigned int *largeint_in, unsigned int *largeint_out, int n); + +#endif diff --git a/src/compression/CMakeLists.txt b/src/compression/CMakeLists.txt new file mode 100644 index 0000000..3cd090b --- /dev/null +++ b/src/compression/CMakeLists.txt @@ -0,0 +1,14 @@ +add_library(tng_io SHARED tng_io.c md5.c) + +if(HAVE_INTTYPES_H) + set_property(TARGET tng_io APPEND PROPERTY COMPILE_DEFINITIONS USE_STD_INTTYPES_H) +endif() + +if(BUILD_FORTRAN) + set_property(TARGET tng_io APPEND PROPERTY COMPILE_DEFINITIONS BUILD_FORTRAN) +endif() + +if(ZLIB_FOUND) + set_property(TARGET tng_io APPEND PROPERTY COMPILE_DEFINITIONS USE_ZLIB) + target_link_libraries(tng_io ${ZLIB_LIBRARIES}) +endif() diff --git a/src/compression/bwlzh.c b/src/compression/bwlzh.c new file mode 100644 index 0000000..b5590e3 --- /dev/null +++ b/src/compression/bwlzh.c @@ -0,0 +1,805 @@ +/* This code is part of the tng compression routines. + * + * Written by Daniel Spangberg + * Copyright (c) 2010, 2013, The GROMACS development team. + * + * + * This program is free software; you can redistribute it and/or + * modify it under the terms of the GNU Lesser General Public License + * as published by the Free Software Foundation; either version 2.1 + * of the License, or (at your option) any later version. + */ + + +#if HAVE_CONFIG_H +#include +#endif + +#include +#include +#include +#include "warnmalloc.h" +#include "tng_compress.h" +#include "bwlzh.h" +#include "dict.h" +#include "vals16.h" +#include "rle.h" +#include "mtf.h" +#include "bwt.h" +#include "lz77.h" + +#if 0 +#define SHOWIT +#endif + +#if 0 +#define SHOWTEST +#endif + +#define MAX_VALS_PER_BLOCK 200000 + +#if 1 +#define PARTIAL_MTF3 +#endif + +#if 0 +#define PARTIAL_MTF +#endif + +int bwlzh_get_buflen(int nvals) +{ + return 132000+nvals*8+12*((nvals+MAX_VALS_PER_BLOCK)/MAX_VALS_PER_BLOCK); +} + +#ifdef SHOWIT +static void printvals(char *name, unsigned int *vals, int nvals) +{ + int i; + int nvalsmax=nvals; + if (nvalsmax>99) + nvalsmax=99; +#if 0 + for (i=0; i>8)&0xFFU; + output[outdata++]=(((unsigned int)nvals)>>16)&0xFFU; + output[outdata++]=(((unsigned int)nvals)>>24)&0xFFU; + + valsleft=nvals; + valstart=0; + while (valsleft) + { + int reducealgo=1; /* Reduce algo is LZ77. */ + if (!enable_lz77) + reducealgo=0; + thisvals=valsleft; + if (thisvals>max_vals_per_block) + thisvals=max_vals_per_block; + valsleft-=thisvals; + if (verbose) + fprintf(stderr,"Creating vals16 block from %d values.\n",thisvals); + +#ifdef SHOWIT + printvals("vals",vals+valstart,thisvals); +#endif + + Ptngc_comp_conv_to_vals16(vals+valstart,thisvals,vals16,&nvals16); + valstart+=thisvals; + +#ifdef SHOWTEST + nvals16=99; +#endif + +#ifdef SHOWIT + printvals("vals16",vals16,nvals16); +#endif + + if (verbose) + { + fprintf(stderr,"Resulting vals16 values: %d\n",nvals16); + } + if (verbose) + { + fprintf(stderr,"BWT\n"); + } + Ptngc_comp_to_bwt(vals16,nvals16,bwt,&bwt_index); + +#ifdef SHOWIT + printvals("bwt",bwt,nvals16); + fprintf(stderr,"BWT INDEX is %d\n",bwt_index); +#endif + + /* Store the number of real values in this block. */ + output[outdata++]=((unsigned int)thisvals)&0xFFU; + output[outdata++]=(((unsigned int)thisvals)>>8)&0xFFU; + output[outdata++]=(((unsigned int)thisvals)>>16)&0xFFU; + output[outdata++]=(((unsigned int)thisvals)>>24)&0xFFU; + + /* Store the number of nvals16 values in this block. */ + output[outdata++]=((unsigned int)nvals16)&0xFFU; + output[outdata++]=(((unsigned int)nvals16)>>8)&0xFFU; + output[outdata++]=(((unsigned int)nvals16)>>16)&0xFFU; + output[outdata++]=(((unsigned int)nvals16)>>24)&0xFFU; + + /* Store the BWT index. */ + output[outdata++]=((unsigned int)bwt_index)&0xFFU; + output[outdata++]=(((unsigned int)bwt_index)>>8)&0xFFU; + output[outdata++]=(((unsigned int)bwt_index)>>16)&0xFFU; + output[outdata++]=(((unsigned int)bwt_index)>>24)&0xFFU; + + if (verbose) + fprintf(stderr,"MTF\n"); +#ifdef PARTIAL_MTF3 + Ptngc_comp_conv_to_mtf_partial3(bwt,nvals16, + mtf3); + for (imtfinner=0; imtfinner<3; imtfinner++) + { + int i; + if (verbose) + fprintf(stderr,"Doing partial MTF: %d\n",imtfinner); + for (i=0; i>=1; + } + coarse[numbits-1]+=thist[jj]; + } +#if 1 + for (jj=0; jj>8)&0xFFU; + output[outdata++]=(((unsigned int)nrle)>>16)&0xFFU; + output[outdata++]=(((unsigned int)nrle)>>24)&0xFFU; + + /* Store the size of the huffman block. */ + output[outdata++]=((unsigned int)bwlzhhufflen)&0xFFU; + output[outdata++]=(((unsigned int)bwlzhhufflen)>>8)&0xFFU; + output[outdata++]=(((unsigned int)bwlzhhufflen)>>16)&0xFFU; + output[outdata++]=(((unsigned int)bwlzhhufflen)>>24)&0xFFU; + + /* Store the huffman block. */ + memcpy(output+outdata,bwlzhhuff,bwlzhhufflen); + outdata+=bwlzhhufflen; + + if (reducealgo==1) + { + /* Store the number of values in this block. */ + output[outdata++]=((unsigned int)noffsets)&0xFFU; + output[outdata++]=(((unsigned int)noffsets)>>8)&0xFFU; + output[outdata++]=(((unsigned int)noffsets)>>16)&0xFFU; + output[outdata++]=(((unsigned int)noffsets)>>24)&0xFFU; + + if (noffsets>0) + { + if (verbose) + fprintf(stderr,"Huffman for offsets\n"); + + huffalgo=-1; + Ptngc_comp_huff_compress_verbose(offsets,noffsets,bwlzhhuff,&bwlzhhufflen,&huffdatalen,nhufflen,&huffalgo,1); + if (verbose) + { + int i; + fprintf(stderr,"Huffman data length is %d B.\n",huffdatalen); + for (i=0; i>8)&0xFFU; + output[outdata++]=(((unsigned int)bwlzhhufflen)>>16)&0xFFU; + output[outdata++]=(((unsigned int)bwlzhhufflen)>>24)&0xFFU; + + /* Store the huffman block. */ + memcpy(output+outdata,bwlzhhuff,bwlzhhufflen); + outdata+=bwlzhhufflen; + } + else + { + int i; + output[outdata++]=1; + for (i=0; i>8)&0xFFU; + } + if (verbose) + fprintf(stderr,"Store raw offsets: %d B\n",noffsets*2); + } + } + +#if 0 + { + int i,ndict; + FILE *f=fopen("len.dict","w"); + Ptngc_comp_make_dict_hist(lens,nlens,dict,&ndict,hist); + for (i=0; i>8)&0xFFU; + output[outdata++]=(((unsigned int)nlens)>>16)&0xFFU; + output[outdata++]=(((unsigned int)nlens)>>24)&0xFFU; + + /* Store the size of the huffman block. */ + output[outdata++]=((unsigned int)bwlzhhufflen)&0xFFU; + output[outdata++]=(((unsigned int)bwlzhhufflen)>>8)&0xFFU; + output[outdata++]=(((unsigned int)bwlzhhufflen)>>16)&0xFFU; + output[outdata++]=(((unsigned int)bwlzhhufflen)>>24)&0xFFU; + + /* Store the huffman block. */ + memcpy(output+outdata,bwlzhhuff,bwlzhhufflen); + outdata+=bwlzhhufflen; + } +#ifdef PARTIAL_MTF3 + } +#endif + } + + *output_len=outdata; + free(hist); + free(dict); + free(bwlzhhuff); +#ifdef PARTIAL_MTF3 + free(mtf3); +#endif + free(tmpmem); +} + + +void DECLSPECDLLEXPORT bwlzh_compress(unsigned int *vals, int nvals, + unsigned char *output, int *output_len) +{ + bwlzh_compress_gen(vals,nvals,output,output_len,1,0); +} + +void DECLSPECDLLEXPORT bwlzh_compress_verbose(unsigned int *vals, int nvals, + unsigned char *output, int *output_len) +{ + bwlzh_compress_gen(vals,nvals,output,output_len,1,1); +} + + +void DECLSPECDLLEXPORT bwlzh_compress_no_lz77(unsigned int *vals, int nvals, + unsigned char *output, int *output_len) +{ + bwlzh_compress_gen(vals,nvals,output,output_len,0,0); +} + +void DECLSPECDLLEXPORT bwlzh_compress_no_lz77_verbose(unsigned int *vals, int nvals, + unsigned char *output, int *output_len) +{ + bwlzh_compress_gen(vals,nvals,output,output_len,0,1); +} + + +static void bwlzh_decompress_gen(unsigned char *input, int nvals, + unsigned int *vals, + int verbose) +{ + unsigned int *vals16; + int nvals16; + int ndict; + int bwt_index; + unsigned int *bwt=NULL; + unsigned int *mtf=NULL; +#ifdef PARTIAL_MTF3 + unsigned char *mtf3=NULL; + int imtfinner; +#endif + unsigned int *rle=NULL; + unsigned int *offsets=NULL; + unsigned int *lens=NULL; + unsigned int *dict=warnmalloc(0x20004*sizeof *dict); + unsigned int *hist=warnmalloc(0x20004*sizeof *hist); + int nrle, noffsets, nlens; + unsigned char *bwlzhhuff=NULL; + int bwlzhhufflen; + int max_vals_per_block=MAX_VALS_PER_BLOCK; + int valsleft; + int thisvals; + int valstart; + int inpdata=0; + int nvalsfile; + + unsigned int *tmpmem=warnmalloc(max_vals_per_block*18*sizeof *tmpmem); + +#if 0 + verbose=1; +#endif + + + bwlzhhuff=warnmalloc(Ptngc_comp_huff_buflen(3*nvals)); + vals16=tmpmem; + bwt=tmpmem+max_vals_per_block*3; + mtf=tmpmem+max_vals_per_block*6; + rle=tmpmem+max_vals_per_block*9; + offsets=tmpmem+max_vals_per_block*12; + lens=tmpmem+max_vals_per_block*15; +#ifdef PARTIAL_MTF3 + mtf3=warnmalloc(max_vals_per_block*3*3*sizeof *mtf3); /* 3 due to expansion of 32 bit to 16 bit, 3 due to up to 3 bytes + per 16 value. */ +#endif + + if (verbose) + { + fprintf(stderr,"Number of input values: %d\n",nvals); + } + + /* Read the number of real values in the whole block. */ + nvalsfile=(int)(((unsigned int)input[inpdata]) | + (((unsigned int)input[inpdata+1])<<8) | + (((unsigned int)input[inpdata+2])<<16) | + (((unsigned int)input[inpdata+3])<<24)); + inpdata+=4; + + if (nvalsfile!=nvals) + { + fprintf(stderr,"BWLZH: The number of values found in the file is different from the number of values expected.\n"); + exit(EXIT_FAILURE); + } + + valsleft=nvals; + valstart=0; + while (valsleft) + { + int valsnew; + int reducealgo; + /* Read the number of real values in this block. */ + thisvals=(int)(((unsigned int)input[inpdata]) | + (((unsigned int)input[inpdata+1])<<8) | + (((unsigned int)input[inpdata+2])<<16) | + (((unsigned int)input[inpdata+3])<<24)); + inpdata+=4; + + valsleft-=thisvals; + + /* Read the number of nvals16 values in this block. */ + nvals16=(int)(((unsigned int)input[inpdata]) | + (((unsigned int)input[inpdata+1])<<8) | + (((unsigned int)input[inpdata+2])<<16) | + (((unsigned int)input[inpdata+3])<<24)); + inpdata+=4; + + /* Read the BWT index. */ + bwt_index=(int)(((unsigned int)input[inpdata]) | + (((unsigned int)input[inpdata+1])<<8) | + (((unsigned int)input[inpdata+2])<<16) | + (((unsigned int)input[inpdata+3])<<24)); + inpdata+=4; + + if (thisvals>max_vals_per_block) + { + /* More memory must be allocated for decompression. */ + max_vals_per_block=thisvals; + if (verbose) + fprintf(stderr,"Allocating more memory: %d B\n",(int)(max_vals_per_block*15*sizeof *tmpmem)); + tmpmem=warnrealloc(tmpmem,max_vals_per_block*18*sizeof *tmpmem); + vals16=tmpmem; + bwt=tmpmem+max_vals_per_block*3; + mtf=tmpmem+max_vals_per_block*6; + rle=tmpmem+max_vals_per_block*9; + offsets=tmpmem+max_vals_per_block*12; + lens=tmpmem+max_vals_per_block*15; +#ifdef PARTIAL_MTF3 + mtf3=warnrealloc(mtf3,max_vals_per_block*3*3*sizeof *mtf3); /* 3 due to expansion of 32 bit to 16 bit, 3 due to up to 3 bytes + per 16 value. */ +#endif + } + +#ifdef PARTIAL_MTF3 + for (imtfinner=0; imtfinner<3; imtfinner++) + { + int i; + if (verbose) + fprintf(stderr,"Doing partial MTF: %d\n",imtfinner); +#endif + + reducealgo=(int)input[inpdata]; + inpdata++; + + /* Read the number of huffman values in this block. */ + nrle=(int)(((unsigned int)input[inpdata]) | + (((unsigned int)input[inpdata+1])<<8) | + (((unsigned int)input[inpdata+2])<<16) | + (((unsigned int)input[inpdata+3])<<24)); + inpdata+=4; + + /* Read the size of the huffman block. */ + bwlzhhufflen=(int)(((unsigned int)input[inpdata]) | + (((unsigned int)input[inpdata+1])<<8) | + (((unsigned int)input[inpdata+2])<<16) | + (((unsigned int)input[inpdata+3])<<24)); + inpdata+=4; + + if (verbose) + fprintf(stderr,"Decompressing huffman block of length %d.\n",bwlzhhufflen); + /* Decompress the huffman block. */ + Ptngc_comp_huff_decompress(input+inpdata,bwlzhhufflen,rle); + inpdata+=bwlzhhufflen; + + if (reducealgo==1) /* LZ77 */ + { + int offstore; + /* Read the number of huffman values in this block. */ + noffsets=(int)(((unsigned int)input[inpdata]) | + (((unsigned int)input[inpdata+1])<<8) | + (((unsigned int)input[inpdata+2])<<16) | + (((unsigned int)input[inpdata+3])<<24)); + inpdata+=4; + + if (noffsets>0) + { + /* How are the offsets stored? */ + offstore=(int)input[inpdata++]; + if (offstore==0) + { + /* Read the size of the huffman block. */ + bwlzhhufflen=(int)(((unsigned int)input[inpdata]) | + (((unsigned int)input[inpdata+1])<<8) | + (((unsigned int)input[inpdata+2])<<16) | + (((unsigned int)input[inpdata+3])<<24)); + inpdata+=4; + + if (verbose) + fprintf(stderr,"Decompressing offset huffman block.\n"); + + /* Decompress the huffman block. */ + Ptngc_comp_huff_decompress(input+inpdata,bwlzhhufflen,offsets); + inpdata+=bwlzhhufflen; + } + else + { + int i; + if (verbose) + fprintf(stderr,"Reading offset block.\n"); + for (i=0; i +#include +#include +#include "warnmalloc.h" +#include "bwt.h" + +#if 0 +#define SHOWIT +#endif +#if 0 +#define SHOWIT2 +#endif + +static int compare_index(int i1,int i2,int nvals,unsigned int *vals,unsigned int *nrepeat) +{ + int i,j; + for (i=0; i>8); + int k1=(int)(nrepeat[i1]&0xFFU); + int repeat2=(int)(nrepeat[i2]>>8); + int k2=(int)(nrepeat[i2]&0xFFU); + + if ((repeat1>1) && (repeat1>1) && (k1==k2)) + { + int maxskip=0; + /* Yes. Compare the repeating patterns. */ + for (j=0; jv2) + return 1; + } + /* The repeating patterns are equal. Skip as far as we can + before continuing. */ + maxskip=repeat1; + if (repeat2vals[i2]) + return 1; + i1++; + if (i1>=nvals) + i1=0; + i2++; + if (i2>=nvals) + i2=0; + } + } + return 0; +} + +void Ptngc_bwt_merge_sort_inner(int *indices, int nvals,unsigned int *vals, + int start, int end, + unsigned int *nrepeat, + int *workarray) +{ + int middle; + if ((end-start)>1) + { + middle=start+(end-start)/2; +#if 0 + printf("For start %d end %d obtained new middle: %d\n",start,end,middle); +#endif + Ptngc_bwt_merge_sort_inner(indices,nvals,vals, + start,middle, + nrepeat, + workarray); + Ptngc_bwt_merge_sort_inner(indices,nvals,vals, + middle,end, + nrepeat, + workarray); +#if 0 + printf("For start %d end %d Before merge: Comparing element %d with %d\n",start,end,middle-1,middle); +#endif + if (compare_index(indices[middle-1],indices[middle],nvals,vals,nrepeat)>0) + { + /* Merge to work array. */ + int i, n=end-start; + int ileft=start; + int iright=middle; + for (i=0; i0) + { + workarray[i]=indices[iright]; + iright++; + } + else + { + workarray[i]=indices[ileft]; + ileft++; + } + } + } + /* Copy result back. */ + memcpy(indices+start,workarray,(end-start)*sizeof(int)); + } + } +} + +/* Burrows-Wheeler transform. */ +void Ptngc_comp_to_bwt(unsigned int *vals, int nvals, + unsigned int *output, int *index) +{ + int i; + int *indices=warnmalloc(2*nvals*sizeof *indices); + unsigned int *nrepeat=warnmalloc(nvals*sizeof *nrepeat); + int *warr=indices+nvals; + + if (nvals>0xFFFFFF) + { + fprintf(stderr,"BWT cannot pack more than %d values.\n",0xFFFFFF); + exit(1); + } + + /* Also note that repeat pattern k (kmax) cannot be larger than 255. */ +#if 0 + printf("Size of arrays is %.2f M\n",4*nvals*sizeof *indices/1024./1024.); +#endif + for (i=0; i=1; k--) + { + try_next_k: + if (k>=1) + { +#ifdef SHOWIT + printf("Trying k=%d at i=%d\n",k,i); +#endif + for (j=k; jmaxrepeat) + new_j=j; + if ((new_j>good_j) || ((new_j==good_j) && (knvals) + repeat=nvals; + nrepeat[i+m]=((unsigned int) (good_k)) | (((unsigned int) (repeat))<<8); + } + /* If no repetition was found for this value signal that here. */ + if (!nrepeat[i]) + nrepeat[i+m]=257U; /* This is 1<<8 | 1 */ + } + } +#ifdef SHOWIT + for (i=0; i>8)!=1) + printf("String repeats at %d: %d %d\n",i,nrepeat[i]>>8,nrepeat[i]&0xFFU); +#endif + + /* Sort cyclic shift matrix. */ + Ptngc_bwt_merge_sort_inner(indices,nvals,vals,0,nvals,nrepeat,warr); + +#if 0 + /* Test that it really is sorted. */ + for (i=0; i=0; i--) + { + vals[i]=input[index]; + index=p[index]+c[input[index]]; + } + free(p); + free(c); +} + diff --git a/src/compression/coder.c b/src/compression/coder.c new file mode 100644 index 0000000..bed73ca --- /dev/null +++ b/src/compression/coder.c @@ -0,0 +1,481 @@ +/* This code is part of the tng compression routines. + * + * Written by Daniel Spangberg + * Copyright (c) 2010, 2013, The GROMACS development team. + * + * + * This program is free software; you can redistribute it and/or + * modify it under the terms of the GNU Lesser General Public License + * as published by the Free Software Foundation; either version 2.1 + * of the License, or (at your option) any later version. + */ + + +#if HAVE_CONFIG_H +#include +#endif + +#include +#include +#include "tng_compress.h" +#include "bwlzh.h" +#include "coder.h" +#include "warnmalloc.h" + +struct coder *Ptngc_coder_init(void) +{ + struct coder *coder=warnmalloc(sizeof *coder); + coder->pack_temporary_bits=0; + return coder; +} + +void Ptngc_coder_deinit(struct coder *coder) +{ + free(coder); +} + +void Ptngc_out8bits(struct coder *coder, unsigned char **output) +{ + while (coder->pack_temporary_bits>=8) + { + unsigned int mask=~(0xFFU<<(coder->pack_temporary_bits-8)); + unsigned char out=(unsigned char)(coder->pack_temporary>>(coder->pack_temporary_bits-8)); + **output=out; + (*output)++; + coder->pack_temporary_bits-=8; + coder->pack_temporary&=mask; + } +} + +void Ptngc_write_pattern(struct coder *coder,unsigned int pattern, int nbits, unsigned char **output) +{ + unsigned int mask1,mask2; + mask1=1; + mask2=1<<(nbits-1); + coder->pack_temporary<<=nbits; /* Make room for new data. */ + coder->pack_temporary_bits+=nbits; + while (nbits) + { + if (pattern & mask1) + coder->pack_temporary|=mask2; + nbits--; + mask1<<=1; + mask2>>=1; + } + Ptngc_out8bits(coder,output); +} + +/* Write up to 24 bits */ +void Ptngc_writebits(struct coder *coder,unsigned int value,int nbits, unsigned char **output_ptr) +{ + /* Make room for the bits. */ + coder->pack_temporary<<=nbits; + coder->pack_temporary_bits+=nbits; + coder->pack_temporary|=value; + Ptngc_out8bits(coder,output_ptr); +} + +/* Write up to 32 bits */ +void Ptngc_write32bits(struct coder *coder,unsigned int value,int nbits, unsigned char **output_ptr) +{ + unsigned int mask; + if (nbits>=8) + mask=0xFFU<<(nbits-8); + else + mask=0xFFU>>(8-nbits); + while (nbits>8) + { + /* Make room for the bits. */ + coder->pack_temporary<<=8; + coder->pack_temporary_bits+=8; + coder->pack_temporary|=(value&mask)>>(nbits-8); + Ptngc_out8bits(coder,output_ptr); + nbits-=8; + mask>>=8; + } + if (nbits) + Ptngc_writebits(coder,value&mask,nbits,output_ptr); +} + +/* Write "arbitrary" number of bits */ +void Ptngc_writemanybits(struct coder *coder,unsigned char *value,int nbits, unsigned char **output_ptr) +{ + int vptr=0; + while (nbits>=24) + { + unsigned int v=((((unsigned int)value[vptr])<<16)| + (((unsigned int)value[vptr+1])<<8)| + (((unsigned int)value[vptr+2]))); + Ptngc_writebits(coder,v,24,output_ptr); + vptr+=3; + nbits-=24; + } + while (nbits>=8) + { + Ptngc_writebits(coder,(unsigned int)value[vptr],8,output_ptr); + vptr++; + nbits-=8; + } + if (nbits) + { + Ptngc_writebits(coder,(unsigned int)value[vptr],nbits,output_ptr); + } +} + +static int write_stop_bit_code(struct coder *coder, unsigned int s,unsigned int coding_parameter, unsigned char **output) +{ + do { + unsigned int extract=~(0xffffffffU<>=coding_parameter; + if (s) + { + this|=1U; + coder->stat_overflow++; + } + coder->pack_temporary<<=(coding_parameter+1); + coder->pack_temporary_bits+=coding_parameter+1; + coder->pack_temporary|=this; + Ptngc_out8bits(coder,output); + if (s) + { + coding_parameter>>=1; + if (coding_parameter<1) + coding_parameter=1; + } + } while (s); + coder->stat_numval++; + return 0; +} + +static int pack_stopbits_item(struct coder *coder,int item, unsigned char **output, int coding_parameter) +{ + /* Find this symbol in table. */ + int s=0; + if (item>0) + s=1+(item-1)*2; + else if (item<0) + s=2+(-item-1)*2; + return write_stop_bit_code(coder,s,coding_parameter,output); + return 0; +} + +static int pack_triplet(struct coder *coder,unsigned int *s, unsigned char **output, int coding_parameter, + unsigned int max_base, int maxbits) +{ + /* Determine base for this triplet. */ + unsigned int min_base=1U<=this_base) + { + this_base*=2; + jbase++; + } + bits_per_value=coding_parameter+jbase; + if (jbase>=3) + { + if (this_base>max_base) + return 1; + this_base=max_base; + bits_per_value=maxbits; + jbase=3; + } + /* 2 bits selects the base */ + coder->pack_temporary<<=2; + coder->pack_temporary_bits+=2; + coder->pack_temporary|=jbase; + Ptngc_out8bits(coder,output); + for (i=0; i<3; i++) + Ptngc_write32bits(coder,s[i],bits_per_value,output); + return 0; +} + +void Ptngc_pack_flush(struct coder *coder,unsigned char **output) +{ + /* Zero-fill just enough. */ + if (coder->pack_temporary_bits>0) + Ptngc_write_pattern(coder,0,8-coder->pack_temporary_bits,output); +} + +unsigned char *Ptngc_pack_array(struct coder *coder,int *input, int *length, int coding, int coding_parameter,int natoms, int speed) +{ + if ((coding==TNG_COMPRESS_ALGO_BWLZH1) || (coding==TNG_COMPRESS_ALGO_BWLZH2)) + { + unsigned char *output=warnmalloc(4+bwlzh_get_buflen(*length)); + int i,j,k,n=*length; + unsigned int *pval=warnmalloc(n*sizeof *pval); + int nframes=n/natoms/3; + int cnt=0; + int most_negative=2147483647; + for (i=0; i>8)&0xFFU; + output[2]=(((unsigned int)most_negative)>>16)&0xFFU; + output[3]=(((unsigned int)most_negative)>>24)&0xFFU; + for (i=0; i=5) + bwlzh_compress(pval,n,output+4,length); + else + bwlzh_compress_no_lz77(pval,n,output+4,length); + (*length)+=4; + free(pval); + return output; + } + else if (coding==TNG_COMPRESS_ALGO_POS_XTC3) + return Ptngc_pack_array_xtc3(input,length,natoms,speed); + else if (coding==TNG_COMPRESS_ALGO_POS_XTC2) + return Ptngc_pack_array_xtc2(coder,input,length); + else + { + unsigned char *output=NULL; + unsigned char *output_ptr=NULL; + int i; + int output_length=0; + + coder->stat_numval=0; + coder->stat_overflow=0; + /* Allocate enough memory for output */ + output=warnmalloc(8* *length*sizeof *output); + output_ptr=output; + if ((coding==TNG_COMPRESS_ALGO_TRIPLET) || + (coding==TNG_COMPRESS_ALGO_POS_TRIPLET_INTRA) || + (coding==TNG_COMPRESS_ALGO_POS_TRIPLET_ONETOONE)) + { + /* Pack triplets. */ + int ntriplets=*length/3; + /* Determine max base and maxbits */ + unsigned int max_base=1U<0) + s=1+(item-1)*2; + else if (item<0) + s=2+(-item-1)*2; + if (s>intmax) + intmax=s; + } + /* Store intmax */ + coder->pack_temporary_bits=32; + coder->pack_temporary=intmax; + Ptngc_out8bits(coder,&output_ptr); + while (intmax>=max_base) + { + max_base*=2; + maxbits++; + } + for (i=0; i0) + s[j]=1+(item-1)*2; + else if (item<0) + s[j]=2+(-item-1)*2; + } + if (pack_triplet(coder,s,&output_ptr,coding_parameter,max_base,maxbits)) + { + free(output); + return NULL; + } + } + } + else + for (i=0; i<*length; i++) + if (pack_stopbits_item(coder,input[i],&output_ptr,coding_parameter)) + { + free(output); + return NULL; + } + Ptngc_pack_flush(coder,&output_ptr); + output_length=(int)(output_ptr-output); + *length=output_length; + return output; + } +} + +static int unpack_array_stop_bits(struct coder *coder,unsigned char *packed,int *output, int length, int coding_parameter) +{ + int i,j; + unsigned int extract_mask=0x80; + unsigned char *ptr=packed; + for (i=0; i>=1; + extract_mask>>=1; + if (!extract_mask) + { + extract_mask=0x80; + ptr++; + } + } + /* Check stop bit */ + bit=*ptr & extract_mask; + extract_mask>>=1; + if (!extract_mask) + { + extract_mask=0x80; + ptr++; + } + if (bit) + { + numbits>>=1; + if (numbits<1) + numbits=1; + inserted_bits+=numbits; + insert_mask=1U<<(inserted_bits-1); + } + } while (bit); + s=(pattern+1)/2; + if ((pattern%2)==0) + s=-s; + output[i]=s; + } + return 0; +} + +static int unpack_array_triplet(struct coder *coder,unsigned char *packed,int *output, int length, int coding_parameter) +{ + int i,j; + unsigned int extract_mask=0x80; + unsigned char *ptr=packed; + /* Determine max base and maxbits */ + unsigned int max_base=1U<=max_base) + { + max_base*=2; + maxbits++; + } + length/=3; + for (i=0; i>=1; + if (!extract_mask) + { + extract_mask=0x80; + ptr++; + } + } + if (jbase==3) + numbits=maxbits; + else + numbits=coding_parameter+jbase; + for (j=0; j<3; j++) + { + int s; + unsigned int jbit; + unsigned int pattern=0; + for (jbit=0; jbit>=1; + if (!extract_mask) + { + extract_mask=0x80; + ptr++; + } + } + s=(pattern+1)/2; + if ((pattern%2)==0) + s=-s; + output[i*3+j]=s; + } + } + return 0; +} + +static int unpack_array_bwlzh(struct coder *coder,unsigned char *packed,int *output, int length, int natoms) +{ + int i,j,k,n=length; + unsigned int *pval=warnmalloc(n*sizeof *pval); + int nframes=n/natoms/3; + int cnt=0; + int most_negative=(int)(((unsigned int)packed[0]) | + (((unsigned int)packed[1])<<8) | + (((unsigned int)packed[2])<<16) | + (((unsigned int)packed[3])<<24)); + bwlzh_decompress(packed+4,length,pval); + for (i=0; i +#include "dict.h" + +void Ptngc_comp_canonical_dict(unsigned int *dict, int *ndict) +{ + int i; + for (i=0; i<0x20004; i++) + dict[i]=i; + *ndict=0x20004; +} + +void Ptngc_comp_make_dict_hist(unsigned int *vals, int nvals, + unsigned int *dict, int *ndict, + unsigned int *hist) +{ + int i; + int j=0; + for (i=0; i<0x20004; i++) + hist[i]=0; + for (i=0; i<0x20004; i++) + dict[i]=i; + for (i=0; i +#endif + +#include +#include +#include "fixpoint.h" + +#define MAX32BIT 4294967295UL +#define MAX31BIT 2147483647UL +#define SIGN32BIT 2147483648UL + +/* Conversion routines from / to double precision */ + +/* Positive double to 32 bit fixed point value */ +fix_t Ptngc_ud_to_fix_t(double d,double max) +{ + fix_t val; + if (d<0.) + d=0.; + if (d>max) + d=max; + val=(fix_t)(MAX32BIT*(d/max)); + if (val>MAX32BIT) + val=MAX32BIT; + return val; +} + +/* double to signed 32 bit fixed point value */ +fix_t Ptngc_d_to_fix_t(double d,double max) +{ + fix_t val; + int sign=0; + if (d<0.) + { + sign=1; + d=-d; + } + if (d>max) + d=max; + val=(fix_t)(MAX31BIT*(d/max)); + if (val>MAX31BIT) + val=MAX31BIT; + if (sign) + val|=SIGN32BIT; + return val; +} + + +/* 32 bit fixed point value to positive double */ +double Ptngc_fix_t_to_ud(fix_t f, double max) +{ + return (double)f*(max/MAX32BIT); +} + +/* signed 32 bit fixed point value to double */ +double Ptngc_fix_t_to_d(fix_t f, double max) +{ + int sign=0; + double d; + if (f&SIGN32BIT) + { + sign=1; + f&=MAX31BIT; + } + d=(double)f*(max/MAX31BIT); + if (sign) + d=-d; + return d; +} + + +/* Convert a floating point variable to two 32 bit integers with range + -2.1e9 to 2.1e9 and precision to somewhere around 1e-9. */ +void Ptngc_d_to_i32x2(double d, fix_t *hi, fix_t *lo) +{ + int sign=0; + double frac; + double ent; + fix_t val,vallo; + if (d<0.) + { + sign=1; + d=-d; + } + /* First the integer part */ + ent=floor(d); + /* Then the fractional part */ + frac=d-ent; + + val=(fix_t)ent; + if (sign) + val|=SIGN32BIT; + + vallo=Ptngc_ud_to_fix_t(frac,1.); + + *hi=val; + *lo=vallo; +} + +/* Convert two 32 bit integers to a floating point variable + -2.1e9 to 2.1e9 and precision to somewhere around 1e-9. */ +double Ptngc_i32x2_to_d(fix_t hi, fix_t lo) +{ + double ent,frac=0.; + double val=0.; + int sign=0; + if (hi&SIGN32BIT) + { + sign=1; + hi&=MAX31BIT; + } + ent=(double)hi; + frac=Ptngc_fix_t_to_ud(lo,1.); + val=ent+frac; + if (sign) + val=-val; + return val; +} + diff --git a/src/compression/huffman.c b/src/compression/huffman.c new file mode 100644 index 0000000..1272e46 --- /dev/null +++ b/src/compression/huffman.c @@ -0,0 +1,579 @@ +/* This code is part of the tng compression routines. + * + * Written by Daniel Spangberg + * Copyright (c) 2010, 2013, The GROMACS development team. + * + * + * This program is free software; you can redistribute it and/or + * modify it under the terms of the GNU Lesser General Public License + * as published by the Free Software Foundation; either version 2.1 + * of the License, or (at your option) any later version. + */ + + +#include +#include +#include +#include "warnmalloc.h" +#include "merge_sort.h" +#include "huffman.h" + +#define MAX_HUFFMAN_LEN 31 + +enum htree_type { htree_leaf, htree_node }; + +struct htree_leaf +{ + enum htree_type nodeleaf; + unsigned int idict; /* Index into input dictionary */ + unsigned int prob; + unsigned int bit; /* One or zero */ +}; + +struct htree_node +{ + enum htree_type nodeleaf; + union htree_nodeleaf *n1; + union htree_nodeleaf *n2; + unsigned int bit; /* One or zero */ + unsigned int prob; +}; + +union htree_nodeleaf +{ + enum htree_type nodeleaf; + struct htree_node node; + struct htree_leaf leaf; +}; + +struct codelength +{ + unsigned int code; + int length; + unsigned int dict; + unsigned int prob; +}; + +static int comp_htree(const void *leafptr1, const void *leafptr2, const void *private) +{ + const union htree_nodeleaf *leaf1=(union htree_nodeleaf *)leafptr1; + const union htree_nodeleaf *leaf2=(union htree_nodeleaf *)leafptr2; + + int rval=0; + if (leaf1->leaf.probleaf.prob) + rval=1; + else if (leaf1->leaf.prob>leaf2->leaf.prob) + rval=-1; + return rval; +} + +static void assign_codes(union htree_nodeleaf *htree, + struct codelength *codelength, + unsigned int code, + int length, + int top) +{ +#if 0 + printf("Assign codes called with code %d length %d\n",code,length); +#endif + if (htree->nodeleaf==htree_leaf) + { + codelength[htree->leaf.idict].length=length+1; + codelength[htree->leaf.idict].code=(code<<1)|htree->leaf.bit; +#if 0 + printf("I am a leaf: %d %d\n", + codelength[htree->leaf.idict].length, + codelength[htree->leaf.idict].code); +#endif + } + else + { + if (!top) + { + code<<=1; + code|=htree->node.bit; + length++; + } +#if 0 + printf("I am a node length: %d\n",length); + printf("I am a node code: %d\n",code); +#endif + assign_codes(htree->node.n1,codelength,code,length,0); + assign_codes(htree->node.n2,codelength,code,length,0); + } +} + +static void free_nodes(union htree_nodeleaf *htree, int top) +{ + if (htree->nodeleaf==htree_leaf) + { + if (!top) + free(htree); + } + else + { + free_nodes(htree->node.n1,0); + free_nodes(htree->node.n2,0); + if (!top) + free(htree); + } +} + +static void flush_8bits(unsigned int *combine, unsigned char **output, int *bitptr) +{ + while ((*bitptr)>=8) + { + unsigned int mask=~(0xFFU<<((*bitptr)-8)); + unsigned char out=(unsigned char)((*combine)>>((*bitptr)-8)); + **output=out; + (*output)++; + (*bitptr)-=8; + (*combine)&=mask; + } +} + +static void writebits(unsigned int value,int length, unsigned char **output, int *bitptr) +{ + unsigned int mask; + unsigned int combine=(unsigned int)**output; + if (length>=8) + mask=0xFFU<<(length-8); + else + mask=0xFFU>>(8-length); + while (length>8) + { + /* Make room for the bits. */ + combine<<=8; + (*bitptr)+=8; + combine|=(value&mask)>>(length-8); + flush_8bits(&combine,output,bitptr); + length-=8; + mask>>=8; + } + if (length) + { + /* Make room for the bits. */ + combine<<=length; + (*bitptr)+=length; + combine|=value; + flush_8bits(&combine,output,bitptr); + } + **output=(unsigned char)combine; +} + +static unsigned int readbits(int length, unsigned char **input, int *bitptr) +{ + unsigned int val=0U; + unsigned int extract_mask=0x80U>>*bitptr; + unsigned char thisval=**input; + while (length--) + { + val<<=1; + val|=((extract_mask & thisval)!=0); + *bitptr=(*bitptr)+1; + extract_mask>>=1; + if (!extract_mask) + { + extract_mask=0x80U; + *input=(*input)+1; + *bitptr=0; + thisval=**input; + } + } + return val; +} + +static int comp_codes(const void *codeptr1, const void *codeptr2, const void *private) +{ + const struct codelength *code1=(struct codelength *)codeptr1; + const struct codelength *code2=(struct codelength *)codeptr2; + + int rval=0; /* It shouldn't be possible to get equal here, though. */ + if (code1->length>code2->length) + rval=1; + else if (code1->lengthlength) + rval=-1; + else if (code1->dict>code2->dict) + rval=1; + else + rval=-1; + return rval; +} + +static int comp_codes_value(const void *codeptr1, const void *codeptr2, const void *private) +{ + const struct codelength *code1=(struct codelength *)codeptr1; + const struct codelength *code2=(struct codelength *)codeptr2; + + int rval=0; /* It shouldn't be possible to get equal here, though. */ + if (code1->dict>code2->dict) + rval=1; + else + rval=-1; + return rval; +} + +/* The huffman_dict array should be 131077 (0x20005) long. The +huffman_dict_unpacked array should be 131077 long (note five longer than +0x20000) */ +void Ptngc_comp_conv_to_huffman(unsigned int *vals, int nvals, + unsigned int *dict, int ndict, + unsigned int *prob, + unsigned char *huffman, + int *huffman_len, + unsigned char *huffman_dict, + int *huffman_dictlen, + unsigned int *huffman_dict_unpacked, + int *huffman_dict_unpackedlen) +{ + int i; + int nleft; + union htree_nodeleaf *htree; + struct codelength *codelength; + int bitptr; + unsigned char *huffman_ptr; + int code; + int longcodes=1; + while (longcodes) + { + /* Create array of leafs (will be array of nodes/trees during + buildup of tree. */ + htree=warnmalloc(ndict*sizeof *htree); + codelength=warnmalloc(ndict*sizeof *codelength); + bitptr=0; + huffman_ptr=huffman; + for (i=0; i1) + { + union htree_nodeleaf *n1=warnmalloc(sizeof *n1); + union htree_nodeleaf *n2=warnmalloc(sizeof *n2); + int new_place; + int p1,p2, new_prob; + *n1=htree[nleft-1]; + *n2=htree[nleft-2]; + if (n1->nodeleaf==htree_leaf) + { + p1=n1->leaf.prob; + n1->leaf.bit=0; + } + else + { + p1=n1->node.prob; + n1->node.bit=0; + } + if (n2->nodeleaf==htree_leaf) + { + p2=n2->leaf.prob; + n2->leaf.bit=1; + } + else + { + p2=n2->node.prob; + n2->node.bit=1; + } + nleft--; + /* Create a new node */ + htree[nleft-1].nodeleaf=htree_node; + htree[nleft-1].node.n1=n1; + htree[nleft-1].node.n2=n2; + new_prob=p1+p2; + htree[nleft-1].node.prob=new_prob; + /* Use insertion sort to place this in the correct place in the + array. */ + /* Where should it be inserted? */ + new_place=nleft; + while (new_place>0) + { + int pc; + if (htree[new_place-1].nodeleaf==htree_node) + pc=htree[new_place-1].node.prob; + else + pc=htree[new_place-1].leaf.prob; + if (new_probMAX_HUFFMAN_LEN) + longcodes=1; + + /* If the codes are too long alter the probabilities. */ + if (longcodes) + { + for (i=0; i>=1; + if (prob[i]==0) + prob[i]=1; + } + + /* Free codelength. We will compute a new one. */ + free(codelength); + } + } + +#if 0 + { + for (i=0; i=0; j--) + { + int bit=c&mask; + if (bit) + printf("1"); + else + printf("0"); + mask>>=1; + } + printf("\n"); + } + } + } +#endif + + /* Simply do compression by writing out the bits. */ + for (i=0; i>8)&0xFFU); + *huffman_ptr++=(unsigned char)((codelength[ndict-1].dict>>16)&0xFFU); + huffman_dict_unpacked[0]=(unsigned char)(codelength[ndict-1].dict&0xFFU); + huffman_dict_unpacked[1]=(unsigned char)((codelength[ndict-1].dict>>8)&0xFFU); + huffman_dict_unpacked[2]=(unsigned char)((codelength[ndict-1].dict>>16)&0xFFU); + for (i=0; i<=(int)codelength[ndict-1].dict; i++) + { + /* Do I have this value? */ + int ihave=0; + int j; + for (j=0; j=0; j--) + { + int bit=c&mask; + if (bit) + printf("1"); + else + printf("0"); + mask>>=1; + } + printf("\n"); + } + } + } +#endif + /* Decompress data. */ + huffman_ptr=huffman; + bitptr=0; + for (i=0; i +#include +#include "warnmalloc.h" +#include "tng_compress.h" +#include "bwlzh.h" +#include "huffman.h" +#include "dict.h" +#include "rle.h" +#include "vals16.h" + +int Ptngc_comp_huff_buflen(int nvals) +{ + return 132000+nvals*8; +} + +/* the value pointed to by chosen_algo should be sent as -1 for autodetect. */ +void Ptngc_comp_huff_compress_verbose(unsigned int *vals, int nvals, + unsigned char *huffman, int *huffman_len, + int *huffdatalen, + int *huffman_lengths,int *chosen_algo, + int isvals16) +{ + unsigned int *dict=warnmalloc(0x20005*sizeof *dict); + unsigned int *hist=warnmalloc(0x20005*sizeof *hist); + unsigned int *vals16=NULL; + unsigned char *huffdict=warnmalloc(0x20005*sizeof *huffdict); + unsigned int *huffdictunpack=warnmalloc(0x20005*sizeof *huffdictunpack); + unsigned char *huffman1=warnmalloc(2*0x20005*sizeof *huffman1); + unsigned char *huffdict1=warnmalloc(0x20005*sizeof *huffdict1); + unsigned int *huffdictunpack1=warnmalloc(0x20005*sizeof *huffdictunpack1); + unsigned int *huffdictrle=warnmalloc((3*0x20005+3)*sizeof *huffdictrle); + unsigned char *huffman2=warnmalloc(6*0x20005*sizeof *huffman2); + unsigned char *huffdict2=warnmalloc(0x20005*sizeof *huffdict2); + unsigned int *huffdictunpack2=warnmalloc(0x20005*sizeof *huffdictunpack2); + int i; + int ndict,ndict1,ndict2; + int nhuff,nhuffdict,nhuffdictunpack; + int nhuff1,nhuffdict1,nhuffdictunpack1; + int nhuffrle,nhuff2,nhuffdict2,nhuffdictunpack2; + int nvals16; + + /* Do I need to convert to vals16? */ + if (!isvals16) + { + vals16=warnmalloc(nvals*3*sizeof *vals16); + Ptngc_comp_conv_to_vals16(vals,nvals,vals16,&nvals16); + nvals=nvals16; + vals=vals16; + } + else + nvals16=nvals; + + /* Determine probabilities. */ + Ptngc_comp_make_dict_hist(vals,nvals,dict,&ndict,hist); + + /* First compress the data using huffman coding (place it ready for output at 14 (code for algorithm+length etc.). */ + Ptngc_comp_conv_to_huffman(vals,nvals,dict,ndict,hist, + huffman+14,&nhuff, + huffdict,&nhuffdict, + huffdictunpack,&nhuffdictunpack); + *huffdatalen=nhuff; + + /* Algorithm 0 stores the huffman dictionary directly (+ a code for + the algorithm) + lengths of the huffman buffer (4) and the huffman dictionary (3). */ + huffman_lengths[0]=nhuff+nhuffdict+1*2+3*4+3+3; + /* Next we try to compress the huffman dictionary using huffman + coding ... (algorithm 1) */ + + /* Determine probabilities. */ + Ptngc_comp_make_dict_hist(huffdictunpack,nhuffdictunpack,dict,&ndict1,hist); + /* Pack huffman dictionary */ + Ptngc_comp_conv_to_huffman(huffdictunpack,nhuffdictunpack, + dict,ndict1,hist, + huffman1,&nhuff1, + huffdict1,&nhuffdict1, + huffdictunpack1,&nhuffdictunpack1); + huffman_lengths[1]=nhuff+nhuff1+nhuffdict1+1*2+3*4+3+3+3+3+3; + + /* ... and rle + huffman coding ... (algorithm 2) Pack any repetetitive patterns. */ + Ptngc_comp_conv_to_rle(huffdictunpack,nhuffdictunpack, + huffdictrle,&nhuffrle,1); + + /* Determine probabilities. */ + Ptngc_comp_make_dict_hist(huffdictrle,nhuffrle,dict,&ndict2,hist); + /* Pack huffman dictionary */ + Ptngc_comp_conv_to_huffman(huffdictrle,nhuffrle, + dict,ndict2,hist, + huffman2,&nhuff2, + huffdict2,&nhuffdict2, + huffdictunpack2,&nhuffdictunpack2); + huffman_lengths[2]=nhuff+nhuff2+nhuffdict2+1*2+3*4+3+3+3+3+3+3; + + /* Choose the best algorithm and output the data. */ + if ((*chosen_algo==0) || ((*chosen_algo==-1) && + (((huffman_lengths[0]>8)&0xFFU; + huffman[4]=(((unsigned int)nvals16)>>16)&0xFFU; + huffman[5]=(((unsigned int)nvals16)>>24)&0xFFU; + huffman[6]=((unsigned int)nvals)&0xFFU; + huffman[7]=(((unsigned int)nvals)>>8)&0xFFU; + huffman[8]=(((unsigned int)nvals)>>16)&0xFFU; + huffman[9]=(((unsigned int)nvals)>>24)&0xFFU; + huffman[10]=((unsigned int)nhuff)&0xFFU; + huffman[11]=(((unsigned int)nhuff)>>8)&0xFFU; + huffman[12]=(((unsigned int)nhuff)>>16)&0xFFU; + huffman[13]=(((unsigned int)nhuff)>>24)&0xFFU; + huffman[14+nhuff]=((unsigned int)nhuffdict)&0xFFU; + huffman[15+nhuff]=(((unsigned int)nhuffdict)>>8)&0xFFU; + huffman[16+nhuff]=(((unsigned int)nhuffdict)>>16)&0xFFU; + huffman[17+nhuff]=((unsigned int)ndict)&0xFFU; + huffman[18+nhuff]=(((unsigned int)ndict)>>8)&0xFFU; + huffman[19+nhuff]=(((unsigned int)ndict)>>16)&0xFFU; + for (i=0; i>8)&0xFFU; + huffman[4]=(((unsigned int)nvals16)>>16)&0xFFU; + huffman[5]=(((unsigned int)nvals16)>>24)&0xFFU; + huffman[6]=((unsigned int)nvals)&0xFFU; + huffman[7]=(((unsigned int)nvals)>>8)&0xFFU; + huffman[8]=(((unsigned int)nvals)>>16)&0xFFU; + huffman[9]=(((unsigned int)nvals)>>24)&0xFFU; + huffman[10]=((unsigned int)nhuff)&0xFFU; + huffman[11]=(((unsigned int)nhuff)>>8)&0xFFU; + huffman[12]=(((unsigned int)nhuff)>>16)&0xFFU; + huffman[13]=(((unsigned int)nhuff)>>24)&0xFFU; + huffman[14+nhuff]=((unsigned int)nhuffdictunpack)&0xFFU; + huffman[15+nhuff]=(((unsigned int)nhuffdictunpack)>>8)&0xFFU; + huffman[16+nhuff]=(((unsigned int)nhuffdictunpack)>>16)&0xFFU; + huffman[17+nhuff]=((unsigned int)ndict)&0xFFU; + huffman[18+nhuff]=(((unsigned int)ndict)>>8)&0xFFU; + huffman[19+nhuff]=(((unsigned int)ndict)>>16)&0xFFU; + huffman[20+nhuff]=((unsigned int)nhuff1)&0xFFU; + huffman[21+nhuff]=(((unsigned int)nhuff1)>>8)&0xFFU; + huffman[22+nhuff]=(((unsigned int)nhuff1)>>16)&0xFFU; + huffman[23+nhuff]=((unsigned int)nhuffdict1)&0xFFU; + huffman[24+nhuff]=(((unsigned int)nhuffdict1)>>8)&0xFFU; + huffman[25+nhuff]=(((unsigned int)nhuffdict1)>>16)&0xFFU; + huffman[26+nhuff]=((unsigned int)ndict1)&0xFFU; + huffman[27+nhuff]=(((unsigned int)ndict1)>>8)&0xFFU; + huffman[28+nhuff]=(((unsigned int)ndict1)>>16)&0xFFU; + for (i=0; i>8)&0xFFU; + huffman[4]=(((unsigned int)nvals16)>>16)&0xFFU; + huffman[5]=(((unsigned int)nvals16)>>24)&0xFFU; + huffman[6]=((unsigned int)nvals)&0xFFU; + huffman[7]=(((unsigned int)nvals)>>8)&0xFFU; + huffman[8]=(((unsigned int)nvals)>>16)&0xFFU; + huffman[9]=(((unsigned int)nvals)>>24)&0xFFU; + huffman[10]=((unsigned int)nhuff)&0xFFU; + huffman[11]=(((unsigned int)nhuff)>>8)&0xFFU; + huffman[12]=(((unsigned int)nhuff)>>16)&0xFFU; + huffman[13]=(((unsigned int)nhuff)>>24)&0xFFU; + huffman[14+nhuff]=((unsigned int)nhuffdictunpack)&0xFFU; + huffman[15+nhuff]=(((unsigned int)nhuffdictunpack)>>8)&0xFFU; + huffman[16+nhuff]=(((unsigned int)nhuffdictunpack)>>16)&0xFFU; + huffman[17+nhuff]=((unsigned int)ndict)&0xFFU; + huffman[18+nhuff]=(((unsigned int)ndict)>>8)&0xFFU; + huffman[19+nhuff]=(((unsigned int)ndict)>>16)&0xFFU; + huffman[20+nhuff]=((unsigned int)nhuffrle)&0xFFU; + huffman[21+nhuff]=(((unsigned int)nhuffrle)>>8)&0xFFU; + huffman[22+nhuff]=(((unsigned int)nhuffrle)>>16)&0xFFU; + huffman[23+nhuff]=((unsigned int)nhuff2)&0xFFU; + huffman[24+nhuff]=(((unsigned int)nhuff2)>>8)&0xFFU; + huffman[25+nhuff]=(((unsigned int)nhuff2)>>16)&0xFFU; + huffman[26+nhuff]=((unsigned int)nhuffdict2)&0xFFU; + huffman[27+nhuff]=(((unsigned int)nhuffdict2)>>8)&0xFFU; + huffman[28+nhuff]=(((unsigned int)nhuffdict2)>>16)&0xFFU; + huffman[29+nhuff]=((unsigned int)ndict2)&0xFFU; + huffman[30+nhuff]=(((unsigned int)ndict2)>>8)&0xFFU; + huffman[31+nhuff]=(((unsigned int)ndict2)>>16)&0xFFU; + for (i=0; i=N_HUFFMAN_ALGO) + return NULL; + return huff_algo_names[algo]; +} diff --git a/src/compression/lz77.c b/src/compression/lz77.c new file mode 100644 index 0000000..705e89d --- /dev/null +++ b/src/compression/lz77.c @@ -0,0 +1,347 @@ +/* This code is part of the tng compression routines. + * + * Written by Daniel Spangberg + * Copyright (c) 2010, 2013, The GROMACS development team. + * + * + * This program is free software; you can redistribute it and/or + * modify it under the terms of the GNU Lesser General Public License + * as published by the Free Software Foundation; either version 2.1 + * of the License, or (at your option) any later version. + */ + + +#include +#include +#include +#include "warnmalloc.h" +#include "bwt.h" +#include "lz77.h" + +/* This is a simple Lempel-Ziv-77 compressor. It has not been set up + to remove every possible repeating pattern, but it might be better + than simple RLE. + + Lempel-Ziv 77 with separate outputs for length, data, and offsets. + */ + +#if 0 +/* Sort the strings (similar to BWT) to find similar patterns in the + input data. The output is an array with two values for each + input value. The sorted index value is the first in each doublet. + The second value is the "inverse" of the first value and can be + used to find the locations of similar strings. */ +static void sort_strings(unsigned int *vals, int nvals, + unsigned int *output) +{ + int i; + int *indices=warnmalloc(2*nvals*sizeof *indices); + unsigned int *nrepeat=warnmalloc(nvals*sizeof *nrepeat); + int *warr=indices+nvals; + + if (nvals>0xFFFFFF) + { + fprintf(stderr,"BWT cannot pack more than %d values.\n",0xFFFFFF); + exit(1); + } + + /* Also note that repeat pattern k (kmax) cannot be larger than 255. */ + for (i=0; i=1; k--) + { + try_next_k: + if (k>=1) + { + for (j=k; jmaxrepeat) + new_j=j; + if ((new_j>good_j) || ((new_j==good_j) && (knvals) + repeat=nvals; + nrepeat[i+m]=((unsigned int) (good_k)) | (((unsigned int) (repeat))<<8); + } + /* If no repetition was found for this value signal that here. */ + if (!nrepeat[i]) + nrepeat[i+m]=257U; /* This is 1<<8 | 1 */ + } + } + + /* Sort cyclic shift matrix. */ + Ptngc_bwt_merge_sort_inner(indices,nvals,vals,0,nvals,nrepeat,warr); + + /* Form output. */ + for (i=0; iNUM_PREVIOUS) + previous[(NUM_PREVIOUS+3)*v]=NUM_PREVIOUS; + previous[(NUM_PREVIOUS+3)*v+3+previous[(NUM_PREVIOUS+3)*v+1]]=i; + previous[(NUM_PREVIOUS+3)*v+1]++; + if (previous[(NUM_PREVIOUS+3)*v+1]>=NUM_PREVIOUS) + previous[(NUM_PREVIOUS+3)*v+1]=0; + } + previous[(NUM_PREVIOUS+3)*v+2]=i; +} + +void Ptngc_comp_to_lz77(unsigned int *vals, int nvals, + unsigned int *data, int *ndata, + unsigned int *len, int *nlens, + unsigned int *offsets, int *noffsets) +{ + int noff=0; + int ndat=0; + int nlen=0; + int i,j; + int *previous=warnmalloc(0x20000*(NUM_PREVIOUS+3)*sizeof *previous); +#if 0 + unsigned int *info=warnmalloc(2*nvals*sizeof *info); + sort_strings(vals,nvals,info); +#endif + for (i=0; i<0x20000; i++) + { + previous[(NUM_PREVIOUS+3)*i]=0; /* Number of items in a circular buffer */ + previous[(NUM_PREVIOUS+3)*i+1]=0; /* Pointer to beginning of circular buffer. */ + previous[(NUM_PREVIOUS+3)*i+2]=-2; /* Last offset that had this value. -2 is really never... */ + } + for (i=0; i=firstoffset) + { + for (k=0; i+klargest_len) && ((k>=(i-j)+16) || ((k>4) && (i-j==1)))) + { + largest_len=k; + largest_offset=j; + } + } + j++; + } + } +#if 0 + /* Search in sorted string buffer too. */ + kmin=info[i*2+1]-MAX_STRING_SEARCH; + kmax=info[i*2+1]+MAX_STRING_SEARCH; + if (kmin<0) + kmin=0; + if (kmax>=nvals) + kmax=nvals; + for (k=kmin; k=i)) + { + for (m=0; i+mlargest_len) && (m>4) && (m+2>=(i-s))) + { + largest_len=m; + largest_offset=s; +#if 0 + fprintf(stderr,"Offset: %d %d\n",m,i-s); +#endif + } + } + } +#endif + /* Check how to write this info. */ + if (largest_len>MAX_LEN) + largest_len=MAX_LEN; + if (largest_len) + { + if (i-largest_offset==1) + { + data[ndat++]=0; + } + else + { + data[ndat++]=1; + offsets[noff++]=i-largest_offset; + } + len[nlen++]=largest_len; +#if 0 + fprintf(stderr,"l:o: %d:%d data=%d i=%d\n",largest_len,i-largest_offset,ndat,i); + fflush(stderr); +#endif + +#if 0 + fprintf(stderr,"Found largest len %d at %d.\n",largest_len,i-largest_offset); +#endif + /* Add these values to the circular buffer. */ + for (k=0; k=nvals) + { + fprintf(stderr,"too many vals.\n"); + exit(EXIT_FAILURE); + } + i++; + } + } + else + vals[i++]=v-2; + } +} diff --git a/src/compression/merge_sort.c b/src/compression/merge_sort.c new file mode 100644 index 0000000..75dd550 --- /dev/null +++ b/src/compression/merge_sort.c @@ -0,0 +1,130 @@ +/* This code is part of the tng compression routines. + * + * Written by Daniel Spangberg + * Copyright (c) 2010, 2013, The GROMACS development team. + * + * + * This program is free software; you can redistribute it and/or + * modify it under the terms of the GNU Lesser General Public License + * as published by the Free Software Foundation; either version 2.1 + * of the License, or (at your option) any later version. + */ + + +#include +#include +#include +#include "warnmalloc.h" +#include "merge_sort.h" + +static void ms_inner(void *base, size_t size, + int start, int end, + int (*compar)(const void *v1,const void *v2,const void *private), + const void *private, + char *workarray) +{ + int middle; + if ((end-start)>1) + { + char *cbase=(char *)base; + middle=start+(end-start)/2; +#if 0 + printf("For start %d end %d obtained new middle: %d\n",start,end,middle); +#endif + ms_inner(base,size, + start,middle, + compar,private,workarray); + ms_inner(base,size, + middle,end, + compar,private,workarray); +#if 0 + printf("For start %d end %d Before merge: Comparing element %d with %d\n",start,end,middle-1,middle); +#endif + if (compar(cbase+(middle-1)*size,cbase+middle*size,private)>0) + { + /* Merge to work array. */ + int i, n=end-start; + int ileft=start; + int iright=middle; + for (i=0; i0) + { + memcpy(workarray+i*size,cbase+iright*size,size); + iright++; + } + else + { + memcpy(workarray+i*size,cbase+ileft*size,size); + ileft++; + } + } + } + /* Copy result back. */ + memcpy(cbase+start*size,workarray,(end-start)*size); + } + } +} + + +void Ptngc_merge_sort(void *base, size_t nmemb, size_t size, + int (*compar)(const void *v1,const void *v2,const void *private), + void *private) +{ + char *warr=warnmalloc(nmemb*size); + ms_inner(base,size,0,nmemb,compar,private,warr); + free(warr); +} + + +#ifdef TEST + +static int compint(const void *v1, const void *v2,const void *private) +{ + const int *i1=(const int *)v1; + const int *i2=(const int *)v2; + if (*i1<*i2) + return -1; + else if (*i1>*i2) + return 1; + else + return 0; +} + +static int qcompint(const void *v1, const void *v2) +{ + return compint(v1,v2,NULL); +} + +#define N 1000000 +int main() +{ + int *arr=warnmalloc(N*sizeof *arr); + int i; + for (i=0; i +#include "warnmalloc.h" +#include "mtf.h" + +/* "Partial" MTF. Byte based. */ +/* Move to front coding. + Acceptable inputs are max 8 bits (0-0xFF) */ +static void comp_conv_to_mtf_byte(unsigned char *vals, int nvals, + unsigned char *valsmtf) +{ + int i; + /* Indices into a linked list */ + int list[256]; + int dict[256]; + /* Head of the linked list */ + int head; + for (i=0; i<256; i++) + dict[i]=i; + for (i=0; i<255; i++) + list[i]=i+1; + list[255]=-1; /* end. */ + head=0; + for (i=0; i>(8*j))&0xFF); + comp_conv_to_mtf_byte(tmp,nvals,tmp+nvals); + for (i=0; i>(8*j))&0xFF); + comp_conv_to_mtf_byte(tmp,nvals,valsmtf+j*nvals); + } + free(tmp); +} + +/* Move to front decoding */ +static void comp_conv_from_mtf_byte(unsigned char *valsmtf, int nvals, + unsigned char *vals) +{ + int i; + /* Indices into a linked list */ + int list[256]; + int dict[256]; + /* Head of the linked list */ + int head; + for (i=0; i<256; i++) + dict[i]=i; + for (i=0; i<255; i++) + list[i]=i+1; + list[255]=-1; /* end. */ + head=0; + for (i=0; i>(8*j))&0xFF); + comp_conv_from_mtf_byte(tmp,nvals,tmp+nvals); + for (i=0; imin_rle) + { + /* Insert run-length */ + unsigned int run=((unsigned int)nsim); + while (run>1) + { + if (run&0x1U) + rle[(*j)++]=1; + else + rle[(*j)++]=0; + run>>=1; + } + nsim=1; + } + while (nsim--) + rle[(*j)++]=v+2; +} + +/* Run length encoding. + Acceptable inputs are about 16 bits (0-0xFFFF) + If input is 0-N output will be be values of 0-(N+2) */ +void Ptngc_comp_conv_to_rle(unsigned int *vals, int nvals, + unsigned int *rle, int *nrle, + int min_rle) +{ + int i; + int j=0; + int nsim=0; + int v=-1; + for (i=0; i +#include +#include +#include "tng_compress.h" +#include "coder.h" +#include "fixpoint.h" + +/* Please see tng_compress.h for info on how to call these routines. */ + +/* This becomes TNGP for positions (little endian) and TNGV for velocities. In ASCII. */ +#define MAGIC_INT_POS 0x50474E54 +#define MAGIC_INT_VEL 0x56474E54 + +#define SPEED_DEFAULT 2 /* Default to relatively fast compression. For very good compression it makes sense to + choose speed=4 or speed=5 */ + +#define PRECISION(hi,lo) (Ptngc_i32x2_to_d(hi,lo)) + +static void quantize(double *x, int natoms, int nframes, + double precision, + int *quant) +{ + int iframe, i, j; + for (iframe=0; iframe>= 8; + c=(unsigned char)(v & 0xFFU); + } +} + +static fix_t readbufferfix(unsigned char *buf, int num) +{ + unsigned char b; + int shift=0; + fix_t f=0UL; + int cnt=0; + do + { + b=buf[cnt++]; + f |= ((fix_t)b & 0xFF)<1) + { + datablock=NULL; + /* Inter-frame compression? */ + if ((coding==TNG_COMPRESS_ALGO_POS_STOPBIT_INTER) || + (coding==TNG_COMPRESS_ALGO_POS_TRIPLET_INTER) || + (coding==TNG_COMPRESS_ALGO_POS_BWLZH_INTER)) + { + struct coder *coder=Ptngc_coder_init(); + length=natoms*3*(nframes-1); + datablock=(char*)Ptngc_pack_array(coder,quant_inter+natoms*3,&length, + coding,coding_parameter,natoms,speed); + Ptngc_coder_deinit(coder); + } + /* One-to-one compression? */ + else if ((coding==TNG_COMPRESS_ALGO_POS_XTC2) || + (coding==TNG_COMPRESS_ALGO_POS_XTC3) || + (coding==TNG_COMPRESS_ALGO_POS_TRIPLET_ONETOONE)) + { + struct coder *coder=Ptngc_coder_init(); + length=natoms*3*(nframes-1); + datablock=(char*)Ptngc_pack_array(coder,quant+natoms*3,&length, + coding,coding_parameter,natoms,speed); + Ptngc_coder_deinit(coder); + } + /* Intra-frame compression? */ + else if ((coding==TNG_COMPRESS_ALGO_POS_TRIPLET_INTRA) || + (coding==TNG_COMPRESS_ALGO_POS_BWLZH_INTRA)) + { + struct coder *coder=Ptngc_coder_init(); + length=natoms*3*(nframes-1); + datablock=(char*)Ptngc_pack_array(coder,quant_intra+natoms*3,&length, + coding,coding_parameter,natoms,speed); + Ptngc_coder_deinit(coder); + } + /* Block length. */ + if (data) + bufferfix((unsigned char*)data+bufloc,(fix_t)length,4); + bufloc+=4; + if (data) + memcpy(data+bufloc,datablock,length); + free(datablock); + bufloc+=length; + } + *nitems=bufloc; +} + +/* Perform velocity compression from vel into the data block */ +static void compress_quantized_vel(int *quant, int *quant_inter, + int natoms, int nframes, + int speed, + int initial_coding, int initial_coding_parameter, + int coding, int coding_parameter, + fix_t prec_hi, fix_t prec_lo, + int *nitems, + char *data) +{ + int bufloc=0; + char *datablock=NULL; + int length; + /* Information needed for decompression. */ + if (data) + bufferfix((unsigned char*)data+bufloc,(fix_t)MAGIC_INT_VEL,4); + bufloc+=4; + /* Number of atoms. */ + if (data) + bufferfix((unsigned char*)data+bufloc,(fix_t)natoms,4); + bufloc+=4; + /* Number of frames. */ + if (data) + bufferfix((unsigned char*)data+bufloc,(fix_t)nframes,4); + bufloc+=4; + /* Initial coding. */ + if (data) + bufferfix((unsigned char*)data+bufloc,(fix_t)initial_coding,4); + bufloc+=4; + /* Initial coding parameter. */ + if (data) + bufferfix((unsigned char*)data+bufloc,(fix_t)initial_coding_parameter,4); + bufloc+=4; + /* Coding. */ + if (data) + bufferfix((unsigned char*)data+bufloc,(fix_t)coding,4); + bufloc+=4; + /* Coding parameter. */ + if (data) + bufferfix((unsigned char*)data+bufloc,(fix_t)coding_parameter,4); + bufloc+=4; + /* Precision. */ + if (data) + bufferfix((unsigned char*)data+bufloc,prec_lo,4); + bufloc+=4; + if (data) + bufferfix((unsigned char*)data+bufloc,prec_hi,4); + bufloc+=4; + /* The initial frame */ + if ((initial_coding==TNG_COMPRESS_ALGO_VEL_STOPBIT_ONETOONE) || + (initial_coding==TNG_COMPRESS_ALGO_VEL_TRIPLET_ONETOONE) || + (initial_coding==TNG_COMPRESS_ALGO_VEL_BWLZH_ONETOONE)) + { + struct coder *coder=Ptngc_coder_init(); + length=natoms*3; + datablock=(char*)Ptngc_pack_array(coder,quant,&length, + initial_coding,initial_coding_parameter,natoms,speed); + Ptngc_coder_deinit(coder); + } + /* Block length. */ + if (data) + bufferfix((unsigned char*)data+bufloc,(fix_t)length,4); + bufloc+=4; + /* The actual data block. */ + if (data) + memcpy(data+bufloc,datablock,length); + free(datablock); + bufloc+=length; + /* The remaining frames */ + if (nframes>1) + { + datablock=NULL; + /* Inter-frame compression? */ + if ((coding==TNG_COMPRESS_ALGO_VEL_TRIPLET_INTER) || + (coding==TNG_COMPRESS_ALGO_VEL_STOPBIT_INTER) || + (coding==TNG_COMPRESS_ALGO_VEL_BWLZH_INTER)) + { + struct coder *coder=Ptngc_coder_init(); + length=natoms*3*(nframes-1); + datablock=(char*)Ptngc_pack_array(coder,quant_inter+natoms*3,&length, + coding,coding_parameter,natoms,speed); + Ptngc_coder_deinit(coder); + } + /* One-to-one compression? */ + else if ((coding==TNG_COMPRESS_ALGO_VEL_STOPBIT_ONETOONE) || + (coding==TNG_COMPRESS_ALGO_VEL_TRIPLET_ONETOONE) || + (coding==TNG_COMPRESS_ALGO_VEL_BWLZH_ONETOONE)) + { + struct coder *coder=Ptngc_coder_init(); + length=natoms*3*(nframes-1); + datablock=(char*)Ptngc_pack_array(coder,quant+natoms*3,&length, + coding,coding_parameter,natoms,speed); + Ptngc_coder_deinit(coder); + } + /* Block length. */ + if (data) + bufferfix((unsigned char*)data+bufloc,(fix_t)length,4); + bufloc+=4; + if (data) + memcpy(data+bufloc,datablock,length); + free(datablock); + bufloc+=length; + } + *nitems=bufloc; +} + +static int determine_best_coding_stop_bits(struct coder *coder,int *input, int *length, + int *coding_parameter, int natoms) +{ + int bits; + unsigned char *packed; + int best_length=0; + int new_parameter=-1; + int io_length; + for (bits=1; bits<20; bits++) + { + io_length=*length; + packed=Ptngc_pack_array(coder,input,&io_length, + TNG_COMPRESS_ALGO_STOPBIT,bits,natoms,0); + if (packed) + { + if ((new_parameter==-1) || (io_length=2) + { + current_coding=TNG_COMPRESS_ALGO_POS_XTC3; + current_coding_parameter=0; + compress_quantized_pos(quant,NULL,quant_intra,natoms,1,speed, + current_coding,current_coding_parameter, + 0,0,prec_hi,prec_lo,¤t_code_size,NULL); + if (current_code_size=6) + { + current_coding=TNG_COMPRESS_ALGO_POS_BWLZH_INTRA; + current_coding_parameter=0; + compress_quantized_pos(quant,NULL,quant_intra,natoms,1,speed, + current_coding,current_coding_parameter, + 0,0,prec_hi,prec_lo,¤t_code_size,NULL); + if (current_code_size=4) + { + current_coding=TNG_COMPRESS_ALGO_POS_BWLZH_INTER; + current_coding_parameter=0; + compress_quantized_pos(quant,quant_inter,quant_intra,natoms,nframes,speed, + TNG_COMPRESS_ALGO_POS_XTC2,0, + current_coding,current_coding_parameter, + prec_hi,prec_lo,¤t_code_size,NULL); + current_code_size-=initial_code_size; /* Correct for the use of XTC2 for the first frame. */ + if (current_code_size=6) + { + current_coding=TNG_COMPRESS_ALGO_POS_BWLZH_INTRA; + current_coding_parameter=0; + compress_quantized_pos(quant,quant_inter,quant_intra,natoms,nframes,speed, + TNG_COMPRESS_ALGO_POS_XTC2,0, + current_coding,current_coding_parameter, + prec_hi,prec_lo,¤t_code_size,NULL); + current_code_size-=initial_code_size; /* Correct for the use of XTC2 for the first frame. */ + if (current_code_size=4) + { + current_coding=TNG_COMPRESS_ALGO_VEL_BWLZH_ONETOONE; + current_coding_parameter=0; + compress_quantized_vel(quant,NULL,natoms,1,speed, + current_coding,current_coding_parameter, + 0,0,prec_hi,prec_lo,¤t_code_size,NULL); + if ((best_coding==-1) || (current_code_size=4) + { + /* Test BWLZH inter */ + current_coding=TNG_COMPRESS_ALGO_VEL_BWLZH_INTER; + current_coding_parameter=0; + compress_quantized_vel(quant,quant_inter,natoms,nframes,speed, + TNG_COMPRESS_ALGO_VEL_STOPBIT_ONETOONE,initial_numbits, + current_coding,current_coding_parameter, + prec_hi,prec_lo,¤t_code_size,NULL); + current_code_size-=initial_code_size; /* Correct for the initial frame */ + if (current_code_size6) + speed=6; + initial_coding=algo[0]; + initial_coding_parameter=algo[1]; + coding=algo[2]; + coding_parameter=algo[3]; + + quantize(pos,natoms,nframes,PRECISION(prec_hi,prec_lo),quant); + quant_inter_differences(quant,natoms,nframes,quant_inter); + quant_intra_differences(quant,natoms,nframes,quant_intra); + + /* If any of the above codings / coding parameters are == -1, the optimal parameters must be found */ + if (initial_coding==-1) + { + initial_coding_parameter=-1; + determine_best_pos_initial_coding(quant,quant_intra,natoms,speed,prec_hi,prec_lo, + &initial_coding,&initial_coding_parameter); + } + else if (initial_coding_parameter==-1) + { + determine_best_pos_initial_coding(quant,quant_intra,natoms,speed,prec_hi,prec_lo, + &initial_coding,&initial_coding_parameter); + } + + if (nframes==1) + { + coding=0; + coding_parameter=0; + } + + if (nframes>1) + { + if (coding==-1) + { + coding_parameter=-1; + determine_best_pos_coding(quant,quant_inter,quant_intra,natoms,nframes,speed,prec_hi,prec_lo, + &coding,&coding_parameter); + } + else if (coding_parameter==-1) + { + determine_best_pos_coding(quant,quant_inter,quant_intra,natoms,nframes,speed,prec_hi,prec_lo, + &coding,&coding_parameter); + } + } + + compress_quantized_pos(quant,quant_inter,quant_intra,natoms,nframes,speed, + initial_coding,initial_coding_parameter, + coding,coding_parameter, + prec_hi,prec_lo,nitems,data); + free(quant_inter); + free(quant_intra); + free(quant); + if (algo[0]==-1) + algo[0]=initial_coding; + if (algo[1]==-1) + algo[1]=initial_coding_parameter; + if (algo[2]==-1) + algo[2]=coding; + if (algo[3]==-1) + algo[3]=coding_parameter; + return data; +} + +char DECLSPECDLLEXPORT *tng_compress_pos_find_algo(double *pos, int natoms, int nframes, + double desired_precision, + int speed, + int *algo, + int *nitems) +{ + algo[0]=-1; + algo[1]=-1; + algo[2]=-1; + algo[3]=-1; + return tng_compress_pos(pos,natoms,nframes,desired_precision,speed,algo,nitems); +} + +int DECLSPECDLLEXPORT tng_compress_nalgo(void) +{ + return 4; /* There are currently four parameters required: + + 1) The compression algorithm for the first frame (initial_coding). + 2) One parameter to the algorithm for the first frame (the initial coding parameter). + 3) The compression algorithm for the remaining frames (coding). + 4) One parameter to the algorithm for the remaining frames (the coding parameter). */ +} + +char DECLSPECDLLEXPORT *tng_compress_vel(double *vel, int natoms, int nframes, + double desired_precision, + int speed, int *algo, + int *nitems) +{ + char *data=malloc(natoms*nframes*14+11*4); /* 12 bytes are required to store 4 32 bit integers + This is 17% extra. The final 11*4 is to store information + needed for decompression. */ + int *quant=malloc(natoms*nframes*3*sizeof *quant); + int *quant_inter=malloc(natoms*nframes*3*sizeof *quant_inter); + + int initial_coding, initial_coding_parameter; + int coding, coding_parameter; + fix_t prec_hi, prec_lo; + Ptngc_d_to_i32x2(desired_precision,&prec_hi,&prec_lo); + if (speed==0) + speed=SPEED_DEFAULT; /* Set the default speed. */ + /* Boundaries of speed. */ + if (speed<1) + speed=1; + if (speed>6) + speed=6; + initial_coding=algo[0]; + initial_coding_parameter=algo[1]; + coding=algo[2]; + coding_parameter=algo[3]; + + quantize(vel,natoms,nframes,PRECISION(prec_hi,prec_lo),quant); + quant_inter_differences(quant,natoms,nframes,quant_inter); + + /* If any of the above codings / coding parameters are == -1, the optimal parameters must be found */ + if (initial_coding==-1) + { + initial_coding_parameter=-1; + determine_best_vel_initial_coding(quant,natoms,speed,prec_hi,prec_lo, + &initial_coding,&initial_coding_parameter); + } + else if (initial_coding_parameter==-1) + { + determine_best_vel_initial_coding(quant,natoms,speed,prec_hi,prec_lo, + &initial_coding,&initial_coding_parameter); + } + + if (nframes==1) + { + coding=0; + coding_parameter=0; + } + + if (nframes>1) + { + if (coding==-1) + { + coding_parameter=-1; + determine_best_vel_coding(quant,quant_inter,natoms,nframes,speed,prec_hi,prec_lo, + &coding,&coding_parameter); + } + else if (coding_parameter==-1) + { + determine_best_vel_coding(quant,quant_inter,natoms,nframes,speed,prec_hi,prec_lo, + &coding,&coding_parameter); + } + } + + compress_quantized_vel(quant,quant_inter,natoms,nframes,speed, + initial_coding,initial_coding_parameter, + coding,coding_parameter, + prec_hi,prec_lo,nitems,data); + free(quant_inter); + free(quant); + if (algo[0]==-1) + algo[0]=initial_coding; + if (algo[1]==-1) + algo[1]=initial_coding_parameter; + if (algo[2]==-1) + algo[2]=coding; + if (algo[3]==-1) + algo[3]=coding_parameter; + return data; +} + +char DECLSPECDLLEXPORT *tng_compress_vel_find_algo(double *vel, int natoms, int nframes, + double desired_precision, + int speed, + int *algo, + int *nitems) +{ + algo[0]=-1; + algo[1]=-1; + algo[2]=-1; + algo[3]=-1; + return tng_compress_vel(vel,natoms,nframes,desired_precision,speed,algo,nitems); +} + +int DECLSPECDLLEXPORT tng_compress_inquire(char *data,int *vel, int *natoms, + int *nframes, double *precision, + int *algo) +{ + int bufloc=0; + fix_t prec_hi, prec_lo; + int initial_coding, initial_coding_parameter; + int coding, coding_parameter; + int magic_int; + magic_int=(int)readbufferfix((unsigned char *)data+bufloc,4); + bufloc+=4; + if (magic_int==MAGIC_INT_POS) + *vel=0; + else if (magic_int==MAGIC_INT_VEL) + *vel=1; + else + return 1; + /* Number of atoms. */ + *natoms=(int)readbufferfix((unsigned char *)data+bufloc,4); + bufloc+=4; + /* Number of frames. */ + *nframes=(int)readbufferfix((unsigned char *)data+bufloc,4); + bufloc+=4; + /* Initial coding. */ + initial_coding=(int)readbufferfix((unsigned char *)data+bufloc,4); + bufloc+=4; + /* Initial coding parameter. */ + initial_coding_parameter=(int)readbufferfix((unsigned char *)data+bufloc,4); + bufloc+=4; + /* Coding. */ + coding=(int)readbufferfix((unsigned char *)data+bufloc,4); + bufloc+=4; + /* Coding parameter. */ + coding_parameter=(int)readbufferfix((unsigned char *)data+bufloc,4); + bufloc+=4; + /* Precision. */ + prec_lo=readbufferfix((unsigned char *)data+bufloc,4); + bufloc+=4; + prec_hi=readbufferfix((unsigned char *)data+bufloc,4); + bufloc+=4; + *precision=PRECISION(prec_hi, prec_lo); + algo[0]=initial_coding; + algo[1]=initial_coding_parameter; + algo[2]=coding; + algo[3]=coding_parameter; + return 0; +} + +static int tng_compress_uncompress_pos(char *data,double *pos) +{ + int bufloc=0; + fix_t prec_hi, prec_lo; + int length; + int natoms, nframes; + int initial_coding, initial_coding_parameter; + int coding, coding_parameter; + int *quant=NULL; + struct coder *coder=NULL; + int rval=0; + int magic_int; + /* Magic integer for positions */ + magic_int=(int)readbufferfix((unsigned char *)data+bufloc,4); + bufloc+=4; + if (magic_int!=MAGIC_INT_POS) + { + rval=1; + goto error; + } + /* Number of atoms. */ + natoms=(int)readbufferfix((unsigned char *)data+bufloc,4); + bufloc+=4; + /* Number of frames. */ + nframes=(int)readbufferfix((unsigned char *)data+bufloc,4); + bufloc+=4; + /* Initial coding. */ + initial_coding=(int)readbufferfix((unsigned char *)data+bufloc,4); + bufloc+=4; + /* Initial coding parameter. */ + initial_coding_parameter=(int)readbufferfix((unsigned char *)data+bufloc,4); + bufloc+=4; + /* Coding. */ + coding=(int)readbufferfix((unsigned char *)data+bufloc,4); + bufloc+=4; + /* Coding parameter. */ + coding_parameter=(int)readbufferfix((unsigned char *)data+bufloc,4); + bufloc+=4; + /* Precision. */ + prec_lo=readbufferfix((unsigned char *)data+bufloc,4); + bufloc+=4; + prec_hi=readbufferfix((unsigned char *)data+bufloc,4); + bufloc+=4; + /* Allocate the memory for the quantized positions */ + quant=malloc(natoms*nframes*3*sizeof *quant); + /* The data block length. */ + length=(int)readbufferfix((unsigned char *)data+bufloc,4); + bufloc+=4; + /* The initial frame */ + coder=Ptngc_coder_init(); + rval=Ptngc_unpack_array(coder,(unsigned char*)data+bufloc,quant,natoms*3, + initial_coding,initial_coding_parameter,natoms); + Ptngc_coder_deinit(coder); + if (rval) + goto error; + /* Skip past the actual data block. */ + bufloc+=length; + /* Obtain the actual positions for the initial block. */ + if ((initial_coding==TNG_COMPRESS_ALGO_POS_XTC2) || + (initial_coding==TNG_COMPRESS_ALGO_POS_TRIPLET_ONETOONE) || + (initial_coding==TNG_COMPRESS_ALGO_POS_XTC3)) + { + unquantize(pos,natoms,1,PRECISION(prec_hi,prec_lo),quant); + } + else if ((initial_coding==TNG_COMPRESS_ALGO_POS_TRIPLET_INTRA) || + (initial_coding==TNG_COMPRESS_ALGO_POS_BWLZH_INTRA)) + { + unquantize_intra_differences(pos,natoms,1,PRECISION(prec_hi,prec_lo),quant); + unquant_intra_differences_first_frame(quant,natoms); + } + /* The remaining frames. */ + if (nframes>1) + { + /* The data block length. */ + length=(int)readbufferfix((unsigned char *)data+bufloc,4); + bufloc+=4; + coder=Ptngc_coder_init(); + rval=Ptngc_unpack_array(coder,(unsigned char *)data+bufloc,quant+natoms*3,(nframes-1)*natoms*3, + coding,coding_parameter,natoms); + Ptngc_coder_deinit(coder); + if (rval) + goto error; + if ((coding==TNG_COMPRESS_ALGO_POS_STOPBIT_INTER) || + (coding==TNG_COMPRESS_ALGO_POS_TRIPLET_INTER) || + (coding==TNG_COMPRESS_ALGO_POS_BWLZH_INTER)) + { + /* This requires that the first frame is already in one-to-one format, even if intra-frame + compression was done there. Therefore the unquant_intra_differences_first_frame should be called + before to convert it correctly. */ + unquantize_inter_differences(pos,natoms,nframes,PRECISION(prec_hi,prec_lo),quant); + } + else if ((coding==TNG_COMPRESS_ALGO_POS_XTC2) || + (coding==TNG_COMPRESS_ALGO_POS_XTC3) || + (coding==TNG_COMPRESS_ALGO_POS_TRIPLET_ONETOONE)) + { + unquantize(pos+natoms*3,natoms,nframes-1,PRECISION(prec_hi,prec_lo),quant+natoms*3); + } + else if ((coding==TNG_COMPRESS_ALGO_POS_TRIPLET_INTRA) || + (coding==TNG_COMPRESS_ALGO_POS_BWLZH_INTRA)) + { + unquantize_intra_differences(pos+natoms*3,natoms,nframes-1,PRECISION(prec_hi,prec_lo),quant+natoms*3); + } + } + error: + free(quant); + return rval; +} + +static int tng_compress_uncompress_vel(char *data,double *vel) +{ + int bufloc=0; + fix_t prec_hi, prec_lo; + int length; + int natoms, nframes; + int initial_coding, initial_coding_parameter; + int coding, coding_parameter; + int *quant=NULL; + struct coder *coder=NULL; + int rval=0; + int magic_int; + /* Magic integer for velocities */ + magic_int=(int)readbufferfix((unsigned char *)data+bufloc,4); + bufloc+=4; + if (magic_int!=MAGIC_INT_VEL) + { + rval=1; + goto error; + } + /* Number of atoms. */ + natoms=(int)readbufferfix((unsigned char *)data+bufloc,4); + bufloc+=4; + /* Number of frames. */ + nframes=(int)readbufferfix((unsigned char *)data+bufloc,4); + bufloc+=4; + /* Initial coding. */ + initial_coding=(int)readbufferfix((unsigned char *)data+bufloc,4); + bufloc+=4; + /* Initial coding parameter. */ + initial_coding_parameter=(int)readbufferfix((unsigned char *)data+bufloc,4); + bufloc+=4; + /* Coding. */ + coding=(int)readbufferfix((unsigned char *)data+bufloc,4); + bufloc+=4; + /* Coding parameter. */ + coding_parameter=(int)readbufferfix((unsigned char *)data+bufloc,4); + bufloc+=4; + /* Precision. */ + prec_lo=readbufferfix((unsigned char *)data+bufloc,4); + bufloc+=4; + prec_hi=readbufferfix((unsigned char *)data+bufloc,4); + bufloc+=4; + /* Allocate the memory for the quantized positions */ + quant=malloc(natoms*nframes*3*sizeof *quant); + /* The data block length. */ + length=(int)readbufferfix((unsigned char *)data+bufloc,4); + bufloc+=4; + /* The initial frame */ + coder=Ptngc_coder_init(); + rval=Ptngc_unpack_array(coder,(unsigned char*)data+bufloc,quant,natoms*3, + initial_coding,initial_coding_parameter,natoms); + Ptngc_coder_deinit(coder); + if (rval) + goto error; + /* Skip past the actual data block. */ + bufloc+=length; + /* Obtain the actual positions for the initial block. */ + if ((initial_coding==TNG_COMPRESS_ALGO_VEL_STOPBIT_ONETOONE) || + (initial_coding==TNG_COMPRESS_ALGO_VEL_TRIPLET_ONETOONE) || + (initial_coding==TNG_COMPRESS_ALGO_VEL_BWLZH_ONETOONE)) + { + unquantize(vel,natoms,1,PRECISION(prec_hi,prec_lo),quant); + } + /* The remaining frames. */ + if (nframes>1) + { + /* The data block length. */ + length=(int)readbufferfix((unsigned char *)data+bufloc,4); + bufloc+=4; + coder=Ptngc_coder_init(); + rval=Ptngc_unpack_array(coder,(unsigned char *)data+bufloc,quant+natoms*3,(nframes-1)*natoms*3, + coding,coding_parameter,natoms); + Ptngc_coder_deinit(coder); + if (rval) + goto error; + /* Inter-frame compression? */ + if ((coding==TNG_COMPRESS_ALGO_VEL_TRIPLET_INTER) || + (coding==TNG_COMPRESS_ALGO_VEL_STOPBIT_INTER) || + (coding==TNG_COMPRESS_ALGO_VEL_BWLZH_INTER)) + { + /* This requires that the first frame is already in one-to-one format. */ + unquantize_inter_differences(vel,natoms,nframes,PRECISION(prec_hi,prec_lo),quant); + } + /* One-to-one compression? */ + else if ((coding==TNG_COMPRESS_ALGO_VEL_STOPBIT_ONETOONE) || + (coding==TNG_COMPRESS_ALGO_VEL_TRIPLET_ONETOONE) || + (coding==TNG_COMPRESS_ALGO_VEL_BWLZH_ONETOONE)) + { + unquantize(vel+natoms*3,natoms,nframes-1,PRECISION(prec_hi,prec_lo),quant+natoms*3); + } + } + error: + free(quant); + return rval; +} + +/* Uncompresses any tng compress block, positions or velocities. It determines whether it is positions or velocities from the data buffer. The return value is 0 if ok, and 1 if not. +*/ +int DECLSPECDLLEXPORT tng_compress_uncompress(char *data,double *posvel) +{ + int magic_int; + magic_int=(int)readbufferfix((unsigned char *)data,4); + if (magic_int==MAGIC_INT_POS) + return tng_compress_uncompress_pos(data,posvel); + else if (magic_int==MAGIC_INT_VEL) + return tng_compress_uncompress_vel(data,posvel); + else + return 1; +} + +static char *compress_algo_pos[TNG_COMPRESS_ALGO_MAX]={ + "Positions invalid algorithm", + "Positions stopbits interframe", + "Positions triplet interframe", + "Positions triplet intraframe", + "Positions invalid algorithm", + "Positions XTC2", + "Positions invalid algorithm", + "Positions triplet one to one", + "Positions BWLZH interframe", + "Positions BWLZH intraframe", + "Positions XTC3" +}; + +static char *compress_algo_vel[TNG_COMPRESS_ALGO_MAX]={ + "Velocities invalid algorithm", + "Velocities stopbits one to one", + "Velocities triplet interframe", + "Velocities triplet one to one", + "Velocities invalid algorithm", + "Velocities invalid algorithm", + "Velocities stopbits interframe", + "Velocities invalid algorithm", + "Velocities BWLZH interframe", + "Velocities BWLZH one to one", + "Velocities invalid algorithm" +}; + +char DECLSPECDLLEXPORT *tng_compress_initial_pos_algo(int *algo) +{ + int i=algo[0]; + if (i<0) + i=0; + if (i>=TNG_COMPRESS_ALGO_MAX) + i=0; + return compress_algo_pos[i]; +} + +char DECLSPECDLLEXPORT *tng_compress_pos_algo(int *algo) +{ + int i=algo[2]; + if (i<0) + i=0; + if (i>=TNG_COMPRESS_ALGO_MAX) + i=0; + return compress_algo_pos[i]; +} + +char DECLSPECDLLEXPORT *tng_compress_initial_vel_algo(int *algo) +{ + int i=algo[0]; + if (i<0) + i=0; + if (i>=TNG_COMPRESS_ALGO_MAX) + i=0; + return compress_algo_vel[i]; +} + +char DECLSPECDLLEXPORT *tng_compress_vel_algo(int *algo) +{ + int i=algo[2]; + if (i<0) + i=0; + if (i>=TNG_COMPRESS_ALGO_MAX) + i=0; + return compress_algo_vel[i]; +} diff --git a/src/compression/vals16.c b/src/compression/vals16.c new file mode 100644 index 0000000..6a8e7e7 --- /dev/null +++ b/src/compression/vals16.c @@ -0,0 +1,73 @@ +/* This code is part of the tng compression routines. + * + * Written by Daniel Spangberg + * Copyright (c) 2010, 2013, The GROMACS development team. + * + * + * This program is free software; you can redistribute it and/or + * modify it under the terms of the GNU Lesser General Public License + * as published by the Free Software Foundation; either version 2.1 + * of the License, or (at your option) any later version. + */ + + +#include "vals16.h" + +/* Coding 32 bit ints in sequences of 16 bit ints. Worst case + the output is 3*nvals long. */ +void Ptngc_comp_conv_to_vals16(unsigned int *vals,int nvals, + unsigned int *vals16, int *nvals16) +{ + int i; + int j=0; + for (i=0; i>15; + vals16[j++]=lo; + if (hi<=0x7FFFU) + vals16[j++]=hi; + else + { + unsigned int lohi=(hi&0x7FFFU)|0x8000U; + unsigned int hihi=hi>>15; + vals16[j++]=lohi; + vals16[j++]=hihi; + } + } + } +#if 0 + /* Test that things that detect that this is bad really works. */ + vals16[0]=0; +#endif + *nvals16=j; +} + +void Ptngc_comp_conv_from_vals16(unsigned int *vals16,int nvals16, + unsigned int *vals, int *nvals) +{ + int i=0; + int j=0; + while (i +#endif + + +#include +#include +#include "tng_compress.h" +#include "warnmalloc.h" + +void DECLSPECDLLEXPORT *Ptngc_warnmalloc_x(size_t size, char *file, int line) +{ + void *mem=malloc(size); + if (!mem) + { + fprintf(stderr,"TRAJNG ERROR: Could not allocate memory of size %lu at %s:%d\n",(unsigned long) size,file,line); + exit(EXIT_FAILURE); + } + return mem; +} + +void DECLSPECDLLEXPORT *Ptngc_warnrealloc_x(void *old, size_t size, char *file, int line) +{ + void *mem=realloc(old,size); + if (!mem) + { + fprintf(stderr,"TRAJNG ERROR: Could not allocate memory of size %lu at %s:%d\n",(unsigned long) size,file,line); + exit(EXIT_FAILURE); + } + return mem; +} diff --git a/src/compression/widemuldiv.c b/src/compression/widemuldiv.c new file mode 100644 index 0000000..98d47cd --- /dev/null +++ b/src/compression/widemuldiv.c @@ -0,0 +1,236 @@ +/* This code is part of the tng compression routines. + * + * Written by Daniel Spangberg + * Copyright (c) 2010, 2013, The GROMACS development team. + * + * + * This program is free software; you can redistribute it and/or + * modify it under the terms of the GNU Lesser General Public License + * as published by the Free Software Foundation; either version 2.1 + * of the License, or (at your option) any later version. + */ + + +#if HAVE_CONFIG_H +#include +#endif + + +#include +#include + +#include "tng_compress.h" + +/* 64 bit integers are not required in this part of the program. But + they improve the running speed if present. If 64 bit integers are + available define the symbol HAVE64BIT. It should get automatically + defined by the defines in my64bit.h */ +#include "my64bit.h" + +#include "widemuldiv.h" + +#ifndef TRAJNG_X86_GCC_INLINE_MULDIV +#if defined(__GNUC__) && defined(__i386__) +#define TRAJNG_X86_GCC_INLINE_MULDIV +#endif /* gcc & i386 */ +#if defined(__GNUC__) && defined(__x86_64__) +#define TRAJNG_X86_GCC_INLINE_MULDIV +#endif /* gcc & x86_64 */ +#endif /* TRAJNG X86 GCC INLINE MULDIV */ + +/* Multiply two 32 bit unsigned integers returning a 64 bit unsigned value (in two integers) */ +void Ptngc_widemul(unsigned int i1, unsigned int i2, unsigned int *ohi, unsigned int *olo) +{ +#if defined(TRAJNG_X86_GCC_INLINE_MULDIV) + __asm__ __volatile__ ("mull %%edx\n\t" + : "=a" (i1), "=d" (i2) + : "a" (i1),"d" (i2) + : "cc"); + *ohi=i2; + *olo=i1; +#else /* TRAJNG X86 GCC INLINE MULDIV */ + +#ifdef HAVE64BIT + my_uint64_t res= ((my_uint64_t)i1) * ((my_uint64_t)i2); + *olo=res & 0xFFFFFFFFU; + *ohi=(res>>32) & 0xFFFFFFFFU; +#else /* HAVE64BIT */ + + unsigned int bits=16; + unsigned int L_m=(1<>bits; + b_L=i2 & L_m; + b_U=i2>>bits; + + /* Do a*a=>2a multiply where a is half number of bits in an int */ + x=a_L*b_L; + x_LL=x & L_m; + x_LU=x>>bits; + + x=a_U*b_L; + x_LU+=x & L_m; + x_UL=x>>bits; + + x=a_L*b_U; + x_LU+=x & L_m; + x_UL+=x>>bits; + + x=a_U*b_U; + x_UL+=x & L_m; + x_UU=x>>bits; + + /* Combine results */ + x_UL+=x_LU>>bits; + x_LU&=L_m; + x_UU+=x_UL>>bits; + x_UL&=L_m; + + x_U=(x_UU<>(bits2-hibit); + s_L=(i<lo) + { + unsigned int t; + hi--; /* Borrow */ + t=allbits-s_L; + lo+=t+1; + } + else + lo-=s_L; + + /* Set bit. */ + res|=rmask; + } + rmask>>=1; + s_L>>=1; + if (s_U & 1U) + s_L|=hibit_mask; + s_U>>=1; + } + *remainder=lo; + *result=res; +#endif /* HAVE64BIT */ +#endif /* TRAJNG X86 GCC INLINE MULDIV */ +} + +/* Add a unsigned int to a largeint. j determines which value in the + largeint to add v1 to. */ +static void largeint_add_gen(unsigned int v1, unsigned int *largeint, int n, int j) +{ + /* Add with carry. unsigned ints in C wrap modulo 2**bits when "overflowed". */ + unsigned int v2=(v1+largeint[j])&0xFFFFFFFFU; /* Add and cap at 32 bits */ + unsigned int carry=0; + if ((((unsigned int)-1)&0xFFFFFFFFU) -v164 mul */ + largeint_add_gen(lo,largeint_out,n,i); + if (i+1 +#endif + +#include +#include +#include +#include +#include "coder.h" +#include "widemuldiv.h" +#include "warnmalloc.h" + +/* Generated by gen_magic.py */ +#define MAX_MAGIC 92 + +static unsigned int magic[MAX_MAGIC]={ +2U, 3U, 4U, 5U, +6U, 8U, 10U, 12U, +16U, 20U, 25U, 32U, +40U, 50U, 64U, 80U, +101U, 128U, 161U, 203U, +256U, 322U, 406U, 512U, +645U, 812U, 1024U, 1290U, +1625U, 2048U, 2580U, 3250U, +4096U, 5160U, 6501U, 8192U, +10321U, 13003U, 16384U, 20642U, +26007U, 32768U, 41285U, 52015U, +65536U, 82570U, 104031U, 131072U, +165140U, 208063U, 262144U, 330280U, +416127U, 524288U, 660561U, 832255U, +1048576U, 1321122U, 1664510U, 2097152U, +2642245U, 3329021U, 4194304U, 5284491U, +6658042U, 8388608U, 10568983U, 13316085U, +16777216U, 21137967U, 26632170U, 33554432U, +42275935U, 53264340U, 67108864U, 84551870U, +106528681U, 134217728U, 169103740U, 213057362U, +268435456U, 338207481U, 426114725U, 536870912U, +676414963U, 852229450U, 1073741824U, 1352829926U, +1704458900U, 2147483648U, 2705659852U, 3408917801U, +}; + +static unsigned int magic_bits[MAX_MAGIC][8]={ +{ 3, 6, 9, 12, 15, 18, 21, 24, }, +{ 5, 10, 15, 20, 24, 29, 34, 39, }, +{ 6, 12, 18, 24, 30, 36, 42, 48, }, +{ 7, 14, 21, 28, 35, 42, 49, 56, }, +{ 8, 16, 24, 32, 39, 47, 55, 63, }, +{ 9, 18, 27, 36, 45, 54, 63, 72, }, +{ 10, 20, 30, 40, 50, 60, 70, 80, }, +{ 11, 22, 33, 44, 54, 65, 76, 87, }, +{ 12, 24, 36, 48, 60, 72, 84, 97, }, +{ 13, 26, 39, 52, 65, 78, 91, 104, }, +{ 14, 28, 42, 56, 70, 84, 98, 112, }, +{ 15, 30, 45, 60, 75, 90, 105, 120, }, +{ 16, 32, 48, 64, 80, 96, 112, 128, }, +{ 17, 34, 51, 68, 85, 102, 119, 136, }, +{ 18, 36, 54, 72, 90, 108, 127, 144, }, +{ 19, 38, 57, 76, 95, 114, 133, 152, }, +{ 20, 40, 60, 80, 100, 120, 140, 160, }, +{ 21, 42, 63, 84, 105, 127, 147, 168, }, +{ 22, 44, 66, 88, 110, 132, 154, 176, }, +{ 23, 46, 69, 92, 115, 138, 161, 184, }, +{ 24, 48, 72, 97, 120, 144, 168, 192, }, +{ 25, 50, 75, 100, 125, 150, 175, 200, }, +{ 26, 52, 78, 104, 130, 156, 182, 208, }, +{ 27, 54, 81, 108, 135, 162, 190, 216, }, +{ 28, 56, 84, 112, 140, 168, 196, 224, }, +{ 29, 58, 87, 116, 145, 174, 203, 232, }, +{ 30, 60, 90, 120, 150, 180, 210, 240, }, +{ 31, 62, 93, 124, 155, 186, 217, 248, }, +{ 32, 64, 96, 128, 160, 192, 224, 256, }, +{ 33, 66, 99, 132, 165, 198, 231, 264, }, +{ 34, 68, 102, 136, 170, 204, 238, 272, }, +{ 35, 70, 105, 140, 175, 210, 245, 280, }, +{ 36, 72, 108, 144, 180, 216, 252, 288, }, +{ 37, 74, 111, 148, 185, 222, 259, 296, }, +{ 38, 76, 114, 152, 190, 228, 266, 304, }, +{ 39, 78, 117, 157, 195, 234, 273, 312, }, +{ 40, 80, 120, 160, 200, 240, 280, 320, }, +{ 41, 82, 123, 164, 205, 246, 287, 328, }, +{ 42, 84, 127, 168, 210, 252, 294, 336, }, +{ 43, 86, 129, 172, 215, 258, 301, 344, }, +{ 44, 88, 132, 176, 220, 264, 308, 352, }, +{ 45, 90, 135, 180, 225, 270, 315, 360, }, +{ 46, 92, 138, 184, 230, 276, 322, 368, }, +{ 47, 94, 141, 188, 235, 282, 329, 376, }, +{ 48, 97, 144, 192, 240, 288, 336, 384, }, +{ 49, 98, 147, 196, 245, 294, 343, 392, }, +{ 50, 100, 150, 200, 250, 300, 350, 400, }, +{ 52, 102, 153, 204, 255, 306, 357, 408, }, +{ 52, 104, 156, 208, 260, 312, 364, 416, }, +{ 53, 106, 159, 212, 265, 318, 371, 424, }, +{ 54, 108, 162, 216, 270, 324, 378, 432, }, +{ 55, 110, 165, 220, 275, 330, 385, 440, }, +{ 56, 112, 168, 224, 280, 336, 392, 448, }, +{ 57, 114, 172, 228, 285, 342, 399, 456, }, +{ 58, 116, 174, 232, 290, 348, 406, 464, }, +{ 59, 118, 177, 236, 295, 354, 413, 472, }, +{ 60, 120, 180, 240, 300, 360, 420, 480, }, +{ 61, 122, 183, 244, 305, 366, 427, 488, }, +{ 62, 124, 186, 248, 310, 372, 434, 496, }, +{ 63, 127, 190, 252, 315, 378, 442, 505, }, +{ 64, 128, 192, 256, 320, 384, 448, 512, }, +{ 65, 130, 195, 260, 325, 390, 455, 520, }, +{ 66, 132, 198, 264, 330, 396, 462, 528, }, +{ 67, 134, 201, 268, 335, 402, 469, 536, }, +{ 68, 136, 204, 272, 340, 408, 476, 544, }, +{ 69, 138, 207, 276, 345, 414, 483, 552, }, +{ 70, 140, 210, 280, 350, 420, 490, 560, }, +{ 71, 142, 213, 284, 355, 426, 497, 568, }, +{ 72, 144, 216, 288, 360, 432, 505, 576, }, +{ 73, 146, 219, 292, 365, 438, 511, 584, }, +{ 74, 148, 222, 296, 370, 444, 518, 592, }, +{ 75, 150, 225, 300, 375, 451, 525, 600, }, +{ 76, 152, 228, 304, 380, 456, 532, 608, }, +{ 77, 154, 231, 308, 385, 462, 539, 616, }, +{ 78, 157, 234, 312, 390, 469, 546, 625, }, +{ 79, 158, 237, 316, 395, 474, 553, 632, }, +{ 80, 160, 240, 320, 400, 480, 560, 640, }, +{ 81, 162, 243, 324, 406, 486, 568, 648, }, +{ 82, 164, 246, 328, 410, 492, 574, 656, }, +{ 83, 166, 249, 332, 415, 498, 581, 664, }, +{ 84, 168, 252, 336, 420, 505, 588, 672, }, +{ 85, 170, 255, 340, 425, 510, 595, 680, }, +{ 86, 172, 258, 344, 430, 516, 602, 688, }, +{ 87, 174, 261, 348, 435, 522, 609, 696, }, +{ 88, 176, 264, 352, 440, 528, 616, 704, }, +{ 89, 178, 267, 356, 445, 534, 623, 712, }, +{ 90, 180, 270, 360, 451, 540, 631, 720, }, +{ 91, 182, 273, 364, 455, 546, 637, 728, }, +{ 92, 184, 276, 368, 460, 552, 644, 736, }, +{ 94, 187, 279, 373, 466, 558, 651, 745, }, +{ 94, 188, 282, 376, 470, 564, 658, 752, }, +{ 95, 190, 285, 380, 475, 570, 665, 760, }, +}; + + +static const double iflipgaincheck=0.89089871814033927; /* 1./(2**(1./6)) */ + + +/* Difference in indices used for determining whether to store as large or small */ +#define QUITE_LARGE 3 +#define IS_LARGE 6 + +#if 0 +#define SHOWIT +#endif + +int Ptngc_magic(unsigned int i) +{ + return magic[i]; +} + +int Ptngc_find_magic_index(unsigned int maxval) +{ + int i=0; + while (magic[i]<=maxval) + i++; + return i; +} + +static unsigned int positive_int(int item) +{ + int s=0; + if (item>0) + s=1+(item-1)*2; + else if (item<0) + s=2+(-item-1)*2; + return s; +} + +static int unpositive_int(int val) +{ + int s=(val+1)/2; + if ((val%2)==0) + s=-s; + return s; +} + + +/* Sequence instructions */ +#define INSTR_DEFAULT 0 +#define INSTR_BASE_RUNLENGTH 1 +#define INSTR_ONLY_LARGE 2 +#define INSTR_ONLY_SMALL 3 +#define INSTR_LARGE_BASE_CHANGE 4 +#define INSTR_FLIP 5 +#define INSTR_LARGE_RLE 6 + +#define MAXINSTR 7 + +#ifdef SHOWIT +static char *instrnames[MAXINSTR]={ + "large+small", + "base+run", + "large", + "small", + "large base change", + "flip", + "large rle", +}; +#endif + +/* Bit patterns in the compressed code stream: */ + +static const int seq_instr[MAXINSTR][2]= + { + { 1,1 }, /* 1 : one large atom + runlength encoded small integers. Use same settings as before. */ + { 0,2 }, /* 00 : set base and runlength in next four bits (x). base (increase/keep/decrease)=x%3-1. runlength=1+x/3. + The special value 1111 in the four bits means runlength=6 and base change=0 */ + { 4,4 }, /* 0100 : next only a large atom comes. */ + { 5,4 }, /* 0101 : next only runlength encoded small integers. Use same settings as before. */ + { 6,4 }, /* 0110 : Large change in base. Change is encoded in the + following 2 bits. change direction (sign) is the first + bit. The next bit +1 is the actual change. This + allows the change of up to +/- 2 indices. */ + { 14,5 }, /* 01110 : flip whether integers should be modified to compress water better */ + { 15,5 }, /* 01111 : Large rle. The next 4 bits encode how many + large atoms are in the following sequence: 3-18. (2 is + more efficiently coded with two large instructions. */ + }; + +static void write_instruction(struct coder *coder,int instr,unsigned char **output_ptr) +{ + Ptngc_writebits(coder,seq_instr[instr][0],seq_instr[instr][1],output_ptr); +#ifdef SHOWIT + fprintf(stderr,"INSTR: %s (%d bits)\n",instrnames[instr],seq_instr[instr][1]); +#endif +} + +static unsigned int readbits(unsigned char **ptr, int *bitptr, int nbits) +{ + unsigned int val=0U; + unsigned int extract_mask=0x80U>>*bitptr; + unsigned char thisval=**ptr; +#ifdef SHOWIT + fprintf(stderr,"Read nbits=%d val=",nbits); +#endif + while (nbits--) + { + val<<=1; + val|=((extract_mask & thisval)!=0); + *bitptr=(*bitptr)+1; + extract_mask>>=1; + if (!extract_mask) + { + extract_mask=0x80U; + *ptr=(*ptr)+1; + *bitptr=0; + if (nbits) + thisval=**ptr; + } + } +#ifdef SHOWIT + fprintf(stderr,"%x\n",val); +#endif + return val; +} + +static void readmanybits(unsigned char **ptr, int *bitptr, int nbits, unsigned char *buffer) +{ + while (nbits>=8) + { + *buffer++=readbits(ptr,bitptr,8); + nbits-=8; +#ifdef SHOWIT + fprintf(stderr,"Read value %02x\n",buffer[-1]); +#endif + } + if (nbits) + { + *buffer++=readbits(ptr,bitptr,nbits); +#ifdef SHOWIT + fprintf(stderr,"Read value %02x\n",buffer[-1]); +#endif + } +} + +static int read_instruction(unsigned char **ptr, int *bitptr) +{ + int instr=-1; + unsigned int bits=readbits(ptr,bitptr,1); + if (bits) + instr=INSTR_DEFAULT; + else + { + bits=readbits(ptr,bitptr,1); + if (!bits) + instr=INSTR_BASE_RUNLENGTH; + else + { + bits=readbits(ptr,bitptr,2); + if (bits==0) + instr=INSTR_ONLY_LARGE; + else if (bits==1) + instr=INSTR_ONLY_SMALL; + else if (bits==2) + instr=INSTR_LARGE_BASE_CHANGE; + else if (bits==3) + { + bits=readbits(ptr,bitptr,1); + if (bits==0) + instr=INSTR_FLIP; + else + instr=INSTR_LARGE_RLE; + } + } + } + return instr; +} + +/* Modifies three integer values for better compression of water */ +static void swap_ints(int *in, int *out) +{ + out[0]=in[0]+in[1]; + out[1]=-in[1]; + out[2]=in[1]+in[2]; +} + +static void swap_is_better(int *input, int *minint, int *sum_normal, int *sum_swapped) +{ + int normal_max=0; + int swapped_max=0; + int i,j; + int normal[3]; + int swapped[3]; + for (i=0; i<3; i++) + { + normal[0]=input[i]-minint[i]; + normal[1]=input[3+i]-input[i]; /* minint[i]-minint[i] cancels out */ + normal[2]=input[6+i]-input[3+i]; /* minint[i]-minint[i] cancels out */ + swap_ints(normal,swapped); + for (j=1; j<3; j++) + { + if (positive_int(normal[j])>normal_max) + normal_max=positive_int(normal[j]); + if (positive_int(swapped[j])>swapped_max) + swapped_max=positive_int(swapped[j]); + } + } + if (normal_max==0) + normal_max=1; + if (swapped_max==0) + swapped_max=1; + *sum_normal=normal_max; + *sum_swapped=swapped_max; +} + +static void swapdecide(struct coder *coder, int *input,int *swapatoms, int *large_index, int *minint, unsigned char **output_ptr) +{ + int didswap=0; + int normal,swapped; + swap_is_better(input,minint,&normal,&swapped); + /* We have to determine if it is worth to change the behaviour. + If diff is positive it means that it is worth something to + swap. But it costs 4 bits to do the change. If we assume that + we gain 0.17 bit by the swap per value, and the runlength>2 + for four molecules in a row, we gain something. So check if we + gain at least 0.17 bits to even attempt the swap. + */ +#ifdef SHOWIT + fprintf(stderr,"Trying Flip: %g %g\n",(double)swapped/normal, (double)normal/swapped); +#endif + if (((swapped>shift) & 0xFFU; + shift+=8; + } + } +} + +/* The opposite of base_compress. */ +static void trajcoder_base_decompress(unsigned char *input, int n, int *index, int *output) +{ + unsigned int largeint[19]; + unsigned int largeint_tmp[19]; + int i,j; + /* Convert the sequence of bytes to a largeint. */ + for (i=0; i<18; i++) + { + int shift=0; + largeint[i]=0U; + for (j=0; j<4; j++) + { + largeint[i]|=((unsigned int)input[i*4+j])<=0; i--) + { + unsigned int remainder=Ptngc_largeint_div(magic[index[i%3]],largeint,largeint_tmp,19); +#if 0 +#ifdef SHOWIT + fprintf(stderr,"Remainder: %u\n",remainder); +#endif +#endif +#if 0 + for (j=0; j<19; j++) + largeint[j]=largeint_tmp[j]; +#endif + memcpy(largeint,largeint_tmp,19*sizeof *largeint); + output[i]=remainder; + } +} + +/* It is "large" if we have to increase the small index quite a + bit. Not so much to be rejected by the not very large check + later. */ +static int is_quite_large(int *input, int small_index, int max_large_index) +{ + int is=0; + int i; + if (small_index+QUITE_LARGE>=max_large_index) + is=1; + else + { + for (i=0; i<3; i++) + if (positive_int(input[i])>magic[small_index+QUITE_LARGE]) + { + is=1; + break; + } + } + return is; +} + +#ifdef SHOWIT +int nbits_sum; +int nvalues_sum; +#endif + +static void write_three_large(struct coder *coder, int *encode_ints, int *large_index, int nbits, unsigned char *compress_buffer, unsigned char **output_ptr) +{ + trajcoder_base_compress(encode_ints,3,large_index,compress_buffer); + Ptngc_writemanybits(coder,compress_buffer,nbits,output_ptr); +#ifdef SHOWIT + fprintf(stderr,"nbits=%d (%g)\n",nbits,nbits/3.); + nbits_sum+=nbits; + nvalues_sum+=3; + fprintf(stderr,"Large: %d %d %d\n",encode_ints[0],encode_ints[1],encode_ints[2]); +#endif +} + +static void insert_batch(int *input_ptr, int ntriplets_left, int *prevcoord,int *minint, int *encode_ints, int startenc, int *nenc) +{ + int nencode=startenc*3; + int tmp_prevcoord[3]; + + tmp_prevcoord[0]=prevcoord[0]; + tmp_prevcoord[1]=prevcoord[1]; + tmp_prevcoord[2]=prevcoord[2]; + + if (startenc) + { + int i; + for (i=0; imaxint[j]) + maxint[j]=input[i*3+j]; + if (input[i*3+j]max_large_index) + max_large_index=large_index[1]; + if (large_index[2]>max_large_index) + max_large_index=large_index[2]; + +#ifdef SHOWIT + for (j=0; j<3; j++) + fprintf(stderr,"minint[%d]=%d. maxint[%d]=%d large_index[%d]=%d value=%d\n",j,minint[j],j,maxint[j], + j,large_index[j],magic[large_index[j]]); + fprintf(stderr,"large_nbits=%d\n",large_nbits); +#endif + + + /* Guess initial small index */ + small_index=max_large_index/2; + + /* Find the largest value that is not large. Not large is half index of + large. */ + max_small=magic[small_index]; + intmax=0; + for (i=0; i<*length; i++) + { + int item=input[i]; + int s=positive_int(item); + if (s>intmax) + if (spack_temporary_bits=32; + coder->pack_temporary=positive_int(minint[0]); + Ptngc_out8bits(coder,&output_ptr); + coder->pack_temporary_bits=32; + coder->pack_temporary=positive_int(minint[1]); + Ptngc_out8bits(coder,&output_ptr); + coder->pack_temporary_bits=32; + coder->pack_temporary=positive_int(minint[2]); + Ptngc_out8bits(coder,&output_ptr); + /* Store max indices */ + coder->pack_temporary_bits=8; + coder->pack_temporary=large_index[0]; + Ptngc_out8bits(coder,&output_ptr); + coder->pack_temporary_bits=8; + coder->pack_temporary=large_index[1]; + Ptngc_out8bits(coder,&output_ptr); + coder->pack_temporary_bits=8; + coder->pack_temporary=large_index[2]; + Ptngc_out8bits(coder,&output_ptr); + /* Store initial small index */ + coder->pack_temporary_bits=8; + coder->pack_temporary=small_index; + Ptngc_out8bits(coder,&output_ptr); + +#if 0 +#ifdef SHOWIT + for (i=0; ilargest_required_base) + largest_required_base=encode_ints[ienc]; + /* Also compute what the largest base is for the current runlength setting! */ + largest_runlength_base=0; + for (ienc=0; (ienclargest_runlength_base) + largest_runlength_base=encode_ints[ienc]; + + largest_required_index=Ptngc_find_magic_index(largest_required_base); + largest_runlength_index=Ptngc_find_magic_index(largest_runlength_base); + + if (largest_required_indexntriplets_left) + new_runlength=ntriplets_left; + + /* We must at least try to get some small integers going. */ + if (new_runlength==0) + { + new_runlength=1; + new_small_index=small_index; + } + + iter_runlength=new_runlength; + iter_small_index=new_small_index; + + /* Iterate to find optimal encoding and runlength */ +#ifdef SHOWIT + fprintf(stderr,"Entering iterative loop.\n"); + fflush(stderr); +#endif + + do { + new_runlength=iter_runlength; + new_small_index=iter_small_index; + +#ifdef SHOWIT + fprintf(stderr,"Test new_small_index=%d Base=%d\n",new_small_index,magic[new_small_index]); +#endif + /* What is the largest runlength + we can do with the currently + selected encoding? Also the max supported runlength is 6 triplets! */ + for (ienc=0; iencnew_small_index) + break; + } + if (ienc/3>new_runlength) + { + iter_runlength=ienc/3; +#ifdef SHOWIT + fprintf(stderr,"I found a new possible runlength: %d\n",iter_runlength); +#endif + } + + /* How large encoding do we have to use? */ + largest_runlength_base=0; + for (ienc=0; ienclargest_runlength_base) + largest_runlength_base=encode_ints[ienc]; + largest_runlength_index=Ptngc_find_magic_index(largest_runlength_base); + if (largest_runlength_index!=new_small_index) + { + iter_small_index=largest_runlength_index; +#ifdef SHOWIT + fprintf(stderr,"I found a new possible small index: %d Base=%d\n",iter_small_index,magic[iter_small_index]); +#endif + } + } while ((new_runlength!=iter_runlength) || + (new_small_index!=iter_small_index)); + +#ifdef SHOWIT + fprintf(stderr,"Exit iterative loop.\n"); + fflush(stderr); +#endif + + /* Verify that we got something good. We may have caught a + substantially larger atom. If so we should just bail + out and let the loop get on another lap. We may have a + minimum runlength though and then we have to fulfill + the request to write out these atoms! */ + rle_index_dep=0; + if (new_runlength<3) + rle_index_dep=IS_LARGE; + else if (new_runlength<6) + rle_index_dep=QUITE_LARGE; + if ((min_runlength) + || ((new_small_index=%g?\n",change,isum,(double)magic[small_index+change]*(double)magic[small_index+change]); +#endif + if (isum>(double)magic[small_index+change]*(double)magic[small_index+change]) + { +#ifdef SHOWIT + fprintf(stderr,"Rejected decrease %d of index due to length of vector: %g>=%g\n",change,isum,(double)magic[small_index+change]*(double)magic[small_index+change]); +#endif + rejected=1; + change++; + } + } while ((change<0) && (rejected)); + if (change==0) + break; + } + } + + /* If the only thing to do is to change the base by + only one -1 it is probably not worth it. */ + if (!((change==-1) && (runlength==new_runlength))) + { + /* If we have a very short runlength we do not + want to do large base changes. It costs 6 + extra bits to do -2. We gain 2/3 + bits per value to decrease the index by -2, + ie 2 bits, so to any changes down we must + have a runlength of 3 or more to do it for + one molecule! If we have several molecules we + will gain of course, so not be so strict. */ + if ((change==-2) && (new_runlength<3)) + { + if (runlength==new_runlength) + change=0; + else + change=-1; +#ifdef SHOWIT + fprintf(stderr,"Rejected change by -2 due to too short runlenght. Change set to %d\n",change); +#endif + } + + /* First adjust base using large base change instruction (if necessary) */ + while ((change>1) || (change<-1) || ((new_runlength==6) && (change))) + { + unsigned int code=0U; + int this_change=change; + if (this_change>2) + this_change=2; + if (this_change<-2) + this_change=-2; + change-=this_change; +#ifdef SHOWIT + fprintf(stderr,"Large base change: %d.\n",this_change); +#endif + small_index+=this_change; + if (this_change<0) + { + code|=2U; + this_change=-this_change; + } + code|=(unsigned int)(this_change-1); + write_instruction(coder,INSTR_LARGE_BASE_CHANGE,&output_ptr); + Ptngc_writebits(coder,code,2,&output_ptr); + } + /* If there is still base change or runlength changes to do, we do them now. */ + if ((new_runlength!=runlength) || (change)) + { + unsigned int ichange=(unsigned int)(change+1); + unsigned int code=0U; + unsigned int irun=(unsigned int)((new_runlength-1)*3); + if (new_runlength==6) + ichange=0; /* Means no change. The change has been taken care of explicitly by a large + base change instruction above. */ + code=ichange+irun; +#ifdef SHOWIT + fprintf(stderr,"Small base change: %d Runlength change: %d\n",change,new_runlength); +#endif + small_index+=change; + write_instruction(coder,INSTR_BASE_RUNLENGTH,&output_ptr); + Ptngc_writebits(coder,code,4,&output_ptr); + runlength=new_runlength; + } + } +#ifdef SHOWIT + else + fprintf(stderr,"Rejected base change due to only change==-1\n"); +#endif +#ifdef SHOWIT + fprintf(stderr,"Current small index: %d Base=%d\n",small_index,magic[small_index]); +#endif + } + /* If we have a large previous integer we can combine it with a sequence of small ints. */ + if (has_large) + { + /* If swapatoms is set to 1 but we did actually not + do any swapping, we must first write out the + large atom and then the small. If swapatoms is 1 + and we did swapping we can use the efficient + encoding. */ + if ((swapatoms) && (!didswap)) + { +#ifdef SHOWIT + fprintf(stderr,"Swapatoms was set to 1 but we did not do swapping!\n"); + fprintf(stderr,"Only one large integer.\n"); +#endif + /* Flush all large atoms. */ + flush_large(coder,&has_large,has_large_ints,has_large,large_index,large_nbits,compress_buffer,&output_ptr); +#ifdef SHOWIT + fprintf(stderr,"Sequence of only small integers.\n"); +#endif + write_instruction(coder,INSTR_ONLY_SMALL,&output_ptr); + } + else + { + +#ifdef SHOWIT + fprintf(stderr,"Sequence of one large and small integers (good compression).\n"); +#endif + /* Flush all large atoms but one! */ + if (has_large>1) + flush_large(coder,&has_large,has_large_ints,has_large-1,large_index,large_nbits,compress_buffer,&output_ptr); + write_instruction(coder,INSTR_DEFAULT,&output_ptr); + write_three_large(coder,has_large_ints,large_index,large_nbits,compress_buffer,&output_ptr); + has_large=0; + } + } + else + { +#ifdef SHOWIT + fprintf(stderr,"Sequence of only small integers.\n"); +#endif + write_instruction(coder,INSTR_ONLY_SMALL,&output_ptr); + } + /* Base compress small integers using the current parameters. */ + nbits=magic_bits[small_index][runlength-1]; + /* The same base is used for the small changes. */ + small_idx[0]=small_index; + small_idx[1]=small_index; + small_idx[2]=small_index; + trajcoder_base_compress(encode_ints,runlength*3,small_idx,compress_buffer); +#ifdef SHOWIT + fprintf(stderr,"nbits=%d (%g)\n",nbits,nbits/(runlength*3.)); + nbits_sum+=nbits; + nvalues_sum+=runlength*3; + fprintf(stderr,"Runlength encoded small integers. runlength=%d\n",runlength); +#endif + /* write out base compressed small integers */ + Ptngc_writemanybits(coder,compress_buffer,nbits,&output_ptr); +#ifdef SHOWIT + for (ienc=0; ienc=0) && (instr +#endif +#include +#include +#include +#include +#include +#include "warnmalloc.h" +#include "widemuldiv.h" +#include "bwlzh.h" + +static const double iflipgaincheck=0.89089871814033927; /* 1./(2**(1./6)) */ + +#define MAX_LARGE_RLE 1024 /* Maximum number of large atoms for large RLE. */ +#define MAX_SMALL_RLE 12 /* Maximum number of small atoms in one group. */ + +#define TRESHOLD_INTRA_INTER_DIRECT 1.5 /* How much larger can the direct + frame deltas for the small + triplets be and be accepted anyway + as better than the intra/inter frame + deltas. For better instructions/RLEs. */ + +#define TRESHOLD_INTER_INTRA 5.0 /* How much larger can the intra + frame deltas for the small + triplets be and be accepted anyway + as better than the inter frame + deltas. */ + +/* Difference in indices used for determining whether to store as + large or small. A fun detail in this compression algorithm is that + if everything works fine, large can often be smaller than small, or + at least not as large as is large in magic.c. This is a key idea of + xtc3. */ +#define QUITE_LARGE 3 +#define IS_LARGE 6 + +#if 0 +#define SHOWIT +#endif + +#if 0 +#define SHOWIT_LIGHT +#endif + +/* These routines are in xtc2.c */ +int Ptngc_magic(unsigned int i); +int Ptngc_find_magic_index(unsigned int maxval); + +static unsigned int positive_int(int item) +{ + int s=0; + if (item>0) + s=1+(item-1)*2; + else if (item<0) + s=2+(-item-1)*2; + return s; +} + +static int unpositive_int(int val) +{ + int s=(val+1)/2; + if ((val%2)==0) + s=-s; + return s; +} + + +/* Sequence instructions */ +#define INSTR_DEFAULT 0U +#define INSTR_SMALL_RUNLENGTH 1U +#define INSTR_ONLY_LARGE 2U +#define INSTR_ONLY_SMALL 3U +#define INSTR_FLIP 4U +#define INSTR_LARGE_RLE 5U +#define INSTR_LARGE_DIRECT 6U +#define INSTR_LARGE_INTRA_DELTA 7U +#define INSTR_LARGE_INTER_DELTA 8U + +#define MAXINSTR 9 + +struct xtc3_context +{ + unsigned int *instructions; + int ninstr, ninstr_alloc; + unsigned int *rle; + int nrle, nrle_alloc; + unsigned int *large_direct; + int nlargedir, nlargedir_alloc; + unsigned int *large_intra_delta; + int nlargeintra, nlargeintra_alloc; + unsigned int *large_inter_delta; + int nlargeinter, nlargeinter_alloc; + unsigned int *smallintra; + int nsmallintra, nsmallintra_alloc; + int minint[3],maxint[3]; + int has_large; + int has_large_ints[MAX_LARGE_RLE*3]; /* Large cache. */ + int has_large_type[MAX_LARGE_RLE]; /* What kind of type this large + int is. */ + int current_large_type; +}; + +static void init_xtc3_context(struct xtc3_context *xtc3_context) +{ + xtc3_context->instructions=NULL; + xtc3_context->ninstr=0; + xtc3_context->ninstr_alloc=0; + xtc3_context->rle=NULL; + xtc3_context->nrle=0; + xtc3_context->nrle_alloc=0; + xtc3_context->large_direct=NULL; + xtc3_context->nlargedir=0; + xtc3_context->nlargedir_alloc=0; + xtc3_context->large_intra_delta=NULL; + xtc3_context->nlargeintra=0; + xtc3_context->nlargeintra_alloc=0; + xtc3_context->large_inter_delta=NULL; + xtc3_context->nlargeinter=0; + xtc3_context->nlargeinter_alloc=0; + xtc3_context->smallintra=NULL; + xtc3_context->nsmallintra=0; + xtc3_context->nsmallintra_alloc=0; + xtc3_context->has_large=0; + xtc3_context->current_large_type=0; +} + +static void free_xtc3_context(struct xtc3_context *xtc3_context) +{ + free(xtc3_context->instructions); + free(xtc3_context->rle); + free(xtc3_context->large_direct); + free(xtc3_context->large_intra_delta); + free(xtc3_context->large_inter_delta); + free(xtc3_context->smallintra); +} + +/* Modifies three integer values for better compression of water */ +static void swap_ints(int *in, int *out) +{ + out[0]=in[0]+in[1]; + out[1]=-in[1]; + out[2]=in[1]+in[2]; +} + +static void swap_is_better(int *input, int *minint, int *sum_normal, int *sum_swapped) +{ + int normal_max=0; + int swapped_max=0; + int i,j; + int normal[3]; + int swapped[3]; + for (i=0; i<3; i++) + { + normal[0]=input[i]-minint[i]; + normal[1]=input[3+i]-input[i]; /* minint[i]-minint[i] cancels out */ + normal[2]=input[6+i]-input[3+i]; /* minint[i]-minint[i] cancels out */ + swap_ints(normal,swapped); + for (j=1; j<3; j++) + { + if (positive_int(normal[j])>normal_max) + normal_max=positive_int(normal[j]); + if (positive_int(swapped[j])>swapped_max) + swapped_max=positive_int(swapped[j]); + } + } + if (normal_max==0) + normal_max=1; + if (swapped_max==0) + swapped_max=1; + *sum_normal=normal_max; + *sum_swapped=swapped_max; +} + +static void allocate_enough_memory(unsigned int **ptr, int *nele, int *nele_alloc) +{ + (*nele)++; + if (*nele>*nele_alloc) + { + *nele_alloc=*nele + *nele/2; + *ptr=warnrealloc(*ptr,*nele_alloc*sizeof **ptr); + } +} + +static void insert_value_in_array(unsigned int **ptr, int *nele, int *nele_alloc, + unsigned int value, + char *arrayname) +{ + allocate_enough_memory(ptr,nele,nele_alloc); +#ifdef SHOWIT + fprintf(stderr,"Inserting value %u into array %s @ %d\n",value,arrayname,(*nele)-1); +#endif + (*ptr)[(*nele)-1]=value; +} + + + +static void swapdecide(struct xtc3_context *xtc3_context, int *input,int *swapatoms, int *large_index, int *minint) +{ + int didswap=0; + int normal,swapped; + swap_is_better(input,minint,&normal,&swapped); + /* We have to determine if it is worth to change the behaviour. + If diff is positive it means that it is worth something to + swap. But it costs 4 bits to do the change. If we assume that + we gain 0.17 bit by the swap per value, and the runlength>2 + for four molecules in a row, we gain something. So check if we + gain at least 0.17 bits to even attempt the swap. + */ +#ifdef SHOWIT + fprintf(stderr,"Trying Flip: %g %g\n",(double)swapped/normal, (double)normal/swapped); +#endif + if (((swappedinstructions, + &xtc3_context->ninstr, + &xtc3_context->ninstr_alloc, + INSTR_FLIP,"instr"); + } +} + +/* It is "large" if we have to increase the small index quite a + bit. Not so much to be rejected by the not very large check + later. */ +static int is_quite_large(int *input, int small_index, int max_large_index) +{ + int is=0; + int i; + if (small_index+QUITE_LARGE>=max_large_index) + is=1; + else + { + for (i=0; i<3; i++) + if (positive_int(input[i])>Ptngc_magic(small_index+QUITE_LARGE)) + { + is=1; + break; + } + } + return is; +} + +#ifdef SHOWIT +int nbits_sum; +int nvalues_sum; +#endif + +static void insert_batch(int *input_ptr, int ntriplets_left, int *prevcoord, int *encode_ints, int startenc, int *nenc) +{ + int nencode=startenc*3; + int tmp_prevcoord[3]; + + tmp_prevcoord[0]=prevcoord[0]; + tmp_prevcoord[1]=prevcoord[1]; + tmp_prevcoord[2]=prevcoord[2]; + + if (startenc) + { + int i; + for (i=0; ihas_large_type[i]!=xtc3_context->current_large_type) + { + unsigned int instr; + xtc3_context->current_large_type=xtc3_context->has_large_type[i]; + if (xtc3_context->current_large_type==0) + instr=INSTR_LARGE_DIRECT; + else if (xtc3_context->current_large_type==1) + instr=INSTR_LARGE_INTRA_DELTA; + else + instr=INSTR_LARGE_INTER_DELTA; + insert_value_in_array(&xtc3_context->instructions, + &xtc3_context->ninstr, + &xtc3_context->ninstr_alloc, + instr,"instr"); + } +} + +static void write_three_large(struct xtc3_context *xtc3_context, + int i) +{ + int m; + if (xtc3_context->current_large_type==0) + { + for (m=0; m<3; m++) + insert_value_in_array(&xtc3_context->large_direct, + &xtc3_context->nlargedir, + &xtc3_context->nlargedir_alloc, + xtc3_context->has_large_ints[i*3+m],"large direct"); + } + else if (xtc3_context->current_large_type==1) + { + for (m=0; m<3; m++) + insert_value_in_array(&xtc3_context->large_intra_delta, + &xtc3_context->nlargeintra, + &xtc3_context->nlargeintra_alloc, + xtc3_context->has_large_ints[i*3+m],"large intra"); + } + else + { + for (m=0; m<3; m++) + insert_value_in_array(&xtc3_context->large_inter_delta, + &xtc3_context->nlargeinter, + &xtc3_context->nlargeinter_alloc, + xtc3_context->has_large_ints[i*3+m],"large inter"); + } +} + +static void flush_large(struct xtc3_context *xtc3_context, + int n) /* How many to flush. */ +{ + int i; + i=0; + while (ihas_large_type[i+j]==xtc3_context->has_large_type[i]); + j++); + if (j<3) + { + for (k=0; kinstructions, + &xtc3_context->ninstr, + &xtc3_context->ninstr_alloc, + INSTR_ONLY_LARGE,"instr"); + write_three_large(xtc3_context,i+k); + } + } + else + { + insert_value_in_array(&xtc3_context->instructions, + &xtc3_context->ninstr, + &xtc3_context->ninstr_alloc, + INSTR_LARGE_RLE,"instr"); + insert_value_in_array(&xtc3_context->rle, + &xtc3_context->nrle, + &xtc3_context->nrle_alloc, + (unsigned int)j,"rle (large)"); + for (k=0; khas_large-n)!=0) + { + int j; + for (i=0; ihas_large-n; i++) + { + xtc3_context->has_large_type[i]=xtc3_context->has_large_type[i+n]; + for (j=0; j<3; j++) + xtc3_context->has_large_ints[i*3+j]=xtc3_context->has_large_ints[(i+n)*3+j]; + } + } + xtc3_context->has_large-=n; /* Number of remaining large atoms in buffer */ +} + +static double compute_intlen(unsigned int *ints) +{ + /* The largest value. */ + unsigned int m=ints[0]; + if (ints[1]>m) + m=ints[1]; + if (ints[2]>m) + m=ints[2]; + return (double)m; +} + +static void buffer_large(struct xtc3_context *xtc3_context, int *input, int inpdata, + int natoms, int intradelta_ok) +{ + unsigned int direct[3], intradelta[3]={0,}, interdelta[3]={0,}; + double minlen; + int best_type; + int frame=inpdata/(natoms*3); + int atomframe=inpdata%(natoms*3); + /* If it is full we must write them all. */ + if (xtc3_context->has_large==MAX_LARGE_RLE) + flush_large(xtc3_context,xtc3_context->has_large); /* Flush all. */ + /* Find out which is the best choice for the large integer. Direct coding, or some + kind of delta coding? */ + /* First create direct coding. */ + direct[0]=(unsigned int)(input[inpdata]-xtc3_context->minint[0]); + direct[1]=(unsigned int)(input[inpdata+1]-xtc3_context->minint[1]); + direct[2]=(unsigned int)(input[inpdata+2]-xtc3_context->minint[2]); + minlen=compute_intlen(direct); + best_type=0; /* Direct. */ +#if 1 + /* Then try intra coding if we can. */ + if ((intradelta_ok) && (atomframe>=3)) + { + double thislen; + intradelta[0]=positive_int(input[inpdata]-input[inpdata-3]); + intradelta[1]=positive_int(input[inpdata+1]-input[inpdata-2]); + intradelta[2]=positive_int(input[inpdata+2]-input[inpdata-1]); + thislen=compute_intlen(intradelta); + if (thislen*TRESHOLD_INTRA_INTER_DIRECT0) + { + double thislen; + interdelta[0]=positive_int(input[inpdata]-input[inpdata-natoms*3]); + interdelta[1]=positive_int(input[inpdata+1]-input[inpdata-natoms*3+1]); + interdelta[2]=positive_int(input[inpdata+2]-input[inpdata-natoms*3+2]); + thislen=compute_intlen(interdelta); + if (thislen*TRESHOLD_INTRA_INTER_DIRECThas_large_type[xtc3_context->has_large]=best_type; + if (best_type==0) + { + xtc3_context->has_large_ints[xtc3_context->has_large*3]=direct[0]; + xtc3_context->has_large_ints[xtc3_context->has_large*3+1]=direct[1]; + xtc3_context->has_large_ints[xtc3_context->has_large*3+2]=direct[2]; + } + else if (best_type==1) + { + xtc3_context->has_large_ints[xtc3_context->has_large*3]=intradelta[0]; + xtc3_context->has_large_ints[xtc3_context->has_large*3+1]=intradelta[1]; + xtc3_context->has_large_ints[xtc3_context->has_large*3+2]=intradelta[2]; + } + else if (best_type==2) + { + xtc3_context->has_large_ints[xtc3_context->has_large*3]=interdelta[0]; + xtc3_context->has_large_ints[xtc3_context->has_large*3+1]=interdelta[1]; + xtc3_context->has_large_ints[xtc3_context->has_large*3+2]=interdelta[2]; + } + xtc3_context->has_large++; +} + +static void output_int(unsigned char *output,int *outdata, unsigned int n) +{ + output[(*outdata)++]=((unsigned int)n)&0xFFU; + output[(*outdata)++]=(((unsigned int)n)>>8)&0xFFU; + output[(*outdata)++]=(((unsigned int)n)>>16)&0xFFU; + output[(*outdata)++]=(((unsigned int)n)>>24)&0xFFU; +} + +#if 0 +static void printarray(unsigned int *a, int n, char *name) +{ + int i; + for (i=0; i>(j*8))&0xFFU) + numbytes=i*4+j+1; + return numbytes; +} + +static void base_compress(unsigned int *data, int len, unsigned char *output, int *outlen) +{ + unsigned int largeint[MAXBASEVALS+1]; + unsigned int largeint_tmp[MAXBASEVALS+1]; + int ixyz, i, j; + int nwrittenout=0; + int numbytes=0; + /* Store the MAXBASEVALS value in the output. */ + output[nwrittenout++]=(unsigned char)(MAXBASEVALS&0xFFU); + output[nwrittenout++]=(unsigned char)((MAXBASEVALS>>8)&0xFFU); + /* Store the BASEINTERVAL value in the output. */ + output[nwrittenout++]=(unsigned char)(BASEINTERVAL&0xFFU); + for (ixyz=0; ixyz<3; ixyz++) + { + unsigned int base=0U; + int nvals=0; + int basegiven=0; + for (j=0; jbase) + base=data[k]; + basecheckvals++; + if (basecheckvals==MAXBASEVALS*BASEINTERVAL) + break; + } + /* The base is one larger than the largest values. */ + base++; + if (base<2) + base=2; + /* Store the base in the output. */ + output[nwrittenout++]=(unsigned char)(base&0xFFU); + output[nwrittenout++]=(unsigned char)((base>>8)&0xFFU); + output[nwrittenout++]=(unsigned char)((base>>16)&0xFFU); + output[nwrittenout++]=(unsigned char)((base>>24)&0xFFU); + basegiven=BASEINTERVAL; + /* How many bytes is needed to store MAXBASEVALS values using this base? */ + numbytes=base_bytes(base,MAXBASEVALS); + } + basegiven--; +#ifdef SHOWIT + fprintf(stderr,"Base for %d is %u. I need %d bytes for %d values.\n",ixyz,base,numbytes,MAXBASEVALS); +#endif + } + if (nvals!=0) + { + Ptngc_largeint_mul(base,largeint,largeint_tmp,MAXBASEVALS+1); + for (j=0; j>(ibyte*8))&(0xFFU)); +#ifdef SHOWIT + fprintf(stderr,"%02x",(unsigned int)output[nwrittenout-1]); +#endif + } +#ifdef SHOWIT + fprintf(stderr,"\n"); +#endif + nvals=0; + for (j=0; j>(ibyte*8))&(0xFFU)); + } + } + } + *outlen=nwrittenout; +} + +static void base_decompress(unsigned char *input, int len, unsigned int *output) +{ + unsigned int largeint[MAXMAXBASEVALS+1]; + unsigned int largeint_tmp[MAXMAXBASEVALS+1]; + int ixyz, i, j; + int maxbasevals=(int)((unsigned int)(input[0])|(((unsigned int)(input[1]))<<8)); + int baseinterval=(int)input[2]; + if (maxbasevals>MAXMAXBASEVALS) + { + fprintf(stderr,"Read a larger maxbasevals value from the file than I can handle. Fix" + " by increasing MAXMAXBASEVALS to at least %d. Although, this is" + " probably a bug in TRAJNG, since MAXMAXBASEVALS should already be insanely large enough.\n",maxbasevals); + exit(EXIT_FAILURE); + } + input+=3; + for (ixyz=0; ixyz<3; ixyz++) + { + int numbytes=0; + int nvals_left=len/3; + int outvals=ixyz; + int basegiven=0; + unsigned int base=0U; +#ifdef SHOWIT + fprintf(stderr,"Base for %d is %u. I need %d bytes for %d values.\n",ixyz,base,numbytes,maxbasevals); +#endif + while (nvals_left) + { + int n; + int numints; + if (basegiven==0) + { + base=(unsigned int)(input[0])| + (((unsigned int)(input[1]))<<8)| + (((unsigned int)(input[2]))<<16)| + (((unsigned int)(input[3]))<<24); + input+=4; + basegiven=baseinterval; + /* How many bytes is needed to store maxbasevals values using this base? */ + numbytes=base_bytes(base,maxbasevals); + } + basegiven--; + if (nvals_leftnvals_left) + n=nvals_left; + numints=(numbytes+3)/4; + for (i=n-1; i>=0; i--) + { + output[outvals+i*3]=Ptngc_largeint_div(base,largeint,largeint_tmp,maxbasevals+1); + for (j=0; j14 bits) we return 0, otherwise 1 */ +static int heuristic_bwlzh(unsigned int *ints, int nints) +{ + int i,num; + num=0; + for (i=0; i=16384) + num++; + if (num>nints/10) + return 0; + else + return 1; +} + +/* Speed selects how careful to try to find the most efficient compression. The BWLZH algo is expensive! + Speed <=2 always avoids BWLZH everywhere it is possible. + Speed 3 and 4 and 5 use heuristics (check proportion of large value). This should mostly be safe. + Speed 5 enables the LZ77 component of BWLZH. + Speed 6 always tests if BWLZH is better and if it is uses it. This can be very slow. + */ +unsigned char *Ptngc_pack_array_xtc3(int *input, int *length, int natoms, int speed) +{ + unsigned char *output=NULL; + int i,ienc,j; + int outdata=0; + /* Pack triplets. */ + int ntriplets=*length/3; + int intmax; + int max_small; + int small_index; + int max_large_index; + int large_index[3]; + int prevcoord[3]; + int runlength=0; /* Initial runlength. "Stupidly" set to zero for + simplicity and explicity */ + int swapatoms=0; /* Initial guess is that we should not swap the + first two atoms in each large+small + transition */ + int didswap; /* Whether swapping was actually done. */ + int inpdata=0; + int encode_ints[3+MAX_SMALL_RLE*3]; /* Up to 3 large + 24 small ints can be encoded at once */ + int nencode; + int ntriplets_left=ntriplets; + int refused=0; + unsigned char *bwlzh_buf=NULL; + int bwlzh_buf_len; + unsigned char *base_buf=NULL; + int base_buf_len; + + struct xtc3_context xtc3_context; + init_xtc3_context(&xtc3_context); + + xtc3_context.maxint[0]=xtc3_context.minint[0]=input[0]; + xtc3_context.maxint[1]=xtc3_context.minint[1]=input[1]; + xtc3_context.maxint[2]=xtc3_context.minint[2]=input[2]; + + /* Values of speed should be sane. */ + if (speed<1) + speed=1; + if (speed>6) + speed=6; + +#ifdef SHOWIT + nbits_sum=0; + nvalues_sum=0; +#endif + /* Allocate enough memory for output */ + output=warnmalloc(8* *length*sizeof *output); + + + for (i=1; ixtc3_context.maxint[j]) + xtc3_context.maxint[j]=input[i*3+j]; + if (input[i*3+j]max_large_index) + max_large_index=large_index[1]; + if (large_index[2]>max_large_index) + max_large_index=large_index[2]; + +#ifdef SHOWIT + for (j=0; j<3; j++) + fprintf(stderr,"minint[%d]=%d. maxint[%d]=%d large_index[%d]=%d value=%d\n",j,xtc3_context.minint[j],j,xtc3_context.maxint[j], + j,large_index[j],Ptngc_magic(large_index[j])); +#endif + + /* Guess initial small index */ + small_index=max_large_index/2; + + /* Find the largest value that is not large. Not large is half index of + large. */ + max_small=Ptngc_magic(small_index); + intmax=0; + for (i=0; i<*length; i++) + { + int item=input[i]; + int s=positive_int(item); + if (s>intmax) + if (s0) + { + unsigned int delta[3], delta2[3]; + delta[0]=positive_int(input[inpdata+3]-input[inpdata-natoms*3+3]); + delta[1]=positive_int(input[inpdata+4]-input[inpdata-natoms*3+4]); + delta[2]=positive_int(input[inpdata+5]-input[inpdata-natoms*3+5]); + delta2[0]=positive_int(encode_ints[3]); + delta2[1]=positive_int(encode_ints[4]); + delta2[2]=positive_int(encode_ints[5]); +#ifdef SHOWIT + fprintf(stderr,"A1: inter delta: %u %u %u. intra delta=%u %u %u\n", + delta[0],delta[1],delta[2], + delta2[0],delta2[1],delta2[2]); +#endif + if (compute_intlen(delta)*TRESHOLD_INTER_INTRAlargest_required_base) + largest_required_base=encode_ints[ienc]; + /* Also compute what the largest base is for the current runlength setting! */ + largest_runlength_base=0; + for (ienc=0; (ienclargest_runlength_base) + largest_runlength_base=encode_ints[ienc]; + + largest_required_index=Ptngc_find_magic_index(largest_required_base); + largest_runlength_index=Ptngc_find_magic_index(largest_runlength_base); + + if (largest_required_indexntriplets_left) + new_runlength=ntriplets_left; + + /* We must at least try to get some small integers going. */ + if (new_runlength==0) + { + new_runlength=1; + new_small_index=small_index; + } + + iter_runlength=new_runlength; + iter_small_index=new_small_index; + + /* Iterate to find optimal encoding and runlength */ +#ifdef SHOWIT + fprintf(stderr,"Entering iterative loop.\n"); + fflush(stderr); +#endif + + do { + new_runlength=iter_runlength; + new_small_index=iter_small_index; + +#ifdef SHOWIT + fprintf(stderr,"Test new_small_index=%d Base=%d\n",new_small_index,Ptngc_magic(new_small_index)); +#endif + /* What is the largest runlength + we can do with the currently + selected encoding? Also the max supported runlength is MAX_SMALL_RLE triplets! */ + for (ienc=0; iencnew_small_index) + break; + } + if (ienc/3>new_runlength) + { + iter_runlength=ienc/3; +#ifdef SHOWIT + fprintf(stderr,"I found a new possible runlength: %d\n",iter_runlength); +#endif + } + + /* How large encoding do we have to use? */ + largest_runlength_base=0; + for (ienc=0; ienclargest_runlength_base) + largest_runlength_base=encode_ints[ienc]; + largest_runlength_index=Ptngc_find_magic_index(largest_runlength_base); + if (largest_runlength_index!=new_small_index) + { + iter_small_index=largest_runlength_index; +#ifdef SHOWIT + fprintf(stderr,"I found a new possible small index: %d Base=%d\n",iter_small_index,Ptngc_magic(iter_small_index)); +#endif + } + } while ((new_runlength!=iter_runlength) || + (new_small_index!=iter_small_index)); + +#ifdef SHOWIT + fprintf(stderr,"Exit iterative loop.\n"); + fflush(stderr); +#endif + + /* Verify that we got something good. We may have caught a + substantially larger atom. If so we should just bail + out and let the loop get on another lap. We may have a + minimum runlength though and then we have to fulfill + the request to write out these atoms! */ + rle_index_dep=0; + if (new_runlength<3) + rle_index_dep=IS_LARGE; + else if (new_runlength<6) + rle_index_dep=QUITE_LARGE; + if ((min_runlength) + || ((new_small_index0)) + { + for (i=0; i=2*new_runlength/3)) + { + /* Put all the values in large arrays, instead of the small array */ + if (new_runlength) + { + for (i=0; i=%g?\n",change,isum,(double)Ptngc_magic(small_index+change)*(double)Ptngc_magic(small_index+change)); +#endif + if (isum>(double)Ptngc_magic(small_index+change)*(double)Ptngc_magic(small_index+change)) + { +#ifdef SHOWIT + fprintf(stderr,"Rejected decrease %d of index due to length of vector: %g>=%g\n",change,isum,(double)Ptngc_magic(small_index+change)*(double)Ptngc_magic(small_index+change)); +#endif + rejected=1; + change++; + } + } while ((change<0) && (rejected)); + if (change==0) + break; + } + } + + /* Always accept the new small indices here. */ + small_index=new_small_index; + /* If we have a new runlength emit it */ + if (runlength!=new_runlength) + { + runlength=new_runlength; + insert_value_in_array(&xtc3_context.instructions, + &xtc3_context.ninstr, + &xtc3_context.ninstr_alloc, + INSTR_SMALL_RUNLENGTH,"instr"); + insert_value_in_array(&xtc3_context.rle, + &xtc3_context.nrle, + &xtc3_context.nrle_alloc, + (unsigned int)runlength,"rle (small)"); + } +#ifdef SHOWIT + fprintf(stderr,"Current small index: %d Base=%d\n",small_index,Ptngc_magic(small_index)); +#endif + } + /* If we have a large previous integer we can combine it with a sequence of small ints. */ + if (xtc3_context.has_large) + { + /* If swapatoms is set to 1 but we did actually not + do any swapping, we must first write out the + large atom and then the small. If swapatoms is 1 + and we did swapping we can use the efficient + encoding. */ + if ((swapatoms) && (!didswap)) + { +#ifdef SHOWIT + fprintf(stderr,"Swapatoms was set to 1 but we did not do swapping!\n"); + fprintf(stderr,"Only one large integer.\n"); +#endif + /* Flush all large atoms. */ + flush_large(&xtc3_context,xtc3_context.has_large); +#ifdef SHOWIT + fprintf(stderr,"Sequence of only small integers.\n"); +#endif + insert_value_in_array(&xtc3_context.instructions, + &xtc3_context.ninstr, + &xtc3_context.ninstr_alloc, + INSTR_ONLY_SMALL,"instr"); + } + else + { + +#ifdef SHOWIT + fprintf(stderr,"Sequence of one large and small integers (good compression).\n"); +#endif + /* Flush all large atoms but one! */ + if (xtc3_context.has_large>1) + flush_large(&xtc3_context,xtc3_context.has_large-1); + + /* Here we must check if we should emit a large + type change instruction. */ + large_instruction_change(&xtc3_context,0); + + insert_value_in_array(&xtc3_context.instructions, + &xtc3_context.ninstr, + &xtc3_context.ninstr_alloc, + INSTR_DEFAULT,"instr"); + + write_three_large(&xtc3_context,0); + xtc3_context.has_large=0; + } + } + else + { +#ifdef SHOWIT + fprintf(stderr,"Sequence of only small integers.\n"); +#endif + insert_value_in_array(&xtc3_context.instructions, + &xtc3_context.ninstr, + &xtc3_context.ninstr_alloc, + INSTR_ONLY_SMALL,"instr"); + } + /* Insert the small integers into the small integer array. */ + for (ienc=0; ienc=5) + bwlzh_compress(xtc3_context.instructions,xtc3_context.ninstr,bwlzh_buf,&bwlzh_buf_len); + else + bwlzh_compress_no_lz77(xtc3_context.instructions,xtc3_context.ninstr,bwlzh_buf,&bwlzh_buf_len); + output_int(output,&outdata,(unsigned int)bwlzh_buf_len); + memcpy(output+outdata,bwlzh_buf,bwlzh_buf_len); + outdata+=bwlzh_buf_len; + free(bwlzh_buf); + } + +#if defined(SHOWIT) || defined(SHOWIT_LIGHT) + fprintf(stderr,"rle: %d\n",xtc3_context.nrle); +#endif + + output_int(output,&outdata,(unsigned int)xtc3_context.nrle); + if (xtc3_context.nrle) + { + bwlzh_buf=warnmalloc(bwlzh_get_buflen(xtc3_context.nrle)); + if (speed>=5) + bwlzh_compress(xtc3_context.rle,xtc3_context.nrle,bwlzh_buf,&bwlzh_buf_len); + else + bwlzh_compress_no_lz77(xtc3_context.rle,xtc3_context.nrle,bwlzh_buf,&bwlzh_buf_len); + output_int(output,&outdata,(unsigned int)bwlzh_buf_len); + memcpy(output+outdata,bwlzh_buf,bwlzh_buf_len); + outdata+=bwlzh_buf_len; + free(bwlzh_buf); + } + +#if defined(SHOWIT) || defined(SHOWIT_LIGHT) + fprintf(stderr,"large direct: %d\n",xtc3_context.nlargedir); +#endif + + output_int(output,&outdata,(unsigned int)xtc3_context.nlargedir); + if (xtc3_context.nlargedir) + { + if ((speed<=2) || ((speed<=5) && (!heuristic_bwlzh(xtc3_context.large_direct,xtc3_context.nlargedir)))) + { + bwlzh_buf=NULL; + bwlzh_buf_len=INT_MAX; + } + else + { + bwlzh_buf=warnmalloc(bwlzh_get_buflen(xtc3_context.nlargedir)); + if (speed>=5) + bwlzh_compress(xtc3_context.large_direct,xtc3_context.nlargedir,bwlzh_buf,&bwlzh_buf_len); + else + bwlzh_compress_no_lz77(xtc3_context.large_direct,xtc3_context.nlargedir,bwlzh_buf,&bwlzh_buf_len); + } + /* If this can be written smaller using base compression we should do that. */ + base_buf=warnmalloc((xtc3_context.nlargedir+3)*sizeof(int)); + base_compress(xtc3_context.large_direct,xtc3_context.nlargedir,base_buf,&base_buf_len); +#if defined(SHOWIT) || defined(SHOWIT_LIGHT) + fprintf(stderr,"Large direct: Base len=%d. BWLZH len=%d\n",base_buf_len,bwlzh_buf_len); +#endif + if (base_buf_len=5) + bwlzh_compress(xtc3_context.large_intra_delta,xtc3_context.nlargeintra,bwlzh_buf,&bwlzh_buf_len); + else + bwlzh_compress_no_lz77(xtc3_context.large_intra_delta,xtc3_context.nlargeintra,bwlzh_buf,&bwlzh_buf_len); + } + /* If this can be written smaller using base compression we should do that. */ + base_buf=warnmalloc((xtc3_context.nlargeintra+3)*sizeof(int)); + base_compress(xtc3_context.large_intra_delta,xtc3_context.nlargeintra,base_buf,&base_buf_len); +#if defined(SHOWIT) || defined(SHOWIT_LIGHT) + fprintf(stderr,"Large intra: Base len=%d. BWLZH len=%d\n",base_buf_len,bwlzh_buf_len); +#endif + if (base_buf_len=5) + bwlzh_compress(xtc3_context.large_inter_delta,xtc3_context.nlargeinter,bwlzh_buf,&bwlzh_buf_len); + else + bwlzh_compress_no_lz77(xtc3_context.large_inter_delta,xtc3_context.nlargeinter,bwlzh_buf,&bwlzh_buf_len); + } + /* If this can be written smaller using base compression we should do that. */ + base_buf=warnmalloc((xtc3_context.nlargeinter+3)*sizeof(int)); + base_compress(xtc3_context.large_inter_delta,xtc3_context.nlargeinter,base_buf,&base_buf_len); +#if defined(SHOWIT) || defined(SHOWIT_LIGHT) + fprintf(stderr,"Large inter: Base len=%d. BWLZH len=%d\n",base_buf_len,bwlzh_buf_len); +#endif + if (base_buf_len=5) + bwlzh_compress(xtc3_context.smallintra,xtc3_context.nsmallintra,bwlzh_buf,&bwlzh_buf_len); + else + bwlzh_compress_no_lz77(xtc3_context.smallintra,xtc3_context.nsmallintra,bwlzh_buf,&bwlzh_buf_len); + } + /* If this can be written smaller using base compression we should do that. */ + base_buf=warnmalloc((xtc3_context.nsmallintra+3)*sizeof(int)); + base_compress(xtc3_context.smallintra,xtc3_context.nsmallintra,base_buf,&base_buf_len); +#if defined(SHOWIT) || defined(SHOWIT_LIGHT) + fprintf(stderr,"Small intra: Base len=%d. BWLZH len=%d\n",base_buf_len,bwlzh_buf_len); +#endif + if (base_buf_lenlarge_direct[(*ilargedir)]+minint[0]; + large_ints[1]=(int)xtc3_context->large_direct[(*ilargedir)+1]+minint[1]; + large_ints[2]=(int)xtc3_context->large_direct[(*ilargedir)+2]+minint[2]; + (*ilargedir)+=3; + } + else if (current_large_type==1) + { + large_ints[0]=unpositive_int(xtc3_context->large_intra_delta[(*ilargeintra)])+prevcoord[0]; + large_ints[1]=unpositive_int(xtc3_context->large_intra_delta[(*ilargeintra)+1])+prevcoord[1]; + large_ints[2]=unpositive_int(xtc3_context->large_intra_delta[(*ilargeintra)+2])+prevcoord[2]; + (*ilargeintra)+=3; + } + else + { + large_ints[0]=unpositive_int(xtc3_context->large_inter_delta[(*ilargeinter)]) + +output[outdata-natoms*3+didswap*3]; + large_ints[1]=unpositive_int(xtc3_context->large_inter_delta[(*ilargeinter)+1]) + +output[outdata-natoms*3+1+didswap*3]; + large_ints[2]=unpositive_int(xtc3_context->large_inter_delta[(*ilargeinter)+2]) + +output[outdata-natoms*3+2+didswap*3]; + (*ilargeinter)+=3; + } + prevcoord[0]=large_ints[0]; + prevcoord[1]=large_ints[1]; + prevcoord[2]=large_ints[2]; + output[outdata]=large_ints[0]; + output[outdata+1]=large_ints[1]; + output[outdata+2]=large_ints[2]; +#ifdef SHOWIT + fprintf(stderr,"Unpack one large: %d %d %d\n",prevcoord[0],prevcoord[1],prevcoord[2]); +#endif +} + + +int Ptngc_unpack_array_xtc3(unsigned char *packed,int *output, int length, int natoms) +{ + int i; + int minint[3]; + unsigned char *ptr=packed; + int prevcoord[3]; + int outdata=0; + int ntriplets_left=length/3; + int swapatoms=0; + int runlength=0; + int current_large_type=0; + int iinstr=0; + int irle=0; + int ilargedir=0; + int ilargeintra=0; + int ilargeinter=0; + int ismallintra=0; + + struct xtc3_context xtc3_context; + init_xtc3_context(&xtc3_context); + + for (i=0; i<3; i++) + { + minint[i]=unpositive_int((int)(((unsigned int)ptr[0]) | + (((unsigned int)ptr[1])<<8) | + (((unsigned int)ptr[2])<<16) | + (((unsigned int)ptr[3])<<24))); + ptr+=4; + } + + xtc3_context.ninstr=(int)(((unsigned int)ptr[0]) | + (((unsigned int)ptr[1])<<8) | + (((unsigned int)ptr[2])<<16) | + (((unsigned int)ptr[3])<<24)); + ptr+=4; + if (xtc3_context.ninstr) + decompress_bwlzh_block(&ptr,xtc3_context.ninstr,&xtc3_context.instructions); + + xtc3_context.nrle=(int)(((unsigned int)ptr[0]) | + (((unsigned int)ptr[1])<<8) | + (((unsigned int)ptr[2])<<16) | + (((unsigned int)ptr[3])<<24)); + ptr+=4; + if (xtc3_context.nrle) + decompress_bwlzh_block(&ptr,xtc3_context.nrle,&xtc3_context.rle); + + xtc3_context.nlargedir=(int)(((unsigned int)ptr[0]) | + (((unsigned int)ptr[1])<<8) | + (((unsigned int)ptr[2])<<16) | + (((unsigned int)ptr[3])<<24)); + ptr+=4; + if (xtc3_context.nlargedir) + { + if (*ptr++==1) + decompress_bwlzh_block(&ptr,xtc3_context.nlargedir,&xtc3_context.large_direct); + else + decompress_base_block(&ptr,xtc3_context.nlargedir,&xtc3_context.large_direct); + } + + xtc3_context.nlargeintra=(int)(((unsigned int)ptr[0]) | + (((unsigned int)ptr[1])<<8) | + (((unsigned int)ptr[2])<<16) | + (((unsigned int)ptr[3])<<24)); + ptr+=4; + if (xtc3_context.nlargeintra) + { + if (*ptr++==1) + decompress_bwlzh_block(&ptr,xtc3_context.nlargeintra,&xtc3_context.large_intra_delta); + else + decompress_base_block(&ptr,xtc3_context.nlargeintra,&xtc3_context.large_intra_delta); + } + + xtc3_context.nlargeinter=(int)(((unsigned int)ptr[0]) | + (((unsigned int)ptr[1])<<8) | + (((unsigned int)ptr[2])<<16) | + (((unsigned int)ptr[3])<<24)); + ptr+=4; + if (xtc3_context.nlargeinter) + { + if (*ptr++==1) + decompress_bwlzh_block(&ptr,xtc3_context.nlargeinter,&xtc3_context.large_inter_delta); + else + decompress_base_block(&ptr,xtc3_context.nlargeinter,&xtc3_context.large_inter_delta); + } + + xtc3_context.nsmallintra=(int)(((unsigned int)ptr[0]) | + (((unsigned int)ptr[1])<<8) | + (((unsigned int)ptr[2])<<16) | + (((unsigned int)ptr[3])<<24)); + ptr+=4; + if (xtc3_context.nsmallintra) + { + if (*ptr++==1) + decompress_bwlzh_block(&ptr,xtc3_context.nsmallintra,&xtc3_context.smallintra); + else + decompress_base_block(&ptr,xtc3_context.nsmallintra,&xtc3_context.smallintra); + } + + /* Initial prevcoord is the minimum integers. */ + prevcoord[0]=minint[0]; + prevcoord[1]=minint[1]; + prevcoord[2]=minint[2]; + + while (ntriplets_left>0) + { + int instr=xtc3_context.instructions[iinstr++]; +#ifdef SHOWIT + fprintf(stderr,"instr=%d @ %d\n",instr,iinstr-1); +#endif +#ifdef SHOWIT + fprintf(stderr,"ntriplets left=%d\n",ntriplets_left); +#endif + if ((instr==INSTR_DEFAULT) /* large+small */ + || (instr==INSTR_ONLY_LARGE) /* only large */ + || (instr==INSTR_ONLY_SMALL)) /* only small */ + { + if (instr!=INSTR_ONLY_SMALL) + { + int didswap=0; + if ((instr==INSTR_DEFAULT) && (swapatoms)) + didswap=1; + unpack_one_large(&xtc3_context,&ilargedir, &ilargeintra, &ilargeinter, + prevcoord, minint, output, outdata, didswap, + natoms, current_large_type); + ntriplets_left--; + outdata+=3; + } + if (instr!=INSTR_ONLY_LARGE) + { + for (i=0; itmp$$.h + echo "#define EXPECTED_FILESIZE" $(ls -l test$testnum.tng |awk '{print $5}'). >>tmp$$.h + mv tmp$$.h test$testnum.h + fi +done diff --git a/src/tests/compression/test1.h b/src/tests/compression/test1.h new file mode 100644 index 0000000..ad802e0 --- /dev/null +++ b/src/tests/compression/test1.h @@ -0,0 +1,23 @@ +#define TESTNAME "Initial coding. Intra frame triple algorithm. Cubic cell." +#define FILENAME "test1.tng" +#define ALGOTEST +#define NATOMS 1000 +#define CHUNKY 1 +#define SCALE 0.1 +#define PRECISION 0.01 +#define WRITEVEL 0 +#define VELPRECISION 0.1 +#define INITIALCODING 3 +#define INITIALCODINGPARAMETER -1 +#define CODING 1 +#define CODINGPARAMETER 0 +#define VELCODING 0 +#define VELCODINGPARAMETER 0 +#define INTMIN1 0 +#define INTMIN2 0 +#define INTMIN3 0 +#define INTMAX1 10000 +#define INTMAX2 10000 +#define INTMAX3 10000 +#define NFRAMES 1000 +#define EXPECTED_FILESIZE 2776230. diff --git a/src/tests/compression/test10.h b/src/tests/compression/test10.h new file mode 100644 index 0000000..471dcab --- /dev/null +++ b/src/tests/compression/test10.h @@ -0,0 +1,23 @@ +#define TESTNAME "Coding. Triple intraframe algorithm. Cubic cell." +#define FILENAME "test10.tng" +#define ALGOTEST +#define NATOMS 1000 +#define CHUNKY 100 +#define SCALE 0.1 +#define PRECISION 0.01 +#define WRITEVEL 0 +#define VELPRECISION 0.1 +#define INITIALCODING 5 +#define INITIALCODINGPARAMETER 0 +#define CODING 3 +#define CODINGPARAMETER -1 +#define VELCODING 0 +#define VELCODINGPARAMETER 0 +#define INTMIN1 0 +#define INTMIN2 0 +#define INTMIN3 0 +#define INTMAX1 10000 +#define INTMAX2 10000 +#define INTMAX3 10000 +#define NFRAMES 1000 +#define EXPECTED_FILESIZE 2728492. diff --git a/src/tests/compression/test11.h b/src/tests/compression/test11.h new file mode 100644 index 0000000..b9d7039 --- /dev/null +++ b/src/tests/compression/test11.h @@ -0,0 +1,23 @@ +#define TESTNAME "Coding. Triple one-to-one algorithm. Cubic cell." +#define FILENAME "test11.tng" +#define ALGOTEST +#define NATOMS 1000 +#define CHUNKY 100 +#define SCALE 0.1 +#define PRECISION 0.01 +#define WRITEVEL 0 +#define VELPRECISION 0.1 +#define INITIALCODING 5 +#define INITIALCODINGPARAMETER 0 +#define CODING 7 +#define CODINGPARAMETER -1 +#define VELCODING 0 +#define VELCODINGPARAMETER 0 +#define INTMIN1 0 +#define INTMIN2 0 +#define INTMIN3 0 +#define INTMAX1 10000 +#define INTMAX2 10000 +#define INTMAX3 10000 +#define NFRAMES 1000 +#define EXPECTED_FILESIZE 4293415. diff --git a/src/tests/compression/test12.h b/src/tests/compression/test12.h new file mode 100644 index 0000000..5d692d8 --- /dev/null +++ b/src/tests/compression/test12.h @@ -0,0 +1,23 @@ +#define TESTNAME "Coding. BWLZH interframe algorithm. Cubic cell." +#define FILENAME "test12.tng" +#define ALGOTEST +#define NATOMS 1000 +#define CHUNKY 100 +#define SCALE 0.1 +#define PRECISION 0.01 +#define WRITEVEL 0 +#define VELPRECISION 0.1 +#define INITIALCODING 5 +#define INITIALCODINGPARAMETER 0 +#define CODING 8 +#define CODINGPARAMETER 0 +#define VELCODING 0 +#define VELCODINGPARAMETER 0 +#define INTMIN1 0 +#define INTMIN2 0 +#define INTMIN3 0 +#define INTMAX1 10000 +#define INTMAX2 10000 +#define INTMAX3 10000 +#define NFRAMES 1000 +#define EXPECTED_FILESIZE 894421. diff --git a/src/tests/compression/test13.h b/src/tests/compression/test13.h new file mode 100644 index 0000000..b650f91 --- /dev/null +++ b/src/tests/compression/test13.h @@ -0,0 +1,23 @@ +#define TESTNAME "Coding. BWLZH intraframe algorithm. Cubic cell." +#define FILENAME "test13.tng" +#define ALGOTEST +#define NATOMS 1000 +#define CHUNKY 100 +#define SCALE 0.1 +#define PRECISION 0.01 +#define WRITEVEL 0 +#define VELPRECISION 0.1 +#define INITIALCODING 5 +#define INITIALCODINGPARAMETER 0 +#define CODING 9 +#define CODINGPARAMETER 0 +#define VELCODING 0 +#define VELCODINGPARAMETER 0 +#define INTMIN1 0 +#define INTMIN2 0 +#define INTMIN3 0 +#define INTMAX1 10000 +#define INTMAX2 10000 +#define INTMAX3 10000 +#define NFRAMES 1000 +#define EXPECTED_FILESIZE 840246. diff --git a/src/tests/compression/test14.h b/src/tests/compression/test14.h new file mode 100644 index 0000000..0db6853 --- /dev/null +++ b/src/tests/compression/test14.h @@ -0,0 +1,23 @@ +#define TESTNAME "Coding. XTC3 algorithm. Cubic cell." +#define FILENAME "test14.tng" +#define ALGOTEST +#define NATOMS 1000 +#define CHUNKY 100 +#define SCALE 0.1 +#define PRECISION 0.01 +#define WRITEVEL 0 +#define VELPRECISION 0.1 +#define INITIALCODING 5 +#define INITIALCODINGPARAMETER 0 +#define CODING 10 +#define CODINGPARAMETER 0 +#define VELCODING 0 +#define VELCODINGPARAMETER 0 +#define INTMIN1 0 +#define INTMIN2 0 +#define INTMIN3 0 +#define INTMAX1 10000 +#define INTMAX2 10000 +#define INTMAX3 10000 +#define NFRAMES 1000 +#define EXPECTED_FILESIZE 1401016. diff --git a/src/tests/compression/test15.h b/src/tests/compression/test15.h new file mode 100644 index 0000000..35f94ba --- /dev/null +++ b/src/tests/compression/test15.h @@ -0,0 +1,23 @@ +#define TESTNAME "Initial coding. Automatic selection of algorithms. Cubic cell." +#define FILENAME "test15.tng" +#define ALGOTEST +#define NATOMS 1000 +#define CHUNKY 1 +#define SCALE 0.1 +#define PRECISION 0.01 +#define WRITEVEL 0 +#define VELPRECISION 0.1 +#define INITIALCODING -1 +#define INITIALCODINGPARAMETER -1 +#define CODING 1 +#define CODINGPARAMETER 0 +#define VELCODING 0 +#define VELCODINGPARAMETER 0 +#define INTMIN1 0 +#define INTMIN2 0 +#define INTMIN3 0 +#define INTMAX1 10000 +#define INTMAX2 10000 +#define INTMAX3 10000 +#define NFRAMES 1000 +#define EXPECTED_FILESIZE 2776230. diff --git a/src/tests/compression/test16.h b/src/tests/compression/test16.h new file mode 100644 index 0000000..5f972cb --- /dev/null +++ b/src/tests/compression/test16.h @@ -0,0 +1,24 @@ +#define TESTNAME "Coding. Automatic selection of algorithms. Cubic cell." +#define FILENAME "test16.tng" +#define ALGOTEST +#define NATOMS 1000 +#define CHUNKY 100 +#define SCALE 0.1 +#define PRECISION 0.01 +#define WRITEVEL 0 +#define VELPRECISION 0.1 +#define INITIALCODING -1 +#define INITIALCODINGPARAMETER -1 +#define CODING -1 +#define CODINGPARAMETER -1 +#define VELCODING 0 +#define VELCODINGPARAMETER 0 +#define INTMIN1 0 +#define INTMIN2 0 +#define INTMIN3 0 +#define INTMAX1 10000 +#define INTMAX2 10000 +#define INTMAX3 10000 +#define NFRAMES 1000 +#define SPEED 6 +#define EXPECTED_FILESIZE 838168. diff --git a/src/tests/compression/test17.h b/src/tests/compression/test17.h new file mode 100644 index 0000000..226e34c --- /dev/null +++ b/src/tests/compression/test17.h @@ -0,0 +1,25 @@ +#define TESTNAME "Initial coding of velocities. Stopbits one-to-one . Cubic cell." +#define FILENAME "test17.tng" +#define ALGOTEST +#define NATOMS 1000 +#define CHUNKY 1 +#define SCALE 0.1 +#define PRECISION 0.01 +#define WRITEVEL 1 +#define VELPRECISION 0.1 +#define INITIALCODING 5 +#define INITIALCODINGPARAMETER 0 +#define CODING 1 +#define CODINGPARAMETER 0 +#define INITIALVELCODING 1 +#define INITIALVELCODINGPARAMETER -1 +#define VELCODING 0 +#define VELCODINGPARAMETER 0 +#define INTMIN1 0 +#define INTMIN2 0 +#define INTMIN3 0 +#define INTMAX1 10000 +#define INTMAX2 10000 +#define INTMAX3 10000 +#define NFRAMES 1000 +#define EXPECTED_FILESIZE 7336171. diff --git a/src/tests/compression/test18.h b/src/tests/compression/test18.h new file mode 100644 index 0000000..1248a6d --- /dev/null +++ b/src/tests/compression/test18.h @@ -0,0 +1,25 @@ +#define TESTNAME "Initial coding of velocities. Triplet one-to-one. Cubic cell." +#define FILENAME "test18.tng" +#define ALGOTEST +#define NATOMS 1000 +#define CHUNKY 1 +#define SCALE 0.1 +#define PRECISION 0.01 +#define WRITEVEL 1 +#define VELPRECISION 0.1 +#define INITIALCODING 5 +#define INITIALCODINGPARAMETER 0 +#define CODING 1 +#define CODINGPARAMETER 0 +#define INITIALVELCODING 3 +#define INITIALVELCODINGPARAMETER -1 +#define VELCODING 0 +#define VELCODINGPARAMETER 0 +#define INTMIN1 0 +#define INTMIN2 0 +#define INTMIN3 0 +#define INTMAX1 10000 +#define INTMAX2 10000 +#define INTMAX3 10000 +#define NFRAMES 1000 +#define EXPECTED_FILESIZE 7089695. diff --git a/src/tests/compression/test19.h b/src/tests/compression/test19.h new file mode 100644 index 0000000..ac09bc2 --- /dev/null +++ b/src/tests/compression/test19.h @@ -0,0 +1,25 @@ +#define TESTNAME "Initial coding of velocities. BWLZH one-to-one. Cubic cell." +#define FILENAME "test19.tng" +#define ALGOTEST +#define NATOMS 1000 +#define CHUNKY 1 +#define SCALE 0.1 +#define PRECISION 0.01 +#define WRITEVEL 1 +#define VELPRECISION 0.1 +#define INITIALCODING 5 +#define INITIALCODINGPARAMETER 0 +#define CODING 1 +#define CODINGPARAMETER 0 +#define INITIALVELCODING 9 +#define INITIALVELCODINGPARAMETER 0 +#define VELCODING 0 +#define VELCODINGPARAMETER 0 +#define INTMIN1 0 +#define INTMIN2 0 +#define INTMIN3 0 +#define INTMAX1 10000 +#define INTMAX2 10000 +#define INTMAX3 10000 +#define NFRAMES 50 +#define EXPECTED_FILESIZE 208809. diff --git a/src/tests/compression/test2.h b/src/tests/compression/test2.h new file mode 100644 index 0000000..f2a0f9a --- /dev/null +++ b/src/tests/compression/test2.h @@ -0,0 +1,23 @@ +#define TESTNAME "Initial coding. XTC2 algorithm. Cubic cell." +#define FILENAME "test2.tng" +#define ALGOTEST +#define NATOMS 1000 +#define CHUNKY 1 +#define SCALE 0.1 +#define PRECISION 0.01 +#define WRITEVEL 0 +#define VELPRECISION 0.1 +#define INITIALCODING 5 +#define INITIALCODINGPARAMETER 0 +#define CODING 1 +#define CODINGPARAMETER 0 +#define VELCODING 0 +#define VELCODINGPARAMETER 0 +#define INTMIN1 0 +#define INTMIN2 0 +#define INTMIN3 0 +#define INTMAX1 10000 +#define INTMAX2 10000 +#define INTMAX3 10000 +#define NFRAMES 1000 +#define EXPECTED_FILESIZE 2796171. diff --git a/src/tests/compression/test20.h b/src/tests/compression/test20.h new file mode 100644 index 0000000..7df8f09 --- /dev/null +++ b/src/tests/compression/test20.h @@ -0,0 +1,25 @@ +#define TESTNAME "Coding of velocities. Stopbit one-to-one. Cubic cell." +#define FILENAME "test20.tng" +#define ALGOTEST +#define NATOMS 1000 +#define CHUNKY 100 +#define SCALE 0.1 +#define PRECISION 0.01 +#define WRITEVEL 1 +#define VELPRECISION 0.1 +#define INITIALCODING 5 +#define INITIALCODINGPARAMETER 0 +#define CODING 5 +#define CODINGPARAMETER 0 +#define INITIALVELCODING 1 +#define INITIALVELCODINGPARAMETER -1 +#define VELCODING 1 +#define VELCODINGPARAMETER -1 +#define INTMIN1 0 +#define INTMIN2 0 +#define INTMIN3 0 +#define INTMAX1 10000 +#define INTMAX2 10000 +#define INTMAX3 10000 +#define NFRAMES 1000 +#define EXPECTED_FILESIZE 7237102. diff --git a/src/tests/compression/test21.h b/src/tests/compression/test21.h new file mode 100644 index 0000000..2ea2353 --- /dev/null +++ b/src/tests/compression/test21.h @@ -0,0 +1,25 @@ +#define TESTNAME "Coding of velocities. Triplet inter. Cubic cell." +#define FILENAME "test21.tng" +#define ALGOTEST +#define NATOMS 1000 +#define CHUNKY 100 +#define SCALE 0.1 +#define PRECISION 0.01 +#define WRITEVEL 1 +#define VELPRECISION 0.1 +#define INITIALCODING 5 +#define INITIALCODINGPARAMETER 0 +#define CODING 5 +#define CODINGPARAMETER 0 +#define INITIALVELCODING 1 +#define INITIALVELCODINGPARAMETER -1 +#define VELCODING 2 +#define VELCODINGPARAMETER -1 +#define INTMIN1 0 +#define INTMIN2 0 +#define INTMIN3 0 +#define INTMAX1 10000 +#define INTMAX2 10000 +#define INTMAX3 10000 +#define NFRAMES 1000 +#define EXPECTED_FILESIZE 6214307. diff --git a/src/tests/compression/test22.h b/src/tests/compression/test22.h new file mode 100644 index 0000000..8b16428 --- /dev/null +++ b/src/tests/compression/test22.h @@ -0,0 +1,25 @@ +#define TESTNAME "Coding of velocities. Triplet one-to-one. Cubic cell." +#define FILENAME "test22.tng" +#define ALGOTEST +#define NATOMS 1000 +#define CHUNKY 100 +#define SCALE 0.1 +#define PRECISION 0.01 +#define WRITEVEL 1 +#define VELPRECISION 0.1 +#define INITIALCODING 5 +#define INITIALCODINGPARAMETER 0 +#define CODING 5 +#define CODINGPARAMETER 0 +#define INITIALVELCODING 1 +#define INITIALVELCODINGPARAMETER -1 +#define VELCODING 3 +#define VELCODINGPARAMETER -1 +#define INTMIN1 0 +#define INTMIN2 0 +#define INTMIN3 0 +#define INTMAX1 10000 +#define INTMAX2 10000 +#define INTMAX3 10000 +#define NFRAMES 1000 +#define EXPECTED_FILESIZE 6988699. diff --git a/src/tests/compression/test23.h b/src/tests/compression/test23.h new file mode 100644 index 0000000..85025d6 --- /dev/null +++ b/src/tests/compression/test23.h @@ -0,0 +1,25 @@ +#define TESTNAME "Coding of velocities. Stopbit interframe. Cubic cell." +#define FILENAME "test23.tng" +#define ALGOTEST +#define NATOMS 1000 +#define CHUNKY 100 +#define SCALE 0.1 +#define PRECISION 0.01 +#define WRITEVEL 1 +#define VELPRECISION 0.1 +#define INITIALCODING 5 +#define INITIALCODINGPARAMETER 0 +#define CODING 5 +#define CODINGPARAMETER 0 +#define INITIALVELCODING 1 +#define INITIALVELCODINGPARAMETER -1 +#define VELCODING 6 +#define VELCODINGPARAMETER -1 +#define INTMIN1 0 +#define INTMIN2 0 +#define INTMIN3 0 +#define INTMAX1 10000 +#define INTMAX2 10000 +#define INTMAX3 10000 +#define NFRAMES 1000 +#define EXPECTED_FILESIZE 6494602. diff --git a/src/tests/compression/test24.h b/src/tests/compression/test24.h new file mode 100644 index 0000000..dae943c --- /dev/null +++ b/src/tests/compression/test24.h @@ -0,0 +1,25 @@ +#define TESTNAME "Coding of velocities. BWLZH interframe. Cubic cell." +#define FILENAME "test24.tng" +#define ALGOTEST +#define NATOMS 1000 +#define CHUNKY 25 +#define SCALE 0.1 +#define PRECISION 0.01 +#define WRITEVEL 1 +#define VELPRECISION 0.1 +#define INITIALCODING 5 +#define INITIALCODINGPARAMETER 0 +#define CODING 5 +#define CODINGPARAMETER 0 +#define INITIALVELCODING 1 +#define INITIALVELCODINGPARAMETER -1 +#define VELCODING 8 +#define VELCODINGPARAMETER 0 +#define INTMIN1 0 +#define INTMIN2 0 +#define INTMIN3 0 +#define INTMAX1 10000 +#define INTMAX2 10000 +#define INTMAX3 10000 +#define NFRAMES 50 +#define EXPECTED_FILESIZE 153520. diff --git a/src/tests/compression/test25.h b/src/tests/compression/test25.h new file mode 100644 index 0000000..6145cb9 --- /dev/null +++ b/src/tests/compression/test25.h @@ -0,0 +1,25 @@ +#define TESTNAME "Coding of velocities. BWLZH one-to-one. Cubic cell." +#define FILENAME "test25.tng" +#define ALGOTEST +#define NATOMS 1000 +#define CHUNKY 25 +#define SCALE 0.1 +#define PRECISION 0.01 +#define WRITEVEL 1 +#define VELPRECISION 0.1 +#define INITIALCODING 5 +#define INITIALCODINGPARAMETER 0 +#define CODING 5 +#define CODINGPARAMETER 0 +#define INITIALVELCODING 1 +#define INITIALVELCODINGPARAMETER -1 +#define VELCODING 9 +#define VELCODINGPARAMETER 0 +#define INTMIN1 0 +#define INTMIN2 0 +#define INTMIN3 0 +#define INTMAX1 10000 +#define INTMAX2 10000 +#define INTMAX3 10000 +#define NFRAMES 50 +#define EXPECTED_FILESIZE 154753. diff --git a/src/tests/compression/test26.h b/src/tests/compression/test26.h new file mode 100644 index 0000000..74523ad --- /dev/null +++ b/src/tests/compression/test26.h @@ -0,0 +1,25 @@ +#define TESTNAME "XTC2 algorithm. Orthorhombic cell." +#define FILENAME "test26.tng" +#define ALGOTEST +#define NATOMS 1000 +#define CHUNKY 100 +#define SCALE 0.1 +#define PRECISION 0.01 +#define WRITEVEL 0 +#define VELPRECISION 0.1 +#define INITIALCODING 5 +#define INITIALCODINGPARAMETER 0 +#define CODING 5 +#define CODINGPARAMETER 0 +#define INITIALVELCODING 1 +#define INITIALVELCODINGPARAMETER -1 +#define VELCODING 9 +#define VELCODINGPARAMETER 0 +#define INTMIN1 0 +#define INTMIN2 0 +#define INTMIN3 0 +#define INTMAX1 20000 +#define INTMAX2 10000 +#define INTMAX3 30000 +#define NFRAMES 1000 +#define EXPECTED_FILESIZE 2861948. diff --git a/src/tests/compression/test27.h b/src/tests/compression/test27.h new file mode 100644 index 0000000..30aea07 --- /dev/null +++ b/src/tests/compression/test27.h @@ -0,0 +1,25 @@ +#define TESTNAME "XTC3 algorithm. Orthorhombic cell." +#define FILENAME "test27.tng" +#define ALGOTEST +#define NATOMS 1000 +#define CHUNKY 100 +#define SCALE 0.1 +#define PRECISION 0.01 +#define WRITEVEL 0 +#define VELPRECISION 0.1 +#define INITIALCODING 10 +#define INITIALCODINGPARAMETER 0 +#define CODING 10 +#define CODINGPARAMETER 0 +#define INITIALVELCODING 1 +#define INITIALVELCODINGPARAMETER -1 +#define VELCODING 9 +#define VELCODINGPARAMETER 0 +#define INTMIN1 0 +#define INTMIN2 0 +#define INTMIN3 0 +#define INTMAX1 20000 +#define INTMAX2 10000 +#define INTMAX3 30000 +#define NFRAMES 200 +#define EXPECTED_FILESIZE 282600. diff --git a/src/tests/compression/test28.h b/src/tests/compression/test28.h new file mode 100644 index 0000000..995b81c --- /dev/null +++ b/src/tests/compression/test28.h @@ -0,0 +1,26 @@ +#define TESTNAME "Initial coding. Autoselect algorithm. Repetitive molecule. Cubic cell." +#define FILENAME "test28.tng" +#define ALGOTEST +#define NATOMS 1000 +#define CHUNKY 1 +#define SCALE 0.1 +#define REGULAR +#define PRECISION 0.01 +#define WRITEVEL 0 +#define VELPRECISION 0.1 +#define INITIALCODING -1 +#define INITIALCODINGPARAMETER -1 +#define CODING 0 +#define CODINGPARAMETER 0 +#define INITIALVELCODING 0 +#define INITIALVELCODINGPARAMETER -1 +#define VELCODING 0 +#define VELCODINGPARAMETER 0 +#define INTMIN1 0 +#define INTMIN2 0 +#define INTMIN3 0 +#define INTMAX1 10000 +#define INTMAX2 10000 +#define INTMAX3 10000 +#define NFRAMES 1000 +#define EXPECTED_FILESIZE 1677619. diff --git a/src/tests/compression/test29.h b/src/tests/compression/test29.h new file mode 100644 index 0000000..f48f909 --- /dev/null +++ b/src/tests/compression/test29.h @@ -0,0 +1,26 @@ +#define TESTNAME "Position coding. Autoselect algorithm. Repetitive molecule. Cubic cell." +#define FILENAME "test29.tng" +#define ALGOTEST +#define NATOMS 1000 +#define CHUNKY 100 +#define SCALE 0.1 +#define REGULAR +#define PRECISION 0.01 +#define WRITEVEL 0 +#define VELPRECISION 0.1 +#define INITIALCODING -1 +#define INITIALCODINGPARAMETER -1 +#define CODING -1 +#define CODINGPARAMETER -1 +#define INITIALVELCODING 0 +#define INITIALVELCODINGPARAMETER -1 +#define VELCODING 0 +#define VELCODINGPARAMETER 0 +#define INTMIN1 0 +#define INTMIN2 0 +#define INTMIN3 0 +#define INTMAX1 10000 +#define INTMAX2 10000 +#define INTMAX3 10000 +#define NFRAMES 1000 +#define EXPECTED_FILESIZE 228148. diff --git a/src/tests/compression/test3.h b/src/tests/compression/test3.h new file mode 100644 index 0000000..c4bbebd --- /dev/null +++ b/src/tests/compression/test3.h @@ -0,0 +1,23 @@ +#define TESTNAME "Initial coding. Triplet one-to-one algorithm. Cubic cell." +#define FILENAME "test3.tng" +#define ALGOTEST +#define NATOMS 1000 +#define CHUNKY 1 +#define SCALE 0.1 +#define PRECISION 0.01 +#define WRITEVEL 0 +#define VELPRECISION 0.1 +#define INITIALCODING 7 +#define INITIALCODINGPARAMETER -1 +#define CODING 1 +#define CODINGPARAMETER 0 +#define VELCODING 0 +#define VELCODINGPARAMETER 0 +#define INTMIN1 0 +#define INTMIN2 0 +#define INTMIN3 0 +#define INTMAX1 10000 +#define INTMAX2 10000 +#define INTMAX3 10000 +#define NFRAMES 1000 +#define EXPECTED_FILESIZE 4356773. diff --git a/src/tests/compression/test30.h b/src/tests/compression/test30.h new file mode 100644 index 0000000..2ea607b --- /dev/null +++ b/src/tests/compression/test30.h @@ -0,0 +1,23 @@ +#define TESTNAME "Initial coding. Intra frame triple algorithm. Large system. Cubic cell." +#define FILENAME "test30.tng" +#define ALGOTEST +#define NATOMS 5000000 +#define CHUNKY 1 +#define SCALE 1. +#define PRECISION 1. +#define WRITEVEL 0 +#define VELPRECISION 0.1 +#define INITIALCODING 3 +#define INITIALCODINGPARAMETER -1 +#define CODING 1 +#define CODINGPARAMETER -1 +#define VELCODING 4 +#define VELCODINGPARAMETER 0 +#define INTMIN1 -536870911 +#define INTMIN2 -536870911 +#define INTMIN3 -536870911 +#define INTMAX1 536870911 +#define INTMAX2 536870911 +#define INTMAX3 536870911 +#define NFRAMES 10 +#define EXPECTED_FILESIZE 280198420. diff --git a/src/tests/compression/test31.h b/src/tests/compression/test31.h new file mode 100644 index 0000000..d5d86d4 --- /dev/null +++ b/src/tests/compression/test31.h @@ -0,0 +1,23 @@ +#define TESTNAME "Initial coding. XTC2 algorithm. Large system. Cubic cell." +#define FILENAME "test31.tng" +#define ALGOTEST +#define NATOMS 5000000 +#define CHUNKY 1 +#define SCALE 1. +#define PRECISION 1. +#define WRITEVEL 0 +#define VELPRECISION 0.1 +#define INITIALCODING 5 +#define INITIALCODINGPARAMETER 0 +#define CODING 1 +#define CODINGPARAMETER -1 +#define VELCODING 4 +#define VELCODINGPARAMETER 0 +#define INTMIN1 -536870911 +#define INTMIN2 -536870911 +#define INTMIN3 -536870911 +#define INTMAX1 536870911 +#define INTMAX2 536870911 +#define INTMAX3 536870911 +#define NFRAMES 10 +#define EXPECTED_FILESIZE 301463456. diff --git a/src/tests/compression/test32.h b/src/tests/compression/test32.h new file mode 100644 index 0000000..7b66438 --- /dev/null +++ b/src/tests/compression/test32.h @@ -0,0 +1,23 @@ +#define TESTNAME "Initial coding. XTC3 algorithm. Large system. Cubic cell." +#define FILENAME "test32.tng" +#define ALGOTEST +#define NATOMS 5000000 +#define CHUNKY 1 +#define SCALE 1. +#define PRECISION 1. +#define WRITEVEL 0 +#define VELPRECISION 0.1 +#define INITIALCODING 10 +#define INITIALCODINGPARAMETER 0 +#define CODING 1 +#define CODINGPARAMETER -1 +#define VELCODING 4 +#define VELCODINGPARAMETER 0 +#define INTMIN1 -536870911 +#define INTMIN2 -536870911 +#define INTMIN3 -536870911 +#define INTMAX1 536870911 +#define INTMAX2 536870911 +#define INTMAX3 536870911 +#define NFRAMES 2 +#define EXPECTED_FILESIZE 31668187. diff --git a/src/tests/compression/test33.h b/src/tests/compression/test33.h new file mode 100644 index 0000000..9463eeb --- /dev/null +++ b/src/tests/compression/test33.h @@ -0,0 +1,23 @@ +#define TESTNAME "Initial coding. Intra frame BWLZH algorithm. Large system. Cubic cell." +#define FILENAME "test33.tng" +#define ALGOTEST +#define NATOMS 5000000 +#define CHUNKY 1 +#define SCALE 1. +#define PRECISION 1. +#define WRITEVEL 0 +#define VELPRECISION 0.1 +#define INITIALCODING 9 +#define INITIALCODINGPARAMETER 0 +#define CODING 1 +#define CODINGPARAMETER -1 +#define VELCODING 4 +#define VELCODINGPARAMETER 0 +#define INTMIN1 -536870911 +#define INTMIN2 -536870911 +#define INTMIN3 -536870911 +#define INTMAX1 536870911 +#define INTMAX2 536870911 +#define INTMAX3 536870911 +#define NFRAMES 2 +#define EXPECTED_FILESIZE 7121047. diff --git a/src/tests/compression/test34.h b/src/tests/compression/test34.h new file mode 100644 index 0000000..14fe88f --- /dev/null +++ b/src/tests/compression/test34.h @@ -0,0 +1,23 @@ +#define TESTNAME "Position coding. Stop bits algorithm. Large system. Cubic cell." +#define FILENAME "test34.tng" +#define ALGOTEST +#define NATOMS 5000000 +#define CHUNKY 2 +#define SCALE 1. +#define PRECISION 1. +#define WRITEVEL 0 +#define VELPRECISION 0.1 +#define INITIALCODING 5 +#define INITIALCODINGPARAMETER 0 +#define CODING 1 +#define CODINGPARAMETER -1 +#define VELCODING 4 +#define VELCODINGPARAMETER 0 +#define INTMIN1 -536870911 +#define INTMIN2 -536870911 +#define INTMIN3 -536870911 +#define INTMAX1 536870911 +#define INTMAX2 536870911 +#define INTMAX3 536870911 +#define NFRAMES 10 +#define EXPECTED_FILESIZE 250247372. diff --git a/src/tests/compression/test35.h b/src/tests/compression/test35.h new file mode 100644 index 0000000..44af4c4 --- /dev/null +++ b/src/tests/compression/test35.h @@ -0,0 +1,23 @@ +#define TESTNAME "Position coding. Inter frame triple algorithm. Large system. Cubic cell." +#define FILENAME "test35.tng" +#define ALGOTEST +#define NATOMS 5000000 +#define CHUNKY 2 +#define SCALE 1. +#define PRECISION 1. +#define WRITEVEL 0 +#define VELPRECISION 0.1 +#define INITIALCODING 5 +#define INITIALCODINGPARAMETER 0 +#define CODING 2 +#define CODINGPARAMETER -1 +#define VELCODING 4 +#define VELCODINGPARAMETER 0 +#define INTMIN1 -536870911 +#define INTMIN2 -536870911 +#define INTMIN3 -536870911 +#define INTMAX1 536870911 +#define INTMAX2 536870911 +#define INTMAX3 536870911 +#define NFRAMES 10 +#define EXPECTED_FILESIZE 243598962. diff --git a/src/tests/compression/test36.h b/src/tests/compression/test36.h new file mode 100644 index 0000000..56f5dcd --- /dev/null +++ b/src/tests/compression/test36.h @@ -0,0 +1,23 @@ +#define TESTNAME "Position coding. Intra frame triple algorithm. Large system. Cubic cell." +#define FILENAME "test36.tng" +#define ALGOTEST +#define NATOMS 5000000 +#define CHUNKY 2 +#define SCALE 1. +#define PRECISION 1. +#define WRITEVEL 0 +#define VELPRECISION 0.1 +#define INITIALCODING 5 +#define INITIALCODINGPARAMETER 0 +#define CODING 3 +#define CODINGPARAMETER -1 +#define VELCODING 4 +#define VELCODINGPARAMETER 0 +#define INTMIN1 -536870911 +#define INTMIN2 -536870911 +#define INTMIN3 -536870911 +#define INTMAX1 536870911 +#define INTMAX2 536870911 +#define INTMAX3 536870911 +#define NFRAMES 10 +#define EXPECTED_FILESIZE 290800607. diff --git a/src/tests/compression/test37.h b/src/tests/compression/test37.h new file mode 100644 index 0000000..68b4872 --- /dev/null +++ b/src/tests/compression/test37.h @@ -0,0 +1,23 @@ +#define TESTNAME "Position coding. XTC2 algorithm. Large system. Cubic cell." +#define FILENAME "test37.tng" +#define ALGOTEST +#define NATOMS 5000000 +#define CHUNKY 2 +#define SCALE 1. +#define PRECISION 1. +#define WRITEVEL 0 +#define VELPRECISION 0.1 +#define INITIALCODING 5 +#define INITIALCODINGPARAMETER 0 +#define CODING 5 +#define CODINGPARAMETER 0 +#define VELCODING 4 +#define VELCODINGPARAMETER 0 +#define INTMIN1 -536870911 +#define INTMIN2 -536870911 +#define INTMIN3 -536870911 +#define INTMAX1 536870911 +#define INTMAX2 536870911 +#define INTMAX3 536870911 +#define NFRAMES 10 +#define EXPECTED_FILESIZE 301463256. diff --git a/src/tests/compression/test38.h b/src/tests/compression/test38.h new file mode 100644 index 0000000..7e43348 --- /dev/null +++ b/src/tests/compression/test38.h @@ -0,0 +1,23 @@ +#define TESTNAME "Position coding. XTC3 algorithm. Large system. Cubic cell." +#define FILENAME "test38.tng" +#define ALGOTEST +#define NATOMS 5000000 +#define CHUNKY 2 +#define SCALE 1. +#define PRECISION 1. +#define WRITEVEL 0 +#define VELPRECISION 0.1 +#define INITIALCODING 10 +#define INITIALCODINGPARAMETER 0 +#define CODING 10 +#define CODINGPARAMETER 0 +#define VELCODING 4 +#define VELCODINGPARAMETER 0 +#define INTMIN1 -536870911 +#define INTMIN2 -536870911 +#define INTMIN3 -536870911 +#define INTMAX1 536870911 +#define INTMAX2 536870911 +#define INTMAX3 536870911 +#define NFRAMES 4 +#define EXPECTED_FILESIZE 63482016. diff --git a/src/tests/compression/test39.h b/src/tests/compression/test39.h new file mode 100644 index 0000000..bb3ecc0 --- /dev/null +++ b/src/tests/compression/test39.h @@ -0,0 +1,23 @@ +#define TESTNAME "Position coding. Intra frame BWLZH algorithm. Large system. Cubic cell." +#define FILENAME "test39.tng" +#define ALGOTEST +#define NATOMS 5000000 +#define CHUNKY 2 +#define SCALE 1. +#define PRECISION 1. +#define WRITEVEL 0 +#define VELPRECISION 0.1 +#define INITIALCODING 5 +#define INITIALCODINGPARAMETER 0 +#define CODING 9 +#define CODINGPARAMETER 0 +#define VELCODING 4 +#define VELCODINGPARAMETER 0 +#define INTMIN1 -536870911 +#define INTMIN2 -536870911 +#define INTMIN3 -536870911 +#define INTMAX1 536870911 +#define INTMAX2 536870911 +#define INTMAX3 536870911 +#define NFRAMES 4 +#define EXPECTED_FILESIZE 67631371. diff --git a/src/tests/compression/test4.h b/src/tests/compression/test4.h new file mode 100644 index 0000000..63ec9da --- /dev/null +++ b/src/tests/compression/test4.h @@ -0,0 +1,23 @@ +#define TESTNAME "Initial coding. BWLZH one-to-one algorithm. Cubic cell." +#define FILENAME "test4.tng" +#define ALGOTEST +#define NATOMS 1000 +#define CHUNKY 1 +#define SCALE 0.1 +#define PRECISION 0.01 +#define WRITEVEL 0 +#define VELPRECISION 0.1 +#define INITIALCODING 9 +#define INITIALCODINGPARAMETER 0 +#define CODING 1 +#define CODINGPARAMETER 0 +#define VELCODING 0 +#define VELCODINGPARAMETER 0 +#define INTMIN1 0 +#define INTMIN2 0 +#define INTMIN3 0 +#define INTMAX1 10000 +#define INTMAX2 10000 +#define INTMAX3 10000 +#define NFRAMES 1000 +#define EXPECTED_FILESIZE 2572043. diff --git a/src/tests/compression/test40.h b/src/tests/compression/test40.h new file mode 100644 index 0000000..5eb9c0d --- /dev/null +++ b/src/tests/compression/test40.h @@ -0,0 +1,23 @@ +#define TESTNAME "Position coding. Inter frame BWLZH algorithm. Large system. Cubic cell." +#define FILENAME "test40.tng" +#define ALGOTEST +#define NATOMS 5000000 +#define CHUNKY 2 +#define SCALE 1. +#define PRECISION 1. +#define WRITEVEL 0 +#define VELPRECISION 0.1 +#define INITIALCODING 5 +#define INITIALCODINGPARAMETER 0 +#define CODING 8 +#define CODINGPARAMETER 0 +#define VELCODING 4 +#define VELCODINGPARAMETER 0 +#define INTMIN1 -536870911 +#define INTMIN2 -536870911 +#define INTMIN3 -536870911 +#define INTMAX1 536870911 +#define INTMAX2 536870911 +#define INTMAX3 536870911 +#define NFRAMES 4 +#define EXPECTED_FILESIZE 63822378. diff --git a/src/tests/compression/test41.h b/src/tests/compression/test41.h new file mode 100644 index 0000000..f673944 --- /dev/null +++ b/src/tests/compression/test41.h @@ -0,0 +1,23 @@ +#define TESTNAME "Initial coding. Intra frame triple algorithm. High accuracy. Cubic cell." +#define FILENAME "test41.tng" +#define ALGOTEST +#define NATOMS 100000 +#define CHUNKY 1 +#define SCALE 0.5 +#define PRECISION 1e-8 +#define WRITEVEL 0 +#define VELPRECISION 0.1 +#define INITIALCODING 3 +#define INITIALCODINGPARAMETER -1 +#define CODING 1 +#define CODINGPARAMETER -1 +#define VELCODING 4 +#define VELCODINGPARAMETER 0 +#define INTMIN1 0 +#define INTMIN2 0 +#define INTMIN3 0 +#define INTMAX1 1610612736 +#define INTMAX2 1610612736 +#define INTMAX3 1610612736 +#define NFRAMES 100 +#define EXPECTED_FILESIZE 53179342. diff --git a/src/tests/compression/test42.h b/src/tests/compression/test42.h new file mode 100644 index 0000000..b40f963 --- /dev/null +++ b/src/tests/compression/test42.h @@ -0,0 +1,23 @@ +#define TESTNAME "Initial coding. XTC2 algorithm. High accuracy. Cubic cell." +#define FILENAME "test42.tng" +#define ALGOTEST +#define NATOMS 100000 +#define CHUNKY 1 +#define SCALE 0.5 +#define PRECISION 1e-8 +#define WRITEVEL 0 +#define VELPRECISION 0.1 +#define INITIALCODING 5 +#define INITIALCODINGPARAMETER 0 +#define CODING 1 +#define CODINGPARAMETER -1 +#define VELCODING 4 +#define VELCODINGPARAMETER 0 +#define INTMIN1 0 +#define INTMIN2 0 +#define INTMIN3 0 +#define INTMAX1 1610612736 +#define INTMAX2 1610612736 +#define INTMAX3 1610612736 +#define NFRAMES 100 +#define EXPECTED_FILESIZE 57283715. diff --git a/src/tests/compression/test43.h b/src/tests/compression/test43.h new file mode 100644 index 0000000..38a59db --- /dev/null +++ b/src/tests/compression/test43.h @@ -0,0 +1,23 @@ +#define TESTNAME "Initial coding. XTC3 algorithm. High accuracy. Cubic cell." +#define FILENAME "test43.tng" +#define ALGOTEST +#define NATOMS 100000 +#define CHUNKY 1 +#define SCALE 0.5 +#define PRECISION 1e-8 +#define WRITEVEL 0 +#define VELPRECISION 0.1 +#define INITIALCODING 10 +#define INITIALCODINGPARAMETER 0 +#define CODING 1 +#define CODINGPARAMETER -1 +#define VELCODING 4 +#define VELCODINGPARAMETER 0 +#define INTMIN1 0 +#define INTMIN2 0 +#define INTMIN3 0 +#define INTMAX1 1610612736 +#define INTMAX2 1610612736 +#define INTMAX3 1610612736 +#define NFRAMES 10 +#define EXPECTED_FILESIZE 3783912. diff --git a/src/tests/compression/test44.h b/src/tests/compression/test44.h new file mode 100644 index 0000000..0bab871 --- /dev/null +++ b/src/tests/compression/test44.h @@ -0,0 +1,23 @@ +#define TESTNAME "Initial coding. Intra frame BWLZH algorithm. High accuracy. Cubic cell." +#define FILENAME "test44.tng" +#define ALGOTEST +#define NATOMS 100000 +#define CHUNKY 1 +#define SCALE 0.5 +#define PRECISION 1e-8 +#define WRITEVEL 0 +#define VELPRECISION 0.1 +#define INITIALCODING 9 +#define INITIALCODINGPARAMETER 0 +#define CODING 1 +#define CODINGPARAMETER -1 +#define VELCODING 4 +#define VELCODINGPARAMETER 0 +#define INTMIN1 0 +#define INTMIN2 0 +#define INTMIN3 0 +#define INTMAX1 1610612736 +#define INTMAX2 1610612736 +#define INTMAX3 1610612736 +#define NFRAMES 10 +#define EXPECTED_FILESIZE 1436901. diff --git a/src/tests/compression/test45.h b/src/tests/compression/test45.h new file mode 100644 index 0000000..095efed --- /dev/null +++ b/src/tests/compression/test45.h @@ -0,0 +1,23 @@ +#define TESTNAME "Position coding. Stop bits algorithm. High accuracy. Cubic cell." +#define FILENAME "test45.tng" +#define ALGOTEST +#define NATOMS 100000 +#define CHUNKY 10 +#define SCALE 0.5 +#define PRECISION 1e-8 +#define WRITEVEL 0 +#define VELPRECISION 0.1 +#define INITIALCODING 3 +#define INITIALCODINGPARAMETER -1 +#define CODING 1 +#define CODINGPARAMETER -1 +#define VELCODING 4 +#define VELCODINGPARAMETER 0 +#define INTMIN1 0 +#define INTMIN2 0 +#define INTMIN3 0 +#define INTMAX1 1610612736 +#define INTMAX2 1610612736 +#define INTMAX3 1610612736 +#define NFRAMES 100 +#define EXPECTED_FILESIZE 36794379. diff --git a/src/tests/compression/test46.h b/src/tests/compression/test46.h new file mode 100644 index 0000000..dc0b8bb --- /dev/null +++ b/src/tests/compression/test46.h @@ -0,0 +1,23 @@ +#define TESTNAME "Position coding. Inter frame triple algorithm. High accuracy. Cubic cell." +#define FILENAME "test46.tng" +#define ALGOTEST +#define NATOMS 100000 +#define CHUNKY 10 +#define SCALE 0.5 +#define PRECISION 1e-8 +#define WRITEVEL 0 +#define VELPRECISION 0.1 +#define INITIALCODING 3 +#define INITIALCODINGPARAMETER -1 +#define CODING 2 +#define CODINGPARAMETER -1 +#define VELCODING 4 +#define VELCODINGPARAMETER 0 +#define INTMIN1 0 +#define INTMIN2 0 +#define INTMIN3 0 +#define INTMAX1 1610612736 +#define INTMAX2 1610612736 +#define INTMAX3 1610612736 +#define NFRAMES 100 +#define EXPECTED_FILESIZE 34508770. diff --git a/src/tests/compression/test47.h b/src/tests/compression/test47.h new file mode 100644 index 0000000..5a5182d --- /dev/null +++ b/src/tests/compression/test47.h @@ -0,0 +1,23 @@ +#define TESTNAME "Position coding. Intra frame triple algorithm. High accuracy. Cubic cell." +#define FILENAME "test47.tng" +#define ALGOTEST +#define NATOMS 100000 +#define CHUNKY 10 +#define SCALE 0.5 +#define PRECISION 1e-8 +#define WRITEVEL 0 +#define VELPRECISION 0.1 +#define INITIALCODING 3 +#define INITIALCODINGPARAMETER -1 +#define CODING 3 +#define CODINGPARAMETER -1 +#define VELCODING 4 +#define VELCODINGPARAMETER 0 +#define INTMIN1 0 +#define INTMIN2 0 +#define INTMIN3 0 +#define INTMAX1 1610612736 +#define INTMAX2 1610612736 +#define INTMAX3 1610612736 +#define NFRAMES 100 +#define EXPECTED_FILESIZE 53174711. diff --git a/src/tests/compression/test48.h b/src/tests/compression/test48.h new file mode 100644 index 0000000..1562779 --- /dev/null +++ b/src/tests/compression/test48.h @@ -0,0 +1,23 @@ +#define TESTNAME "Position coding. XTC2 algorithm. High accuracy. Cubic cell." +#define FILENAME "test48.tng" +#define ALGOTEST +#define NATOMS 100000 +#define CHUNKY 10 +#define SCALE 0.5 +#define PRECISION 1e-8 +#define WRITEVEL 0 +#define VELPRECISION 0.1 +#define INITIALCODING 3 +#define INITIALCODINGPARAMETER -1 +#define CODING 5 +#define CODINGPARAMETER 0 +#define VELCODING 4 +#define VELCODINGPARAMETER 0 +#define INTMIN1 0 +#define INTMIN2 0 +#define INTMIN3 0 +#define INTMAX1 805306368 +#define INTMAX2 805306368 +#define INTMAX3 805306368 +#define NFRAMES 100 +#define EXPECTED_FILESIZE 55638414. diff --git a/src/tests/compression/test49.h b/src/tests/compression/test49.h new file mode 100644 index 0000000..dcad4ec --- /dev/null +++ b/src/tests/compression/test49.h @@ -0,0 +1,23 @@ +#define TESTNAME "Position coding. XTC3 algorithm. High accuracy. Cubic cell." +#define FILENAME "test49.tng" +#define ALGOTEST +#define NATOMS 100000 +#define CHUNKY 10 +#define SCALE 0.5 +#define PRECISION 1e-8 +#define WRITEVEL 0 +#define VELPRECISION 0.1 +#define INITIALCODING 3 +#define INITIALCODINGPARAMETER -1 +#define CODING 10 +#define CODINGPARAMETER 0 +#define VELCODING 4 +#define VELCODINGPARAMETER 0 +#define INTMIN1 0 +#define INTMIN2 0 +#define INTMIN3 0 +#define INTMAX1 805306368 +#define INTMAX2 805306368 +#define INTMAX3 805306368 +#define NFRAMES 20 +#define EXPECTED_FILESIZE 3585605. diff --git a/src/tests/compression/test5.h b/src/tests/compression/test5.h new file mode 100644 index 0000000..a933044 --- /dev/null +++ b/src/tests/compression/test5.h @@ -0,0 +1,23 @@ +#define TESTNAME "Initial coding. XTC3 algorithm. Cubic cell." +#define FILENAME "test5.tng" +#define ALGOTEST +#define NATOMS 1000 +#define CHUNKY 1 +#define SCALE 0.1 +#define PRECISION 0.01 +#define WRITEVEL 0 +#define VELPRECISION 0.1 +#define INITIALCODING 10 +#define INITIALCODINGPARAMETER 0 +#define CODING 1 +#define CODINGPARAMETER 0 +#define VELCODING 0 +#define VELCODINGPARAMETER 0 +#define INTMIN1 0 +#define INTMIN2 0 +#define INTMIN3 0 +#define INTMAX1 10000 +#define INTMAX2 10000 +#define INTMAX3 10000 +#define NFRAMES 1000 +#define EXPECTED_FILESIZE 3346179. diff --git a/src/tests/compression/test50.h b/src/tests/compression/test50.h new file mode 100644 index 0000000..a00e209 --- /dev/null +++ b/src/tests/compression/test50.h @@ -0,0 +1,23 @@ +#define TESTNAME "Position coding. Intra frame BWLZH algorithm. High accuracy. Cubic cell." +#define FILENAME "test50.tng" +#define ALGOTEST +#define NATOMS 100000 +#define CHUNKY 10 +#define SCALE 0.5 +#define PRECISION 1e-8 +#define WRITEVEL 0 +#define VELPRECISION 0.1 +#define INITIALCODING 3 +#define INITIALCODINGPARAMETER -1 +#define CODING 9 +#define CODINGPARAMETER 0 +#define VELCODING 4 +#define VELCODINGPARAMETER 0 +#define INTMIN1 0 +#define INTMIN2 0 +#define INTMIN3 0 +#define INTMAX1 805306368 +#define INTMAX2 805306368 +#define INTMAX3 805306368 +#define NFRAMES 20 +#define EXPECTED_FILESIZE 3143379. diff --git a/src/tests/compression/test51.h b/src/tests/compression/test51.h new file mode 100644 index 0000000..4641b61 --- /dev/null +++ b/src/tests/compression/test51.h @@ -0,0 +1,23 @@ +#define TESTNAME "Position coding. Inter frame BWLZH algorithm. High accuracy. Cubic cell." +#define FILENAME "test51.tng" +#define ALGOTEST +#define NATOMS 100000 +#define CHUNKY 10 +#define SCALE 0.5 +#define PRECISION 1e-8 +#define WRITEVEL 0 +#define VELPRECISION 0.1 +#define INITIALCODING 3 +#define INITIALCODINGPARAMETER -1 +#define CODING 8 +#define CODINGPARAMETER 0 +#define VELCODING 4 +#define VELCODINGPARAMETER 0 +#define INTMIN1 0 +#define INTMIN2 0 +#define INTMIN3 0 +#define INTMAX1 805306368 +#define INTMAX2 805306368 +#define INTMAX3 805306368 +#define NFRAMES 20 +#define EXPECTED_FILESIZE 2897696. diff --git a/src/tests/compression/test52.h b/src/tests/compression/test52.h new file mode 100644 index 0000000..2f49a8b --- /dev/null +++ b/src/tests/compression/test52.h @@ -0,0 +1,24 @@ +#define TESTNAME "Velocity coding. Stop bits algorithm. High accuracy. Cubic cell." +#define FILENAME "test52.tng" +#define ALGOTEST +#define NATOMS 100000 +#define CHUNKY 10 +#define SCALE 0.5 +#define PRECISION 1e-8 +#define WRITEVEL 1 +#define VELINTMUL 100000 +#define VELPRECISION 1e-8 +#define INITIALCODING 3 +#define INITIALCODINGPARAMETER -1 +#define CODING 5 +#define CODINGPARAMETER 0 +#define VELCODING 1 +#define VELCODINGPARAMETER -1 +#define INTMIN1 0 +#define INTMIN2 0 +#define INTMIN3 0 +#define INTMAX1 805306368 +#define INTMAX2 805306368 +#define INTMAX3 805306368 +#define NFRAMES 100 +#define EXPECTED_FILESIZE 173083705. diff --git a/src/tests/compression/test53.h b/src/tests/compression/test53.h new file mode 100644 index 0000000..679b680 --- /dev/null +++ b/src/tests/compression/test53.h @@ -0,0 +1,24 @@ +#define TESTNAME "Velocity coding. Triple algorithm. High accuracy. Cubic cell." +#define FILENAME "test53.tng" +#define ALGOTEST +#define NATOMS 100000 +#define CHUNKY 10 +#define SCALE 0.5 +#define PRECISION 1e-8 +#define WRITEVEL 1 +#define VELINTMUL 100000 +#define VELPRECISION 1e-8 +#define INITIALCODING 3 +#define INITIALCODINGPARAMETER -1 +#define CODING 5 +#define CODINGPARAMETER 0 +#define VELCODING 3 +#define VELCODINGPARAMETER -1 +#define INTMIN1 0 +#define INTMIN2 0 +#define INTMIN3 0 +#define INTMAX1 805306368 +#define INTMAX2 805306368 +#define INTMAX3 805306368 +#define NFRAMES 100 +#define EXPECTED_FILESIZE 168548573. diff --git a/src/tests/compression/test54.h b/src/tests/compression/test54.h new file mode 100644 index 0000000..c2c2869 --- /dev/null +++ b/src/tests/compression/test54.h @@ -0,0 +1,24 @@ +#define TESTNAME "Velocity coding. Interframe triple algorithm. High accuracy. Cubic cell." +#define FILENAME "test54.tng" +#define ALGOTEST +#define NATOMS 100000 +#define CHUNKY 10 +#define SCALE 0.5 +#define PRECISION 1e-8 +#define WRITEVEL 1 +#define VELINTMUL 100000 +#define VELPRECISION 1e-8 +#define INITIALCODING 3 +#define INITIALCODINGPARAMETER -1 +#define CODING 5 +#define CODINGPARAMETER 0 +#define VELCODING 2 +#define VELCODINGPARAMETER -1 +#define INTMIN1 0 +#define INTMIN2 0 +#define INTMIN3 0 +#define INTMAX1 805306368 +#define INTMAX2 805306368 +#define INTMAX3 805306368 +#define NFRAMES 100 +#define EXPECTED_FILESIZE 161798573. diff --git a/src/tests/compression/test55.h b/src/tests/compression/test55.h new file mode 100644 index 0000000..a094d97 --- /dev/null +++ b/src/tests/compression/test55.h @@ -0,0 +1,24 @@ +#define TESTNAME "Velocity coding. Interframe stop-bits algorithm. High accuracy. Cubic cell." +#define FILENAME "test55.tng" +#define ALGOTEST +#define NATOMS 100000 +#define CHUNKY 10 +#define SCALE 0.5 +#define PRECISION 1e-8 +#define WRITEVEL 1 +#define VELINTMUL 100000 +#define VELPRECISION 1e-8 +#define INITIALCODING 3 +#define INITIALCODINGPARAMETER -1 +#define CODING 5 +#define CODINGPARAMETER 0 +#define VELCODING 6 +#define VELCODINGPARAMETER -1 +#define INTMIN1 0 +#define INTMIN2 0 +#define INTMIN3 0 +#define INTMAX1 805306368 +#define INTMAX2 805306368 +#define INTMAX3 805306368 +#define NFRAMES 100 +#define EXPECTED_FILESIZE 166298533. diff --git a/src/tests/compression/test56.h b/src/tests/compression/test56.h new file mode 100644 index 0000000..6655ad0 --- /dev/null +++ b/src/tests/compression/test56.h @@ -0,0 +1,24 @@ +#define TESTNAME "Velocity coding. Intraframe BWLZH algorithm. High accuracy. Cubic cell." +#define FILENAME "test56.tng" +#define ALGOTEST +#define NATOMS 100000 +#define CHUNKY 10 +#define SCALE 0.5 +#define PRECISION 1e-8 +#define WRITEVEL 1 +#define VELINTMUL 100000 +#define VELPRECISION 1e-8 +#define INITIALCODING 3 +#define INITIALCODINGPARAMETER -1 +#define CODING 5 +#define CODINGPARAMETER 0 +#define VELCODING 9 +#define VELCODINGPARAMETER 0 +#define INTMIN1 0 +#define INTMIN2 0 +#define INTMIN3 0 +#define INTMAX1 805306368 +#define INTMAX2 805306368 +#define INTMAX3 805306368 +#define NFRAMES 20 +#define EXPECTED_FILESIZE 23390767. diff --git a/src/tests/compression/test57.h b/src/tests/compression/test57.h new file mode 100644 index 0000000..2788e98 --- /dev/null +++ b/src/tests/compression/test57.h @@ -0,0 +1,24 @@ +#define TESTNAME "Velocity coding. Interframe BWLZH algorithm. High accuracy. Cubic cell." +#define FILENAME "test57.tng" +#define ALGOTEST +#define NATOMS 100000 +#define CHUNKY 10 +#define SCALE 0.5 +#define PRECISION 1e-8 +#define WRITEVEL 1 +#define VELINTMUL 100000 +#define VELPRECISION 1e-8 +#define INITIALCODING 3 +#define INITIALCODINGPARAMETER -1 +#define CODING 5 +#define CODINGPARAMETER 0 +#define VELCODING 8 +#define VELCODINGPARAMETER 0 +#define INTMIN1 0 +#define INTMIN2 0 +#define INTMIN3 0 +#define INTMAX1 805306368 +#define INTMAX2 805306368 +#define INTMAX3 805306368 +#define NFRAMES 20 +#define EXPECTED_FILESIZE 13817974. diff --git a/src/tests/compression/test6.h b/src/tests/compression/test6.h new file mode 100644 index 0000000..cd3cfee --- /dev/null +++ b/src/tests/compression/test6.h @@ -0,0 +1,23 @@ +#define TESTNAME "Coding. XTC2 algorithm. Cubic cell." +#define FILENAME "test6.tng" +#define ALGOTEST +#define NATOMS 1000 +#define CHUNKY 100 +#define SCALE 0.1 +#define PRECISION 0.01 +#define WRITEVEL 0 +#define VELPRECISION 0.1 +#define INITIALCODING 5 +#define INITIALCODINGPARAMETER 0 +#define CODING 5 +#define CODINGPARAMETER 0 +#define VELCODING 0 +#define VELCODINGPARAMETER 0 +#define INTMIN1 0 +#define INTMIN2 0 +#define INTMIN3 0 +#define INTMAX1 10000 +#define INTMAX2 10000 +#define INTMAX3 10000 +#define NFRAMES 1000 +#define EXPECTED_FILESIZE 2736662. diff --git a/src/tests/compression/test7.h b/src/tests/compression/test7.h new file mode 100644 index 0000000..90a2657 --- /dev/null +++ b/src/tests/compression/test7.h @@ -0,0 +1,23 @@ +#define TESTNAME "Coding. Stopbit interframe algorithm. Cubic cell." +#define FILENAME "test7.tng" +#define ALGOTEST +#define NATOMS 1000 +#define CHUNKY 100 +#define SCALE 0.1 +#define PRECISION 0.01 +#define WRITEVEL 0 +#define VELPRECISION 0.1 +#define INITIALCODING 5 +#define INITIALCODINGPARAMETER 0 +#define CODING 1 +#define CODINGPARAMETER -1 +#define VELCODING 0 +#define VELCODINGPARAMETER 0 +#define INTMIN1 0 +#define INTMIN2 0 +#define INTMIN3 0 +#define INTMAX1 10000 +#define INTMAX2 10000 +#define INTMAX3 10000 +#define NFRAMES 1000 +#define EXPECTED_FILESIZE 2545049. diff --git a/src/tests/compression/test8.h b/src/tests/compression/test8.h new file mode 100644 index 0000000..5a21053 --- /dev/null +++ b/src/tests/compression/test8.h @@ -0,0 +1,23 @@ +#define TESTNAME "Coding. Stopbit interframe algorithm with intraframe compression as initial. Cubic cell." +#define FILENAME "test8.tng" +#define ALGOTEST +#define NATOMS 1000 +#define CHUNKY 100 +#define SCALE 0.1 +#define PRECISION 0.01 +#define WRITEVEL 0 +#define VELPRECISION 0.1 +#define INITIALCODING 3 +#define INITIALCODINGPARAMETER -1 +#define CODING 1 +#define CODINGPARAMETER -1 +#define VELCODING 0 +#define VELCODINGPARAMETER 0 +#define INTMIN1 0 +#define INTMIN2 0 +#define INTMIN3 0 +#define INTMAX1 10000 +#define INTMAX2 10000 +#define INTMAX3 10000 +#define NFRAMES 1000 +#define EXPECTED_FILESIZE 2544876. diff --git a/src/tests/compression/test9.h b/src/tests/compression/test9.h new file mode 100644 index 0000000..1247beb --- /dev/null +++ b/src/tests/compression/test9.h @@ -0,0 +1,23 @@ +#define TESTNAME "Coding. Triple interframe algorithm. Cubic cell." +#define FILENAME "test9.tng" +#define ALGOTEST +#define NATOMS 1000 +#define CHUNKY 100 +#define SCALE 0.1 +#define PRECISION 0.01 +#define WRITEVEL 0 +#define VELPRECISION 0.1 +#define INITIALCODING 5 +#define INITIALCODINGPARAMETER 0 +#define CODING 2 +#define CODINGPARAMETER -1 +#define VELCODING 0 +#define VELCODINGPARAMETER 0 +#define INTMIN1 0 +#define INTMIN2 0 +#define INTMIN3 0 +#define INTMAX1 10000 +#define INTMAX2 10000 +#define INTMAX3 10000 +#define NFRAMES 1000 +#define EXPECTED_FILESIZE 2418212. diff --git a/src/tests/compression/testsuite.c b/src/tests/compression/testsuite.c new file mode 100644 index 0000000..1c5b613 --- /dev/null +++ b/src/tests/compression/testsuite.c @@ -0,0 +1,657 @@ +/* tng compression routines */ + +/* Only modify testsuite.c + *Then* run testsuite.sh to perform the test. + */ + +#include +#include +#include +#include +#include +#include +#include TESTPARAM + +#define FUDGE 1.1 /* 10% off target precision is acceptable */ + +static void keepinbox(int *val) +{ + while (val[0]>INTMAX1) + val[0]-=(INTMAX1-INTMIN1+1); + while (val[0]INTMAX2) + val[1]-=(INTMAX2-INTMIN2+1); + while (val[1]INTMAX3) + val[2]-=(INTMAX3-INTMIN3+1); + while (val[2]1)) + ifl=1; + else if ((i==1) && (flip) && (length>1)) + ifl=0; + target[ifl*3]=base[0]+(intsin((i+iframe)*direction[0])*scale)/256; + target[ifl*3+1]=base[1]+(intcos((i+iframe)*direction[1])*scale)/256; + target[ifl*3+2]=base[2]+(intcos((i+iframe)*direction[2])*scale)/256; + keepinbox(target+ifl*3); + } +} + +#ifndef FRAMESCALE +#define FRAMESCALE 1 +#endif + +static void genibox(int *intbox, int iframe) +{ + int molecule_length=1; + int molpos[3]; + int direction[3]={1,1,1}; + int scale=1; + int flip=0; + int i=0; + molpos[0]=intsin(iframe*FRAMESCALE)/32; + molpos[1]=1+intcos(iframe*FRAMESCALE)/32; + molpos[2]=2+intsin(iframe*FRAMESCALE)/16; + keepinbox(molpos); + while (iNATOMS) + this_mol_length=NATOMS-i; + /* We must test the large rle as well. This requires special + sequencies to get triggered. So insert these from time to + time */ +#ifndef REGULAR + if ((i%10)==0) + { + int j; + intbox[i*3]=molpos[0]; + intbox[i*3+1]=molpos[1]; + intbox[i*3+2]=molpos[2]; + for (j=1; j5) + scale=1; + + molecule_length++; + if (molecule_length>30) + molecule_length=1; + if (i%9) + flip=1-flip; + } +} + +static void genivelbox(int *intvelbox, int iframe) +{ + int i; + for (i=0; imaxdiff) + maxdiff=(double)fabs(arr1[i*stride1+j]-arr2[i*stride2+j]); + } +#if 0 + for (i=0; iprec*0.5*FUDGE) + { + return 0; + } + else + return 1; +} + +struct tng_file +{ + FILE *f; + int natoms; + int chunky; + double precision; + double velprecision; + int initial_coding; + int initial_coding_parameter; + int coding; + int coding_parameter; + int initial_velcoding; + int initial_velcoding_parameter; + int velcoding; + int velcoding_parameter; + int speed; + int nframes; + int nframes_delivered; + int writevel; + double *pos; + double *vel; +}; + +static size_t fwrite_int_le(int *x,FILE *f) +{ + unsigned char c[4]; + unsigned int i=(unsigned int)*x; + c[0]=(unsigned char)(i&0xFFU); + c[1]=(unsigned char)((i>>8)&0xFFU); + c[2]=(unsigned char)((i>>16)&0xFFU); + c[3]=(unsigned char)((i>>24)&0xFFU); + return fwrite(c,1,4,f); +} + +static size_t fread_int_le(int *x,FILE *f) +{ + unsigned char c[4]; + unsigned int i; + size_t n=fread(c,1,4,f); + if (n) + { + i=(((unsigned int)c[3])<<24)|(((unsigned int)c[2])<<16)|(((unsigned int)c[1])<<8)|((unsigned int)c[0]); + *x=(int)i; + } + return n; +} + +static struct tng_file *open_tng_file_write(char *filename, + int natoms,int chunky, + double precision, + int writevel, + double velprecision, + int initial_coding, + int initial_coding_parameter, + int coding, + int coding_parameter, + int initial_velcoding, + int initial_velcoding_parameter, + int velcoding, + int velcoding_parameter, + int speed) +{ + struct tng_file *tng_file=malloc(sizeof *tng_file); + tng_file->pos=NULL; + tng_file->vel=NULL; + tng_file->nframes=0; + tng_file->chunky=chunky; + tng_file->precision=precision; + tng_file->natoms=natoms; + tng_file->writevel=writevel; + tng_file->velprecision=velprecision; + tng_file->initial_coding=initial_coding; + tng_file->initial_coding_parameter=initial_coding_parameter; + tng_file->coding=coding; + tng_file->coding_parameter=coding_parameter; + tng_file->initial_velcoding=initial_velcoding; + tng_file->initial_velcoding_parameter=initial_velcoding_parameter; + tng_file->velcoding=velcoding; + tng_file->velcoding_parameter=velcoding_parameter; + tng_file->speed=speed; + tng_file->pos=malloc(natoms*chunky*3*sizeof *tng_file->pos); + tng_file->f=fopen(filename,"wb"); + if (writevel) + tng_file->vel=malloc(natoms*chunky*3*sizeof *tng_file->vel); + fwrite_int_le(&natoms,tng_file->f); + return tng_file; +} + +static void flush_tng_frames(struct tng_file *tng_file) +{ + int algo[4]; + char *buf; + int nitems; + fwrite_int_le(&tng_file->nframes,tng_file->f); + algo[0]=tng_file->initial_coding; + algo[1]=tng_file->initial_coding_parameter; + algo[2]=tng_file->coding; + algo[3]=tng_file->coding_parameter; + buf=tng_compress_pos(tng_file->pos, + tng_file->natoms, + tng_file->nframes, + tng_file->precision, + tng_file->speed,algo,&nitems); + tng_file->initial_coding=algo[0]; + tng_file->initial_coding_parameter=algo[1]; + tng_file->coding=algo[2]; + tng_file->coding_parameter=algo[3]; + fwrite_int_le(&nitems,tng_file->f); + fwrite(buf,1,nitems,tng_file->f); + free(buf); + if (tng_file->writevel) + { + algo[0]=tng_file->initial_velcoding; + algo[1]=tng_file->initial_velcoding_parameter; + algo[2]=tng_file->velcoding; + algo[3]=tng_file->velcoding_parameter; + buf=tng_compress_vel(tng_file->vel, + tng_file->natoms, + tng_file->nframes, + tng_file->velprecision, + tng_file->speed,algo,&nitems); + tng_file->initial_velcoding=algo[0]; + tng_file->initial_velcoding_parameter=algo[1]; + tng_file->velcoding=algo[2]; + tng_file->velcoding_parameter=algo[3]; + fwrite_int_le(&nitems,tng_file->f); + fwrite(buf,1,nitems,tng_file->f); + free(buf); + } + tng_file->nframes=0; +} + +static void write_tng_file(struct tng_file *tng_file, + double *pos,double *vel) +{ + memcpy(tng_file->pos+tng_file->nframes*tng_file->natoms*3,pos,tng_file->natoms*3*sizeof *tng_file->pos); + if (tng_file->writevel) + memcpy(tng_file->vel+tng_file->nframes*tng_file->natoms*3,vel,tng_file->natoms*3*sizeof *tng_file->vel); + tng_file->nframes++; + if (tng_file->nframes==tng_file->chunky) + flush_tng_frames(tng_file); +} + +static void close_tng_file_write(struct tng_file *tng_file) +{ + if (tng_file->nframes) + flush_tng_frames(tng_file); + fclose(tng_file->f); + free(tng_file->pos); + free(tng_file->vel); + free(tng_file); +} + +static struct tng_file *open_tng_file_read(char *filename, int writevel) +{ + struct tng_file *tng_file=malloc(sizeof *tng_file); + tng_file->pos=NULL; + tng_file->vel=NULL; + tng_file->f=fopen(filename,"rb"); + tng_file->nframes=0; + tng_file->nframes_delivered=0; + tng_file->writevel=writevel; + fread_int_le(&tng_file->natoms,tng_file->f); + return tng_file; +} + +static int read_tng_file(struct tng_file *tng_file, + double *pos, + double *vel) +{ + if (tng_file->nframes==tng_file->nframes_delivered) + { + int nitems; + char *buf; + free(tng_file->pos); + free(tng_file->vel); + if (!fread_int_le(&tng_file->nframes,tng_file->f)) + return 1; + if (!fread_int_le(&nitems,tng_file->f)) + return 1; + buf=malloc(nitems); + if (!fread(buf,1,nitems,tng_file->f)) + return 1; + tng_file->pos=malloc(tng_file->natoms*tng_file->nframes*3*sizeof *tng_file->pos); + if (tng_file->writevel) + tng_file->vel=malloc(tng_file->natoms*tng_file->nframes*3*sizeof *tng_file->vel); +#if 0 + { + int natoms, nframes, algo[4]; + double precision; + int ivel; + char *initial_coding, *coding; + tng_compress_inquire(buf,&ivel,&natoms,&nframes,&precision,algo); + initial_coding=tng_compress_initial_pos_algo(algo); + coding=tng_compress_pos_algo(algo); + printf("ivel=%d natoms=%d nframes=%d precision=%g initial pos=%s pos=%s\n",ivel,natoms,nframes,precision,initial_coding,coding); + } +#endif + tng_compress_uncompress(buf,tng_file->pos); + free(buf); + if (tng_file->writevel) + { + if (!fread_int_le(&nitems,tng_file->f)) + return 1; + buf=malloc(nitems); + if (!fread(buf,1,nitems,tng_file->f)) + return 1; +#if 0 + { + int natoms, nframes, algo[4]; + double precision; + int ivel; + char *initial_coding, *coding; + tng_compress_inquire(buf,&ivel,&natoms,&nframes,&precision,algo); + initial_coding=tng_compress_initial_vel_algo(algo); + coding=tng_compress_vel_algo(algo); + printf("ivel=%d natoms=%d nframes=%d precision=%g initial vel=%s vel=%s\n",ivel,natoms,nframes,precision,initial_coding,coding); + } +#endif + tng_compress_uncompress(buf,tng_file->vel); + free(buf); + } + tng_file->nframes_delivered=0; + } + memcpy(pos,tng_file->pos+tng_file->nframes_delivered*tng_file->natoms*3,tng_file->natoms*3*sizeof *pos); + if (tng_file->writevel) + memcpy(vel,tng_file->vel+tng_file->nframes_delivered*tng_file->natoms*3,tng_file->natoms*3*sizeof *vel); + tng_file->nframes_delivered++; + return 0; +} + +static void close_tng_file_read(struct tng_file *tng_file) +{ + free(tng_file->vel); + free(tng_file->pos); + fclose(tng_file->f); + free(tng_file); +} + + + +#ifndef EXPECTED_FILESIZE +#define EXPECTED_FILESIZE 1 +#endif + +#ifndef INITIALVELCODING +#define INITIALVELCODING -1 +#endif +#ifndef INITIALVELCODINGPARAMETER +#define INITIALVELCODINGPARAMETER -1 +#endif + +#ifndef SPEED +#define SPEED 5 +#endif + +/* Return value 1 means file error. + Return value 4 means coding error in coordinates. + Return value 5 means coding error in velocities. + Return value 9 means filesize seems too off. + + Return value 100+ means test specific error. + */ +static int algotest() +{ + int i; + int *intbox=warnmalloc(NATOMS*3*sizeof *intbox); + int *intvelbox=warnmalloc(NATOMS*3*sizeof *intvelbox); + double *box1=warnmalloc(NATOMS*STRIDE1*sizeof *box1); + double *velbox1=warnmalloc(NATOMS*STRIDE1*sizeof *velbox1); + double time1, lambda1; + double H1[9]; + int startframe=0; + int endframe=NFRAMES; +#ifdef GEN + FILE *file; + double filesize; +#else + int i2; + int readreturn; + double H2[9]; + double time2, lambda2; + double *box2=warnmalloc(NATOMS*STRIDE2*sizeof *box2); + double *velbox2=warnmalloc(NATOMS*STRIDE2*sizeof *velbox2); +#endif +#ifdef GEN + void *dumpfile=open_tng_file_write(FILENAME,NATOMS,CHUNKY, + PRECISION,WRITEVEL,VELPRECISION, + INITIALCODING, + INITIALCODINGPARAMETER,CODING,CODINGPARAMETER, + INITIALVELCODING,INITIALVELCODINGPARAMETER, + VELCODING,VELCODINGPARAMETER,SPEED); +#else + void *dumpfile=open_tng_file_read(FILENAME,WRITEVEL); +#endif + if (!dumpfile) + return 1; + for (i=0; i<9; i++) + H1[i]=0.; + H1[0]=INTMAX1*PRECISION*SCALE; + H1[4]=INTMAX2*PRECISION*SCALE; + H1[8]=INTMAX3*PRECISION*SCALE; + for (i=startframe; i0) + { + if ((fabs(filesize-EXPECTED_FILESIZE)/EXPECTED_FILESIZE)>0.05) + return 9; + } +#endif + return 0; +} + +int main() +{ + int testval; + if (sizeof(int)<4) + { + fprintf(stderr,"ERROR: sizeof(int) is too small: %d<4\n",(int)sizeof(int)); + exit(EXIT_FAILURE); + } +#ifdef GEN + printf("Tng compress testsuite generating test: %s\n",TESTNAME); +#else + printf("Tng compress testsuite running test: %s\n",TESTNAME); +#endif + testval=algotest(); + if (testval==0) + printf("Passed.\n"); + else if (testval==1) + { + printf("ERROR: File error.\n"); + exit(EXIT_FAILURE); + } + else if (testval==4) + { + printf("ERROR: Read coding error in coordinates.\n"); + exit(EXIT_FAILURE); + } + else if (testval==5) + { + printf("ERROR: Read coding error in velocities.\n"); + exit(EXIT_FAILURE); + } + else if (testval==9) + { + printf("ERROR: Generated filesize differs too much.\n"); + exit(EXIT_FAILURE); + } + else + { + printf("ERROR: Unknown error.\n"); + exit(EXIT_FAILURE); + } + return 0; +} diff --git a/src/tests/compression/testsuite.sh b/src/tests/compression/testsuite.sh new file mode 100755 index 0000000..da21de3 --- /dev/null +++ b/src/tests/compression/testsuite.sh @@ -0,0 +1,33 @@ +#!/bin/sh +do_write_test="Yes" +if [ -n "$1" ]; then + do_write_test="" +fi +STARTTEST=1 +ENDTEST=57 +#CFLAGS="-O2 -Wall" +CFLAGS="-O2 -g" +LIBS="-lm" +#CFLAGS="-O0 -Wall -g" +#LIBS="-lm -lefence" +CC="gcc" +# 32 bit +#CC="gcc -m32" +for testnum in $(seq $STARTTEST $ENDTEST); do + testname=$(grep "TESTNAME" test$testnum.h|sed 's/#define TESTNAME//') + sed "s/TESTPARAM/\"test$testnum.h\"/" test$testnum.c + if [ -n "$do_write_test" ]; then + echo Write test $testnum: $testname + $CC -DGEN $CFLAGS -I../ -L../ -o gen$testnum test$testnum.c -ltng_compress $LIBS + ./gen$testnum + rm -f gen$testnum + fi + echo Read test $testnum: $testname + $CC $CFLAGS -I../ -L../ -o read$testnum test$testnum.c -ltng_compress $LIBS + ./read$testnum + rm -f read$testnum + rm -f test$testnum.c +done + + + -- cgit v0.10.1