1 /* This code is part of the tng compression routines.
3 * Written by Daniel Spangberg and Magnus Lundborg
4 * Copyright (c) 2010, 2013-2014 The GROMACS development team.
7 * This program is free software; you can redistribute it and/or
8 * modify it under the terms of the Revised BSD License.
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
18 /* The cost estimates are ripped right out of xtc2.c, so take these
19 with a grain (truckload) of salt. */
26 #include "../../include/compression/warnmalloc.h"
27 #include "../../include/compression/widemuldiv.h"
28 #include "../../include/compression/bwlzh.h"
30 static const double iflipgaincheck=0.89089871814033927; /* 1./(2**(1./6)) */
32 #define MAX_LARGE_RLE 1024 /* Maximum number of large atoms for large RLE. */
33 #define MAX_SMALL_RLE 12 /* Maximum number of small atoms in one group. */
35 #define TRESHOLD_INTRA_INTER_DIRECT 1.5 /* How much larger can the direct
36 frame deltas for the small
37 triplets be and be accepted anyway
38 as better than the intra/inter frame
39 deltas. For better instructions/RLEs. */
41 #define TRESHOLD_INTER_INTRA 5.0 /* How much larger can the intra
42 frame deltas for the small
43 triplets be and be accepted anyway
44 as better than the inter frame
47 /* Difference in indices used for determining whether to store as
48 large or small. A fun detail in this compression algorithm is that
49 if everything works fine, large can often be smaller than small, or
50 at least not as large as is large in magic.c. This is a key idea of
64 #define TNG_INLINE __inline
66 #define TNG_INLINE inline
69 /* These routines are in xtc2.c */
70 int Ptngc_magic(unsigned int i);
71 int Ptngc_find_magic_index(const unsigned int maxval);
73 static TNG_INLINE unsigned int positive_int(const int item)
83 static TNG_INLINE int unpositive_int(const int val)
92 /* Sequence instructions */
93 #define INSTR_DEFAULT 0U
94 #define INSTR_SMALL_RUNLENGTH 1U
95 #define INSTR_ONLY_LARGE 2U
96 #define INSTR_ONLY_SMALL 3U
98 #define INSTR_LARGE_RLE 5U
99 #define INSTR_LARGE_DIRECT 6U
100 #define INSTR_LARGE_INTRA_DELTA 7U
101 #define INSTR_LARGE_INTER_DELTA 8U
107 unsigned int *instructions;
108 int ninstr, ninstr_alloc;
110 int nrle, nrle_alloc;
111 unsigned int *large_direct;
112 int nlargedir, nlargedir_alloc;
113 unsigned int *large_intra_delta;
114 int nlargeintra, nlargeintra_alloc;
115 unsigned int *large_inter_delta;
116 int nlargeinter, nlargeinter_alloc;
117 unsigned int *smallintra;
118 int nsmallintra, nsmallintra_alloc;
119 int minint[3],maxint[3];
121 int has_large_ints[MAX_LARGE_RLE*3]; /* Large cache. */
122 int has_large_type[MAX_LARGE_RLE]; /* What kind of type this large
124 int current_large_type;
127 static void init_xtc3_context(struct xtc3_context *xtc3_context)
129 xtc3_context->instructions=NULL;
130 xtc3_context->ninstr=0;
131 xtc3_context->ninstr_alloc=0;
132 xtc3_context->rle=NULL;
133 xtc3_context->nrle=0;
134 xtc3_context->nrle_alloc=0;
135 xtc3_context->large_direct=NULL;
136 xtc3_context->nlargedir=0;
137 xtc3_context->nlargedir_alloc=0;
138 xtc3_context->large_intra_delta=NULL;
139 xtc3_context->nlargeintra=0;
140 xtc3_context->nlargeintra_alloc=0;
141 xtc3_context->large_inter_delta=NULL;
142 xtc3_context->nlargeinter=0;
143 xtc3_context->nlargeinter_alloc=0;
144 xtc3_context->smallintra=NULL;
145 xtc3_context->nsmallintra=0;
146 xtc3_context->nsmallintra_alloc=0;
147 xtc3_context->has_large=0;
148 xtc3_context->current_large_type=0;
151 static void free_xtc3_context(struct xtc3_context *xtc3_context)
153 free(xtc3_context->instructions);
154 free(xtc3_context->rle);
155 free(xtc3_context->large_direct);
156 free(xtc3_context->large_intra_delta);
157 free(xtc3_context->large_inter_delta);
158 free(xtc3_context->smallintra);
161 /* Modifies three integer values for better compression of water */
162 static void swap_ints(const int *in, int *out)
169 static void swap_is_better(const int *input, const int *minint, int *sum_normal, int *sum_swapped)
178 normal[0]=input[i]-minint[i];
179 normal[1]=input[3+i]-input[i]; /* minint[i]-minint[i] cancels out */
180 normal[2]=input[6+i]-input[3+i]; /* minint[i]-minint[i] cancels out */
181 swap_ints(normal,swapped);
184 if (positive_int(normal[j])>(unsigned int)normal_max)
185 normal_max=positive_int(normal[j]);
186 if (positive_int(swapped[j])>(unsigned int)swapped_max)
187 swapped_max=positive_int(swapped[j]);
194 *sum_normal=normal_max;
195 *sum_swapped=swapped_max;
198 static void allocate_enough_memory(unsigned int **ptr, int *nele, int *nele_alloc)
201 if (*nele>*nele_alloc)
203 *nele_alloc=*nele + *nele/2;
204 *ptr=warnrealloc(*ptr,*nele_alloc*sizeof **ptr);
208 static void insert_value_in_array(unsigned int **ptr, int *nele, int *nele_alloc,
209 const unsigned int value,
215 allocate_enough_memory(ptr,nele,nele_alloc);
217 fprintf(stderr,"Inserting value %u into array %s @ %d\n",value,arrayname,(*nele)-1);
219 (*ptr)[(*nele)-1]=value;
224 static void swapdecide(struct xtc3_context *xtc3_context, int *input,int *swapatoms, int *large_index, int *minint)
229 swap_is_better(input,minint,&normal,&swapped);
230 /* We have to determine if it is worth to change the behaviour.
231 If diff is positive it means that it is worth something to
232 swap. But it costs 4 bits to do the change. If we assume that
233 we gain 0.17 bit by the swap per value, and the runlength>2
234 for four molecules in a row, we gain something. So check if we
235 gain at least 0.17 bits to even attempt the swap.
238 fprintf(stderr,"Trying Flip: %g %g\n",(double)swapped/normal, (double)normal/swapped);
240 if (((swapped<normal) && (fabs((double)swapped/normal)<iflipgaincheck)) ||
241 ((normal<swapped) && (fabs((double)normal/swapped)<iflipgaincheck)))
263 fprintf(stderr,"Flip. Swapatoms is now %d\n",*swapatoms);
265 insert_value_in_array(&xtc3_context->instructions,
266 &xtc3_context->ninstr,
267 &xtc3_context->ninstr_alloc,
272 /* It is "large" if we have to increase the small index quite a
273 bit. Not so much to be rejected by the not very large check
275 static int is_quite_large(const int *input, const int small_index, const int max_large_index)
279 if (small_index+QUITE_LARGE>=max_large_index)
284 if (positive_int(input[i])>(unsigned int)Ptngc_magic(small_index+QUITE_LARGE))
298 static void insert_batch(const int *input_ptr, const int ntriplets_left, const int *prevcoord, int *encode_ints, const int startenc, int *nenc)
300 int nencode=startenc*3;
301 int tmp_prevcoord[3];
303 tmp_prevcoord[0]=prevcoord[0];
304 tmp_prevcoord[1]=prevcoord[1];
305 tmp_prevcoord[2]=prevcoord[2];
310 for (i=0; i<startenc; i++)
312 tmp_prevcoord[0]+=encode_ints[i*3];
313 tmp_prevcoord[1]+=encode_ints[i*3+1];
314 tmp_prevcoord[2]+=encode_ints[i*3+2];
316 fprintf(stderr,"%6d: %6d %6d %6d\t\t%6d %6d %6d\t\t%6d %6d %6d\n",i*3,
317 tmp_prevcoord[0],tmp_prevcoord[1],tmp_prevcoord[2],
321 positive_int(encode_ints[i*3]),
322 positive_int(encode_ints[i*3+1]),
323 positive_int(encode_ints[i*3+2]));
329 fprintf(stderr,"New batch\n");
331 while ((nencode<3+MAX_SMALL_RLE*3) && (nencode<ntriplets_left*3))
333 encode_ints[nencode]=input_ptr[nencode]-tmp_prevcoord[0];
334 encode_ints[nencode+1]=input_ptr[nencode+1]-tmp_prevcoord[1];
335 encode_ints[nencode+2]=input_ptr[nencode+2]-tmp_prevcoord[2];
337 fprintf(stderr,"%6d: %6d %6d %6d\t\t%6d %6d %6d\t\t%6d %6d %6d\n",nencode,
339 input_ptr[nencode+1],
340 input_ptr[nencode+2],
341 encode_ints[nencode],
342 encode_ints[nencode+1],
343 encode_ints[nencode+2],
344 positive_int(encode_ints[nencode]),
345 positive_int(encode_ints[nencode+1]),
346 positive_int(encode_ints[nencode+2]));
348 tmp_prevcoord[0]=input_ptr[nencode];
349 tmp_prevcoord[1]=input_ptr[nencode+1];
350 tmp_prevcoord[2]=input_ptr[nencode+2];
356 static void large_instruction_change(struct xtc3_context *xtc3_context, const int i)
358 /* If the first large is of a different kind than the currently used we must
359 emit an "instruction" to change the large type. */
360 if (xtc3_context->has_large_type[i]!=xtc3_context->current_large_type)
363 xtc3_context->current_large_type=xtc3_context->has_large_type[i];
364 if (xtc3_context->current_large_type==0)
365 instr=INSTR_LARGE_DIRECT;
366 else if (xtc3_context->current_large_type==1)
367 instr=INSTR_LARGE_INTRA_DELTA;
369 instr=INSTR_LARGE_INTER_DELTA;
370 insert_value_in_array(&xtc3_context->instructions,
371 &xtc3_context->ninstr,
372 &xtc3_context->ninstr_alloc,
377 static void write_three_large(struct xtc3_context *xtc3_context,
381 if (xtc3_context->current_large_type==0)
384 insert_value_in_array(&xtc3_context->large_direct,
385 &xtc3_context->nlargedir,
386 &xtc3_context->nlargedir_alloc,
387 xtc3_context->has_large_ints[i*3+m],"large direct");
389 else if (xtc3_context->current_large_type==1)
392 insert_value_in_array(&xtc3_context->large_intra_delta,
393 &xtc3_context->nlargeintra,
394 &xtc3_context->nlargeintra_alloc,
395 xtc3_context->has_large_ints[i*3+m],"large intra");
400 insert_value_in_array(&xtc3_context->large_inter_delta,
401 &xtc3_context->nlargeinter,
402 &xtc3_context->nlargeinter_alloc,
403 xtc3_context->has_large_ints[i*3+m],"large inter");
407 static void flush_large(struct xtc3_context *xtc3_context,
408 const int n) /* How many to flush. */
415 /* If the first large is of a different kind than the currently used we must
416 emit an "instruction" to change the large type. */
417 large_instruction_change(xtc3_context,i);
418 /* How many large of the same kind in a row? */
421 (xtc3_context->has_large_type[i+j]==xtc3_context->has_large_type[i]);
427 insert_value_in_array(&xtc3_context->instructions,
428 &xtc3_context->ninstr,
429 &xtc3_context->ninstr_alloc,
430 INSTR_ONLY_LARGE,"instr");
431 write_three_large(xtc3_context,i+k);
436 insert_value_in_array(&xtc3_context->instructions,
437 &xtc3_context->ninstr,
438 &xtc3_context->ninstr_alloc,
439 INSTR_LARGE_RLE,"instr");
440 insert_value_in_array(&xtc3_context->rle,
442 &xtc3_context->nrle_alloc,
443 (unsigned int)j,"rle (large)");
445 write_three_large(xtc3_context,i+k);
449 if ((xtc3_context->has_large-n)!=0)
452 for (i=0; i<xtc3_context->has_large-n; i++)
454 xtc3_context->has_large_type[i]=xtc3_context->has_large_type[i+n];
456 xtc3_context->has_large_ints[i*3+j]=xtc3_context->has_large_ints[(i+n)*3+j];
459 xtc3_context->has_large-=n; /* Number of remaining large atoms in buffer */
462 static double compute_intlen(unsigned int *ints)
464 /* The largest value. */
465 unsigned int m=ints[0];
473 static void buffer_large(struct xtc3_context *xtc3_context, int *input, const int inpdata,
474 const int natoms, const int intradelta_ok)
476 unsigned int direct[3], intradelta[3]={0,}, interdelta[3]={0,};
479 int frame=inpdata/(natoms*3);
480 int atomframe=inpdata%(natoms*3);
481 /* If it is full we must write them all. */
482 if (xtc3_context->has_large==MAX_LARGE_RLE)
483 flush_large(xtc3_context,xtc3_context->has_large); /* Flush all. */
484 /* Find out which is the best choice for the large integer. Direct coding, or some
485 kind of delta coding? */
486 /* First create direct coding. */
487 direct[0]=(unsigned int)(input[inpdata]-xtc3_context->minint[0]);
488 direct[1]=(unsigned int)(input[inpdata+1]-xtc3_context->minint[1]);
489 direct[2]=(unsigned int)(input[inpdata+2]-xtc3_context->minint[2]);
490 minlen=compute_intlen(direct);
491 best_type=0; /* Direct. */
493 /* Then try intra coding if we can. */
494 if ((intradelta_ok) && (atomframe>=3))
497 intradelta[0]=positive_int(input[inpdata]-input[inpdata-3]);
498 intradelta[1]=positive_int(input[inpdata+1]-input[inpdata-2]);
499 intradelta[2]=positive_int(input[inpdata+2]-input[inpdata-1]);
500 thislen=compute_intlen(intradelta);
501 if (thislen*TRESHOLD_INTRA_INTER_DIRECT<minlen)
504 best_type=1; /* Intra delta */
509 /* Then try inter coding if we can. */
513 interdelta[0]=positive_int(input[inpdata]-input[inpdata-natoms*3]);
514 interdelta[1]=positive_int(input[inpdata+1]-input[inpdata-natoms*3+1]);
515 interdelta[2]=positive_int(input[inpdata+2]-input[inpdata-natoms*3+2]);
516 thislen=compute_intlen(interdelta);
517 if (thislen*TRESHOLD_INTRA_INTER_DIRECT<minlen)
519 best_type=2; /* Inter delta */
523 xtc3_context->has_large_type[xtc3_context->has_large]=best_type;
526 xtc3_context->has_large_ints[xtc3_context->has_large*3]=direct[0];
527 xtc3_context->has_large_ints[xtc3_context->has_large*3+1]=direct[1];
528 xtc3_context->has_large_ints[xtc3_context->has_large*3+2]=direct[2];
530 else if (best_type==1)
532 xtc3_context->has_large_ints[xtc3_context->has_large*3]=intradelta[0];
533 xtc3_context->has_large_ints[xtc3_context->has_large*3+1]=intradelta[1];
534 xtc3_context->has_large_ints[xtc3_context->has_large*3+2]=intradelta[2];
536 else if (best_type==2)
538 xtc3_context->has_large_ints[xtc3_context->has_large*3]=interdelta[0];
539 xtc3_context->has_large_ints[xtc3_context->has_large*3+1]=interdelta[1];
540 xtc3_context->has_large_ints[xtc3_context->has_large*3+2]=interdelta[2];
542 xtc3_context->has_large++;
545 static void output_int(unsigned char *output,int *outdata, const unsigned int n)
547 output[(*outdata)++]=((unsigned int)n)&0xFFU;
548 output[(*outdata)++]=(((unsigned int)n)>>8)&0xFFU;
549 output[(*outdata)++]=(((unsigned int)n)>>16)&0xFFU;
550 output[(*outdata)++]=(((unsigned int)n)>>24)&0xFFU;
554 static void printarray(unsigned int *a, int n, char *name)
558 fprintf(stderr,"%u %s\n",a[i],name);
562 /* The base_compress routine first compresses all x coordinates, then
563 y and finally z. The bases used for each can be different. The
564 MAXBASEVALS value determines how many coordinates are compressed
565 into a single number. Only resulting whole bytes are dealt with for
566 simplicity. MAXMAXBASEVALS is the insanely large value to accept
567 files written with that value. BASEINTERVAL determines how often a
568 new base is actually computed and stored in the output
569 file. MAXBASEVALS*BASEINTERVAL values are stored using the same
570 base in BASEINTERVAL different integers. Note that the primarily
571 the decompression using a large MAXBASEVALS becomes very slow. */
572 #define MAXMAXBASEVALS 16384U
573 #define MAXBASEVALS 24U
574 #define BASEINTERVAL 8
576 /* How many bytes are needed to store n values in base base */
577 static int base_bytes(const unsigned int base, const int n)
580 unsigned int largeint[MAXMAXBASEVALS+1];
581 unsigned int largeint_tmp[MAXMAXBASEVALS+1];
584 memset(largeint, 0U, sizeof(unsigned int) * (n+1));
590 Ptngc_largeint_mul(base,largeint,largeint_tmp,n+1);
591 memcpy(largeint, largeint_tmp, (n+1)*sizeof *largeint);
593 Ptngc_largeint_add(base-1U,largeint,n+1);
598 if ((largeint[i]>>(j*8))&0xFFU)
603 static void base_compress(unsigned int *data, const int len, unsigned char *output, int *outlen)
605 unsigned int largeint[MAXBASEVALS+1];
606 unsigned int largeint_tmp[MAXBASEVALS+1];
610 unsigned int numbytes=0;
611 /* Store the MAXBASEVALS value in the output. */
612 output[nwrittenout++]=(unsigned char)(MAXBASEVALS&0xFFU);
613 output[nwrittenout++]=(unsigned char)((MAXBASEVALS>>8)&0xFFU);
614 /* Store the BASEINTERVAL value in the output. */
615 output[nwrittenout++]=(unsigned char)(BASEINTERVAL&0xFFU);
616 for (ixyz=0; ixyz<3; ixyz++)
618 unsigned int base=0U;
622 memset(largeint, 0U, sizeof(unsigned int) * (MAXBASEVALS+1));
624 for (i=ixyz; i<len; i+=3)
633 /* Find the largest value for this particular coordinate. */
634 for (k=i; k<len; k+=3)
639 if (basecheckvals==MAXBASEVALS*BASEINTERVAL)
642 /* The base is one larger than the largest values. */
646 /* Store the base in the output. */
647 output[nwrittenout++]=(unsigned char)(base&0xFFU);
648 output[nwrittenout++]=(unsigned char)((base>>8)&0xFFU);
649 output[nwrittenout++]=(unsigned char)((base>>16)&0xFFU);
650 output[nwrittenout++]=(unsigned char)((base>>24)&0xFFU);
651 basegiven=BASEINTERVAL;
652 /* How many bytes is needed to store MAXBASEVALS values using this base? */
653 numbytes=base_bytes(base,MAXBASEVALS);
657 fprintf(stderr,"Base for %d is %u. I need %d bytes for %d values.\n",ixyz,base,numbytes,MAXBASEVALS);
662 Ptngc_largeint_mul(base,largeint,largeint_tmp,MAXBASEVALS+1);
663 for (j=0; j<MAXBASEVALS+1; j++)
664 largeint[j]=largeint_tmp[j];
666 Ptngc_largeint_add(data[i],largeint,MAXBASEVALS+1);
668 fprintf(stderr,"outputting value %u\n",data[i]);
671 if (nvals==MAXBASEVALS)
674 fprintf(stderr,"Writing largeint: ");
676 for (j=0; j<numbytes; j++)
680 output[nwrittenout++]=(unsigned char)((largeint[ilarge]>>(ibyte*8))&(0xFFU));
682 fprintf(stderr,"%02x",(unsigned int)output[nwrittenout-1]);
686 fprintf(stderr,"\n");
690 memset(largeint, 0U, sizeof(unsigned int) * (MAXBASEVALS+1));
695 numbytes=base_bytes(base,nvals);
697 fprintf(stderr,"Base for %d is %u. I need %d bytes for %d values.\n",ixyz,base,numbytes,nvals);
699 for (j=0; j<numbytes; j++)
703 output[nwrittenout++]=(unsigned char)((largeint[ilarge]>>(ibyte*8))&(0xFFU));
710 static void base_decompress(unsigned char *input, const int len, unsigned int *output)
712 unsigned int largeint[MAXMAXBASEVALS+1];
713 unsigned int largeint_tmp[MAXMAXBASEVALS+1];
715 int maxbasevals=(int)((unsigned int)(input[0])|(((unsigned int)(input[1]))<<8));
716 int baseinterval=(int)input[2];
717 if (maxbasevals>(int)MAXMAXBASEVALS)
719 fprintf(stderr,"Read a larger maxbasevals value from the file than I can handle. Fix"
720 " by increasing MAXMAXBASEVALS to at least %d. Although, this is"
721 " probably a bug in TRAJNG, since MAXMAXBASEVALS should already be insanely large enough.\n",maxbasevals);
725 for (ixyz=0; ixyz<3; ixyz++)
728 int nvals_left=len/3;
731 unsigned int base=0U;
733 fprintf(stderr,"Base for %d is %u. I need %d bytes for %d values.\n",ixyz,base,numbytes,maxbasevals);
740 base=(unsigned int)(input[0])|
741 (((unsigned int)(input[1]))<<8)|
742 (((unsigned int)(input[2]))<<16)|
743 (((unsigned int)(input[3]))<<24);
745 basegiven=baseinterval;
746 /* How many bytes is needed to store maxbasevals values using this base? */
747 numbytes=base_bytes(base,maxbasevals);
750 if (nvals_left<maxbasevals)
752 numbytes=base_bytes(base,nvals_left);
754 fprintf(stderr,"Base for %d is %u. I need %d bytes for %d values.\n",ixyz,base,numbytes,nvals_left);
757 memset(largeint, 0U, sizeof(unsigned int) * (maxbasevals+1));
759 fprintf(stderr,"Reading largeint: ");
761 if (numbytes/4 < maxbasevals+1)
763 for (j=0; j<numbytes; j++)
767 largeint[ilarge]|=((unsigned int)input[j])<<(ibyte*8);
769 fprintf(stderr,"%02x",(unsigned int)input[j]);
774 fprintf(stderr,"\n");
777 /* Do the long division required to get the output values. */
781 for (i=n-1; i>=0; i--)
783 output[outvals+i*3]=Ptngc_largeint_div(base,largeint,largeint_tmp,maxbasevals+1);
784 for (j=0; j<maxbasevals+1; j++)
785 largeint[j]=largeint_tmp[j];
789 fprintf(stderr,"outputting value %u\n",output[outvals+i*3]);
797 /* If a large proportion of the integers are large (More than 10\% are >14 bits) we return 0, otherwise 1 */
798 static int heuristic_bwlzh(unsigned int *ints, const int nints)
802 for (i=0; i<nints; i++)
811 /* Speed selects how careful to try to find the most efficient compression. The BWLZH algo is expensive!
812 Speed <=2 always avoids BWLZH everywhere it is possible.
813 Speed 3 and 4 and 5 use heuristics (check proportion of large value). This should mostly be safe.
814 Speed 5 enables the LZ77 component of BWLZH.
815 Speed 6 always tests if BWLZH is better and if it is uses it. This can be very slow.
817 unsigned char *Ptngc_pack_array_xtc3(int *input, int *length, const int natoms, int speed)
819 unsigned char *output=NULL;
823 int ntriplets=*length/3;
830 int runlength=0; /* Initial runlength. "Stupidly" set to zero for
831 simplicity and explicity */
832 int swapatoms=0; /* Initial guess is that we should not swap the
833 first two atoms in each large+small
835 int didswap; /* Whether swapping was actually done. */
837 int encode_ints[3+MAX_SMALL_RLE*3]; /* Up to 3 large + 24 small ints can be encoded at once */
839 int ntriplets_left=ntriplets;
841 unsigned char *bwlzh_buf=NULL;
843 unsigned char *base_buf=NULL;
846 struct xtc3_context xtc3_context;
847 init_xtc3_context(&xtc3_context);
849 memcpy(xtc3_context.maxint, input, 3*sizeof *xtc3_context.maxint);
850 memcpy(xtc3_context.minint, input, 3*sizeof *xtc3_context.maxint);
852 /* Values of speed should be sane. */
862 /* Allocate enough memory for output */
864 output=warnmalloc(8*48*sizeof *output);
866 output=warnmalloc(8* *length*sizeof *output);
869 for (i=1; i<ntriplets; i++)
872 if (input[i*3+j]>xtc3_context.maxint[j])
873 xtc3_context.maxint[j]=input[i*3+j];
874 if (input[i*3+j]<xtc3_context.minint[j])
875 xtc3_context.minint[j]=input[i*3+j];
878 large_index[0]=Ptngc_find_magic_index(xtc3_context.maxint[0]-xtc3_context.minint[0]+1);
879 large_index[1]=Ptngc_find_magic_index(xtc3_context.maxint[1]-xtc3_context.minint[1]+1);
880 large_index[2]=Ptngc_find_magic_index(xtc3_context.maxint[2]-xtc3_context.minint[2]+1);
881 max_large_index=large_index[0];
882 if (large_index[1]>max_large_index)
883 max_large_index=large_index[1];
884 if (large_index[2]>max_large_index)
885 max_large_index=large_index[2];
889 fprintf(stderr,"minint[%d]=%d. maxint[%d]=%d large_index[%d]=%d value=%d\n",j,xtc3_context.minint[j],j,xtc3_context.maxint[j],
890 j,large_index[j],Ptngc_magic(large_index[j]));
893 /* Guess initial small index */
894 small_index=max_large_index/2;
896 /* Find the largest value that is not large. Not large is half index of
898 max_small=Ptngc_magic(small_index);
900 for (i=0; i<*length; i++)
903 int s=positive_int(item);
908 /* This value is not critical, since if I guess wrong, the code will
909 just insert instructions to increase this value. */
910 small_index=Ptngc_find_magic_index(intmax);
912 fprintf(stderr,"initial small_index=%d value=%d\n",small_index,Ptngc_magic(small_index));
915 output_int(output,&outdata,positive_int(xtc3_context.minint[0]));
916 output_int(output,&outdata,positive_int(xtc3_context.minint[1]));
917 output_int(output,&outdata,positive_int(xtc3_context.minint[2]));
921 for (i=0; i<ntriplets_left; i++)
922 fprintf(stderr,"VALUE:%d %6d %6d %6d\n",
925 input[inpdata+i*3+1],
926 input[inpdata+i*3+2]);
930 /* Initial prevcoord is the minimum integers. */
931 memcpy(prevcoord, xtc3_context.minint, 3*sizeof *prevcoord);
932 prevcoord[0]=xtc3_context.minint[0];
933 prevcoord[1]=xtc3_context.minint[1];
934 prevcoord[2]=xtc3_context.minint[2];
936 while (ntriplets_left)
938 if (ntriplets_left<0)
940 fprintf(stderr,"TRAJNG: BUG! ntriplets_left<0!\n");
943 /* If only less than three atoms left we just write them all as large integers. Here no swapping is done! */
944 if (ntriplets_left<3)
946 for (ienc=0; ienc<ntriplets_left; ienc++)
948 buffer_large(&xtc3_context,input,inpdata,natoms,1);
952 flush_large(&xtc3_context,xtc3_context.has_large); /* Flush all */
957 int largest_required_base;
958 int largest_runlength_base;
959 int largest_required_index;
960 int largest_runlength_index;
964 int iter_small_index;
967 /* Insert the next batch of integers to be encoded into the buffer */
969 fprintf(stderr,"Initial batch\n");
971 insert_batch(input+inpdata,ntriplets_left,prevcoord,encode_ints,0,&nencode);
973 /* First we must decide if the next value is large (does not reasonably fit in current small encoding)
974 Also, if we have not written any values yet, we must begin by writing a large atom. */
975 if ((inpdata==0) || (is_quite_large(encode_ints,small_index,max_large_index)) || (refused))
977 /* If any of the next two atoms are large we should probably write them as large and not swap them */
979 if ((is_quite_large(encode_ints+3,small_index,max_large_index)) || (is_quite_large(encode_ints+6,small_index,max_large_index)))
984 /* If doing inter-frame coding results in smaller values we should not do any swapping either. */
985 int frame=inpdata/(natoms*3);
988 unsigned int delta[3], delta2[3];
989 delta[0]=positive_int(input[inpdata+3]-input[inpdata-natoms*3+3]);
990 delta[1]=positive_int(input[inpdata+4]-input[inpdata-natoms*3+4]);
991 delta[2]=positive_int(input[inpdata+5]-input[inpdata-natoms*3+5]);
992 delta2[0]=positive_int(encode_ints[3]);
993 delta2[1]=positive_int(encode_ints[4]);
994 delta2[2]=positive_int(encode_ints[5]);
996 fprintf(stderr,"A1: inter delta: %u %u %u. intra delta=%u %u %u\n",
997 delta[0],delta[1],delta[2],
998 delta2[0],delta2[1],delta2[2]);
1000 if (compute_intlen(delta)*TRESHOLD_INTER_INTRA<compute_intlen(delta2))
1002 delta[0]=positive_int(input[inpdata+6]-input[inpdata-natoms*3+6]);
1003 delta[1]=positive_int(input[inpdata+7]-input[inpdata-natoms*3+7]);
1004 delta[2]=positive_int(input[inpdata+8]-input[inpdata-natoms*3+8]);
1005 delta2[0]=positive_int(encode_ints[6]);
1006 delta2[1]=positive_int(encode_ints[7]);
1007 delta2[2]=positive_int(encode_ints[8]);
1009 fprintf(stderr,"A2: inter delta: %u %u %u. intra delta=%u %u %u\n",
1010 delta[0],delta[1],delta[2],
1011 delta2[0],delta2[1],delta2[2]);
1013 if (compute_intlen(delta)*TRESHOLD_INTER_INTRA<compute_intlen(delta2))
1017 fprintf(stderr,"SETTING NO SWAP!\n");
1026 /* Next we must decide if we should swap the first
1029 swapdecide(&xtc3_context,input+inpdata,&swapatoms,large_index,xtc3_context.minint);
1033 /* If we should do the integer swapping manipulation we should do it now */
1040 in[0]=input[inpdata+i];
1041 in[1]=input[inpdata+3+i]-input[inpdata+i];
1042 in[2]=input[inpdata+6+i]-input[inpdata+3+i];
1044 encode_ints[i]=out[0];
1045 encode_ints[3+i]=out[1];
1046 encode_ints[6+i]=out[2];
1048 /* We have swapped atoms, so the minimum run-length is 2 */
1050 fprintf(stderr,"Swap atoms results in:\n");
1052 fprintf(stderr,"%d: %6d %6d %6d\t\t%6d %6d %6d\n",i*3,
1056 positive_int(encode_ints[i*3]),
1057 positive_int(encode_ints[i*3+1]),
1058 positive_int(encode_ints[i*3+2]));
1064 /* Cache large value for later possible combination with
1065 a sequence of small integers. */
1066 if ((swapatoms) && (didswap))
1068 buffer_large(&xtc3_context,input,inpdata+3,natoms,0); /* This is a swapped integer, so inpdata is one atom later and
1069 intra coding is not ok. */
1070 for (ienc=0; ienc<3; ienc++)
1071 prevcoord[ienc]=input[inpdata+3+ienc];
1075 buffer_large(&xtc3_context,input,inpdata,natoms,1);
1076 for (ienc=0; ienc<3; ienc++)
1077 prevcoord[ienc]=input[inpdata+ienc];
1082 fprintf(stderr,"Prevcoord after packing of large: %d %d %d\n",
1083 prevcoord[0],prevcoord[1],prevcoord[2]);
1086 /* We have written a large integer so we have one less atoms to worry about */
1092 /* Insert the next batch of integers to be encoded into the buffer */
1094 fprintf(stderr,"Update batch due to large int.\n");
1096 if ((swapatoms) && (didswap))
1098 /* Keep swapped values. */
1100 for (ienc=0; ienc<3; ienc++)
1101 encode_ints[i*3+ienc]=encode_ints[(i+1)*3+ienc];
1103 insert_batch(input+inpdata,ntriplets_left,prevcoord,encode_ints,min_runlength,&nencode);
1105 /* Here we should only have differences for the atom coordinates. */
1106 /* Convert the ints to positive ints */
1107 for (ienc=0; ienc<nencode; ienc++)
1109 int pint=positive_int(encode_ints[ienc]);
1110 encode_ints[ienc]=pint;
1112 /* Now we must decide what base and runlength to do. If we have swapped atoms it will be at least 2.
1113 If even the next atom is large, we will not do anything. */
1114 largest_required_base=0;
1115 /* Determine required base */
1116 for (ienc=0; ienc<min_runlength*3; ienc++)
1117 if (encode_ints[ienc]>largest_required_base)
1118 largest_required_base=encode_ints[ienc];
1119 /* Also compute what the largest base is for the current runlength setting! */
1120 largest_runlength_base=0;
1121 for (ienc=0; (ienc<runlength*3) && (ienc<nencode); ienc++)
1122 if (encode_ints[ienc]>largest_runlength_base)
1123 largest_runlength_base=encode_ints[ienc];
1125 largest_required_index=Ptngc_find_magic_index(largest_required_base);
1126 largest_runlength_index=Ptngc_find_magic_index(largest_runlength_base);
1128 if (largest_required_index<largest_runlength_index)
1130 new_runlength=min_runlength;
1131 new_small_index=largest_required_index;
1135 new_runlength=runlength;
1136 new_small_index=largest_runlength_index;
1139 /* Only allow increase of runlength wrt min_runlength */
1140 if (new_runlength<min_runlength)
1141 new_runlength=min_runlength;
1143 /* If the current runlength is longer than the number of
1144 triplets left stop it from being so. */
1145 if (new_runlength>ntriplets_left)
1146 new_runlength=ntriplets_left;
1148 /* We must at least try to get some small integers going. */
1149 if (new_runlength==0)
1152 new_small_index=small_index;
1155 iter_runlength=new_runlength;
1156 iter_small_index=new_small_index;
1158 /* Iterate to find optimal encoding and runlength */
1160 fprintf(stderr,"Entering iterative loop.\n");
1165 new_runlength=iter_runlength;
1166 new_small_index=iter_small_index;
1169 fprintf(stderr,"Test new_small_index=%d Base=%d\n",new_small_index,Ptngc_magic(new_small_index));
1171 /* What is the largest runlength
1172 we can do with the currently
1173 selected encoding? Also the max supported runlength is MAX_SMALL_RLE triplets! */
1174 for (ienc=0; ienc<nencode && ienc<MAX_SMALL_RLE*3; ienc++)
1176 int test_index=Ptngc_find_magic_index(encode_ints[ienc]);
1177 if (test_index>new_small_index)
1180 if (ienc/3>new_runlength)
1182 iter_runlength=ienc/3;
1184 fprintf(stderr,"I found a new possible runlength: %d\n",iter_runlength);
1188 /* How large encoding do we have to use? */
1189 largest_runlength_base=0;
1190 for (ienc=0; ienc<iter_runlength*3; ienc++)
1191 if (encode_ints[ienc]>largest_runlength_base)
1192 largest_runlength_base=encode_ints[ienc];
1193 largest_runlength_index=Ptngc_find_magic_index(largest_runlength_base);
1194 if (largest_runlength_index!=new_small_index)
1196 iter_small_index=largest_runlength_index;
1198 fprintf(stderr,"I found a new possible small index: %d Base=%d\n",iter_small_index,Ptngc_magic(iter_small_index));
1201 } while ((new_runlength!=iter_runlength) ||
1202 (new_small_index!=iter_small_index));
1205 fprintf(stderr,"Exit iterative loop.\n");
1209 /* Verify that we got something good. We may have caught a
1210 substantially larger atom. If so we should just bail
1211 out and let the loop get on another lap. We may have a
1212 minimum runlength though and then we have to fulfill
1213 the request to write out these atoms! */
1215 if (new_runlength<3)
1216 rle_index_dep=IS_LARGE;
1217 else if (new_runlength<6)
1218 rle_index_dep=QUITE_LARGE;
1220 || ((new_small_index<small_index+IS_LARGE) && (new_small_index+rle_index_dep<max_large_index))
1222 || (new_small_index+IS_LARGE<max_large_index)
1226 /* If doing inter-frame coding of large integers results
1227 in smaller values than the small value we should not
1228 produce a sequence of small values here. */
1229 int frame=inpdata/(natoms*3);
1232 if ((!swapatoms) && (frame>0))
1234 for (i=0; i<new_runlength; i++)
1236 unsigned int delta[3];
1237 unsigned int delta2[3];
1238 delta[0]=positive_int(input[inpdata+i*3]-input[inpdata-natoms*3+i*3]);
1239 delta[1]=positive_int(input[inpdata+i*3+1]-input[inpdata-natoms*3+i*3+1]);
1240 delta[2]=positive_int(input[inpdata+i*3+2]-input[inpdata-natoms*3+i*3+2]);
1241 delta2[0]=positive_int(encode_ints[i*3]);
1242 delta2[1]=positive_int(encode_ints[i*3+1]);
1243 delta2[2]=positive_int(encode_ints[i*3+2]);
1244 if (compute_intlen(delta)*TRESHOLD_INTER_INTRA<compute_intlen(delta2))
1249 /* Most of the values should become smaller, otherwise
1250 we should encode them with intra coding. */
1251 if ((!swapatoms) && (numsmaller>=2*new_runlength/3))
1253 /* Put all the values in large arrays, instead of the small array */
1256 for (i=0; i<new_runlength; i++)
1257 buffer_large(&xtc3_context,input,inpdata+i*3,natoms,1);
1259 prevcoord[i]=input[inpdata+(new_runlength-1)*3+i];
1260 inpdata+=3*new_runlength;
1261 ntriplets_left-=new_runlength;
1266 if ((new_runlength!=runlength) || (new_small_index!=small_index))
1268 int change=new_small_index-small_index;
1270 if (new_small_index<=0)
1276 for (ixx=0; ixx<new_runlength; ixx++)
1281 double isum=0.; /* ints can be almost 32 bit so multiplication will overflow. So do doubles. */
1282 for (ixyz=0; ixyz<3; ixyz++)
1284 /* encode_ints is already positive (and multiplied by 2 versus the original, just as magic ints) */
1285 double id=encode_ints[ixx*3+ixyz];
1290 fprintf(stderr,"Tested decrease %d of index: %g>=%g?\n",change,isum,(double)Ptngc_magic(small_index+change)*(double)Ptngc_magic(small_index+change));
1292 if (isum>(double)Ptngc_magic(small_index+change)*(double)Ptngc_magic(small_index+change))
1295 fprintf(stderr,"Rejected decrease %d of index due to length of vector: %g>=%g\n",change,isum,(double)Ptngc_magic(small_index+change)*(double)Ptngc_magic(small_index+change));
1300 } while ((change<0) && (rejected));
1306 /* Always accept the new small indices here. */
1307 small_index=new_small_index;
1308 /* If we have a new runlength emit it */
1309 if (runlength!=new_runlength)
1311 runlength=new_runlength;
1312 insert_value_in_array(&xtc3_context.instructions,
1313 &xtc3_context.ninstr,
1314 &xtc3_context.ninstr_alloc,
1315 INSTR_SMALL_RUNLENGTH,"instr");
1316 insert_value_in_array(&xtc3_context.rle,
1318 &xtc3_context.nrle_alloc,
1319 (unsigned int)runlength,"rle (small)");
1322 fprintf(stderr,"Current small index: %d Base=%d\n",small_index,Ptngc_magic(small_index));
1325 /* If we have a large previous integer we can combine it with a sequence of small ints. */
1326 if (xtc3_context.has_large)
1328 /* If swapatoms is set to 1 but we did actually not
1329 do any swapping, we must first write out the
1330 large atom and then the small. If swapatoms is 1
1331 and we did swapping we can use the efficient
1333 if ((swapatoms) && (!didswap))
1336 fprintf(stderr,"Swapatoms was set to 1 but we did not do swapping!\n");
1337 fprintf(stderr,"Only one large integer.\n");
1339 /* Flush all large atoms. */
1340 flush_large(&xtc3_context,xtc3_context.has_large);
1342 fprintf(stderr,"Sequence of only small integers.\n");
1344 insert_value_in_array(&xtc3_context.instructions,
1345 &xtc3_context.ninstr,
1346 &xtc3_context.ninstr_alloc,
1347 INSTR_ONLY_SMALL,"instr");
1353 fprintf(stderr,"Sequence of one large and small integers (good compression).\n");
1355 /* Flush all large atoms but one! */
1356 if (xtc3_context.has_large>1)
1357 flush_large(&xtc3_context,xtc3_context.has_large-1);
1359 /* Here we must check if we should emit a large
1360 type change instruction. */
1361 large_instruction_change(&xtc3_context,0);
1363 insert_value_in_array(&xtc3_context.instructions,
1364 &xtc3_context.ninstr,
1365 &xtc3_context.ninstr_alloc,
1366 INSTR_DEFAULT,"instr");
1368 write_three_large(&xtc3_context,0);
1369 xtc3_context.has_large=0;
1375 fprintf(stderr,"Sequence of only small integers.\n");
1377 insert_value_in_array(&xtc3_context.instructions,
1378 &xtc3_context.ninstr,
1379 &xtc3_context.ninstr_alloc,
1380 INSTR_ONLY_SMALL,"instr");
1382 /* Insert the small integers into the small integer array. */
1383 for (ienc=0; ienc<runlength*3; ienc++)
1384 insert_value_in_array(&xtc3_context.smallintra,
1385 &xtc3_context.nsmallintra,
1386 &xtc3_context.nsmallintra_alloc,
1387 (unsigned int)encode_ints[ienc],"smallintra");
1390 for (ienc=0; ienc<runlength; ienc++)
1391 fprintf(stderr,"Small: %d %d %d\n",
1392 encode_ints[ienc*3],
1393 encode_ints[ienc*3+1],
1394 encode_ints[ienc*3+2]);
1396 /* Update prevcoord. */
1397 for (ienc=0; ienc<runlength; ienc++)
1400 fprintf(stderr,"Prevcoord in packing: %d %d %d\n",
1401 prevcoord[0],prevcoord[1],prevcoord[2]);
1403 prevcoord[0]+=unpositive_int(encode_ints[ienc*3]);
1404 prevcoord[1]+=unpositive_int(encode_ints[ienc*3+1]);
1405 prevcoord[2]+=unpositive_int(encode_ints[ienc*3+2]);
1408 fprintf(stderr,"Prevcoord in packing: %d %d %d\n",
1409 prevcoord[0],prevcoord[1],prevcoord[2]);
1412 inpdata+=3*runlength;
1413 ntriplets_left-=runlength;
1421 fprintf(stderr,"Refused value: %d old is %d max is %d\n",new_small_index,small_index,max_large_index);
1428 fprintf(stderr,"Number of triplets left is %d\n",ntriplets_left);
1432 /* If we have large previous integers we must flush them now. */
1433 if (xtc3_context.has_large)
1434 flush_large(&xtc3_context,xtc3_context.has_large);
1436 /* Now it is time to compress all the data in the buffers with the bwlzh or base algo. */
1439 /* Inspect the data. */
1440 printarray(xtc3_context.instructions,xtc3_context.ninstr,"A instr");
1441 printarray(xtc3_context.rle,xtc3_context.nrle,"A rle");
1442 printarray(xtc3_context.large_direct,xtc3_context.nlargedir,"A largedir");
1443 printarray(xtc3_context.large_intra_delta,xtc3_context.nlargeintra,"A largeintra");
1444 printarray(xtc3_context.large_inter_delta,xtc3_context.nlargeinter,"A largeinter");
1445 printarray(xtc3_context.smallintra,xtc3_context.nsmallintra,"A smallintra");
1449 #if defined(SHOWIT) || defined(SHOWIT_LIGHT)
1450 fprintf(stderr,"instructions: %d\n",xtc3_context.ninstr);
1453 #if defined(SHOWIT) || defined(SHOWIT_LIGHT)
1454 #define bwlzh_compress bwlzh_compress_verbose
1455 #define bwlzh_compress_no_lz77 bwlzh_compress_no_lz77_verbose
1458 output_int(output,&outdata,(unsigned int)xtc3_context.ninstr);
1459 if (xtc3_context.ninstr)
1461 bwlzh_buf=warnmalloc(bwlzh_get_buflen(xtc3_context.ninstr));
1463 bwlzh_compress(xtc3_context.instructions,xtc3_context.ninstr,bwlzh_buf,&bwlzh_buf_len);
1465 bwlzh_compress_no_lz77(xtc3_context.instructions,xtc3_context.ninstr,bwlzh_buf,&bwlzh_buf_len);
1466 output_int(output,&outdata,(unsigned int)bwlzh_buf_len);
1467 memcpy(output+outdata,bwlzh_buf,bwlzh_buf_len);
1468 outdata+=bwlzh_buf_len;
1472 #if defined(SHOWIT) || defined(SHOWIT_LIGHT)
1473 fprintf(stderr,"rle: %d\n",xtc3_context.nrle);
1476 output_int(output,&outdata,(unsigned int)xtc3_context.nrle);
1477 if (xtc3_context.nrle)
1479 bwlzh_buf=warnmalloc(bwlzh_get_buflen(xtc3_context.nrle));
1481 bwlzh_compress(xtc3_context.rle,xtc3_context.nrle,bwlzh_buf,&bwlzh_buf_len);
1483 bwlzh_compress_no_lz77(xtc3_context.rle,xtc3_context.nrle,bwlzh_buf,&bwlzh_buf_len);
1484 output_int(output,&outdata,(unsigned int)bwlzh_buf_len);
1485 memcpy(output+outdata,bwlzh_buf,bwlzh_buf_len);
1486 outdata+=bwlzh_buf_len;
1490 #if defined(SHOWIT) || defined(SHOWIT_LIGHT)
1491 fprintf(stderr,"large direct: %d\n",xtc3_context.nlargedir);
1494 output_int(output,&outdata,(unsigned int)xtc3_context.nlargedir);
1495 if (xtc3_context.nlargedir)
1497 if ((speed<=2) || ((speed<=5) && (!heuristic_bwlzh(xtc3_context.large_direct,xtc3_context.nlargedir))))
1500 bwlzh_buf_len=INT_MAX;
1504 bwlzh_buf=warnmalloc(bwlzh_get_buflen(xtc3_context.nlargedir));
1506 bwlzh_compress(xtc3_context.large_direct,xtc3_context.nlargedir,bwlzh_buf,&bwlzh_buf_len);
1508 bwlzh_compress_no_lz77(xtc3_context.large_direct,xtc3_context.nlargedir,bwlzh_buf,&bwlzh_buf_len);
1510 /* If this can be written smaller using base compression we should do that. */
1511 base_buf=warnmalloc((xtc3_context.nlargedir+3)*sizeof(int));
1512 base_compress(xtc3_context.large_direct,xtc3_context.nlargedir,base_buf,&base_buf_len);
1513 #if defined(SHOWIT) || defined(SHOWIT_LIGHT)
1514 fprintf(stderr,"Large direct: Base len=%d. BWLZH len=%d\n",base_buf_len,bwlzh_buf_len);
1516 if (base_buf_len<bwlzh_buf_len)
1518 output[outdata++]=0U;
1519 output_int(output,&outdata,(unsigned int)base_buf_len);
1520 memcpy(output+outdata,base_buf,base_buf_len);
1521 outdata+=base_buf_len;
1525 output[outdata++]=1U;
1526 output_int(output,&outdata,(unsigned int)bwlzh_buf_len);
1527 memcpy(output+outdata,bwlzh_buf,bwlzh_buf_len);
1528 outdata+=bwlzh_buf_len;
1534 #if defined(SHOWIT) || defined(SHOWIT_LIGHT)
1535 fprintf(stderr,"large intra: %d\n",xtc3_context.nlargeintra);
1538 output_int(output,&outdata,(unsigned int)xtc3_context.nlargeintra);
1539 if (xtc3_context.nlargeintra)
1541 if ((speed<=2) || ((speed<=5) && (!heuristic_bwlzh(xtc3_context.large_intra_delta,xtc3_context.nlargeintra))))
1544 bwlzh_buf_len=INT_MAX;
1548 bwlzh_buf=warnmalloc(bwlzh_get_buflen(xtc3_context.nlargeintra));
1550 bwlzh_compress(xtc3_context.large_intra_delta,xtc3_context.nlargeintra,bwlzh_buf,&bwlzh_buf_len);
1552 bwlzh_compress_no_lz77(xtc3_context.large_intra_delta,xtc3_context.nlargeintra,bwlzh_buf,&bwlzh_buf_len);
1554 /* If this can be written smaller using base compression we should do that. */
1555 base_buf=warnmalloc((xtc3_context.nlargeintra+3)*sizeof(int));
1556 base_compress(xtc3_context.large_intra_delta,xtc3_context.nlargeintra,base_buf,&base_buf_len);
1557 #if defined(SHOWIT) || defined(SHOWIT_LIGHT)
1558 fprintf(stderr,"Large intra: Base len=%d. BWLZH len=%d\n",base_buf_len,bwlzh_buf_len);
1560 if (base_buf_len<bwlzh_buf_len)
1562 output[outdata++]=0U;
1563 output_int(output,&outdata,(unsigned int)base_buf_len);
1564 memcpy(output+outdata,base_buf,base_buf_len);
1565 outdata+=base_buf_len;
1569 output[outdata++]=1U;
1570 output_int(output,&outdata,(unsigned int)bwlzh_buf_len);
1571 memcpy(output+outdata,bwlzh_buf,bwlzh_buf_len);
1572 outdata+=bwlzh_buf_len;
1578 #if defined(SHOWIT) || defined(SHOWIT_LIGHT)
1579 fprintf(stderr,"large inter: %d\n",xtc3_context.nlargeinter);
1582 output_int(output,&outdata,(unsigned int)xtc3_context.nlargeinter);
1583 if (xtc3_context.nlargeinter)
1585 if ((speed<=2) || ((speed<=5) && (!heuristic_bwlzh(xtc3_context.large_inter_delta,xtc3_context.nlargeinter))))
1588 bwlzh_buf_len=INT_MAX;
1592 bwlzh_buf=warnmalloc(bwlzh_get_buflen(xtc3_context.nlargeinter));
1594 bwlzh_compress(xtc3_context.large_inter_delta,xtc3_context.nlargeinter,bwlzh_buf,&bwlzh_buf_len);
1596 bwlzh_compress_no_lz77(xtc3_context.large_inter_delta,xtc3_context.nlargeinter,bwlzh_buf,&bwlzh_buf_len);
1598 /* If this can be written smaller using base compression we should do that. */
1599 base_buf=warnmalloc((xtc3_context.nlargeinter+3)*sizeof(int));
1600 base_compress(xtc3_context.large_inter_delta,xtc3_context.nlargeinter,base_buf,&base_buf_len);
1601 #if defined(SHOWIT) || defined(SHOWIT_LIGHT)
1602 fprintf(stderr,"Large inter: Base len=%d. BWLZH len=%d\n",base_buf_len,bwlzh_buf_len);
1604 if (base_buf_len<bwlzh_buf_len)
1606 output[outdata++]=0U;
1607 output_int(output,&outdata,(unsigned int)base_buf_len);
1608 memcpy(output+outdata,base_buf,base_buf_len);
1609 outdata+=base_buf_len;
1613 output[outdata++]=1U;
1614 output_int(output,&outdata,(unsigned int)bwlzh_buf_len);
1615 memcpy(output+outdata,bwlzh_buf,bwlzh_buf_len);
1616 outdata+=bwlzh_buf_len;
1622 #if defined(SHOWIT) || defined(SHOWIT_LIGHT)
1623 fprintf(stderr,"small intra: %d\n",xtc3_context.nsmallintra);
1626 output_int(output,&outdata,(unsigned int)xtc3_context.nsmallintra);
1627 if (xtc3_context.nsmallintra)
1629 if ((speed<=2) || ((speed<=5) && (!heuristic_bwlzh(xtc3_context.smallintra,xtc3_context.nsmallintra))))
1632 bwlzh_buf_len=INT_MAX;
1636 bwlzh_buf=warnmalloc(bwlzh_get_buflen(xtc3_context.nsmallintra));
1638 bwlzh_compress(xtc3_context.smallintra,xtc3_context.nsmallintra,bwlzh_buf,&bwlzh_buf_len);
1640 bwlzh_compress_no_lz77(xtc3_context.smallintra,xtc3_context.nsmallintra,bwlzh_buf,&bwlzh_buf_len);
1642 /* If this can be written smaller using base compression we should do that. */
1643 base_buf=warnmalloc((xtc3_context.nsmallintra+3)*sizeof(int));
1644 base_compress(xtc3_context.smallintra,xtc3_context.nsmallintra,base_buf,&base_buf_len);
1645 #if defined(SHOWIT) || defined(SHOWIT_LIGHT)
1646 fprintf(stderr,"Small intra: Base len=%d. BWLZH len=%d\n",base_buf_len,bwlzh_buf_len);
1648 if (base_buf_len<bwlzh_buf_len)
1650 output[outdata++]=0U;
1651 output_int(output,&outdata,(unsigned int)base_buf_len);
1652 memcpy(output+outdata,base_buf,base_buf_len);
1653 outdata+=base_buf_len;
1657 output[outdata++]=1U;
1658 output_int(output,&outdata,(unsigned int)bwlzh_buf_len);
1659 memcpy(output+outdata,bwlzh_buf,bwlzh_buf_len);
1660 outdata+=bwlzh_buf_len;
1667 free_xtc3_context(&xtc3_context);
1671 static void decompress_bwlzh_block(unsigned char **ptr,
1673 unsigned int **vals)
1675 int bwlzh_buf_len=(int)(((unsigned int)(*ptr)[0]) |
1676 (((unsigned int)(*ptr)[1])<<8) |
1677 (((unsigned int)(*ptr)[2])<<16) |
1678 (((unsigned int)(*ptr)[3])<<24));
1680 *vals=warnmalloc(nvals*sizeof (**vals));
1681 bwlzh_decompress(*ptr,nvals,*vals);
1682 (*ptr)+=bwlzh_buf_len;
1685 static void decompress_base_block(unsigned char **ptr,
1687 unsigned int **vals)
1689 int base_buf_len=(int)(((unsigned int)(*ptr)[0]) |
1690 (((unsigned int)(*ptr)[1])<<8) |
1691 (((unsigned int)(*ptr)[2])<<16) |
1692 (((unsigned int)(*ptr)[3])<<24));
1694 *vals=warnmalloc(nvals*sizeof (**vals));
1695 base_decompress(*ptr,nvals,*vals);
1696 (*ptr)+=base_buf_len;
1699 static void unpack_one_large(struct xtc3_context *xtc3_context,
1700 int *ilargedir, int *ilargeintra,
1701 int *ilargeinter, int *prevcoord,
1702 int *minint, int *output,
1703 const int outdata, const int didswap,
1704 const int natoms, const int current_large_type)
1706 int large_ints[3]={0,0,0};
1707 if (current_large_type==0 && xtc3_context->large_direct)
1709 large_ints[0]=(int)xtc3_context->large_direct[(*ilargedir)]+minint[0];
1710 large_ints[1]=(int)xtc3_context->large_direct[(*ilargedir)+1]+minint[1];
1711 large_ints[2]=(int)xtc3_context->large_direct[(*ilargedir)+2]+minint[2];
1714 else if (current_large_type==1 && xtc3_context->large_intra_delta)
1716 large_ints[0]=unpositive_int(xtc3_context->large_intra_delta[(*ilargeintra)])+prevcoord[0];
1717 large_ints[1]=unpositive_int(xtc3_context->large_intra_delta[(*ilargeintra)+1])+prevcoord[1];
1718 large_ints[2]=unpositive_int(xtc3_context->large_intra_delta[(*ilargeintra)+2])+prevcoord[2];
1721 else if (xtc3_context->large_inter_delta)
1723 large_ints[0]=unpositive_int(xtc3_context->large_inter_delta[(*ilargeinter)])
1724 +output[outdata-natoms*3+didswap*3];
1725 large_ints[1]=unpositive_int(xtc3_context->large_inter_delta[(*ilargeinter)+1])
1726 +output[outdata-natoms*3+1+didswap*3];
1727 large_ints[2]=unpositive_int(xtc3_context->large_inter_delta[(*ilargeinter)+2])
1728 +output[outdata-natoms*3+2+didswap*3];
1731 memcpy(prevcoord, large_ints, 3*sizeof *prevcoord);
1732 output[outdata]=large_ints[0];
1733 output[outdata+1]=large_ints[1];
1734 output[outdata+2]=large_ints[2];
1736 fprintf(stderr,"Unpack one large: %d %d %d\n",prevcoord[0],prevcoord[1],prevcoord[2]);
1741 int Ptngc_unpack_array_xtc3(unsigned char *packed,int *output, const int length, const int natoms)
1745 unsigned char *ptr=packed;
1748 int ntriplets_left=length/3;
1751 int current_large_type=0;
1759 struct xtc3_context xtc3_context;
1760 init_xtc3_context(&xtc3_context);
1764 minint[i]=unpositive_int((int)(((unsigned int)ptr[0]) |
1765 (((unsigned int)ptr[1])<<8) |
1766 (((unsigned int)ptr[2])<<16) |
1767 (((unsigned int)ptr[3])<<24)));
1771 xtc3_context.ninstr=(int)(((unsigned int)ptr[0]) |
1772 (((unsigned int)ptr[1])<<8) |
1773 (((unsigned int)ptr[2])<<16) |
1774 (((unsigned int)ptr[3])<<24));
1776 if (xtc3_context.ninstr)
1777 decompress_bwlzh_block(&ptr,xtc3_context.ninstr,&xtc3_context.instructions);
1779 xtc3_context.nrle=(int)(((unsigned int)ptr[0]) |
1780 (((unsigned int)ptr[1])<<8) |
1781 (((unsigned int)ptr[2])<<16) |
1782 (((unsigned int)ptr[3])<<24));
1784 if (xtc3_context.nrle)
1785 decompress_bwlzh_block(&ptr,xtc3_context.nrle,&xtc3_context.rle);
1787 xtc3_context.nlargedir=(int)(((unsigned int)ptr[0]) |
1788 (((unsigned int)ptr[1])<<8) |
1789 (((unsigned int)ptr[2])<<16) |
1790 (((unsigned int)ptr[3])<<24));
1792 if (xtc3_context.nlargedir)
1795 decompress_bwlzh_block(&ptr,xtc3_context.nlargedir,&xtc3_context.large_direct);
1797 decompress_base_block(&ptr,xtc3_context.nlargedir,&xtc3_context.large_direct);
1800 xtc3_context.nlargeintra=(int)(((unsigned int)ptr[0]) |
1801 (((unsigned int)ptr[1])<<8) |
1802 (((unsigned int)ptr[2])<<16) |
1803 (((unsigned int)ptr[3])<<24));
1805 if (xtc3_context.nlargeintra)
1808 decompress_bwlzh_block(&ptr,xtc3_context.nlargeintra,&xtc3_context.large_intra_delta);
1810 decompress_base_block(&ptr,xtc3_context.nlargeintra,&xtc3_context.large_intra_delta);
1813 xtc3_context.nlargeinter=(int)(((unsigned int)ptr[0]) |
1814 (((unsigned int)ptr[1])<<8) |
1815 (((unsigned int)ptr[2])<<16) |
1816 (((unsigned int)ptr[3])<<24));
1818 if (xtc3_context.nlargeinter)
1821 decompress_bwlzh_block(&ptr,xtc3_context.nlargeinter,&xtc3_context.large_inter_delta);
1823 decompress_base_block(&ptr,xtc3_context.nlargeinter,&xtc3_context.large_inter_delta);
1826 xtc3_context.nsmallintra=(int)(((unsigned int)ptr[0]) |
1827 (((unsigned int)ptr[1])<<8) |
1828 (((unsigned int)ptr[2])<<16) |
1829 (((unsigned int)ptr[3])<<24));
1831 if (xtc3_context.nsmallintra)
1834 decompress_bwlzh_block(&ptr,xtc3_context.nsmallintra,&xtc3_context.smallintra);
1836 decompress_base_block(&ptr,xtc3_context.nsmallintra,&xtc3_context.smallintra);
1839 /* Initial prevcoord is the minimum integers. */
1840 memcpy(prevcoord, minint, 3*sizeof *prevcoord);
1842 while (ntriplets_left>0 && iinstr<xtc3_context.ninstr)
1844 int instr=xtc3_context.instructions[iinstr++];
1846 fprintf(stderr,"instr=%d @ %d\n",instr,iinstr-1);
1849 fprintf(stderr,"ntriplets left=%d\n",ntriplets_left);
1851 if ((instr==INSTR_DEFAULT) /* large+small */
1852 || (instr==INSTR_ONLY_LARGE) /* only large */
1853 || (instr==INSTR_ONLY_SMALL)) /* only small */
1855 if (instr!=INSTR_ONLY_SMALL)
1858 if ((instr==INSTR_DEFAULT) && (swapatoms))
1860 unpack_one_large(&xtc3_context,&ilargedir, &ilargeintra, &ilargeinter,
1861 prevcoord, minint, output, outdata, didswap,
1862 natoms, current_large_type);
1866 if (instr!=INSTR_ONLY_LARGE)
1868 for (i=0; i<runlength; i++)
1870 prevcoord[0]+=unpositive_int(xtc3_context.smallintra[ismallintra]);
1871 prevcoord[1]+=unpositive_int(xtc3_context.smallintra[ismallintra+1]);
1872 prevcoord[2]+=unpositive_int(xtc3_context.smallintra[ismallintra+2]);
1874 output[outdata+i*3]=prevcoord[0];
1875 output[outdata+i*3+1]=prevcoord[1];
1876 output[outdata+i*3+2]=prevcoord[2];
1878 fprintf(stderr,"Unpack small: %d %d %d\n",prevcoord[0],prevcoord[1],prevcoord[2]);
1881 if ((instr==INSTR_DEFAULT) && (swapatoms))
1885 int tmp=output[outdata-3+i];
1886 output[outdata-3+i]=output[outdata+i];
1887 output[outdata+i]=tmp;
1890 fprintf(stderr,"Unswap results in\n");
1892 fprintf(stderr,"%d %d %d\n",output[outdata-3+i*3],output[outdata-2+i*3],output[outdata-1+i*3]);
1895 ntriplets_left-=runlength;
1896 outdata+=runlength*3;
1899 else if (instr==INSTR_LARGE_RLE && irle<xtc3_context.nrle)
1901 int large_rle=xtc3_context.rle[irle++];
1903 fprintf(stderr,"large_rle=%d @ %d\n",large_rle,irle-1);
1905 for (i=0; i<large_rle; i++)
1907 unpack_one_large(&xtc3_context,&ilargedir, &ilargeintra, &ilargeinter,
1908 prevcoord, minint, output, outdata, 0,
1909 natoms, current_large_type);
1914 else if (instr==INSTR_SMALL_RUNLENGTH && irle<xtc3_context.nrle)
1916 runlength=xtc3_context.rle[irle++];
1918 fprintf(stderr,"small_rle=%d @ %d\n",runlength,irle-1);
1921 else if (instr==INSTR_FLIP)
1923 swapatoms=1-swapatoms;
1925 fprintf(stderr,"new flip=%d\n",swapatoms);
1928 else if (instr==INSTR_LARGE_DIRECT)
1930 current_large_type=0;
1932 fprintf(stderr,"large direct\n");
1935 else if (instr==INSTR_LARGE_INTRA_DELTA)
1937 current_large_type=1;
1939 fprintf(stderr,"large intra delta\n");
1942 else if (instr==INSTR_LARGE_INTER_DELTA)
1944 current_large_type=2;
1946 fprintf(stderr,"large inter delta\n");
1950 if (ntriplets_left<0)
1952 fprintf(stderr,"TRAJNG XTC3: A bug has been found. At end ntriplets_left<0\n");
1955 free_xtc3_context(&xtc3_context);