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

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