src/external/tng_io/src/compression/coder.c

   1 /* This code is part of the tng compression routines.
   2  *
   3  * Written by Daniel Spangberg
   4  * Copyright (c) 2010, 2013, The GROMACS development team.
   5  *
   6  *
   7  * This program is free software; you can redistribute it and/or
   8  * modify it under the terms of the Revised BSD License.
   9  */
  10
  11
  12 #include <stdio.h>
  13 #include <stdlib.h>
  14 #include "../../include/compression/tng_compress.h"
  15 #include "../../include/compression/bwlzh.h"
  16 #include "../../include/compression/coder.h"
  17 #include "../../include/compression/warnmalloc.h"
  18
  19 #ifndef USE_WINDOWS
  20 #if defined(WIN32) || defined(_WIN32) || defined(WIN64) || defined(_WIN64)
  21 #define USE_WINDOWS
  22 #endif /* win32... */
  23 #endif /* not defined USE_WINDOWS */
  24
  25 #ifdef USE_WINDOWS
  26 #define TNG_INLINE __inline
  27 #define TNG_SNPRINTF _snprintf
  28 #else
  29 #define TNG_INLINE inline
  30 #define TNG_SNPRINTF snprintf
  31 #endif
  32
  33 struct coder DECLSPECDLLEXPORT *Ptngc_coder_init(void)
  34 {
  35     struct coder *coder_inst=warnmalloc(sizeof *coder_inst);
  36     coder_inst->pack_temporary_bits=0;
  37     return coder_inst;
  38 }
  39
  40 void DECLSPECDLLEXPORT Ptngc_coder_deinit(struct coder *coder_inst)
  41 {
  42     free(coder_inst);
  43 }
  44
  45 TNG_INLINE void DECLSPECDLLEXPORT Ptngc_out8bits(struct coder *coder_inst, unsigned char **output)
  46 {
  47   int pack_temporary_bits=coder_inst->pack_temporary_bits;
  48   unsigned int pack_temporary=coder_inst->pack_temporary;
  49   while (pack_temporary_bits>=8)
  50     {
  51       unsigned int mask=~(0xFFU<<(pack_temporary_bits-8));
  52       unsigned char out=(unsigned char)(pack_temporary>>(pack_temporary_bits-8));
  53       **output=out;
  54       (*output)++;
  55       pack_temporary_bits-=8;
  56       pack_temporary&=mask;
  57     }
  58   coder_inst->pack_temporary_bits=pack_temporary_bits;
  59   coder_inst->pack_temporary=pack_temporary;
  60 }
  61
  62 void DECLSPECDLLEXPORT Ptngc_write_pattern(struct coder *coder_inst, unsigned int pattern,
  63                          int nbits, unsigned char **output)
  64 {
  65     unsigned int mask1,mask2;
  66     mask1=1;
  67     mask2=1<<(nbits-1);
  68     coder_inst->pack_temporary<<=nbits; /* Make room for new data. */
  69     coder_inst->pack_temporary_bits+=nbits;
  70     while (nbits)
  71     {
  72         if (pattern & mask1)
  73             coder_inst->pack_temporary|=mask2;
  74         nbits--;
  75         mask1<<=1;
  76         mask2>>=1;
  77     }
  78     Ptngc_out8bits(coder_inst,output);
  79 }
  80
  81 /* Write up to 24 bits */
  82 TNG_INLINE void DECLSPECDLLEXPORT Ptngc_writebits(struct coder *coder_inst,
  83                                 unsigned int value, int nbits,
  84                                 unsigned char **output_ptr)
  85 {
  86   /* Make room for the bits. */
  87   coder_inst->pack_temporary<<=nbits;
  88   coder_inst->pack_temporary_bits+=nbits;
  89   coder_inst->pack_temporary|=value;
  90   Ptngc_out8bits(coder_inst,output_ptr);
  91 }
  92
  93 /* Write up to 32 bits */
  94 void DECLSPECDLLEXPORT Ptngc_write32bits(struct coder *coder_inst,unsigned int value,
  95                        int nbits, unsigned char **output_ptr)
  96 {
  97   unsigned int mask;
  98   if (nbits>=8)
  99     mask=0xFFU<<(nbits-8);
 100   else
 101     mask=0xFFU>>(8-nbits);
 102   while (nbits>8)
 103     {
 104       /* Make room for the bits. */
 105       coder_inst->pack_temporary<<=8;
 106       coder_inst->pack_temporary_bits+=8;
 107       coder_inst->pack_temporary|=(value&mask)>>(nbits-8);
 108       Ptngc_out8bits(coder_inst,output_ptr);
 109       nbits-=8;
 110       mask>>=8;
 111     }
 112   if (nbits)
 113     Ptngc_writebits(coder_inst,value&mask,nbits,output_ptr);
 114 }
 115
 116 /* Write "arbitrary" number of bits */
 117 void DECLSPECDLLEXPORT Ptngc_writemanybits(struct coder *coder_inst, unsigned char *value,
 118                          int nbits, unsigned char **output_ptr)
 119 {
 120   int vptr=0;
 121   while (nbits>=24)
 122     {
 123       unsigned int v=((((unsigned int)value[vptr])<<16)|
 124                       (((unsigned int)value[vptr+1])<<8)|
 125                       (((unsigned int)value[vptr+2])));
 126       Ptngc_writebits(coder_inst,v,24,output_ptr);
 127       vptr+=3;
 128       nbits-=24;
 129     }
 130   while (nbits>=8)
 131     {
 132       Ptngc_writebits(coder_inst,(unsigned int)value[vptr],8,output_ptr);
 133       vptr++;
 134       nbits-=8;
 135     }
 136   if (nbits)
 137     {
 138       Ptngc_writebits(coder_inst,(unsigned int)value[vptr],nbits,output_ptr);
 139     }
 140 }
 141
 142 static int write_stop_bit_code(struct coder *coder_inst, unsigned int s,
 143                                unsigned int coding_parameter,
 144                                unsigned char **output)
 145 {
 146   do {
 147     unsigned int extract=~(0xffffffffU<<coding_parameter);
 148     unsigned int this=(s&extract)<<1;
 149     s>>=coding_parameter;
 150     if (s)
 151       {
 152         this|=1U;
 153         coder_inst->stat_overflow++;
 154       }
 155     coder_inst->pack_temporary<<=(coding_parameter+1);
 156     coder_inst->pack_temporary_bits+=coding_parameter+1;
 157     coder_inst->pack_temporary|=this;
 158     Ptngc_out8bits(coder_inst,output);
 159     if (s)
 160       {
 161         coding_parameter>>=1;
 162         if (coding_parameter<1)
 163           coding_parameter=1;
 164       }
 165   } while (s);
 166   coder_inst->stat_numval++;
 167   return 0;
 168 }
 169
 170 static int pack_stopbits_item(struct coder *coder_inst,int item,
 171                               unsigned char **output, int coding_parameter)
 172 {
 173     /* Find this symbol in table. */
 174     int s=0;
 175     if (item>0)
 176         s=1+(item-1)*2;
 177     else if (item<0)
 178         s=2+(-item-1)*2;
 179     return write_stop_bit_code(coder_inst,s,coding_parameter,output);
 180 }
 181
 182 static int pack_triplet(struct coder *coder_inst, unsigned int *s,
 183                         unsigned char **output, int coding_parameter,
 184                         unsigned int max_base, int maxbits)
 185 {
 186   /* Determine base for this triplet. */
 187   unsigned int min_base=1U<<coding_parameter;
 188   unsigned int this_base=min_base;
 189   int i;
 190   unsigned int jbase=0;
 191   unsigned int bits_per_value;
 192   for (i=0; i<3; i++)
 193     while (s[i]>=this_base)
 194       {
 195         this_base*=2;
 196         jbase++;
 197       }
 198   bits_per_value=coding_parameter+jbase;
 199   if (jbase>=3)
 200     {
 201       if (this_base>max_base)
 202         return 1;
 203       bits_per_value=maxbits;
 204       jbase=3;
 205     }
 206   /* 2 bits selects the base */
 207   coder_inst->pack_temporary<<=2;
 208   coder_inst->pack_temporary_bits+=2;
 209   coder_inst->pack_temporary|=jbase;
 210   Ptngc_out8bits(coder_inst,output);
 211   for (i=0; i<3; i++)
 212     Ptngc_write32bits(coder_inst,s[i],bits_per_value,output);
 213   return 0;
 214 }
 215
 216 void DECLSPECDLLEXPORT Ptngc_pack_flush(struct coder *coder_inst,unsigned char **output)
 217 {
 218   /* Zero-fill just enough. */
 219   if (coder_inst->pack_temporary_bits>0)
 220     Ptngc_write_pattern(coder_inst,0,8-coder_inst->pack_temporary_bits,output);
 221 }
 222
 223 unsigned char DECLSPECDLLEXPORT *Ptngc_pack_array(struct coder *coder_inst,
 224                                 int *input, int *length, int coding,
 225                                 int coding_parameter, int natoms, int speed)
 226 {
 227   if ((coding==TNG_COMPRESS_ALGO_BWLZH1) || (coding==TNG_COMPRESS_ALGO_BWLZH2))
 228     {
 229       unsigned char *output=warnmalloc(4+bwlzh_get_buflen(*length));
 230       int i,j,k,n=*length;
 231       unsigned int *pval=warnmalloc(n*sizeof *pval);
 232       int nframes=n/natoms/3;
 233       int cnt=0;
 234       int most_negative=2147483647;
 235       for (i=0; i<n; i++)
 236         if (input[i]<most_negative)
 237           most_negative=input[i];
 238       most_negative=-most_negative;
 239       output[0]=((unsigned int)most_negative)&0xFFU;
 240       output[1]=(((unsigned int)most_negative)>>8)&0xFFU;
 241       output[2]=(((unsigned int)most_negative)>>16)&0xFFU;
 242       output[3]=(((unsigned int)most_negative)>>24)&0xFFU;
 243       for (i=0; i<natoms; i++)
 244         for (j=0; j<3; j++)
 245           for (k=0; k<nframes; k++)
 246             {
 247               int item=input[k*3*natoms+i*3+j];
 248               pval[cnt++]=(unsigned int)(item+most_negative);
 249
 250             }
 251       if (speed>=5)
 252         bwlzh_compress(pval,n,output+4,length);
 253       else
 254         bwlzh_compress_no_lz77(pval,n,output+4,length);
 255       (*length)+=4;
 256       free(pval);
 257       return output;
 258     }
 259   else if (coding==TNG_COMPRESS_ALGO_POS_XTC3)
 260     return Ptngc_pack_array_xtc3(input,length,natoms,speed);
 261   else if (coding==TNG_COMPRESS_ALGO_POS_XTC2)
 262     return Ptngc_pack_array_xtc2(coder_inst,input,length);
 263   else
 264     {
 265       unsigned char *output=NULL;
 266       unsigned char *output_ptr=NULL;
 267       int i;
 268       int output_length=0;
 269
 270       coder_inst->stat_numval=0;
 271       coder_inst->stat_overflow=0;
 272       /* Allocate enough memory for output */
 273       output=warnmalloc(8* *length*sizeof *output);
 274       output_ptr=output;
 275       if ((coding==TNG_COMPRESS_ALGO_TRIPLET) ||
 276           (coding==TNG_COMPRESS_ALGO_POS_TRIPLET_INTRA) ||
 277           (coding==TNG_COMPRESS_ALGO_POS_TRIPLET_ONETOONE))
 278         {
 279           /* Pack triplets. */
 280           int ntriplets=*length/3;
 281           /* Determine max base and maxbits */
 282           unsigned int max_base=1U<<coding_parameter;
 283           unsigned int maxbits=coding_parameter;
 284           unsigned int intmax=0;
 285           for (i=0; i<*length; i++)
 286             {
 287               int item=input[i];
 288               unsigned int s=0;
 289               if (item>0)
 290                 s=1+(item-1)*2;
 291               else if (item<0)
 292                 s=2+(-item-1)*2;
 293               if (s>intmax)
 294                 intmax=s;
 295             }
 296           /* Store intmax */
 297           coder_inst->pack_temporary_bits=32;
 298           coder_inst->pack_temporary=intmax;
 299           Ptngc_out8bits(coder_inst,&output_ptr);
 300           while (intmax>=max_base)
 301             {
 302               max_base*=2;
 303               maxbits++;
 304             }
 305           for (i=0; i<ntriplets; i++)
 306             {
 307               int j;
 308               unsigned int s[3];
 309               for (j=0; j<3; j++)
 310                 {
 311                   int item=input[i*3+j];
 312                   /* Find this symbol in table. */
 313                   s[j]=0;
 314                   if (item>0)
 315                     s[j]=1+(item-1)*2;
 316                   else if (item<0)
 317                     s[j]=2+(-item-1)*2;
 318                 }
 319               if (pack_triplet(coder_inst, s, &output_ptr,
 320                                coding_parameter, max_base,maxbits))
 321                 {
 322                   free(output);
 323                   return NULL;
 324                 }
 325             }
 326         }
 327       else
 328         for (i=0; i<*length; i++)
 329           if (pack_stopbits_item(coder_inst,input[i],&output_ptr,coding_parameter))
 330             {
 331               free(output);
 332               return NULL;
 333             }
 334       Ptngc_pack_flush(coder_inst,&output_ptr);
 335       output_length=(int)(output_ptr-output);
 336       *length=output_length;
 337       return output;
 338     }
 339 }
 340
 341 static int unpack_array_stop_bits(struct coder *coder_inst,
 342                                   unsigned char *packed,int *output,
 343                                   int length, int coding_parameter)
 344 {
 345   int i,j;
 346   unsigned int extract_mask=0x80;
 347   unsigned char *ptr=packed;
 348   (void) coder_inst;
 349   for (i=0; i<length; i++)
 350     {
 351       unsigned int pattern=0;
 352       int numbits=coding_parameter;
 353       unsigned int bit;
 354       int s;
 355       unsigned int insert_mask=1U<<(numbits-1);
 356       int inserted_bits=numbits;
 357       do {
 358         for (j=0; j<numbits; j++)
 359           {
 360             bit=*ptr & extract_mask;
 361             if (bit)
 362               pattern|=insert_mask;
 363             insert_mask>>=1;
 364             extract_mask>>=1;
 365             if (!extract_mask)
 366               {
 367                 extract_mask=0x80;
 368                 ptr++;
 369               }
 370           }
 371         /* Check stop bit */
 372         bit=*ptr & extract_mask;
 373         extract_mask>>=1;
 374         if (!extract_mask)
 375           {
 376             extract_mask=0x80;
 377             ptr++;
 378           }
 379         if (bit)
 380           {
 381             numbits>>=1;
 382             if (numbits<1)
 383               numbits=1;
 384             inserted_bits+=numbits;
 385             insert_mask=1U<<(inserted_bits-1);
 386           }
 387       } while (bit);
 388       s=(pattern+1)/2;
 389       if ((pattern%2)==0)
 390         s=-s;
 391       output[i]=s;
 392     }
 393   return 0;
 394 }
 395
 396 static int unpack_array_triplet(struct coder *coder_inst,
 397                                 unsigned char *packed, int *output,
 398                                 int length, int coding_parameter)
 399 {
 400   int i,j;
 401   unsigned int extract_mask=0x80;
 402   unsigned char *ptr=packed;
 403   /* Determine max base and maxbits */
 404   unsigned int max_base=1U<<coding_parameter;
 405   unsigned int maxbits=coding_parameter;
 406   unsigned int intmax;
 407   /* Get intmax */
 408   (void) coder_inst;
 409   intmax=((unsigned int)ptr[0])<<24|
 410     ((unsigned int)ptr[1])<<16|
 411     ((unsigned int)ptr[2])<<8|
 412     ((unsigned int)ptr[3]);
 413   ptr+=4;
 414   while (intmax>=max_base)
 415     {
 416       max_base*=2;
 417       maxbits++;
 418     }
 419   length/=3;
 420   for (i=0; i<length; i++)
 421     {
 422       /* Find base */
 423       unsigned int jbase=0;
 424       unsigned int numbits;
 425       unsigned int bit;
 426       for (j=0; j<2; j++)
 427         {
 428           bit=*ptr & extract_mask;
 429           jbase<<=1;
 430           if (bit)
 431             jbase|=1U;
 432           extract_mask>>=1;
 433           if (!extract_mask)
 434             {
 435               extract_mask=0x80;
 436               ptr++;
 437             }
 438         }
 439       if (jbase==3)
 440         numbits=maxbits;
 441       else
 442         numbits=coding_parameter+jbase;
 443       for (j=0; j<3; j++)
 444         {
 445           int s;
 446           unsigned int jbit;
 447           unsigned int pattern=0;
 448           for (jbit=0; jbit<numbits; jbit++)
 449             {
 450               bit=*ptr & extract_mask;
 451               pattern<<=1;
 452               if (bit)
 453                 pattern|=1U;
 454               extract_mask>>=1;
 455               if (!extract_mask)
 456                 {
 457                   extract_mask=0x80;
 458                   ptr++;
 459                 }
 460             }
 461           s=(pattern+1)/2;
 462           if ((pattern%2)==0)
 463             s=-s;
 464           output[i*3+j]=s;
 465         }
 466     }
 467   return 0;
 468 }
 469
 470 static int unpack_array_bwlzh(struct coder *coder_inst,
 471                               unsigned char *packed, int *output,
 472                               int length, int natoms)
 473 {
 474   int i,j,k,n=length;
 475   unsigned int *pval=warnmalloc(n*sizeof *pval);
 476   int nframes=n/natoms/3;
 477   int cnt=0;
 478   int most_negative=(int)(((unsigned int)packed[0]) |
 479                           (((unsigned int)packed[1])<<8) |
 480                           (((unsigned int)packed[2])<<16) |
 481                           (((unsigned int)packed[3])<<24));
 482   (void) coder_inst;
 483   bwlzh_decompress(packed+4,length,pval);
 484   for (i=0; i<natoms; i++)
 485     for (j=0; j<3; j++)
 486       for (k=0; k<nframes; k++)
 487         {
 488           unsigned int s=pval[cnt++];
 489           output[k*3*natoms+i*3+j]=(int)s-most_negative;
 490         }
 491   free(pval);
 492   return 0;
 493 }
 494
 495 int DECLSPECDLLEXPORT Ptngc_unpack_array(struct coder *coder_inst,
 496                        unsigned char *packed, int *output,
 497                        int length, int coding, int coding_parameter,
 498                        int natoms)
 499 {
 500   if ((coding==TNG_COMPRESS_ALGO_STOPBIT) ||
 501       (coding==TNG_COMPRESS_ALGO_VEL_STOPBIT_INTER))
 502     return unpack_array_stop_bits(coder_inst, packed, output, length, coding_parameter);
 503   else if ((coding==TNG_COMPRESS_ALGO_TRIPLET) ||
 504            (coding==TNG_COMPRESS_ALGO_POS_TRIPLET_INTRA) ||
 505            (coding==TNG_COMPRESS_ALGO_POS_TRIPLET_ONETOONE))
 506     return unpack_array_triplet(coder_inst, packed, output, length, coding_parameter);
 507   else if (coding==TNG_COMPRESS_ALGO_POS_XTC2)
 508     return Ptngc_unpack_array_xtc2(coder_inst, packed, output, length);
 509   else if ((coding==TNG_COMPRESS_ALGO_BWLZH1) || (coding==TNG_COMPRESS_ALGO_BWLZH2))
 510     return unpack_array_bwlzh(coder_inst, packed, output, length,natoms);
 511   else if (coding==TNG_COMPRESS_ALGO_POS_XTC3)
 512     return Ptngc_unpack_array_xtc3(packed, output, length,natoms);
 513   return 1;
 514 }
 515