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
Ptngc_coder_init(void)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
Ptngc_coder_deinit(struct coder * coder_inst)40 void DECLSPECDLLEXPORT Ptngc_coder_deinit(struct coder *coder_inst)
41 {
42 free(coder_inst);
43 }
44
Ptngc_out8bits(struct coder * coder_inst,unsigned char ** output)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
Ptngc_write_pattern(struct coder * coder_inst,unsigned int pattern,int nbits,unsigned char ** output)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 */
Ptngc_writebits(struct coder * coder_inst,unsigned int value,const int nbits,unsigned char ** output_ptr)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 */
Ptngc_write32bits(struct coder * coder_inst,unsigned int value,int nbits,unsigned char ** output_ptr)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 */
Ptngc_writemanybits(struct coder * coder_inst,unsigned char * value,int nbits,unsigned char ** output_ptr)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
write_stop_bit_code(struct coder * coder_inst,unsigned int s,unsigned int coding_parameter,unsigned char ** output)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
pack_stopbits_item(struct coder * coder_inst,const int item,unsigned char ** output,const int coding_parameter)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
pack_triplet(struct coder * coder_inst,unsigned int * s,unsigned char ** output,const int coding_parameter,const unsigned int max_base,const int maxbits)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
Ptngc_pack_flush(struct coder * coder_inst,unsigned char ** output)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
Ptngc_pack_array(struct coder * coder_inst,int * input,int * length,const int coding,const int coding_parameter,const int natoms,const int speed)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
unpack_array_stop_bits(struct coder * coder_inst,unsigned char * packed,int * output,const int length,const int coding_parameter)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
unpack_array_triplet(struct coder * coder_inst,unsigned char * packed,int * output,int length,const int coding_parameter)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
unpack_array_bwlzh(struct coder * coder_inst,unsigned char * packed,int * output,const int length,const int natoms)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
Ptngc_unpack_array(struct coder * coder_inst,unsigned char * packed,int * output,const int length,const int coding,const int coding_parameter,const int natoms)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
513