?  search:
summary | log | graphiclog1 | graphiclog2 | commit | commitdiff | tree | refs | edit | fork
blame | history | raw | HEAD
TNG version 1.7.3
[gromacs.git] / src / external / tng_io / src / compression / xtc2.c
blob7332f5ffe20eace47bb9834ccfed9ee81e10f525
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 /* This code is heavily influenced by
12 http://hpcv100.rc.rug.nl/xdrf.html
13 Based on coordinate compression (c) by Frans van Hoesel.
14 and GROMACS xtc files (http://www.gromacs.org)
15 (c) Copyright (c) Erik Lindahl, David van der Spoel
16 */
17
18 #include <stdio.h>
19 #include <stdlib.h>
20 #include <string.h>
21 #include <math.h>
22 #include "../../include/compression/coder.h"
23 #include "../../include/compression/widemuldiv.h"
24 #include "../../include/compression/warnmalloc.h"
25
26 /* Generated by gen_magic.py */
27 #define MAX_MAGIC 92
28
29 #ifdef USE_WINDOWS
30 #define TNG_INLINE __inline
31 #else
32 #define TNG_INLINE inline
33 #endif
34
35 static unsigned int magic[MAX_MAGIC]={
36 2U, 3U, 4U, 5U,
37 6U, 8U, 10U, 12U,
38 16U, 20U, 25U, 32U,
39 40U, 50U, 64U, 80U,
40 101U, 128U, 161U, 203U,
41 256U, 322U, 406U, 512U,
42 645U, 812U, 1024U, 1290U,
43 1625U, 2048U, 2580U, 3250U,
44 4096U, 5160U, 6501U, 8192U,
45 10321U, 13003U, 16384U, 20642U,
46 26007U, 32768U, 41285U, 52015U,
47 65536U, 82570U, 104031U, 131072U,
48 165140U, 208063U, 262144U, 330280U,
49 416127U, 524288U, 660561U, 832255U,
50 1048576U, 1321122U, 1664510U, 2097152U,
51 2642245U, 3329021U, 4194304U, 5284491U,
52 6658042U, 8388608U, 10568983U, 13316085U,
53 16777216U, 21137967U, 26632170U, 33554432U,
54 42275935U, 53264340U, 67108864U, 84551870U,
55 106528681U, 134217728U, 169103740U, 213057362U,
56 268435456U, 338207481U, 426114725U, 536870912U,
57 676414963U, 852229450U, 1073741824U, 1352829926U,
58 1704458900U, 2147483648U, 2705659852U, 3408917801U,
59 };
60
61 static unsigned int magic_bits[MAX_MAGIC][8]={
62 { 3, 6, 9, 12, 15, 18, 21, 24, },
63 { 5, 10, 15, 20, 24, 29, 34, 39, },
64 { 6, 12, 18, 24, 30, 36, 42, 48, },
65 { 7, 14, 21, 28, 35, 42, 49, 56, },
66 { 8, 16, 24, 32, 39, 47, 55, 63, },
67 { 9, 18, 27, 36, 45, 54, 63, 72, },
68 { 10, 20, 30, 40, 50, 60, 70, 80, },
69 { 11, 22, 33, 44, 54, 65, 76, 87, },
70 { 12, 24, 36, 48, 60, 72, 84, 97, },
71 { 13, 26, 39, 52, 65, 78, 91, 104, },
72 { 14, 28, 42, 56, 70, 84, 98, 112, },
73 { 15, 30, 45, 60, 75, 90, 105, 120, },
74 { 16, 32, 48, 64, 80, 96, 112, 128, },
75 { 17, 34, 51, 68, 85, 102, 119, 136, },
76 { 18, 36, 54, 72, 90, 108, 127, 144, },
77 { 19, 38, 57, 76, 95, 114, 133, 152, },
78 { 20, 40, 60, 80, 100, 120, 140, 160, },
79 { 21, 42, 63, 84, 105, 127, 147, 168, },
80 { 22, 44, 66, 88, 110, 132, 154, 176, },
81 { 23, 46, 69, 92, 115, 138, 161, 184, },
82 { 24, 48, 72, 97, 120, 144, 168, 192, },
83 { 25, 50, 75, 100, 125, 150, 175, 200, },
84 { 26, 52, 78, 104, 130, 156, 182, 208, },
85 { 27, 54, 81, 108, 135, 162, 190, 216, },
86 { 28, 56, 84, 112, 140, 168, 196, 224, },
87 { 29, 58, 87, 116, 145, 174, 203, 232, },
88 { 30, 60, 90, 120, 150, 180, 210, 240, },
89 { 31, 62, 93, 124, 155, 186, 217, 248, },
90 { 32, 64, 96, 128, 160, 192, 224, 256, },
91 { 33, 66, 99, 132, 165, 198, 231, 264, },
92 { 34, 68, 102, 136, 170, 204, 238, 272, },
93 { 35, 70, 105, 140, 175, 210, 245, 280, },
94 { 36, 72, 108, 144, 180, 216, 252, 288, },
95 { 37, 74, 111, 148, 185, 222, 259, 296, },
96 { 38, 76, 114, 152, 190, 228, 266, 304, },
97 { 39, 78, 117, 157, 195, 234, 273, 312, },
98 { 40, 80, 120, 160, 200, 240, 280, 320, },
99 { 41, 82, 123, 164, 205, 246, 287, 328, },
100 { 42, 84, 127, 168, 210, 252, 294, 336, },
101 { 43, 86, 129, 172, 215, 258, 301, 344, },
102 { 44, 88, 132, 176, 220, 264, 308, 352, },
103 { 45, 90, 135, 180, 225, 270, 315, 360, },
104 { 46, 92, 138, 184, 230, 276, 322, 368, },
105 { 47, 94, 141, 188, 235, 282, 329, 376, },
106 { 48, 97, 144, 192, 240, 288, 336, 384, },
107 { 49, 98, 147, 196, 245, 294, 343, 392, },
108 { 50, 100, 150, 200, 250, 300, 350, 400, },
109 { 52, 102, 153, 204, 255, 306, 357, 408, },
110 { 52, 104, 156, 208, 260, 312, 364, 416, },
111 { 53, 106, 159, 212, 265, 318, 371, 424, },
112 { 54, 108, 162, 216, 270, 324, 378, 432, },
113 { 55, 110, 165, 220, 275, 330, 385, 440, },
114 { 56, 112, 168, 224, 280, 336, 392, 448, },
115 { 57, 114, 172, 228, 285, 342, 399, 456, },
116 { 58, 116, 174, 232, 290, 348, 406, 464, },
117 { 59, 118, 177, 236, 295, 354, 413, 472, },
118 { 60, 120, 180, 240, 300, 360, 420, 480, },
119 { 61, 122, 183, 244, 305, 366, 427, 488, },
120 { 62, 124, 186, 248, 310, 372, 434, 496, },
121 { 63, 127, 190, 252, 315, 378, 442, 505, },
122 { 64, 128, 192, 256, 320, 384, 448, 512, },
123 { 65, 130, 195, 260, 325, 390, 455, 520, },
124 { 66, 132, 198, 264, 330, 396, 462, 528, },
125 { 67, 134, 201, 268, 335, 402, 469, 536, },
126 { 68, 136, 204, 272, 340, 408, 476, 544, },
127 { 69, 138, 207, 276, 345, 414, 483, 552, },
128 { 70, 140, 210, 280, 350, 420, 490, 560, },
129 { 71, 142, 213, 284, 355, 426, 497, 568, },
130 { 72, 144, 216, 288, 360, 432, 505, 576, },
131 { 73, 146, 219, 292, 365, 438, 511, 584, },
132 { 74, 148, 222, 296, 370, 444, 518, 592, },
133 { 75, 150, 225, 300, 375, 451, 525, 600, },
134 { 76, 152, 228, 304, 380, 456, 532, 608, },
135 { 77, 154, 231, 308, 385, 462, 539, 616, },
136 { 78, 157, 234, 312, 390, 469, 546, 625, },
137 { 79, 158, 237, 316, 395, 474, 553, 632, },
138 { 80, 160, 240, 320, 400, 480, 560, 640, },
139 { 81, 162, 243, 324, 406, 486, 568, 648, },
140 { 82, 164, 246, 328, 410, 492, 574, 656, },
141 { 83, 166, 249, 332, 415, 498, 581, 664, },
142 { 84, 168, 252, 336, 420, 505, 588, 672, },
143 { 85, 170, 255, 340, 425, 510, 595, 680, },
144 { 86, 172, 258, 344, 430, 516, 602, 688, },
145 { 87, 174, 261, 348, 435, 522, 609, 696, },
146 { 88, 176, 264, 352, 440, 528, 616, 704, },
147 { 89, 178, 267, 356, 445, 534, 623, 712, },
148 { 90, 180, 270, 360, 451, 540, 631, 720, },
149 { 91, 182, 273, 364, 455, 546, 637, 728, },
150 { 92, 184, 276, 368, 460, 552, 644, 736, },
151 { 94, 187, 279, 373, 466, 558, 651, 745, },
152 { 94, 188, 282, 376, 470, 564, 658, 752, },
153 { 95, 190, 285, 380, 475, 570, 665, 760, },
154 };
155
156
157 static const double iflipgaincheck=0.89089871814033927; /* 1./(2**(1./6)) */
158
159
160 /* Difference in indices used for determining whether to store as large or small */
161 #define QUITE_LARGE 3
162 #define IS_LARGE 6
163
164 #if 0
165 #define SHOWIT
166 #endif
167
168 #ifdef USE_WINDOWS
169 #define TNG_INLINE __inline
170 #else
171 #define TNG_INLINE inline
172 #endif
173
174 int Ptngc_magic(const unsigned int i)
175 {
176 return magic[i];
177 }
178
179 int Ptngc_find_magic_index(const unsigned int maxval)
180 {
181 int i;
182
183 if(maxval > magic[MAX_MAGIC/4])
184 {
185 if(maxval > magic[MAX_MAGIC/2])
186 {
187 i = MAX_MAGIC/2 + 1;
188 }
189 else
190 {
191 i = MAX_MAGIC/4 + 1;
192 }
193 }
194 else
195 {
196 i = 0;
197 }
198
199 while (magic[i]<=maxval)
200 i++;
201 return i;
202 }
203
204 static TNG_INLINE unsigned int positive_int(const int item)
205 {
206 int s=0;
207 if (item>0)
208 s=1+(item-1)*2;
209 else if (item<0)
210 s=2+(-item-1)*2;
211 return s;
212 }
213
214 static TNG_INLINE int unpositive_int(const int val)
215 {
216 int s=(val+1)/2;
217 if ((val%2)==0)
218 s=-s;
219 return s;
220 }
221
222
223 /* Sequence instructions */
224 #define INSTR_DEFAULT 0
225 #define INSTR_BASE_RUNLENGTH 1
226 #define INSTR_ONLY_LARGE 2
227 #define INSTR_ONLY_SMALL 3
228 #define INSTR_LARGE_BASE_CHANGE 4
229 #define INSTR_FLIP 5
230 #define INSTR_LARGE_RLE 6
231
232 #define MAXINSTR 7
233
234 #ifdef SHOWIT
235 static char *instrnames[MAXINSTR]={
236 "large+small",
237 "base+run",
238 "large",
239 "small",
240 "large base change",
241 "flip",
242 "large rle",
243 };
244 #endif
245
246 /* Bit patterns in the compressed code stream: */
247
248 static const int seq_instr[MAXINSTR][2]=
249 {
250 { 1,1 }, /* 1 : one large atom + runlength encoded small integers. Use same settings as before. */
251 { 0,2 }, /* 00 : set base and runlength in next four bits (x). base (increase/keep/decrease)=x%3-1. runlength=1+x/3.
252 The special value 1111 in the four bits means runlength=6 and base change=0 */
253 { 4,4 }, /* 0100 : next only a large atom comes. */
254 { 5,4 }, /* 0101 : next only runlength encoded small integers. Use same settings as before. */
255 { 6,4 }, /* 0110 : Large change in base. Change is encoded in the
256 following 2 bits. change direction (sign) is the first
257 bit. The next bit +1 is the actual change. This
258 allows the change of up to +/- 2 indices. */
259 { 14,5 }, /* 01110 : flip whether integers should be modified to compress water better */
260 { 15,5 }, /* 01111 : Large rle. The next 4 bits encode how many
261 large atoms are in the following sequence: 3-18. (2 is
262 more efficiently coded with two large instructions. */
263 };
264
265 static void write_instruction(struct coder *coder, const int instr, unsigned char **output_ptr)
266 {
267 Ptngc_writebits(coder,seq_instr[instr][0],seq_instr[instr][1],output_ptr);
268 #ifdef SHOWIT
269 fprintf(stderr,"INSTR: %s (%d bits)\n",instrnames[instr],seq_instr[instr][1]);
270 #endif
271 }
272
273 static unsigned int readbits(unsigned char **ptr, int *bitptr, int nbits)
274 {
275 unsigned int val=0U;
276 unsigned int extract_mask=0x80U>>*bitptr;
277 unsigned char thisval=**ptr;
278 #ifdef SHOWIT
279 fprintf(stderr,"Read nbits=%d val=",nbits);
280 #endif
281 while (nbits--)
282 {
283 val<<=1;
284 val|=((extract_mask & thisval)!=0);
285 *bitptr=(*bitptr)+1;
286 extract_mask>>=1;
287 if (!extract_mask)
288 {
289 extract_mask=0x80U;
290 *ptr=(*ptr)+1;
291 *bitptr=0;
292 if (nbits)
293 thisval=**ptr;
294 }
295 }
296 #ifdef SHOWIT
297 fprintf(stderr,"%x\n",val);
298 #endif
299 return val;
300 }
301
302 static void readmanybits(unsigned char **ptr, int *bitptr, int nbits, unsigned char *buffer)
303 {
304 while (nbits>=8)
305 {
306 *buffer++=readbits(ptr,bitptr,8);
307 nbits-=8;
308 #ifdef SHOWIT
309 fprintf(stderr,"Read value %02x\n",buffer[-1]);
310 #endif
311 }
312 if (nbits)
313 {
314 *buffer++=readbits(ptr,bitptr,nbits);
315 #ifdef SHOWIT
316 fprintf(stderr,"Read value %02x\n",buffer[-1]);
317 #endif
318 }
319 }
320
321 static int read_instruction(unsigned char **ptr, int *bitptr)
322 {
323 int instr=-1;
324 unsigned int bits=readbits(ptr,bitptr,1);
325 if (bits)
326 instr=INSTR_DEFAULT;
327 else
328 {
329 bits=readbits(ptr,bitptr,1);
330 if (!bits)
331 instr=INSTR_BASE_RUNLENGTH;
332 else
333 {
334 bits=readbits(ptr,bitptr,2);
335 if (bits==0)
336 instr=INSTR_ONLY_LARGE;
337 else if (bits==1)
338 instr=INSTR_ONLY_SMALL;
339 else if (bits==2)
340 instr=INSTR_LARGE_BASE_CHANGE;
341 else if (bits==3)
342 {
343 bits=readbits(ptr,bitptr,1);
344 if (bits==0)
345 instr=INSTR_FLIP;
346 else
347 instr=INSTR_LARGE_RLE;
348 }
349 }
350 }
351 return instr;
352 }
353
354 /* Modifies three integer values for better compression of water */
355 static void swap_ints(int *in, int *out)
356 {
357 out[0]=in[0]+in[1];
358 out[1]=-in[1];
359 out[2]=in[1]+in[2];
360 }
361
362 static void swap_is_better(int *input, int *minint, int *sum_normal, int *sum_swapped)
363 {
364 int normal_max=0;
365 int swapped_max=0;
366 int i,j;
367 int normal[3];
368 int swapped[3];
369 for (i=0; i<3; i++)
370 {
371 normal[0]=input[i]-minint[i];
372 normal[1]=input[3+i]-input[i]; /* minint[i]-minint[i] cancels out */
373 normal[2]=input[6+i]-input[3+i]; /* minint[i]-minint[i] cancels out */
374 swap_ints(normal,swapped);
375 for (j=1; j<3; j++)
376 {
377 if (positive_int(normal[j])>(unsigned int)normal_max)
378 normal_max=positive_int(normal[j]);
379 if (positive_int(swapped[j])>(unsigned int)swapped_max)
380 swapped_max=positive_int(swapped[j]);
381 }
382 }
383 if (normal_max==0)
384 normal_max=1;
385 if (swapped_max==0)
386 swapped_max=1;
387 *sum_normal=normal_max;
388 *sum_swapped=swapped_max;
389 }
390
391 static void swapdecide(struct coder *coder, int *input, int *swapatoms, int *large_index, int *minint, unsigned char **output_ptr)
392 {
393 int didswap=0;
394 int normal,swapped;
395 (void)large_index;
396 swap_is_better(input,minint,&normal,&swapped);
397 /* We have to determine if it is worth to change the behaviour.
398 If diff is positive it means that it is worth something to
399 swap. But it costs 4 bits to do the change. If we assume that
400 we gain 0.17 bit by the swap per value, and the runlength>2
401 for four molecules in a row, we gain something. So check if we
402 gain at least 0.17 bits to even attempt the swap.
403 */
404 #ifdef SHOWIT
405 fprintf(stderr,"Trying Flip: %g %g\n",(double)swapped/normal, (double)normal/swapped);
406 #endif
407 if (((swapped<normal) && (fabs((double)swapped/normal)<iflipgaincheck)) ||
408 ((normal<swapped) && (fabs((double)normal/swapped)<iflipgaincheck)))
409 {
410 if (swapped<normal)
411 {
412 if (!*swapatoms)
413 {
414 *swapatoms=1;
415 didswap=1;
416 }
417 }
418 else
419 {
420 if (*swapatoms)
421 {
422 *swapatoms=0;
423 didswap=1;
424 }
425 }
426 }
427 if (didswap)
428 {
429 #ifdef SHOWIT
430 fprintf(stderr,"Flip. Swapatoms is now %d\n",*swapatoms);
431 #endif
432 write_instruction(coder,INSTR_FLIP,output_ptr);
433 }
434 }
435
436 /* Compute number of bits required to store values using three different bases in the index array */
437 static int compute_magic_bits(int *index)
438 {
439 unsigned int largeint[4];
440 unsigned int largeint_tmp[4];
441 int i,j,onebit;
442 for (i=0; i<4; i++)
443 largeint[i]=0U;
444 for (i=0; i<3; i++)
445 {
446 if (i!=0)
447 {
448 /* We must do the multiplication of the largeint with the integer base */
449 Ptngc_largeint_mul(magic[index[i]],largeint,largeint_tmp,4);
450 for (j=0; j<4; j++)
451 largeint[j]=largeint_tmp[j];
452 }
453 Ptngc_largeint_add(magic[index[i]]-1,largeint,4);
454 }
455 /* Find last bit. */
456 #if 0
457 printf("Largeint is %u %u %u\n",largeint[0],largeint[1],largeint[2]);
458 #endif
459 onebit=0;
460 for (i=0; i<3; i++)
461 for (j=0; j<32; j++)
462 if (largeint[i]&(1U<<j))
463 onebit=i*32+j+1;
464 return onebit;
465 }
466
467 /* Convert a sequence of (hopefully) small positive integers
468 using the base pointed to by index (x base, y base and z base can be different).
469 The largest number of integers supported is 18 (29 to handle/detect overflow) */
470 static void trajcoder_base_compress(int *input, const int n, int *index, unsigned char *result)
471 {
472 unsigned int largeint[19];
473 unsigned int largeint_tmp[19];
474 int i, j;
475
476 memset(largeint, 0U, sizeof(unsigned int) * 19);
477
478 if(n > 0)
479 {
480 Ptngc_largeint_add(input[0],largeint,19);
481 }
482
483 for (i=1; i<n; i++)
484 {
485 /* We must do the multiplication of the largeint with the integer base */
486 Ptngc_largeint_mul(magic[index[i%3]],largeint,largeint_tmp,19);
487
488 memcpy(largeint,largeint_tmp,19*sizeof *largeint);
489
490 Ptngc_largeint_add(input[i],largeint,19);
491 }
492 if (largeint[18])
493 {
494 fprintf(stderr,"TRAJNG: BUG! Overflow in compression detected.\n");
495 exit(EXIT_FAILURE);
496 }
497 #if 0
498 #ifdef SHOWIT
499 for (i=0; i<19; i++)
500 fprintf(stderr,"Largeint[%d]=0x%x\n",i,largeint[i]);
501 #endif
502 #endif
503 /* Convert the largeint to a sequence of bytes. */
504 for (i=0; i<18; i++)
505 {
506 int shift=0;
507 for (j=0; j<4; j++)
508 {
509 result[i*4+j]=(largeint[i]>>shift) & 0xFFU;
510 shift+=8;
511 }
512 }
513 }
514
515 /* The opposite of base_compress. */
516 static void trajcoder_base_decompress(unsigned char *input, const int n, int *index, int *output)
517 {
518 unsigned int largeint[19];
519 unsigned int largeint_tmp[19];
520 int i,j;
521 /* Convert the sequence of bytes to a largeint. */
522 for (i=0; i<18; i++)
523 {
524 int shift=0;
525 largeint[i]=0U;
526 for (j=0; j<4; j++)
527 {
528 largeint[i]|=((unsigned int)input[i*4+j])<<shift;
529 shift+=8;
530 }
531 }
532 largeint[18]=0U;
533 #if 0
534 #ifdef SHOWIT
535 for (i=0; i<19; i++)
536 fprintf(stderr,"Largeint[%d]=0x%x\n",i,largeint[i]);
537 #endif
538 #endif
539 for (i=n-1; i>=0; i--)
540 {
541 unsigned int remainder=Ptngc_largeint_div(magic[index[i%3]],largeint,largeint_tmp,19);
542 #if 0
543 #ifdef SHOWIT
544 fprintf(stderr,"Remainder: %u\n",remainder);
545 #endif
546 #endif
547 #if 0
548 for (j=0; j<19; j++)
549 largeint[j]=largeint_tmp[j];
550 #endif
551 memcpy(largeint,largeint_tmp,19*sizeof *largeint);
552 output[i]=remainder;
553 }
554 }
555
556 /* It is "large" if we have to increase the small index quite a
557 bit. Not so much to be rejected by the not very large check
558 later. */
559 static int is_quite_large(int *input, const int small_index, const int max_large_index)
560 {
561 int is=0;
562 int i;
563 if (small_index+QUITE_LARGE>=max_large_index)
564 is=1;
565 else
566 {
567 for (i=0; i<3; i++)
568 if (positive_int(input[i])>magic[small_index+QUITE_LARGE])
569 {
570 is=1;
571 break;
572 }
573 }
574 return is;
575 }
576
577 #ifdef SHOWIT
578 int nbits_sum;
579 int nvalues_sum;
580 #endif
581
582 static void write_three_large(struct coder *coder, int *encode_ints, int *large_index, const int nbits, unsigned char *compress_buffer, unsigned char **output_ptr)
583 {
584 trajcoder_base_compress(encode_ints,3,large_index,compress_buffer);
585 Ptngc_writemanybits(coder,compress_buffer,nbits,output_ptr);
586 #ifdef SHOWIT
587 fprintf(stderr,"nbits=%d (%g)\n",nbits,nbits/3.);
588 nbits_sum+=nbits;
589 nvalues_sum+=3;
590 fprintf(stderr,"Large: %d %d %d\n",encode_ints[0],encode_ints[1],encode_ints[2]);
591 #endif
592 }
593
594 static void insert_batch(int *input_ptr, int ntriplets_left, const int *prevcoord,int *minint, int *encode_ints, const int startenc, int *nenc)
595 {
596 int nencode=startenc*3;
597 int tmp_prevcoord[3];
598
599 tmp_prevcoord[0]=prevcoord[0];
600 tmp_prevcoord[1]=prevcoord[1];
601 tmp_prevcoord[2]=prevcoord[2];
602
603 if (startenc)
604 {
605 int i;
606 for (i=0; i<startenc; i++)
607 {
608 tmp_prevcoord[0]+=encode_ints[i*3];
609 tmp_prevcoord[1]+=encode_ints[i*3+1];
610 tmp_prevcoord[2]+=encode_ints[i*3+2];
611 #ifdef SHOWIT
612 fprintf(stderr,"%6d: %6d %6d %6d\t\t%6d %6d %6d\t\t%6d %6d %6d\n",i*3,
613 tmp_prevcoord[0],tmp_prevcoord[1],tmp_prevcoord[2],
614 encode_ints[i*3],
615 encode_ints[i*3+1],
616 encode_ints[i*3+2],
617 positive_int(encode_ints[i*3]),
618 positive_int(encode_ints[i*3+1]),
619 positive_int(encode_ints[i*3+2]));
620 #endif
621 }
622 }
623
624 #ifdef SHOWIT
625 fprintf(stderr,"New batch\n");
626 #endif
627 while ((nencode<21) && (nencode<ntriplets_left*3))
628 {
629 encode_ints[nencode]=input_ptr[nencode]-minint[0]-tmp_prevcoord[0];
630 encode_ints[nencode+1]=input_ptr[nencode+1]-minint[1]-tmp_prevcoord[1];
631 encode_ints[nencode+2]=input_ptr[nencode+2]-minint[2]-tmp_prevcoord[2];
632 #ifdef SHOWIT
633 fprintf(stderr,"%6d: %6d %6d %6d\t\t%6d %6d %6d\t\t%6d %6d %6d\n",nencode,
634 input_ptr[nencode]-minint[0],
635 input_ptr[nencode+1]-minint[1],
636 input_ptr[nencode+2]-minint[2],
637 encode_ints[nencode],
638 encode_ints[nencode+1],
639 encode_ints[nencode+2],
640 positive_int(encode_ints[nencode]),
641 positive_int(encode_ints[nencode+1]),
642 positive_int(encode_ints[nencode+2]));
643 #endif
644 tmp_prevcoord[0]=input_ptr[nencode]-minint[0];
645 tmp_prevcoord[1]=input_ptr[nencode+1]-minint[1];
646 tmp_prevcoord[2]=input_ptr[nencode+2]-minint[2];
647 nencode+=3;
648 }
649 *nenc=nencode;
650 }
651
652 static void flush_large(struct coder *coder, int *has_large, int *has_large_ints, const int n,
653 int *large_index, const int large_nbits, unsigned char *compress_buffer,
654 unsigned char **output_ptr)
655 {
656 int i;
657 if (n<3)
658 {
659 for (i=0; i<n; i++)
660 {
661 write_instruction(coder,INSTR_ONLY_LARGE,output_ptr);
662 write_three_large(coder,has_large_ints+i*3,large_index,large_nbits,compress_buffer,output_ptr);
663 }
664 }
665 else
666 {
667 #if 0
668 printf("CODING RLE: %d\n",n);
669 #endif
670 write_instruction(coder,INSTR_LARGE_RLE,output_ptr);
671 Ptngc_writebits(coder,n-3,4,output_ptr); /* Number of large atoms in this sequence: 3 to 18 */
672 for (i=0; i<n; i++)
673 write_three_large(coder,has_large_ints+i*3,large_index,large_nbits,compress_buffer,output_ptr);
674 }
675 if (((*has_large)-n)!=0)
676 {
677 int j;
678 for (i=0; i<(*has_large)-n; i++)
679 for (j=0; j<3; j++)
680 has_large_ints[i*3+j]=has_large_ints[(i+n)*3+j];
681 }
682 *has_large=(*has_large)-n; /* Number of remaining large atoms in buffer */
683 }
684
685 static void buffer_large(struct coder *coder, int *has_large, int *has_large_ints, int *new_large_ints,
686 int *large_index, const int large_nbits, unsigned char *compress_buffer,
687 unsigned char **output_ptr)
688 {
689 /* If it is full we must write them all. */
690 if (*has_large==18)
691 flush_large(coder,has_large,has_large_ints, *has_large,
692 large_index,large_nbits,compress_buffer,output_ptr); /* Flush them all. */
693 has_large_ints[(*has_large)*3]=new_large_ints[0];
694 has_large_ints[(*has_large)*3+1]=new_large_ints[1];
695 has_large_ints[(*has_large)*3+2]=new_large_ints[2];
696 *has_large=(*has_large)+1;
697 }
698
699
700 unsigned char *Ptngc_pack_array_xtc2(struct coder *coder,int *input, int *length)
701 {
702 unsigned char *output=NULL;
703 unsigned char *output_ptr=NULL;
704 int i,ienc,j;
705 int output_length=0;
706 /* Pack triplets. */
707 int ntriplets=*length/3;
708 int intmax;
709 int max_small;
710 int large_index[3];
711 int max_large_index;
712 int large_nbits;
713 int small_index;
714 int small_idx[3];
715 int minint[3],maxint[3];
716 int has_large=0;
717 int has_large_ints[54]; /* Large cache. Up to 18 large atoms. */
718 int prevcoord[3];
719 int runlength=0; /* Initial runlength. "Stupidly" set to zero for
720 simplicity and explicity */
721 int swapatoms=0; /* Initial guess is that we should not swap the
722 first two atoms in each large+small
723 transition */
724 int didswap; /* Whether swapping was actually done. */
725 int *input_ptr=input;
726 int encode_ints[21]; /* Up to 3 large + 18 small ints can be encoded at once */
727 int nencode;
728 unsigned char compress_buffer[18*4]; /* Holds compressed result for 3 large ints or up to 18 small ints. */
729 int ntriplets_left=ntriplets;
730 int refused=0;
731 #ifdef SHOWIT
732 nbits_sum=0;
733 nvalues_sum=0;
734 #endif
735 /* Allocate enough memory for output */
736 output=warnmalloc(8* *length*sizeof *output);
737 output_ptr=output;
738
739 maxint[0]=minint[0]=input[0];
740 maxint[1]=minint[1]=input[1];
741 maxint[2]=minint[2]=input[2];
742
743 for (i=1; i<ntriplets; i++)
744 {
745 for (j=0; j<3; j++)
746 {
747 if (input[i*3+j]>maxint[j])
748 maxint[j]=input[i*3+j];
749 if (input[i*3+j]<minint[j])
750 minint[j]=input[i*3+j];
751 }
752 }
753
754 large_index[0]=Ptngc_find_magic_index(maxint[0]-minint[0]+1);
755 large_index[1]=Ptngc_find_magic_index(maxint[1]-minint[1]+1);
756 large_index[2]=Ptngc_find_magic_index(maxint[2]-minint[2]+1);
757 large_nbits=compute_magic_bits(large_index);
758 max_large_index=large_index[0];
759 if (large_index[1]>max_large_index)
760 max_large_index=large_index[1];
761 if (large_index[2]>max_large_index)
762 max_large_index=large_index[2];
763
764 #ifdef SHOWIT
765 for (j=0; j<3; j++)
766 fprintf(stderr,"minint[%d]=%d. maxint[%d]=%d large_index[%d]=%d value=%d\n",j,minint[j],j,maxint[j],
767 j,large_index[j],magic[large_index[j]]);
768 fprintf(stderr,"large_nbits=%d\n",large_nbits);
769 #endif
770
771
772 /* Guess initial small index */
773 small_index=max_large_index/2;
774
775 /* Find the largest value that is not large. Not large is half index of
776 large. */
777 max_small=magic[small_index];
778 intmax=0;
779 for (i=0; i<*length; i++)
780 {
781 int item=input[i];
782 int s=positive_int(item);
783 if (s>intmax)
784 if (s<max_small)
785 intmax=s;
786 }
787 /* This value is not critical, since if I guess wrong, the code will
788 just insert instructions to increase this value. */
789 small_index=Ptngc_find_magic_index(intmax);
790 #ifdef SHOWIT
791 fprintf(stderr,"initial_small intmax=%d. small_index=%d value=%d\n",intmax,small_index,magic[small_index]);
792 #endif
793
794 /* Store min integers */
795 coder->pack_temporary_bits=32;
796 coder->pack_temporary=positive_int(minint[0]);
797 Ptngc_out8bits(coder,&output_ptr);
798 coder->pack_temporary_bits=32;
799 coder->pack_temporary=positive_int(minint[1]);
800 Ptngc_out8bits(coder,&output_ptr);
801 coder->pack_temporary_bits=32;
802 coder->pack_temporary=positive_int(minint[2]);
803 Ptngc_out8bits(coder,&output_ptr);
804 /* Store max indices */
805 coder->pack_temporary_bits=8;
806 coder->pack_temporary=large_index[0];
807 Ptngc_out8bits(coder,&output_ptr);
808 coder->pack_temporary_bits=8;
809 coder->pack_temporary=large_index[1];
810 Ptngc_out8bits(coder,&output_ptr);
811 coder->pack_temporary_bits=8;
812 coder->pack_temporary=large_index[2];
813 Ptngc_out8bits(coder,&output_ptr);
814 /* Store initial small index */
815 coder->pack_temporary_bits=8;
816 coder->pack_temporary=small_index;
817 Ptngc_out8bits(coder,&output_ptr);
818
819 #if 0
820 #ifdef SHOWIT
821 for (i=0; i<ntriplets_left; i++)
822 fprintf(stderr,"VALUE:%d %6d %6d %6d\n",
823 i,
824 input_ptr[i*3],
825 input_ptr[i*3+1],
826 input_ptr[i*3+2]);
827 #endif
828 #endif
829
830 /* Initial prevcoord is the minimum integers. */
831 memcpy(prevcoord, minint, 3*sizeof *prevcoord);
832
833 while (ntriplets_left)
834 {
835 if (ntriplets_left<0)
836 {
837 fprintf(stderr,"TRAJNG: BUG! ntriplets_left<0!\n");
838 exit(EXIT_FAILURE);
839 }
840 /* If only less than three atoms left we just write them all as large integers. Here no swapping is done! */
841 if (ntriplets_left<3)
842 {
843 for (ienc=0; ienc<ntriplets_left; ienc++)
844 {
845 int jenc;
846 for (jenc=0; jenc<3; jenc++)
847 encode_ints[jenc]=input_ptr[ienc*3+jenc]-minint[jenc];
848 buffer_large(coder,&has_large,has_large_ints,encode_ints,large_index,large_nbits,compress_buffer,&output_ptr);
849 input_ptr+=3;
850 ntriplets_left--;
851 }
852 flush_large(coder,&has_large,has_large_ints,has_large,large_index,large_nbits,compress_buffer,&output_ptr);
853 }
854 else
855 {
856 int min_runlength=0;
857 int largest_required_base;
858 int largest_runlength_base;
859 int largest_required_index;
860 int largest_runlength_index;
861 int new_runlength;
862 int new_small_index;
863 int iter_runlength;
864 int iter_small_index;
865 int rle_index_dep;
866 didswap=0;
867 /* Insert the next batch of integers to be encoded into the buffer */
868 #ifdef SHOWIT
869 fprintf(stderr,"Initial batch\n");
870 #endif
871 insert_batch(input_ptr,ntriplets_left,prevcoord,minint,encode_ints,0,&nencode);
872
873 /* First we must decide if the next value is large (does not reasonably fit in current small encoding)
874 Also, if we have not written any values yet, we must begin by writing a large atom. */
875 if ((input_ptr==input) || (is_quite_large(encode_ints,small_index,max_large_index)) || (refused))
876 {
877 /* If any of the next two atoms are large we should probably write them as large and not swap them */
878 int no_swap=0;
879 if ((is_quite_large(encode_ints+3,small_index,max_large_index)) || (is_quite_large(encode_ints+6,small_index,max_large_index)))
880 no_swap=1;
881 if (!no_swap)
882 {
883 /* Next we must decide if we should swap the first
884 two values. */
885 #if 1
886 swapdecide(coder,input_ptr,&swapatoms,large_index,minint,&output_ptr);
887 #else
888 swapatoms=0;
889 #endif
890 /* If we should do the integer swapping manipulation we should do it now */
891 if (swapatoms)
892 {
893 didswap=1;
894 for (i=0; i<3; i++)
895 {
896 int in[3], out[3];
897 in[0]=input_ptr[i]-minint[i];
898 in[1]=input_ptr[3+i]-input_ptr[i]; /* minint[i]-minint[i] cancels out */
899 in[2]=input_ptr[6+i]-input_ptr[3+i]; /* minint[i]-minint[i] cancels out */
900 swap_ints(in,out);
901 encode_ints[i]=out[0];
902 encode_ints[3+i]=out[1];
903 encode_ints[6+i]=out[2];
904 }
905 /* We have swapped atoms, so the minimum run-length is 2 */
906 #ifdef SHOWIT
907 fprintf(stderr,"Swap atoms results in:\n");
908 for (i=0; i<3; i++)
909 fprintf(stderr,"%d: %6d %6d %6d\t\t%6d %6d %6d\n",i*3,
910 encode_ints[i*3],
911 encode_ints[i*3+1],
912 encode_ints[i*3+2],
913 positive_int(encode_ints[i*3]),
914 positive_int(encode_ints[i*3+1]),
915 positive_int(encode_ints[i*3+2]));
916
917 #endif
918 min_runlength=2;
919 }
920 }
921 if ((swapatoms) && (didswap))
922 {
923 for (ienc=0; ienc<3; ienc++)
924 prevcoord[ienc]=encode_ints[ienc];
925 }
926 else
927 {
928 for (ienc=0; ienc<3; ienc++)
929 prevcoord[ienc]=input_ptr[ienc]-minint[ienc];
930 }
931 /* Cache large value for later possible combination with
932 a sequence of small integers. */
933 buffer_large(coder,&has_large,has_large_ints,prevcoord,large_index,large_nbits,compress_buffer,&output_ptr);
934 #ifdef SHOWIT
935 fprintf(stderr,"Prevcoord after packing of large: %d %d %d\n",
936 prevcoord[0],prevcoord[1],prevcoord[2]);
937 #endif
938
939 /* We have written a large integer so we have one less atoms to worry about */
940 input_ptr+=3;
941 ntriplets_left--;
942
943 refused=0;
944
945 /* Insert the next batch of integers to be encoded into the buffer */
946 #ifdef SHOWIT
947 fprintf(stderr,"Update batch due to large int.\n");
948 #endif
949 if ((swapatoms) && (didswap))
950 {
951 /* Keep swapped values. */
952 for (i=0; i<2; i++)
953 for (ienc=0; ienc<3; ienc++)
954 encode_ints[i*3+ienc]=encode_ints[(i+1)*3+ienc];
955 }
956 insert_batch(input_ptr,ntriplets_left,prevcoord,minint,encode_ints,min_runlength,&nencode);
957 }
958 /* Here we should only have differences for the atom coordinates. */
959 /* Convert the ints to positive ints */
960 for (ienc=0; ienc<nencode; ienc++)
961 {
962 int pint=positive_int(encode_ints[ienc]);
963 encode_ints[ienc]=pint;
964 }
965 /* Now we must decide what base and runlength to do. If we have swapped atoms it will be at least 2.
966 If even the next atom is large, we will not do anything. */
967 largest_required_base=0;
968 /* Determine required base */
969 for (ienc=0; ienc<min_runlength*3; ienc++)
970 if (encode_ints[ienc]>largest_required_base)
971 largest_required_base=encode_ints[ienc];
972 /* Also compute what the largest base is for the current runlength setting! */
973 largest_runlength_base=0;
974 for (ienc=0; (ienc<runlength*3) && (ienc<nencode); ienc++)
975 if (encode_ints[ienc]>largest_runlength_base)
976 largest_runlength_base=encode_ints[ienc];
977
978 largest_required_index=Ptngc_find_magic_index(largest_required_base);
979 largest_runlength_index=Ptngc_find_magic_index(largest_runlength_base);
980
981 if (largest_required_index<largest_runlength_index)
982 {
983 new_runlength=min_runlength;
984 new_small_index=largest_required_index;
985 }
986 else
987 {
988 new_runlength=runlength;
989 new_small_index=largest_runlength_index;
990 }
991
992 /* Only allow increase of runlength wrt min_runlength */
993 if (new_runlength<min_runlength)
994 new_runlength=min_runlength;
995
996 /* If the current runlength is longer than the number of
997 triplets left stop it from being so. */
998 if (new_runlength>ntriplets_left)
999 new_runlength=ntriplets_left;
1000
1001 /* We must at least try to get some small integers going. */
1002 if (new_runlength==0)
1003 {
1004 new_runlength=1;
1005 new_small_index=small_index;
1006 }
1007
1008 iter_runlength=new_runlength;
1009 iter_small_index=new_small_index;
1010
1011 /* Iterate to find optimal encoding and runlength */
1012 #ifdef SHOWIT
1013 fprintf(stderr,"Entering iterative loop.\n");
1014 fflush(stderr);
1015 #endif
1016
1017 do {
1018 new_runlength=iter_runlength;
1019 new_small_index=iter_small_index;
1020
1021 #ifdef SHOWIT
1022 fprintf(stderr,"Test new_small_index=%d Base=%d\n",new_small_index,magic[new_small_index]);
1023 #endif
1024 /* What is the largest runlength
1025 we can do with the currently
1026 selected encoding? Also the max supported runlength is 6 triplets! */
1027 for (ienc=0; ienc<nencode && ienc<18; ienc++)
1028 {
1029 int test_index=Ptngc_find_magic_index(encode_ints[ienc]);
1030 if (test_index>new_small_index)
1031 break;
1032 }
1033 if (ienc/3>new_runlength)
1034 {
1035 iter_runlength=ienc/3;
1036 #ifdef SHOWIT
1037 fprintf(stderr,"I found a new possible runlength: %d\n",iter_runlength);
1038 #endif
1039 }
1040
1041 /* How large encoding do we have to use? */
1042 largest_runlength_base=0;
1043 for (ienc=0; ienc<iter_runlength*3; ienc++)
1044 if (encode_ints[ienc]>largest_runlength_base)
1045 largest_runlength_base=encode_ints[ienc];
1046 largest_runlength_index=Ptngc_find_magic_index(largest_runlength_base);
1047 if (largest_runlength_index!=new_small_index)
1048 {
1049 iter_small_index=largest_runlength_index;
1050 #ifdef SHOWIT
1051 fprintf(stderr,"I found a new possible small index: %d Base=%d\n",iter_small_index,magic[iter_small_index]);
1052 #endif
1053 }
1054 } while ((new_runlength!=iter_runlength) ||
1055 (new_small_index!=iter_small_index));
1056
1057 #ifdef SHOWIT
1058 fprintf(stderr,"Exit iterative loop.\n");
1059 fflush(stderr);
1060 #endif
1061
1062 /* Verify that we got something good. We may have caught a
1063 substantially larger atom. If so we should just bail
1064 out and let the loop get on another lap. We may have a
1065 minimum runlength though and then we have to fulfill
1066 the request to write out these atoms! */
1067 rle_index_dep=0;
1068 if (new_runlength<3)
1069 rle_index_dep=IS_LARGE;
1070 else if (new_runlength<6)
1071 rle_index_dep=QUITE_LARGE;
1072 if ((min_runlength)
1073 || ((new_small_index<small_index+IS_LARGE) && (new_small_index+rle_index_dep<max_large_index))
1074 #if 1
1075 || (new_small_index+IS_LARGE<max_large_index)
1076 #endif
1077 )
1078 {
1079 int nbits;
1080 if ((new_runlength!=runlength) || (new_small_index!=small_index))
1081 {
1082 int change=new_small_index-small_index;
1083
1084 if (new_small_index<=0)
1085 change=0;
1086
1087 if (change<0)
1088 {
1089 int ixx;
1090 for (ixx=0; ixx<new_runlength; ixx++)
1091 {
1092 int rejected;
1093 do {
1094 int ixyz;
1095 double isum=0.; /* ints can be almost 32 bit so multiplication will overflow. So do doubles. */
1096 for (ixyz=0; ixyz<3; ixyz++)
1097 {
1098 /* encode_ints is already positive (and multiplied by 2 versus the original, just as magic ints) */
1099 double id=encode_ints[ixx*3+ixyz];
1100 isum+=id*id;
1101 }
1102 rejected=0;
1103 #ifdef SHOWIT
1104 fprintf(stderr,"Tested decrease %d of index: %g>=%g?\n",change,isum,(double)magic[small_index+change]*(double)magic[small_index+change]);
1105 #endif
1106 if (isum>(double)magic[small_index+change]*(double)magic[small_index+change])
1107 {
1108 #ifdef SHOWIT
1109 fprintf(stderr,"Rejected decrease %d of index due to length of vector: %g>=%g\n",change,isum,(double)magic[small_index+change]*(double)magic[small_index+change]);
1110 #endif
1111 rejected=1;
1112 change++;
1113 }
1114 } while ((change<0) && (rejected));
1115 if (change==0)
1116 break;
1117 }
1118 }
1119
1120 /* If the only thing to do is to change the base by
1121 only one -1 it is probably not worth it. */
1122 if (!((change==-1) && (runlength==new_runlength)))
1123 {
1124 /* If we have a very short runlength we do not
1125 want to do large base changes. It costs 6
1126 extra bits to do -2. We gain 2/3
1127 bits per value to decrease the index by -2,
1128 ie 2 bits, so to any changes down we must
1129 have a runlength of 3 or more to do it for
1130 one molecule! If we have several molecules we
1131 will gain of course, so not be so strict. */
1132 if ((change==-2) && (new_runlength<3))
1133 {
1134 if (runlength==new_runlength)
1135 change=0;
1136 else
1137 change=-1;
1138 #ifdef SHOWIT
1139 fprintf(stderr,"Rejected change by -2 due to too short runlenght. Change set to %d\n",change);
1140 #endif
1141 }
1142
1143 /* First adjust base using large base change instruction (if necessary) */
1144 while ((change>1) || (change<-1) || ((new_runlength==6) && (change)))
1145 {
1146 unsigned int code=0U;
1147 int this_change=change;
1148 if (this_change>2)
1149 this_change=2;
1150 if (this_change<-2)
1151 this_change=-2;
1152 change-=this_change;
1153 #ifdef SHOWIT
1154 fprintf(stderr,"Large base change: %d.\n",this_change);
1155 #endif
1156 small_index+=this_change;
1157 if (this_change<0)
1158 {
1159 code|=2U;
1160 this_change=-this_change;
1161 }
1162 code|=(unsigned int)(this_change-1);
1163 write_instruction(coder,INSTR_LARGE_BASE_CHANGE,&output_ptr);
1164 Ptngc_writebits(coder,code,2,&output_ptr);
1165 }
1166 /* If there is still base change or runlength changes to do, we do them now. */
1167 if ((new_runlength!=runlength) || (change))
1168 {
1169 unsigned int ichange=(unsigned int)(change+1);
1170 unsigned int code=0U;
1171 unsigned int irun=(unsigned int)((new_runlength-1)*3);
1172 if (new_runlength==6)
1173 ichange=0; /* Means no change. The change has been taken care of explicitly by a large
1174 base change instruction above. */
1175 code=ichange+irun;
1176 #ifdef SHOWIT
1177 fprintf(stderr,"Small base change: %d Runlength change: %d\n",change,new_runlength);
1178 #endif
1179 small_index+=change;
1180 write_instruction(coder,INSTR_BASE_RUNLENGTH,&output_ptr);
1181 Ptngc_writebits(coder,code,4,&output_ptr);
1182 runlength=new_runlength;
1183 }
1184 }
1185 #ifdef SHOWIT
1186 else
1187 fprintf(stderr,"Rejected base change due to only change==-1\n");
1188 #endif
1189 #ifdef SHOWIT
1190 fprintf(stderr,"Current small index: %d Base=%d\n",small_index,magic[small_index]);
1191 #endif
1192 }
1193 /* If we have a large previous integer we can combine it with a sequence of small ints. */
1194 if (has_large)
1195 {
1196 /* If swapatoms is set to 1 but we did actually not
1197 do any swapping, we must first write out the
1198 large atom and then the small. If swapatoms is 1
1199 and we did swapping we can use the efficient
1200 encoding. */
1201 if ((swapatoms) && (!didswap))
1202 {
1203 #ifdef SHOWIT
1204 fprintf(stderr,"Swapatoms was set to 1 but we did not do swapping!\n");
1205 fprintf(stderr,"Only one large integer.\n");
1206 #endif
1207 /* Flush all large atoms. */
1208 flush_large(coder,&has_large,has_large_ints,has_large,large_index,large_nbits,compress_buffer,&output_ptr);
1209 #ifdef SHOWIT
1210 fprintf(stderr,"Sequence of only small integers.\n");
1211 #endif
1212 write_instruction(coder,INSTR_ONLY_SMALL,&output_ptr);
1213 }
1214 else
1215 {
1216
1217 #ifdef SHOWIT
1218 fprintf(stderr,"Sequence of one large and small integers (good compression).\n");
1219 #endif
1220 /* Flush all large atoms but one! */
1221 if (has_large>1)
1222 flush_large(coder,&has_large,has_large_ints,has_large-1,large_index,large_nbits,compress_buffer,&output_ptr);
1223 write_instruction(coder,INSTR_DEFAULT,&output_ptr);
1224 write_three_large(coder,has_large_ints,large_index,large_nbits,compress_buffer,&output_ptr);
1225 has_large=0;
1226 }
1227 }
1228 else
1229 {
1230 #ifdef SHOWIT
1231 fprintf(stderr,"Sequence of only small integers.\n");
1232 #endif
1233 write_instruction(coder,INSTR_ONLY_SMALL,&output_ptr);
1234 }
1235 /* Base compress small integers using the current parameters. */
1236 nbits=magic_bits[small_index][runlength-1];
1237 /* The same base is used for the small changes. */
1238 small_idx[0]=small_index;
1239 small_idx[1]=small_index;
1240 small_idx[2]=small_index;
1241 trajcoder_base_compress(encode_ints,runlength*3,small_idx,compress_buffer);
1242 #ifdef SHOWIT
1243 fprintf(stderr,"nbits=%d (%g)\n",nbits,nbits/(runlength*3.));
1244 nbits_sum+=nbits;
1245 nvalues_sum+=runlength*3;
1246 fprintf(stderr,"Runlength encoded small integers. runlength=%d\n",runlength);
1247 #endif
1248 /* write out base compressed small integers */
1249 Ptngc_writemanybits(coder,compress_buffer,nbits,&output_ptr);
1250 #ifdef SHOWIT
1251 for (ienc=0; ienc<runlength; ienc++)
1252 fprintf(stderr,"Small: %d %d %d\n",
1253 encode_ints[ienc*3],
1254 encode_ints[ienc*3+1],
1255 encode_ints[ienc*3+2]);
1256 #endif
1257 /* Update prevcoord. */
1258 for (ienc=0; ienc<runlength; ienc++)
1259 {
1260 #ifdef SHOWIT
1261 fprintf(stderr,"Prevcoord in packing: %d %d %d\n",
1262 prevcoord[0],prevcoord[1],prevcoord[2]);
1263 #endif
1264 prevcoord[0]+=unpositive_int(encode_ints[ienc*3]);
1265 prevcoord[1]+=unpositive_int(encode_ints[ienc*3+1]);
1266 prevcoord[2]+=unpositive_int(encode_ints[ienc*3+2]);
1267 }
1268 #ifdef SHOWIT
1269 fprintf(stderr,"Prevcoord in packing: %d %d %d\n",
1270 prevcoord[0],prevcoord[1],prevcoord[2]);
1271 #endif
1272
1273 input_ptr+=3*runlength;
1274 ntriplets_left-=runlength;
1275 }
1276 else
1277 {
1278 #ifdef SHOWIT
1279 fprintf(stderr,"Refused value: %d old is %d max is %d\n",new_small_index,small_index,max_large_index);
1280 fflush(stderr);
1281 #endif
1282 refused=1;
1283 }
1284 }
1285 #ifdef SHOWIT
1286 fprintf(stderr,"Number of triplets left is %d\n",ntriplets_left);
1287 #endif
1288 }
1289 /* If we have large previous integers we must flush them now. */
1290 if (has_large)
1291 flush_large(coder,&has_large,has_large_ints,has_large,large_index,large_nbits,compress_buffer,&output_ptr);
1292 Ptngc_pack_flush(coder,&output_ptr);
1293 output_length=(int)(output_ptr-output);
1294 #ifdef SHOWIT
1295 fprintf(stderr,"Done block: nbits=%d nvalues=%d (%g)\n",nbits_sum,nvalues_sum,(double)nbits_sum/nvalues_sum);
1296 #endif
1297 *length=output_length;
1298 return output;
1299 }
1300
1301
1302 int Ptngc_unpack_array_xtc2(struct coder *coder, unsigned char *packed, int *output, const int length)
1303 {
1304 unsigned char *ptr=packed;
1305 int bitptr=0;
1306 int minint[3];
1307 int large_index[3];
1308 int small_index;
1309 int prevcoord[3];
1310 int ntriplets_left=length/3;
1311 int swapatoms=0;
1312 int runlength=0;
1313 int large_nbits;
1314 unsigned char compress_buffer[18*4]; /* Holds compressed result for 3 large ints or up to 18 small ints. */
1315 int encode_ints[21]; /* Up to 3 large + 18 small ints can be encoded at once */
1316 (void)coder;
1317
1318 /* Read min integers. */
1319 minint[0]=unpositive_int(readbits(&ptr,&bitptr,32));
1320 minint[1]=unpositive_int(readbits(&ptr,&bitptr,32));
1321 minint[2]=unpositive_int(readbits(&ptr,&bitptr,32));
1322 /* Read large indices */
1323 large_index[0]=readbits(&ptr,&bitptr,8);
1324 large_index[1]=readbits(&ptr,&bitptr,8);
1325 large_index[2]=readbits(&ptr,&bitptr,8);
1326 /* Read small index */
1327 small_index=readbits(&ptr,&bitptr,8);
1328
1329 large_nbits=compute_magic_bits(large_index);
1330
1331 #ifdef SHOWIT
1332 fprintf(stderr,"Minimum integers: %d %d %d\n",minint[0],minint[1],minint[2]);
1333 fprintf(stderr,"Large indices: %d %d %d\n",large_index[0],large_index[1],large_index[2]);
1334 fprintf(stderr,"Small index: %d\n",small_index);
1335 fprintf(stderr,"large_nbits=%d\n",large_nbits);
1336 #endif
1337
1338 /* Initial prevcoord is the minimum integers. */
1339 memcpy(prevcoord, minint, 3*sizeof *prevcoord);
1340
1341 while (ntriplets_left)
1342 {
1343 int instr=read_instruction(&ptr,&bitptr);
1344 #ifdef SHOWIT
1345 if ((instr>=0) && (instr<MAXINSTR))
1346 fprintf(stderr,"Decoded instruction %s\n",instrnames[instr]);
1347 #endif
1348 if ((instr==INSTR_DEFAULT) /* large+small */
1349 || (instr==INSTR_ONLY_LARGE) /* only large */
1350 || (instr==INSTR_ONLY_SMALL)) /* only small */
1351 {
1352 int large_ints[3]={0,0,0};
1353 if (instr!=INSTR_ONLY_SMALL)
1354 {
1355 /* Clear the compress buffer. */
1356 int i;
1357 for (i=0; i<18*4; i++)
1358 compress_buffer[i]=0;
1359 /* Get the large value. */
1360 readmanybits(&ptr,&bitptr,large_nbits,compress_buffer);
1361 trajcoder_base_decompress(compress_buffer,3,large_index,encode_ints);
1362 memcpy(large_ints, encode_ints, 3*sizeof *large_ints);
1363 #ifdef SHOWIT
1364 fprintf(stderr,"large ints: %d %d %d\n",large_ints[0],large_ints[1],large_ints[2]);
1365 #endif
1366 }
1367 if (instr!=INSTR_ONLY_LARGE)
1368 {
1369 int small_idx[3];
1370 int i;
1371 /* The same base is used for the small changes. */
1372 small_idx[0]=small_index;
1373 small_idx[1]=small_index;
1374 small_idx[2]=small_index;
1375 /* Clear the compress buffer. */
1376 for (i=0; i<18*4; i++)
1377 compress_buffer[i]=0;
1378 /* Get the small values. */
1379 readmanybits(&ptr,&bitptr,magic_bits[small_index][runlength-1],compress_buffer);
1380 trajcoder_base_decompress(compress_buffer,3*runlength,small_idx,encode_ints);
1381 #ifdef SHOWIT
1382 for (i=0; i<runlength; i++)
1383 fprintf(stderr,"small ints: %d %d %d\n",encode_ints[i*3+0],encode_ints[i*3+1],encode_ints[i*3+2]);
1384 #endif
1385 }
1386 if (instr==INSTR_DEFAULT)
1387 {
1388 /* Check for swapped atoms */
1389 if (swapatoms)
1390 {
1391 /* Unswap the atoms. */
1392 int i;
1393 for (i=0; i<3; i++)
1394 {
1395 int in[3], out[3];
1396 in[0]=large_ints[i];
1397 in[1]=unpositive_int(encode_ints[i]);
1398 in[2]=unpositive_int(encode_ints[3+i]);
1399 swap_ints(in,out);
1400 large_ints[i]=out[0];
1401 encode_ints[i]=positive_int(out[1]);
1402 encode_ints[3+i]=positive_int(out[2]);
1403 }
1404 }
1405 }
1406 /* Output result. */
1407 if (instr!=INSTR_ONLY_SMALL)
1408 {
1409 /* Output large value */
1410 *output++=large_ints[0]+minint[0];
1411 *output++=large_ints[1]+minint[1];
1412 *output++=large_ints[2]+minint[2];
1413 prevcoord[0]=large_ints[0];
1414 prevcoord[1]=large_ints[1];
1415 prevcoord[2]=large_ints[2];
1416 #ifdef SHOWIT
1417 fprintf(stderr,"Prevcoord after unpacking of large: %d %d %d\n",
1418 prevcoord[0],prevcoord[1],prevcoord[2]);
1419 fprintf(stderr,"VALUE:%d %6d %6d %6d\n",
1420 length/3-ntriplets_left,
1421 prevcoord[0]+minint[0],
1422 prevcoord[1]+minint[1],
1423 prevcoord[2]+minint[2]);
1424 #endif
1425 ntriplets_left--;
1426 }
1427 if (instr!=INSTR_ONLY_LARGE)
1428 {
1429 /* Output small values */
1430 int i;
1431 #ifdef SHOWIT
1432 fprintf(stderr,"Prevcoord before unpacking of small: %d %d %d\n",
1433 prevcoord[0],prevcoord[1],prevcoord[2]);
1434 #endif
1435 for (i=0; i<runlength; i++)
1436 {
1437 int v[3];
1438 v[0]=unpositive_int(encode_ints[i*3]);
1439 v[1]=unpositive_int(encode_ints[i*3+1]);
1440 v[2]=unpositive_int(encode_ints[i*3+2]);
1441 prevcoord[0]+=v[0];
1442 prevcoord[1]+=v[1];
1443 prevcoord[2]+=v[2];
1444 #ifdef SHOWIT
1445 fprintf(stderr,"Prevcoord after unpacking of small: %d %d %d\n",
1446 prevcoord[0],prevcoord[1],prevcoord[2]);
1447 fprintf(stderr,"Unpacked small values: %6d %6d %6d\t\t%6d %6d %6d\n",v[0],v[1],v[2],prevcoord[0],prevcoord[1],prevcoord[2]);
1448 fprintf(stderr,"VALUE:%d %6d %6d %6d\n",
1449 length/3-(ntriplets_left-i),
1450 prevcoord[0]+minint[0],
1451 prevcoord[1]+minint[1],
1452 prevcoord[2]+minint[2]);
1453 #endif
1454 *output++=prevcoord[0]+minint[0];
1455 *output++=prevcoord[1]+minint[1];
1456 *output++=prevcoord[2]+minint[2];
1457 }
1458 ntriplets_left-=runlength;
1459 }
1460 }
1461 else if (instr==INSTR_LARGE_RLE)
1462 {
1463 int i,j;
1464 int large_ints[3];
1465 /* How many large atoms in this sequence? */
1466 int n=(int)readbits(&ptr,&bitptr,4)+3; /* 3-18 large atoms */
1467 for (i=0; i<n; i++)
1468 {
1469 /* Clear the compress buffer. */
1470 for (j=0; j<18*4; j++)
1471 compress_buffer[j]=0;
1472 /* Get the large value. */
1473 readmanybits(&ptr,&bitptr,large_nbits,compress_buffer);
1474 trajcoder_base_decompress(compress_buffer,3,large_index,encode_ints);
1475 memcpy(large_ints, encode_ints, 3*sizeof *large_ints);
1476 /* Output large value */
1477 *output++=large_ints[0]+minint[0];
1478 *output++=large_ints[1]+minint[1];
1479 *output++=large_ints[2]+minint[2];
1480 memcpy(prevcoord, large_ints, 3*sizeof *prevcoord);
1481 }
1482 ntriplets_left-=n;
1483 }
1484 else if (instr==INSTR_BASE_RUNLENGTH)
1485 {
1486 unsigned int code=readbits(&ptr,&bitptr,4);
1487 int change;
1488 if (code==15)
1489 {
1490 change=0;
1491 runlength=6;
1492 }
1493 else
1494 {
1495 int ichange=code%3;
1496 runlength=code/3+1;
1497 change=ichange-1;
1498 }
1499 small_index+=change;
1500 }
1501 else if (instr==INSTR_FLIP)
1502 {
1503 swapatoms=1-swapatoms;
1504 }
1505 else if (instr==INSTR_LARGE_BASE_CHANGE)
1506 {
1507 unsigned int ichange=readbits(&ptr,&bitptr,2);
1508 int change=(int)(ichange&0x1U)+1;
1509 if (ichange&0x2U)
1510 change=-change;
1511 small_index+=change;
1512 }
1513 else
1514 {
1515 fprintf(stderr,"TRAJNG: BUG! Encoded unknown instruction.\n");
1516 exit(EXIT_FAILURE);
1517 }
1518 #ifdef SHOWIT
1519 fprintf(stderr,"Number of triplets left is %d\n",ntriplets_left);
1520 #endif
1521 }
1522 return 0;
1523 }
The ultimate molecular dynamics simulation package