diff options
Diffstat (limited to 'src/compression/tng_compress.c')
| -rw-r--r-- | src/compression/tng_compress.c | 1545 |
1 files changed, 1545 insertions, 0 deletions
diff --git a/src/compression/tng_compress.c b/src/compression/tng_compress.c new file mode 100644 index 0000000..e92f24d --- /dev/null +++ b/src/compression/tng_compress.c @@ -0,0 +1,1545 @@ +/* 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 <stdlib.h> +#include <string.h> +#include <math.h> +#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<nframes; iframe++) + for (i=0; i<natoms; i++) + for (j=0; j<3; j++) + quant[iframe*natoms*3+i*3+j]=(int)floor((x[iframe*natoms*3+i*3+j]/precision)+0.5); +#if 0 + printf("Q precision=%g\n",precision); + for (iframe=0; iframe<nframes; iframe++) + for (i=0; i<natoms; i++) + for (j=0; j<3; j++) + printf("Q: %d %d %d: %d %g\n",iframe,i,j,quant[iframe*natoms*3+i*3+j],x[iframe*natoms*3+i*3+j]); +#endif +} + +static void quant_inter_differences(int *quant, int natoms, int nframes, + int *quant_inter) +{ + int iframe, i, j; + /* The first frame is used for absolute positions. */ + for (i=0; i<natoms; i++) + for (j=0; j<3; j++) + quant_inter[i*3+j]=quant[i*3+j]; + /* For all other frames, the difference to the previous frame is used. */ + for (iframe=1; iframe<nframes; iframe++) + for (i=0; i<natoms; i++) + for (j=0; j<3; j++) + quant_inter[iframe*natoms*3+i*3+j]=quant[iframe*natoms*3+i*3+j]-quant[(iframe-1)*natoms*3+i*3+j]; +} + +static void quant_intra_differences(int *quant, int natoms, int nframes, + int *quant_intra) +{ + int iframe, i, j; + for (iframe=0; iframe<nframes; iframe++) + { + /* The first atom is used with its absolute position. */ + for (j=0; j<3; j++) + quant_intra[iframe*natoms*3+j]=quant[iframe*natoms*3+j]; + /* For all other atoms the intraframe differences are computed. */ + for (i=1; i<natoms; i++) + for (j=0; j<3; j++) + quant_intra[iframe*natoms*3+i*3+j]=quant[iframe*natoms*3+i*3+j]-quant[iframe*natoms*3+(i-1)*3+j]; + } +} + +static void unquantize(double *x, int natoms, int nframes, + double precision, + int *quant) +{ + int iframe, i, j; + for (iframe=0; iframe<nframes; iframe++) + for (i=0; i<natoms; i++) + for (j=0; j<3; j++) + x[iframe*natoms*3+i*3+j]=(double)quant[iframe*natoms*3+i*3+j]*precision; +} + +static void unquantize_inter_differences(double *x, int natoms, int nframes, + double precision, + int *quant) +{ + int iframe, i, j; + for (i=0; i<natoms; i++) + for (j=0; j<3; j++) + { + int q=quant[i*3+j]; /* First value. */ + x[i*3+j]=(double)q*precision; + for (iframe=1; iframe<nframes; iframe++) + { + q+=quant[iframe*natoms*3+i*3+j]; + x[iframe*natoms*3+i*3+j]=(double)q*precision; + } + } +} + +/* In frame update required for the initial frame if intra-frame + compression was used. */ +static void unquant_intra_differences_first_frame(int *quant, int natoms) +{ + int i,j; + for (j=0; j<3; j++) + { + int q=quant[j]; + for (i=1; i<natoms; i++) + { + q+=quant[i*3+j]; + quant[i*3+j]=q; + } + } +#if 0 + for (j=0; j<3; j++) + for (i=0; i<natoms; i++) + { + printf("UQ: %d %d %d: %d\n",0,j,i,quant[i*3+j]); + } +#endif +} + +static void unquantize_intra_differences(double *x, int natoms, int nframes, + double precision, + int *quant) +{ + int iframe, i, j; +#if 0 + printf("UQ precision=%g\n",precision); +#endif + for (iframe=0; iframe<nframes; iframe++) + for (j=0; j<3; j++) + { + int q=quant[iframe*natoms*3+j]; + x[iframe*natoms*3+j]=(double)q*precision; + for (i=1; i<natoms; i++) + { + q+=quant[iframe*natoms*3+i*3+j]; + x[iframe*natoms*3+i*3+j]=(double)q*precision; + } + } +} + +/* Buffer num 8 bit bytes into buffer location buf */ +static void bufferfix(unsigned char *buf, fix_t v, int num) +{ + /* Store in little endian format. */ + unsigned char c; /* at least 8 bits. */ + c=(unsigned char)(v & 0xFFU); + while (num--) + { + *buf++=c; + v >>= 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)<<shift; + shift+=8; + } while (--num); + return f; +} + +/* Perform position compression from the quantized data. */ +static void compress_quantized_pos(int *quant, int *quant_inter, int *quant_intra, + 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_POS,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_POS_XTC2) || + (initial_coding==TNG_COMPRESS_ALGO_POS_TRIPLET_ONETOONE) || + (initial_coding==TNG_COMPRESS_ALGO_POS_XTC3)) + { + 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); + } + else if ((initial_coding==TNG_COMPRESS_ALGO_POS_TRIPLET_INTRA) || + (initial_coding==TNG_COMPRESS_ALGO_POS_BWLZH_INTRA)) + { + struct coder *coder=Ptngc_coder_init(); + length=natoms*3; + datablock=(char*)Ptngc_pack_array(coder,quant_intra,&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_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<best_length)) + { + new_parameter=bits; + best_length=io_length; + } + free(packed); + } + } + if (new_parameter==-1) + return 1; + + *coding_parameter=new_parameter; + *length=best_length; + return 0; +} + +static int determine_best_coding_triple(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_TRIPLET,bits,natoms,0); + if (packed) + { + if ((new_parameter==-1) || (io_length<best_length)) + { + new_parameter=bits; + best_length=io_length; + } + free(packed); + } + } + if (new_parameter==-1) + return 1; + + *coding_parameter=new_parameter; + *length=best_length; + return 0; +} + +static void determine_best_pos_initial_coding(int *quant, int *quant_intra, int natoms, int speed, + fix_t prec_hi, fix_t prec_lo, + int *initial_coding, int *initial_coding_parameter) +{ + if (*initial_coding==-1) + { + /* Determine all parameters automatically */ + int best_coding; + int best_coding_parameter; + int best_code_size; + int current_coding; + int current_coding_parameter; + int current_code_size; + struct coder *coder; + /* Start with XTC2, it should always work. */ + current_coding=TNG_COMPRESS_ALGO_POS_XTC2; + 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); + best_coding=current_coding; + best_coding_parameter=current_coding_parameter; + best_code_size=current_code_size; + + /* Determine best parameter for triplet intra. */ + current_coding=TNG_COMPRESS_ALGO_POS_TRIPLET_INTRA; + coder=Ptngc_coder_init(); + current_code_size=natoms*3; + current_coding_parameter=0; + if (!determine_best_coding_triple(coder,quant_intra,¤t_code_size,¤t_coding_parameter,natoms)) + { + if (current_code_size<best_code_size) + { + best_coding=current_coding; + best_coding_parameter=current_coding_parameter; + best_code_size=current_code_size; + } + } + Ptngc_coder_deinit(coder); + + /* Determine best parameter for triplet one-to-one. */ + current_coding=TNG_COMPRESS_ALGO_POS_TRIPLET_ONETOONE; + coder=Ptngc_coder_init(); + current_code_size=natoms*3; + current_coding_parameter=0; + if (!determine_best_coding_triple(coder,quant,¤t_code_size,¤t_coding_parameter,natoms)) + { + if (current_code_size<best_code_size) + { + best_coding=current_coding; + best_coding_parameter=current_coding_parameter; + best_code_size=current_code_size; + } + } + Ptngc_coder_deinit(coder); + + if (speed>=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<best_code_size) + { + best_coding=current_coding; + best_coding_parameter=current_coding_parameter; + best_code_size=current_code_size; + } + } + /* Test BWLZH intra */ + if (speed>=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<best_code_size) + { + best_coding=current_coding; + best_coding_parameter=current_coding_parameter; + best_code_size=current_code_size; + } + } + *initial_coding=best_coding; + *initial_coding_parameter=best_coding_parameter; + } + else + { + if (*initial_coding_parameter==-1) + { + if ((*initial_coding==TNG_COMPRESS_ALGO_POS_XTC2) || + (*initial_coding==TNG_COMPRESS_ALGO_POS_XTC3) || + (*initial_coding==TNG_COMPRESS_ALGO_POS_BWLZH_INTRA)) + *initial_coding_parameter=0; + else if (*initial_coding==TNG_COMPRESS_ALGO_POS_TRIPLET_INTRA) + { + struct coder *coder=Ptngc_coder_init(); + int current_code_size=natoms*3; + determine_best_coding_triple(coder,quant_intra,¤t_code_size,initial_coding_parameter,natoms); + Ptngc_coder_deinit(coder); + } + else if (*initial_coding==TNG_COMPRESS_ALGO_POS_TRIPLET_ONETOONE) + { + struct coder *coder=Ptngc_coder_init(); + int current_code_size=natoms*3; + determine_best_coding_triple(coder,quant,¤t_code_size,initial_coding_parameter,natoms); + Ptngc_coder_deinit(coder); + } + } + } +} + +static void determine_best_pos_coding(int *quant, int *quant_inter, int *quant_intra, int natoms, int nframes, int speed, + fix_t prec_hi, fix_t prec_lo, + int *coding, int *coding_parameter) +{ + if (*coding==-1) + { + /* Determine all parameters automatically */ + int best_coding; + int best_coding_parameter; + int best_code_size; + int current_coding; + int current_coding_parameter; + int current_code_size; + int initial_code_size; + struct coder *coder; + /* Always use XTC2 for the initial coding. */ + compress_quantized_pos(quant,quant_inter,quant_intra,natoms,1,speed, + TNG_COMPRESS_ALGO_POS_XTC2,0, + 0,0, + prec_hi,prec_lo,&initial_code_size,NULL); + /* Start with XTC2, it should always work. */ + current_coding=TNG_COMPRESS_ALGO_POS_XTC2; + 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); + best_coding=current_coding; + best_coding_parameter=current_coding_parameter; + best_code_size=current_code_size-initial_code_size; /* Correct for the use of XTC2 for the first frame. */ + + /* Determine best parameter for stopbit interframe coding. */ + current_coding=TNG_COMPRESS_ALGO_POS_STOPBIT_INTER; + coder=Ptngc_coder_init(); + current_code_size=natoms*3*(nframes-1); + current_coding_parameter=0; + if (!determine_best_coding_stop_bits(coder,quant_inter+natoms*3,¤t_code_size, + ¤t_coding_parameter,natoms)) + { + if (current_code_size<best_code_size) + { + best_coding=current_coding; + best_coding_parameter=current_coding_parameter; + best_code_size=current_code_size; + } + } + Ptngc_coder_deinit(coder); + + /* Determine best parameter for triplet interframe coding. */ + current_coding=TNG_COMPRESS_ALGO_POS_TRIPLET_INTER; + coder=Ptngc_coder_init(); + current_code_size=natoms*3*(nframes-1); + current_coding_parameter=0; + if (!determine_best_coding_triple(coder,quant_inter+natoms*3,¤t_code_size, + ¤t_coding_parameter,natoms)) + { + if (current_code_size<best_code_size) + { + best_coding=current_coding; + best_coding_parameter=current_coding_parameter; + best_code_size=current_code_size; + } + } + Ptngc_coder_deinit(coder); + + /* Determine best parameter for triplet intraframe coding. */ + current_coding=TNG_COMPRESS_ALGO_POS_TRIPLET_INTRA; + coder=Ptngc_coder_init(); + current_code_size=natoms*3*(nframes-1); + current_coding_parameter=0; + if (!determine_best_coding_triple(coder,quant_intra+natoms*3,¤t_code_size, + ¤t_coding_parameter,natoms)) + { + if (current_code_size<best_code_size) + { + best_coding=current_coding; + best_coding_parameter=current_coding_parameter; + best_code_size=current_code_size; + } + } + Ptngc_coder_deinit(coder); + + /* Determine best parameter for triplet one-to-one coding. */ + current_coding=TNG_COMPRESS_ALGO_POS_TRIPLET_ONETOONE; + coder=Ptngc_coder_init(); + current_code_size=natoms*3*(nframes-1); + current_coding_parameter=0; + if (!determine_best_coding_triple(coder,quant+natoms*3,¤t_code_size, + ¤t_coding_parameter,natoms)) + { + if (current_code_size<best_code_size) + { + best_coding=current_coding; + best_coding_parameter=current_coding_parameter; + best_code_size=current_code_size; + } + } + Ptngc_coder_deinit(coder); + + /* Test BWLZH inter */ + if (speed>=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<best_code_size) + { + best_coding=current_coding; + best_coding_parameter=current_coding_parameter; + best_code_size=current_code_size; + } + } + + /* Test BWLZH intra */ + if (speed>=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<best_code_size) + { + best_coding=current_coding; + best_coding_parameter=current_coding_parameter; + best_code_size=current_code_size; + } + } + *coding=best_coding; + *coding_parameter=best_coding_parameter; + } + else if (*coding_parameter==-1) + { + if ((*coding==TNG_COMPRESS_ALGO_POS_XTC2) || + (*coding==TNG_COMPRESS_ALGO_POS_XTC3) || + (*coding==TNG_COMPRESS_ALGO_POS_BWLZH_INTER) || + (*coding==TNG_COMPRESS_ALGO_POS_BWLZH_INTRA)) + *coding_parameter=0; + else if (*coding==TNG_COMPRESS_ALGO_POS_STOPBIT_INTER) + { + struct coder *coder=Ptngc_coder_init(); + int current_code_size=natoms*3*(nframes-1); + determine_best_coding_stop_bits(coder,quant_inter+natoms*3,¤t_code_size, + coding_parameter,natoms); + Ptngc_coder_deinit(coder); + } + else if (*coding==TNG_COMPRESS_ALGO_POS_TRIPLET_INTER) + { + struct coder *coder=Ptngc_coder_init(); + int current_code_size=natoms*3*(nframes-1); + determine_best_coding_triple(coder,quant_inter+natoms*3,¤t_code_size, + coding_parameter,natoms); + Ptngc_coder_deinit(coder); + } + else if (*coding==TNG_COMPRESS_ALGO_POS_TRIPLET_INTRA) + { + struct coder *coder=Ptngc_coder_init(); + int current_code_size=natoms*3*(nframes-1); + determine_best_coding_triple(coder,quant_intra+natoms*3,¤t_code_size, + coding_parameter,natoms); + Ptngc_coder_deinit(coder); + } + else if (*coding==TNG_COMPRESS_ALGO_POS_TRIPLET_ONETOONE) + { + struct coder *coder=Ptngc_coder_init(); + int current_code_size=natoms*3*(nframes-1); + determine_best_coding_triple(coder,quant+natoms*3,¤t_code_size, + coding_parameter,natoms); + Ptngc_coder_deinit(coder); + } + } +} + +static void determine_best_vel_initial_coding(int *quant, int natoms, int speed, + fix_t prec_hi, fix_t prec_lo, + int *initial_coding, int *initial_coding_parameter) +{ + if (*initial_coding==-1) + { + /* Determine all parameters automatically */ + int best_coding=-1; + int best_coding_parameter=-1; + int best_code_size=-1; + int current_coding; + int current_coding_parameter; + int current_code_size; + struct coder *coder; + /* Start to determine best parameter for stopbit one-to-one. */ + current_coding=TNG_COMPRESS_ALGO_VEL_STOPBIT_ONETOONE; + current_code_size=natoms*3; + current_coding_parameter=0; + coder=Ptngc_coder_init(); + if (!determine_best_coding_stop_bits(coder,quant,¤t_code_size, + ¤t_coding_parameter,natoms)) + { + best_coding=current_coding; + best_coding_parameter=current_coding_parameter; + best_code_size=current_code_size; + } + Ptngc_coder_deinit(coder); + + /* Determine best parameter for triplet one-to-one. */ + current_coding=TNG_COMPRESS_ALGO_VEL_TRIPLET_ONETOONE; + coder=Ptngc_coder_init(); + current_code_size=natoms*3; + current_coding_parameter=0; + if (!determine_best_coding_triple(coder,quant,¤t_code_size,¤t_coding_parameter,natoms)) + { + if ((best_coding==-1) || (current_code_size<best_code_size)) + { + best_coding=current_coding; + best_coding_parameter=current_coding_parameter; + best_code_size=current_code_size; + } + } + Ptngc_coder_deinit(coder); + + /* Test BWLZH one-to-one */ + if (speed>=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<best_code_size)) + { + best_coding=current_coding; + best_coding_parameter=current_coding_parameter; + best_code_size=current_code_size; + } + } + *initial_coding=best_coding; + *initial_coding_parameter=best_coding_parameter; + } + else if (*initial_coding_parameter==-1) + { + if (*initial_coding==TNG_COMPRESS_ALGO_VEL_BWLZH_ONETOONE) + *initial_coding_parameter=0; + else if (*initial_coding==TNG_COMPRESS_ALGO_VEL_STOPBIT_ONETOONE) + { + struct coder *coder=Ptngc_coder_init(); + int current_code_size=natoms*3; + determine_best_coding_stop_bits(coder,quant,¤t_code_size, + initial_coding_parameter,natoms); + Ptngc_coder_deinit(coder); + } + else if (*initial_coding==TNG_COMPRESS_ALGO_VEL_TRIPLET_ONETOONE) + { + struct coder *coder=Ptngc_coder_init(); + int current_code_size=natoms*3; + determine_best_coding_triple(coder,quant,¤t_code_size,initial_coding_parameter,natoms); + Ptngc_coder_deinit(coder); + } + } +} + +static void determine_best_vel_coding(int *quant, int *quant_inter, int natoms, int nframes, int speed, + fix_t prec_hi, fix_t prec_lo, + int *coding, int *coding_parameter) +{ + if (*coding==-1) + { + /* Determine all parameters automatically */ + int best_coding; + int best_coding_parameter; + int best_code_size; + int current_coding; + int current_coding_parameter; + int current_code_size; + int initial_code_size; + int initial_numbits=5; + struct coder *coder; + /* Use stopbits one-to-one coding for the initial coding. */ + compress_quantized_vel(quant,NULL,natoms,1,speed, + TNG_COMPRESS_ALGO_VEL_STOPBIT_ONETOONE,initial_numbits, + 0,0,prec_hi,prec_lo,&initial_code_size,NULL); + + /* Test stopbit one-to-one */ + current_coding=TNG_COMPRESS_ALGO_VEL_STOPBIT_ONETOONE; + current_code_size=natoms*3*(nframes-1); + current_coding_parameter=0; + coder=Ptngc_coder_init(); + determine_best_coding_stop_bits(coder,quant+natoms*3,¤t_code_size, + ¤t_coding_parameter,natoms); + Ptngc_coder_deinit(coder); + best_coding=current_coding; + best_code_size=current_code_size; + best_coding_parameter=current_coding_parameter; + + /* Test triplet interframe */ + current_coding=TNG_COMPRESS_ALGO_VEL_TRIPLET_INTER; + current_code_size=natoms*3*(nframes-1); + current_coding_parameter=0; + coder=Ptngc_coder_init(); + if (!determine_best_coding_triple(coder,quant_inter+natoms*3,¤t_code_size, + ¤t_coding_parameter,natoms)) + { + if (current_code_size<best_code_size) + { + best_coding=current_coding; + best_code_size=current_code_size; + best_coding_parameter=current_coding_parameter; + } + } + Ptngc_coder_deinit(coder); + + /* Test triplet one-to-one */ + current_coding=TNG_COMPRESS_ALGO_VEL_TRIPLET_ONETOONE; + current_code_size=natoms*3*(nframes-1); + current_coding_parameter=0; + coder=Ptngc_coder_init(); + if (!determine_best_coding_triple(coder,quant+natoms*3,¤t_code_size, + ¤t_coding_parameter,natoms)) + { + if (current_code_size<best_code_size) + { + best_coding=current_coding; + best_code_size=current_code_size; + best_coding_parameter=current_coding_parameter; + } + } + Ptngc_coder_deinit(coder); + + /* Test stopbit interframe */ + current_coding=TNG_COMPRESS_ALGO_VEL_STOPBIT_INTER; + current_code_size=natoms*3*(nframes-1); + current_coding_parameter=0; + coder=Ptngc_coder_init(); + if (!determine_best_coding_stop_bits(coder,quant_inter+natoms*3,¤t_code_size, + ¤t_coding_parameter,natoms)) + { + if (current_code_size<best_code_size) + { + best_coding=current_coding; + best_code_size=current_code_size; + best_coding_parameter=current_coding_parameter; + } + } + Ptngc_coder_deinit(coder); + + if (speed>=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_size<best_code_size) + { + best_coding=current_coding; + best_code_size=current_code_size; + best_coding_parameter=current_coding_parameter; + } + + /* Test BWLZH one-to-one */ + current_coding=TNG_COMPRESS_ALGO_VEL_BWLZH_ONETOONE; + 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_size<best_code_size) + { + best_coding=current_coding; + best_code_size=current_code_size; + best_coding_parameter=current_coding_parameter; + } + } + *coding=best_coding; + *coding_parameter=best_coding_parameter; + } + else if (*coding_parameter==-1) + { + if ((*coding==TNG_COMPRESS_ALGO_VEL_BWLZH_INTER) || + (*coding==TNG_COMPRESS_ALGO_VEL_BWLZH_ONETOONE)) + *coding_parameter=0; + else if (*coding==TNG_COMPRESS_ALGO_VEL_STOPBIT_ONETOONE) + { + struct coder *coder=Ptngc_coder_init(); + int current_code_size=natoms*3*(nframes-1); + determine_best_coding_stop_bits(coder,quant+natoms*3,¤t_code_size, + coding_parameter,natoms); + Ptngc_coder_deinit(coder); + } + else if (*coding==TNG_COMPRESS_ALGO_VEL_TRIPLET_INTER) + { + struct coder *coder=Ptngc_coder_init(); + int current_code_size=natoms*3*(nframes-1); + determine_best_coding_triple(coder,quant_inter+natoms*3,¤t_code_size, + coding_parameter,natoms); + Ptngc_coder_deinit(coder); + } + else if (*coding==TNG_COMPRESS_ALGO_VEL_TRIPLET_ONETOONE) + { + struct coder *coder=Ptngc_coder_init(); + int current_code_size=natoms*3*(nframes-1); + determine_best_coding_triple(coder,quant+natoms*3,¤t_code_size, + coding_parameter,natoms); + Ptngc_coder_deinit(coder); + } + else if (*coding==TNG_COMPRESS_ALGO_VEL_STOPBIT_INTER) + { + struct coder *coder=Ptngc_coder_init(); + int current_code_size=natoms*3*(nframes-1); + determine_best_coding_stop_bits(coder,quant_inter+natoms*3,¤t_code_size, + coding_parameter,natoms); + Ptngc_coder_deinit(coder); + } + } +} + +char DECLSPECDLLEXPORT *tng_compress_pos(double *pos, 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_intra=malloc(natoms*nframes*3*sizeof *quant_intra); + 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(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]; +} |