1 /* This code is part of the tng compression routines.
3 * Written by Daniel Spangberg
4 * Copyright (c) 2010, 2013, 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.
14 #include "../../include/compression/tng_compress.h"
15 #include "../../include/compression/coder.h"
16 #include "../../include/compression/fixpoint.h"
18 /* Please see tng_compress.h for info on how to call these routines. */
20 /* This becomes TNGP for positions (little endian) and TNGV for velocities. In ASCII. */
21 #define MAGIC_INT_POS 0x50474E54
22 #define MAGIC_INT_VEL 0x56474E54
24 #define SPEED_DEFAULT 2 /* Default to relatively fast compression. For very good compression it makes sense to
25 choose speed=4 or speed=5 */
27 #define PRECISION(hi,lo) (Ptngc_i32x2_to_d(hi,lo))
29 #define MAX_FVAL 2147483647.
31 static int verify_input_data(double *x, int natoms, int nframes, double precision)
34 for (iframe=0; iframe<nframes; iframe++)
35 for (i=0; i<natoms; i++)
37 if (fabs(x[iframe*natoms*3+i*3+j]/precision+0.5)>=MAX_FVAL)
42 for (iframe=0; iframe<nframes; iframe++)
43 for (i=0; i<natoms; i++)
45 if (fabs(x[iframe*natoms*3+i*3+j]/precision+0.5)>=MAX_FVAL)
46 printf("ERROR. Too large value: %d %d %d: %g %g %g\n",iframe,i,j,x[iframe*natoms*3+i*3+j],precision,x[iframe*natoms*3+i*3+j]/precision/MAX_FVAL);
51 static int verify_input_data_float(float *x, int natoms, int nframes, float precision)
54 for (iframe=0; iframe<nframes; iframe++)
55 for (i=0; i<natoms; i++)
57 if (fabs(x[iframe*natoms*3+i*3+j]/precision+0.5)>=MAX_FVAL)
62 for (iframe=0; iframe<nframes; iframe++)
63 for (i=0; i<natoms; i++)
65 if (fabs(x[iframe*natoms*3+i*3+j]/precision+0.5)>=MAX_FVAL)
66 printf("ERROR. Too large value: %d %d %d: %g %g %g\n",iframe,i,j,x[iframe*natoms*3+i*3+j],precision,x[iframe*natoms*3+i*3+j]/precision/MAX_FVAL);
71 static int quantize(double *x, int natoms, int nframes,
76 for (iframe=0; iframe<nframes; iframe++)
77 for (i=0; i<natoms; i++)
79 quant[iframe*natoms*3+i*3+j]=(int)floor((x[iframe*natoms*3+i*3+j]/precision)+0.5);
80 return verify_input_data(x,natoms,nframes,precision);
83 static int quantize_float(float *x, int natoms, int nframes,
88 for (iframe=0; iframe<nframes; iframe++)
89 for (i=0; i<natoms; i++)
91 quant[iframe*natoms*3+i*3+j]=(int)floor((x[iframe*natoms*3+i*3+j]/precision)+0.5);
92 return verify_input_data_float(x,natoms,nframes,precision);
95 static void quant_inter_differences(int *quant, int natoms, int nframes,
99 /* The first frame is used for absolute positions. */
100 for (i=0; i<natoms; i++)
102 quant_inter[i*3+j]=quant[i*3+j];
103 /* For all other frames, the difference to the previous frame is used. */
104 for (iframe=1; iframe<nframes; iframe++)
105 for (i=0; i<natoms; i++)
107 quant_inter[iframe*natoms*3+i*3+j]=quant[iframe*natoms*3+i*3+j]-quant[(iframe-1)*natoms*3+i*3+j];
110 static void quant_intra_differences(int *quant, int natoms, int nframes,
114 for (iframe=0; iframe<nframes; iframe++)
116 /* The first atom is used with its absolute position. */
118 quant_intra[iframe*natoms*3+j]=quant[iframe*natoms*3+j];
119 /* For all other atoms the intraframe differences are computed. */
120 for (i=1; i<natoms; i++)
122 quant_intra[iframe*natoms*3+i*3+j]=quant[iframe*natoms*3+i*3+j]-quant[iframe*natoms*3+(i-1)*3+j];
126 static void unquantize(double *x, int natoms, int nframes,
131 for (iframe=0; iframe<nframes; iframe++)
132 for (i=0; i<natoms; i++)
134 x[iframe*natoms*3+i*3+j]=(double)quant[iframe*natoms*3+i*3+j]*precision;
137 static void unquantize_float(float *x, int natoms, int nframes,
142 for (iframe=0; iframe<nframes; iframe++)
143 for (i=0; i<natoms; i++)
145 x[iframe*natoms*3+i*3+j]=(float)quant[iframe*natoms*3+i*3+j]*precision;
148 static void unquantize_inter_differences(double *x, int natoms, int nframes,
153 for (i=0; i<natoms; i++)
156 int q=quant[i*3+j]; /* First value. */
157 x[i*3+j]=(double)q*precision;
158 for (iframe=1; iframe<nframes; iframe++)
160 q+=quant[iframe*natoms*3+i*3+j];
161 x[iframe*natoms*3+i*3+j]=(double)q*precision;
166 static void unquantize_inter_differences_float(float *x, int natoms, int nframes,
171 for (i=0; i<natoms; i++)
174 int q=quant[i*3+j]; /* First value. */
175 x[i*3+j]=(float)q*precision;
176 for (iframe=1; iframe<nframes; iframe++)
178 q+=quant[iframe*natoms*3+i*3+j];
179 x[iframe*natoms*3+i*3+j]=(float)q*precision;
184 static void unquantize_inter_differences_int(int *x, int natoms, int nframes,
188 for (i=0; i<natoms; i++)
191 int q=quant[i*3+j]; /* First value. */
193 for (iframe=1; iframe<nframes; iframe++)
195 q+=quant[iframe*natoms*3+i*3+j];
196 x[iframe*natoms*3+i*3+j]=q;
201 /* In frame update required for the initial frame if intra-frame
202 compression was used. */
203 static void unquant_intra_differences_first_frame(int *quant, int natoms)
209 for (i=1; i<natoms; i++)
217 for (i=0; i<natoms; i++)
219 printf("UQ: %d %d %d: %d\n",0,j,i,quant[i*3+j]);
224 static void unquantize_intra_differences(double *x, int natoms, int nframes,
230 printf("UQ precision=%g\n",precision);
232 for (iframe=0; iframe<nframes; iframe++)
235 int q=quant[iframe*natoms*3+j];
236 x[iframe*natoms*3+j]=(double)q*precision;
237 for (i=1; i<natoms; i++)
239 q+=quant[iframe*natoms*3+i*3+j];
240 x[iframe*natoms*3+i*3+j]=(double)q*precision;
245 static void unquantize_intra_differences_float(float *x, int natoms, int nframes,
250 for (iframe=0; iframe<nframes; iframe++)
253 int q=quant[iframe*natoms*3+j];
254 x[iframe*natoms*3+j]=(float)q*precision;
255 for (i=1; i<natoms; i++)
257 q+=quant[iframe*natoms*3+i*3+j];
258 x[iframe*natoms*3+i*3+j]=(float)q*precision;
263 static void unquantize_intra_differences_int(int *x, int natoms, int nframes,
267 for (iframe=0; iframe<nframes; iframe++)
270 int q=quant[iframe*natoms*3+j];
271 x[iframe*natoms*3+j]=q;
272 for (i=1; i<natoms; i++)
274 q+=quant[iframe*natoms*3+i*3+j];
275 x[iframe*natoms*3+i*3+j]=q;
280 /* Buffer num 8 bit bytes into buffer location buf */
281 static void bufferfix(unsigned char *buf, fix_t v, int num)
283 /* Store in little endian format. */
284 unsigned char c; /* at least 8 bits. */
285 c=(unsigned char)(v & 0xFFU);
290 c=(unsigned char)(v & 0xFFU);
294 static fix_t readbufferfix(unsigned char *buf, int num)
303 f |= ((fix_t)b & 0xFF)<<shift;
309 /* Perform position compression from the quantized data. */
310 static void compress_quantized_pos(int *quant, int *quant_inter, int *quant_intra,
311 int natoms, int nframes,
313 int initial_coding, int initial_coding_parameter,
314 int coding, int coding_parameter,
315 fix_t prec_hi, fix_t prec_lo,
320 char *datablock=NULL;
322 /* Information needed for decompression. */
324 bufferfix((unsigned char*)data+bufloc,(fix_t)MAGIC_INT_POS,4);
326 /* Number of atoms. */
328 bufferfix((unsigned char*)data+bufloc,(fix_t)natoms,4);
330 /* Number of frames. */
332 bufferfix((unsigned char*)data+bufloc,(fix_t)nframes,4);
334 /* Initial coding. */
336 bufferfix((unsigned char*)data+bufloc,(fix_t)initial_coding,4);
338 /* Initial coding parameter. */
340 bufferfix((unsigned char*)data+bufloc,(fix_t)initial_coding_parameter,4);
344 bufferfix((unsigned char*)data+bufloc,(fix_t)coding,4);
346 /* Coding parameter. */
348 bufferfix((unsigned char*)data+bufloc,(fix_t)coding_parameter,4);
352 bufferfix((unsigned char*)data+bufloc,prec_lo,4);
355 bufferfix((unsigned char*)data+bufloc,prec_hi,4);
357 /* The initial frame */
358 if ((initial_coding==TNG_COMPRESS_ALGO_POS_XTC2) ||
359 (initial_coding==TNG_COMPRESS_ALGO_POS_TRIPLET_ONETOONE) ||
360 (initial_coding==TNG_COMPRESS_ALGO_POS_XTC3))
362 struct coder *coder=Ptngc_coder_init();
364 datablock=(char*)Ptngc_pack_array(coder,quant,&length,
365 initial_coding,initial_coding_parameter,natoms,speed);
366 Ptngc_coder_deinit(coder);
368 else if ((initial_coding==TNG_COMPRESS_ALGO_POS_TRIPLET_INTRA) ||
369 (initial_coding==TNG_COMPRESS_ALGO_POS_BWLZH_INTRA))
371 struct coder *coder=Ptngc_coder_init();
373 datablock=(char*)Ptngc_pack_array(coder,quant_intra,&length,
374 initial_coding,initial_coding_parameter,natoms,speed);
375 Ptngc_coder_deinit(coder);
379 bufferfix((unsigned char*)data+bufloc,(fix_t)length,4);
381 /* The actual data block. */
383 memcpy(data+bufloc,datablock,length);
386 /* The remaining frames */
390 /* Inter-frame compression? */
391 if ((coding==TNG_COMPRESS_ALGO_POS_STOPBIT_INTER) ||
392 (coding==TNG_COMPRESS_ALGO_POS_TRIPLET_INTER) ||
393 (coding==TNG_COMPRESS_ALGO_POS_BWLZH_INTER))
395 struct coder *coder=Ptngc_coder_init();
396 length=natoms*3*(nframes-1);
397 datablock=(char*)Ptngc_pack_array(coder,quant_inter+natoms*3,&length,
398 coding,coding_parameter,natoms,speed);
399 Ptngc_coder_deinit(coder);
401 /* One-to-one compression? */
402 else if ((coding==TNG_COMPRESS_ALGO_POS_XTC2) ||
403 (coding==TNG_COMPRESS_ALGO_POS_XTC3) ||
404 (coding==TNG_COMPRESS_ALGO_POS_TRIPLET_ONETOONE))
406 struct coder *coder=Ptngc_coder_init();
407 length=natoms*3*(nframes-1);
408 datablock=(char*)Ptngc_pack_array(coder,quant+natoms*3,&length,
409 coding,coding_parameter,natoms,speed);
410 Ptngc_coder_deinit(coder);
412 /* Intra-frame compression? */
413 else if ((coding==TNG_COMPRESS_ALGO_POS_TRIPLET_INTRA) ||
414 (coding==TNG_COMPRESS_ALGO_POS_BWLZH_INTRA))
416 struct coder *coder=Ptngc_coder_init();
417 length=natoms*3*(nframes-1);
418 datablock=(char*)Ptngc_pack_array(coder,quant_intra+natoms*3,&length,
419 coding,coding_parameter,natoms,speed);
420 Ptngc_coder_deinit(coder);
424 bufferfix((unsigned char*)data+bufloc,(fix_t)length,4);
429 memcpy(data+bufloc,datablock,length);
437 /* Perform velocity compression from vel into the data block */
438 static void compress_quantized_vel(int *quant, int *quant_inter,
439 int natoms, int nframes,
441 int initial_coding, int initial_coding_parameter,
442 int coding, int coding_parameter,
443 fix_t prec_hi, fix_t prec_lo,
448 char *datablock=NULL;
450 /* Information needed for decompression. */
452 bufferfix((unsigned char*)data+bufloc,(fix_t)MAGIC_INT_VEL,4);
454 /* Number of atoms. */
456 bufferfix((unsigned char*)data+bufloc,(fix_t)natoms,4);
458 /* Number of frames. */
460 bufferfix((unsigned char*)data+bufloc,(fix_t)nframes,4);
462 /* Initial coding. */
464 bufferfix((unsigned char*)data+bufloc,(fix_t)initial_coding,4);
466 /* Initial coding parameter. */
468 bufferfix((unsigned char*)data+bufloc,(fix_t)initial_coding_parameter,4);
472 bufferfix((unsigned char*)data+bufloc,(fix_t)coding,4);
474 /* Coding parameter. */
476 bufferfix((unsigned char*)data+bufloc,(fix_t)coding_parameter,4);
480 bufferfix((unsigned char*)data+bufloc,prec_lo,4);
483 bufferfix((unsigned char*)data+bufloc,prec_hi,4);
487 /* The initial frame */
488 if ((initial_coding==TNG_COMPRESS_ALGO_VEL_STOPBIT_ONETOONE) ||
489 (initial_coding==TNG_COMPRESS_ALGO_VEL_TRIPLET_ONETOONE) ||
490 (initial_coding==TNG_COMPRESS_ALGO_VEL_BWLZH_ONETOONE))
492 struct coder *coder=Ptngc_coder_init();
493 datablock=(char*)Ptngc_pack_array(coder,quant,&length,
494 initial_coding,initial_coding_parameter,natoms,speed);
495 Ptngc_coder_deinit(coder);
499 bufferfix((unsigned char*)data+bufloc,(fix_t)length,4);
501 /* The actual data block. */
502 if (data && datablock)
504 memcpy(data+bufloc,datablock,length);
508 /* The remaining frames */
512 /* Inter-frame compression? */
513 if ((coding==TNG_COMPRESS_ALGO_VEL_TRIPLET_INTER) ||
514 (coding==TNG_COMPRESS_ALGO_VEL_STOPBIT_INTER) ||
515 (coding==TNG_COMPRESS_ALGO_VEL_BWLZH_INTER))
517 struct coder *coder=Ptngc_coder_init();
518 length=natoms*3*(nframes-1);
519 datablock=(char*)Ptngc_pack_array(coder,quant_inter+natoms*3,&length,
520 coding,coding_parameter,natoms,speed);
521 Ptngc_coder_deinit(coder);
523 /* One-to-one compression? */
524 else if ((coding==TNG_COMPRESS_ALGO_VEL_STOPBIT_ONETOONE) ||
525 (coding==TNG_COMPRESS_ALGO_VEL_TRIPLET_ONETOONE) ||
526 (coding==TNG_COMPRESS_ALGO_VEL_BWLZH_ONETOONE))
528 struct coder *coder=Ptngc_coder_init();
529 length=natoms*3*(nframes-1);
530 datablock=(char*)Ptngc_pack_array(coder,quant+natoms*3,&length,
531 coding,coding_parameter,natoms,speed);
532 Ptngc_coder_deinit(coder);
536 bufferfix((unsigned char*)data+bufloc,(fix_t)length,4);
539 memcpy(data+bufloc,datablock,length);
546 static int determine_best_coding_stop_bits(struct coder *coder,int *input, int *length,
547 int *coding_parameter, int natoms)
550 unsigned char *packed;
552 int new_parameter=-1;
554 for (bits=1; bits<20; bits++)
557 packed=Ptngc_pack_array(coder,input,&io_length,
558 TNG_COMPRESS_ALGO_STOPBIT,bits,natoms,0);
561 if ((new_parameter==-1) || (io_length<best_length))
564 best_length=io_length;
569 if (new_parameter==-1)
572 *coding_parameter=new_parameter;
577 static int determine_best_coding_triple(struct coder *coder,int *input, int *length,
578 int *coding_parameter, int natoms)
581 unsigned char *packed;
583 int new_parameter=-1;
585 for (bits=1; bits<20; bits++)
588 packed=Ptngc_pack_array(coder,input,&io_length,
589 TNG_COMPRESS_ALGO_TRIPLET,bits,natoms,0);
592 if ((new_parameter==-1) || (io_length<best_length))
595 best_length=io_length;
600 if (new_parameter==-1)
603 *coding_parameter=new_parameter;
608 static void determine_best_pos_initial_coding(int *quant, int *quant_intra, int natoms, int speed,
609 fix_t prec_hi, fix_t prec_lo,
610 int *initial_coding, int *initial_coding_parameter)
612 if (*initial_coding==-1)
614 /* Determine all parameters automatically */
616 int best_coding_parameter;
619 int current_coding_parameter;
620 int current_code_size;
622 /* Start with XTC2, it should always work. */
623 current_coding=TNG_COMPRESS_ALGO_POS_XTC2;
624 current_coding_parameter=0;
625 compress_quantized_pos(quant,NULL,quant_intra,natoms,1,speed,
626 current_coding,current_coding_parameter,
627 0,0,prec_hi,prec_lo,¤t_code_size,NULL);
628 best_coding=current_coding;
629 best_coding_parameter=current_coding_parameter;
630 best_code_size=current_code_size;
632 /* Determine best parameter for triplet intra. */
633 current_coding=TNG_COMPRESS_ALGO_POS_TRIPLET_INTRA;
634 coder=Ptngc_coder_init();
635 current_code_size=natoms*3;
636 current_coding_parameter=0;
637 if (!determine_best_coding_triple(coder,quant_intra,¤t_code_size,¤t_coding_parameter,natoms))
639 if (current_code_size<best_code_size)
641 best_coding=current_coding;
642 best_coding_parameter=current_coding_parameter;
643 best_code_size=current_code_size;
646 Ptngc_coder_deinit(coder);
648 /* Determine best parameter for triplet one-to-one. */
649 current_coding=TNG_COMPRESS_ALGO_POS_TRIPLET_ONETOONE;
650 coder=Ptngc_coder_init();
651 current_code_size=natoms*3;
652 current_coding_parameter=0;
653 if (!determine_best_coding_triple(coder,quant,¤t_code_size,¤t_coding_parameter,natoms))
655 if (current_code_size<best_code_size)
657 best_coding=current_coding;
658 best_coding_parameter=current_coding_parameter;
659 best_code_size=current_code_size;
662 Ptngc_coder_deinit(coder);
666 current_coding=TNG_COMPRESS_ALGO_POS_XTC3;
667 current_coding_parameter=0;
668 compress_quantized_pos(quant,NULL,quant_intra,natoms,1,speed,
669 current_coding,current_coding_parameter,
670 0,0,prec_hi,prec_lo,¤t_code_size,NULL);
671 if (current_code_size<best_code_size)
673 best_coding=current_coding;
674 best_coding_parameter=current_coding_parameter;
675 best_code_size=current_code_size;
678 /* Test BWLZH intra */
681 current_coding=TNG_COMPRESS_ALGO_POS_BWLZH_INTRA;
682 current_coding_parameter=0;
683 compress_quantized_pos(quant,NULL,quant_intra,natoms,1,speed,
684 current_coding,current_coding_parameter,
685 0,0,prec_hi,prec_lo,¤t_code_size,NULL);
686 if (current_code_size<best_code_size)
688 best_coding=current_coding;
689 best_coding_parameter=current_coding_parameter;
692 *initial_coding=best_coding;
693 *initial_coding_parameter=best_coding_parameter;
697 if (*initial_coding_parameter==-1)
699 if ((*initial_coding==TNG_COMPRESS_ALGO_POS_XTC2) ||
700 (*initial_coding==TNG_COMPRESS_ALGO_POS_XTC3) ||
701 (*initial_coding==TNG_COMPRESS_ALGO_POS_BWLZH_INTRA))
702 *initial_coding_parameter=0;
703 else if (*initial_coding==TNG_COMPRESS_ALGO_POS_TRIPLET_INTRA)
705 struct coder *coder=Ptngc_coder_init();
706 int current_code_size=natoms*3;
707 determine_best_coding_triple(coder,quant_intra,¤t_code_size,initial_coding_parameter,natoms);
708 Ptngc_coder_deinit(coder);
710 else if (*initial_coding==TNG_COMPRESS_ALGO_POS_TRIPLET_ONETOONE)
712 struct coder *coder=Ptngc_coder_init();
713 int current_code_size=natoms*3;
714 determine_best_coding_triple(coder,quant,¤t_code_size,initial_coding_parameter,natoms);
715 Ptngc_coder_deinit(coder);
721 static void determine_best_pos_coding(int *quant, int *quant_inter, int *quant_intra, int natoms, int nframes, int speed,
722 fix_t prec_hi, fix_t prec_lo,
723 int *coding, int *coding_parameter)
727 /* Determine all parameters automatically */
729 int best_coding_parameter;
732 int current_coding_parameter;
733 int current_code_size;
734 int initial_code_size;
736 /* Always use XTC2 for the initial coding. */
737 compress_quantized_pos(quant,quant_inter,quant_intra,natoms,1,speed,
738 TNG_COMPRESS_ALGO_POS_XTC2,0,
740 prec_hi,prec_lo,&initial_code_size,NULL);
741 /* Start with XTC2, it should always work. */
742 current_coding=TNG_COMPRESS_ALGO_POS_XTC2;
743 current_coding_parameter=0;
744 compress_quantized_pos(quant,quant_inter,quant_intra,natoms,nframes,speed,
745 TNG_COMPRESS_ALGO_POS_XTC2,0,
746 current_coding,current_coding_parameter,
747 prec_hi,prec_lo,¤t_code_size,NULL);
748 best_coding=current_coding;
749 best_coding_parameter=current_coding_parameter;
750 best_code_size=current_code_size-initial_code_size; /* Correct for the use of XTC2 for the first frame. */
752 /* Determine best parameter for stopbit interframe coding. */
753 current_coding=TNG_COMPRESS_ALGO_POS_STOPBIT_INTER;
754 coder=Ptngc_coder_init();
755 current_code_size=natoms*3*(nframes-1);
756 current_coding_parameter=0;
757 if (!determine_best_coding_stop_bits(coder,quant_inter+natoms*3,¤t_code_size,
758 ¤t_coding_parameter,natoms))
760 if (current_code_size<best_code_size)
762 best_coding=current_coding;
763 best_coding_parameter=current_coding_parameter;
764 best_code_size=current_code_size;
767 Ptngc_coder_deinit(coder);
769 /* Determine best parameter for triplet interframe coding. */
770 current_coding=TNG_COMPRESS_ALGO_POS_TRIPLET_INTER;
771 coder=Ptngc_coder_init();
772 current_code_size=natoms*3*(nframes-1);
773 current_coding_parameter=0;
774 if (!determine_best_coding_triple(coder,quant_inter+natoms*3,¤t_code_size,
775 ¤t_coding_parameter,natoms))
777 if (current_code_size<best_code_size)
779 best_coding=current_coding;
780 best_coding_parameter=current_coding_parameter;
781 best_code_size=current_code_size;
784 Ptngc_coder_deinit(coder);
786 /* Determine best parameter for triplet intraframe coding. */
787 current_coding=TNG_COMPRESS_ALGO_POS_TRIPLET_INTRA;
788 coder=Ptngc_coder_init();
789 current_code_size=natoms*3*(nframes-1);
790 current_coding_parameter=0;
791 if (!determine_best_coding_triple(coder,quant_intra+natoms*3,¤t_code_size,
792 ¤t_coding_parameter,natoms))
794 if (current_code_size<best_code_size)
796 best_coding=current_coding;
797 best_coding_parameter=current_coding_parameter;
798 best_code_size=current_code_size;
801 Ptngc_coder_deinit(coder);
803 /* Determine best parameter for triplet one-to-one coding. */
804 current_coding=TNG_COMPRESS_ALGO_POS_TRIPLET_ONETOONE;
805 coder=Ptngc_coder_init();
806 current_code_size=natoms*3*(nframes-1);
807 current_coding_parameter=0;
808 if (!determine_best_coding_triple(coder,quant+natoms*3,¤t_code_size,
809 ¤t_coding_parameter,natoms))
811 if (current_code_size<best_code_size)
813 best_coding=current_coding;
814 best_coding_parameter=current_coding_parameter;
815 best_code_size=current_code_size;
818 Ptngc_coder_deinit(coder);
820 /* Test BWLZH inter */
823 current_coding=TNG_COMPRESS_ALGO_POS_BWLZH_INTER;
824 current_coding_parameter=0;
825 compress_quantized_pos(quant,quant_inter,quant_intra,natoms,nframes,speed,
826 TNG_COMPRESS_ALGO_POS_XTC2,0,
827 current_coding,current_coding_parameter,
828 prec_hi,prec_lo,¤t_code_size,NULL);
829 current_code_size-=initial_code_size; /* Correct for the use of XTC2 for the first frame. */
830 if (current_code_size<best_code_size)
832 best_coding=current_coding;
833 best_coding_parameter=current_coding_parameter;
834 best_code_size=current_code_size;
838 /* Test BWLZH intra */
841 current_coding=TNG_COMPRESS_ALGO_POS_BWLZH_INTRA;
842 current_coding_parameter=0;
843 compress_quantized_pos(quant,quant_inter,quant_intra,natoms,nframes,speed,
844 TNG_COMPRESS_ALGO_POS_XTC2,0,
845 current_coding,current_coding_parameter,
846 prec_hi,prec_lo,¤t_code_size,NULL);
847 current_code_size-=initial_code_size; /* Correct for the use of XTC2 for the first frame. */
848 if (current_code_size<best_code_size)
850 best_coding=current_coding;
851 best_coding_parameter=current_coding_parameter;
855 *coding_parameter=best_coding_parameter;
857 else if (*coding_parameter==-1)
859 if ((*coding==TNG_COMPRESS_ALGO_POS_XTC2) ||
860 (*coding==TNG_COMPRESS_ALGO_POS_XTC3) ||
861 (*coding==TNG_COMPRESS_ALGO_POS_BWLZH_INTER) ||
862 (*coding==TNG_COMPRESS_ALGO_POS_BWLZH_INTRA))
864 else if (*coding==TNG_COMPRESS_ALGO_POS_STOPBIT_INTER)
866 struct coder *coder=Ptngc_coder_init();
867 int current_code_size=natoms*3*(nframes-1);
868 determine_best_coding_stop_bits(coder,quant_inter+natoms*3,¤t_code_size,
869 coding_parameter,natoms);
870 Ptngc_coder_deinit(coder);
872 else if (*coding==TNG_COMPRESS_ALGO_POS_TRIPLET_INTER)
874 struct coder *coder=Ptngc_coder_init();
875 int current_code_size=natoms*3*(nframes-1);
876 determine_best_coding_triple(coder,quant_inter+natoms*3,¤t_code_size,
877 coding_parameter,natoms);
878 Ptngc_coder_deinit(coder);
880 else if (*coding==TNG_COMPRESS_ALGO_POS_TRIPLET_INTRA)
882 struct coder *coder=Ptngc_coder_init();
883 int current_code_size=natoms*3*(nframes-1);
884 determine_best_coding_triple(coder,quant_intra+natoms*3,¤t_code_size,
885 coding_parameter,natoms);
886 Ptngc_coder_deinit(coder);
888 else if (*coding==TNG_COMPRESS_ALGO_POS_TRIPLET_ONETOONE)
890 struct coder *coder=Ptngc_coder_init();
891 int current_code_size=natoms*3*(nframes-1);
892 determine_best_coding_triple(coder,quant+natoms*3,¤t_code_size,
893 coding_parameter,natoms);
894 Ptngc_coder_deinit(coder);
899 static void determine_best_vel_initial_coding(int *quant, int natoms, int speed,
900 fix_t prec_hi, fix_t prec_lo,
901 int *initial_coding, int *initial_coding_parameter)
903 if (*initial_coding==-1)
905 /* Determine all parameters automatically */
907 int best_coding_parameter=-1;
908 int best_code_size=-1;
910 int current_coding_parameter;
911 int current_code_size;
913 /* Start to determine best parameter for stopbit one-to-one. */
914 current_coding=TNG_COMPRESS_ALGO_VEL_STOPBIT_ONETOONE;
915 current_code_size=natoms*3;
916 current_coding_parameter=0;
917 coder=Ptngc_coder_init();
918 if (!determine_best_coding_stop_bits(coder,quant,¤t_code_size,
919 ¤t_coding_parameter,natoms))
921 best_coding=current_coding;
922 best_coding_parameter=current_coding_parameter;
923 best_code_size=current_code_size;
925 Ptngc_coder_deinit(coder);
927 /* Determine best parameter for triplet one-to-one. */
928 current_coding=TNG_COMPRESS_ALGO_VEL_TRIPLET_ONETOONE;
929 coder=Ptngc_coder_init();
930 current_code_size=natoms*3;
931 current_coding_parameter=0;
932 if (!determine_best_coding_triple(coder,quant,¤t_code_size,¤t_coding_parameter,natoms))
934 if ((best_coding==-1) || (current_code_size<best_code_size))
936 best_coding=current_coding;
937 best_coding_parameter=current_coding_parameter;
938 best_code_size=current_code_size;
941 Ptngc_coder_deinit(coder);
943 /* Test BWLZH one-to-one */
946 current_coding=TNG_COMPRESS_ALGO_VEL_BWLZH_ONETOONE;
947 current_coding_parameter=0;
948 compress_quantized_vel(quant,NULL,natoms,1,speed,
949 current_coding,current_coding_parameter,
950 0,0,prec_hi,prec_lo,¤t_code_size,NULL);
951 if ((best_coding==-1) || (current_code_size<best_code_size))
953 best_coding=current_coding;
954 best_coding_parameter=current_coding_parameter;
957 *initial_coding=best_coding;
958 *initial_coding_parameter=best_coding_parameter;
960 else if (*initial_coding_parameter==-1)
962 if (*initial_coding==TNG_COMPRESS_ALGO_VEL_BWLZH_ONETOONE)
963 *initial_coding_parameter=0;
964 else if (*initial_coding==TNG_COMPRESS_ALGO_VEL_STOPBIT_ONETOONE)
966 struct coder *coder=Ptngc_coder_init();
967 int current_code_size=natoms*3;
968 determine_best_coding_stop_bits(coder,quant,¤t_code_size,
969 initial_coding_parameter,natoms);
970 Ptngc_coder_deinit(coder);
972 else if (*initial_coding==TNG_COMPRESS_ALGO_VEL_TRIPLET_ONETOONE)
974 struct coder *coder=Ptngc_coder_init();
975 int current_code_size=natoms*3;
976 determine_best_coding_triple(coder,quant,¤t_code_size,initial_coding_parameter,natoms);
977 Ptngc_coder_deinit(coder);
982 static void determine_best_vel_coding(int *quant, int *quant_inter, int natoms, int nframes, int speed,
983 fix_t prec_hi, fix_t prec_lo,
984 int *coding, int *coding_parameter)
988 /* Determine all parameters automatically */
990 int best_coding_parameter;
993 int current_coding_parameter;
994 int current_code_size;
995 int initial_code_size;
996 int initial_numbits=5;
998 /* Use stopbits one-to-one coding for the initial coding. */
999 compress_quantized_vel(quant,NULL,natoms,1,speed,
1000 TNG_COMPRESS_ALGO_VEL_STOPBIT_ONETOONE,initial_numbits,
1001 0,0,prec_hi,prec_lo,&initial_code_size,NULL);
1003 /* Test stopbit one-to-one */
1004 current_coding=TNG_COMPRESS_ALGO_VEL_STOPBIT_ONETOONE;
1005 current_code_size=natoms*3*(nframes-1);
1006 current_coding_parameter=0;
1007 coder=Ptngc_coder_init();
1008 determine_best_coding_stop_bits(coder,quant+natoms*3,¤t_code_size,
1009 ¤t_coding_parameter,natoms);
1010 Ptngc_coder_deinit(coder);
1011 best_coding=current_coding;
1012 best_code_size=current_code_size;
1013 best_coding_parameter=current_coding_parameter;
1015 /* Test triplet interframe */
1016 current_coding=TNG_COMPRESS_ALGO_VEL_TRIPLET_INTER;
1017 current_code_size=natoms*3*(nframes-1);
1018 current_coding_parameter=0;
1019 coder=Ptngc_coder_init();
1020 if (!determine_best_coding_triple(coder,quant_inter+natoms*3,¤t_code_size,
1021 ¤t_coding_parameter,natoms))
1023 if (current_code_size<best_code_size)
1025 best_coding=current_coding;
1026 best_code_size=current_code_size;
1027 best_coding_parameter=current_coding_parameter;
1030 Ptngc_coder_deinit(coder);
1032 /* Test triplet one-to-one */
1033 current_coding=TNG_COMPRESS_ALGO_VEL_TRIPLET_ONETOONE;
1034 current_code_size=natoms*3*(nframes-1);
1035 current_coding_parameter=0;
1036 coder=Ptngc_coder_init();
1037 if (!determine_best_coding_triple(coder,quant+natoms*3,¤t_code_size,
1038 ¤t_coding_parameter,natoms))
1040 if (current_code_size<best_code_size)
1042 best_coding=current_coding;
1043 best_code_size=current_code_size;
1044 best_coding_parameter=current_coding_parameter;
1047 Ptngc_coder_deinit(coder);
1049 /* Test stopbit interframe */
1050 current_coding=TNG_COMPRESS_ALGO_VEL_STOPBIT_INTER;
1051 current_code_size=natoms*3*(nframes-1);
1052 current_coding_parameter=0;
1053 coder=Ptngc_coder_init();
1054 if (!determine_best_coding_stop_bits(coder,quant_inter+natoms*3,¤t_code_size,
1055 ¤t_coding_parameter,natoms))
1057 if (current_code_size<best_code_size)
1059 best_coding=current_coding;
1060 best_code_size=current_code_size;
1061 best_coding_parameter=current_coding_parameter;
1064 Ptngc_coder_deinit(coder);
1068 /* Test BWLZH inter */
1069 current_coding=TNG_COMPRESS_ALGO_VEL_BWLZH_INTER;
1070 current_coding_parameter=0;
1071 compress_quantized_vel(quant,quant_inter,natoms,nframes,speed,
1072 TNG_COMPRESS_ALGO_VEL_STOPBIT_ONETOONE,initial_numbits,
1073 current_coding,current_coding_parameter,
1074 prec_hi,prec_lo,¤t_code_size,NULL);
1075 current_code_size-=initial_code_size; /* Correct for the initial frame */
1076 if (current_code_size<best_code_size)
1078 best_coding=current_coding;
1079 best_code_size=current_code_size;
1080 best_coding_parameter=current_coding_parameter;
1083 /* Test BWLZH one-to-one */
1084 current_coding=TNG_COMPRESS_ALGO_VEL_BWLZH_ONETOONE;
1085 current_coding_parameter=0;
1086 compress_quantized_vel(quant,quant_inter,natoms,nframes,speed,
1087 TNG_COMPRESS_ALGO_VEL_STOPBIT_ONETOONE,initial_numbits,
1088 current_coding,current_coding_parameter,
1089 prec_hi,prec_lo,¤t_code_size,NULL);
1090 current_code_size-=initial_code_size; /* Correct for the initial frame */
1091 if (current_code_size<best_code_size)
1093 best_coding=current_coding;
1094 best_coding_parameter=current_coding_parameter;
1097 *coding=best_coding;
1098 *coding_parameter=best_coding_parameter;
1100 else if (*coding_parameter==-1)
1102 if ((*coding==TNG_COMPRESS_ALGO_VEL_BWLZH_INTER) ||
1103 (*coding==TNG_COMPRESS_ALGO_VEL_BWLZH_ONETOONE))
1104 *coding_parameter=0;
1105 else if (*coding==TNG_COMPRESS_ALGO_VEL_STOPBIT_ONETOONE)
1107 struct coder *coder=Ptngc_coder_init();
1108 int current_code_size=natoms*3*(nframes-1);
1109 determine_best_coding_stop_bits(coder,quant+natoms*3,¤t_code_size,
1110 coding_parameter,natoms);
1111 Ptngc_coder_deinit(coder);
1113 else if (*coding==TNG_COMPRESS_ALGO_VEL_TRIPLET_INTER)
1115 struct coder *coder=Ptngc_coder_init();
1116 int current_code_size=natoms*3*(nframes-1);
1117 determine_best_coding_triple(coder,quant_inter+natoms*3,¤t_code_size,
1118 coding_parameter,natoms);
1119 Ptngc_coder_deinit(coder);
1121 else if (*coding==TNG_COMPRESS_ALGO_VEL_TRIPLET_ONETOONE)
1123 struct coder *coder=Ptngc_coder_init();
1124 int current_code_size=natoms*3*(nframes-1);
1125 determine_best_coding_triple(coder,quant+natoms*3,¤t_code_size,
1126 coding_parameter,natoms);
1127 Ptngc_coder_deinit(coder);
1129 else if (*coding==TNG_COMPRESS_ALGO_VEL_STOPBIT_INTER)
1131 struct coder *coder=Ptngc_coder_init();
1132 int current_code_size=natoms*3*(nframes-1);
1133 determine_best_coding_stop_bits(coder,quant_inter+natoms*3,¤t_code_size,
1134 coding_parameter,natoms);
1135 Ptngc_coder_deinit(coder);
1140 char DECLSPECDLLEXPORT *tng_compress_pos_int(int *pos, int natoms, int nframes,
1141 unsigned long prec_hi, unsigned long prec_lo,
1142 int speed,int *algo,
1145 char *data=malloc(natoms*nframes*14+11*4); /* 12 bytes are required to store 4 32 bit integers
1146 This is 17% extra. The final 11*4 is to store information
1147 needed for decompression. */
1148 int *quant=pos; /* Already quantized positions. */
1149 int *quant_intra=malloc(natoms*nframes*3*sizeof *quant_intra);
1150 int *quant_inter=malloc(natoms*nframes*3*sizeof *quant_inter);
1152 int initial_coding, initial_coding_parameter;
1153 int coding, coding_parameter;
1155 speed=SPEED_DEFAULT; /* Set the default speed. */
1156 /* Boundaries of speed. */
1161 initial_coding=algo[0];
1162 initial_coding_parameter=algo[1];
1164 coding_parameter=algo[3];
1166 quant_inter_differences(quant,natoms,nframes,quant_inter);
1167 quant_intra_differences(quant,natoms,nframes,quant_intra);
1169 /* If any of the above codings / coding parameters are == -1, the optimal parameters must be found */
1170 if (initial_coding==-1)
1172 initial_coding_parameter=-1;
1173 determine_best_pos_initial_coding(quant,quant_intra,natoms,speed,prec_hi,prec_lo,
1174 &initial_coding,&initial_coding_parameter);
1176 else if (initial_coding_parameter==-1)
1178 determine_best_pos_initial_coding(quant,quant_intra,natoms,speed,prec_hi,prec_lo,
1179 &initial_coding,&initial_coding_parameter);
1192 coding_parameter=-1;
1193 determine_best_pos_coding(quant,quant_inter,quant_intra,natoms,nframes,speed,prec_hi,prec_lo,
1194 &coding,&coding_parameter);
1196 else if (coding_parameter==-1)
1198 determine_best_pos_coding(quant,quant_inter,quant_intra,natoms,nframes,speed,prec_hi,prec_lo,
1199 &coding,&coding_parameter);
1203 compress_quantized_pos(quant,quant_inter,quant_intra,natoms,nframes,speed,
1204 initial_coding,initial_coding_parameter,
1205 coding,coding_parameter,
1206 prec_hi,prec_lo,nitems,data);
1210 algo[0]=initial_coding;
1212 algo[1]=initial_coding_parameter;
1216 algo[3]=coding_parameter;
1220 char DECLSPECDLLEXPORT *tng_compress_pos(double *pos, int natoms, int nframes,
1221 double desired_precision,
1222 int speed,int *algo,
1225 int *quant=malloc(natoms*nframes*3*sizeof *quant);
1227 fix_t prec_hi, prec_lo;
1228 Ptngc_d_to_i32x2(desired_precision,&prec_hi,&prec_lo);
1230 if (quantize(pos,natoms,nframes,PRECISION(prec_hi,prec_lo),quant))
1231 data=NULL; /* Error occured. Too large input values. */
1233 data=tng_compress_pos_int(quant,natoms,nframes,prec_hi,prec_lo,speed,algo,nitems);
1238 char DECLSPECDLLEXPORT *tng_compress_pos_float(float *pos, int natoms, int nframes,
1239 float desired_precision,
1240 int speed,int *algo,
1243 int *quant=malloc(natoms*nframes*3*sizeof *quant);
1245 fix_t prec_hi, prec_lo;
1246 Ptngc_d_to_i32x2((double)desired_precision,&prec_hi,&prec_lo);
1248 if (quantize_float(pos,natoms,nframes,(float)PRECISION(prec_hi,prec_lo),quant))
1249 data=NULL; /* Error occured. Too large input values. */
1251 data=tng_compress_pos_int(quant,natoms,nframes,prec_hi,prec_lo,speed,algo,nitems);
1256 char DECLSPECDLLEXPORT *tng_compress_pos_find_algo(double *pos, int natoms, int nframes,
1257 double desired_precision,
1266 return tng_compress_pos(pos,natoms,nframes,desired_precision,speed,algo,nitems);
1269 char DECLSPECDLLEXPORT *tng_compress_pos_float_find_algo(float *pos, int natoms, int nframes,
1270 float desired_precision,
1279 return tng_compress_pos_float(pos,natoms,nframes,desired_precision,speed,algo,nitems);
1282 char DECLSPECDLLEXPORT *tng_compress_pos_int_find_algo(int *pos, int natoms, int nframes,
1283 unsigned long prec_hi, unsigned long prec_lo,
1284 int speed,int *algo,
1291 return tng_compress_pos_int(pos,natoms,nframes,prec_hi,prec_lo,speed,algo,nitems);
1296 int DECLSPECDLLEXPORT tng_compress_nalgo(void)
1298 return 4; /* There are currently four parameters required:
1300 1) The compression algorithm for the first frame (initial_coding).
1301 2) One parameter to the algorithm for the first frame (the initial coding parameter).
1302 3) The compression algorithm for the remaining frames (coding).
1303 4) One parameter to the algorithm for the remaining frames (the coding parameter). */
1306 char DECLSPECDLLEXPORT *tng_compress_vel_int(int *vel, int natoms, int nframes,
1307 unsigned long prec_hi, unsigned long prec_lo,
1308 int speed, int *algo,
1311 char *data=malloc(natoms*nframes*14+11*4); /* 12 bytes are required to store 4 32 bit integers
1312 This is 17% extra. The final 11*4 is to store information
1313 needed for decompression. */
1315 int *quant_inter=malloc(natoms*nframes*3*sizeof *quant_inter);
1317 int initial_coding, initial_coding_parameter;
1318 int coding, coding_parameter;
1320 speed=SPEED_DEFAULT; /* Set the default speed. */
1321 /* Boundaries of speed. */
1326 initial_coding=algo[0];
1327 initial_coding_parameter=algo[1];
1329 coding_parameter=algo[3];
1331 quant_inter_differences(quant,natoms,nframes,quant_inter);
1333 /* If any of the above codings / coding parameters are == -1, the optimal parameters must be found */
1334 if (initial_coding==-1)
1336 initial_coding_parameter=-1;
1337 determine_best_vel_initial_coding(quant,natoms,speed,prec_hi,prec_lo,
1338 &initial_coding,&initial_coding_parameter);
1340 else if (initial_coding_parameter==-1)
1342 determine_best_vel_initial_coding(quant,natoms,speed,prec_hi,prec_lo,
1343 &initial_coding,&initial_coding_parameter);
1356 coding_parameter=-1;
1357 determine_best_vel_coding(quant,quant_inter,natoms,nframes,speed,prec_hi,prec_lo,
1358 &coding,&coding_parameter);
1360 else if (coding_parameter==-1)
1362 determine_best_vel_coding(quant,quant_inter,natoms,nframes,speed,prec_hi,prec_lo,
1363 &coding,&coding_parameter);
1367 compress_quantized_vel(quant,quant_inter,natoms,nframes,speed,
1368 initial_coding,initial_coding_parameter,
1369 coding,coding_parameter,
1370 prec_hi,prec_lo,nitems,data);
1373 algo[0]=initial_coding;
1375 algo[1]=initial_coding_parameter;
1379 algo[3]=coding_parameter;
1383 char DECLSPECDLLEXPORT *tng_compress_vel(double *vel, int natoms, int nframes,
1384 double desired_precision,
1385 int speed, int *algo,
1388 int *quant=malloc(natoms*nframes*3*sizeof *quant);
1390 fix_t prec_hi, prec_lo;
1391 Ptngc_d_to_i32x2(desired_precision,&prec_hi,&prec_lo);
1392 if (quantize(vel,natoms,nframes,PRECISION(prec_hi,prec_lo),quant))
1393 data=NULL; /* Error occured. Too large input values. */
1395 data=tng_compress_vel_int(quant,natoms,nframes,prec_hi,prec_lo,speed,algo,nitems);
1400 char DECLSPECDLLEXPORT *tng_compress_vel_float(float *vel, int natoms, int nframes,
1401 float desired_precision,
1402 int speed, int *algo,
1405 int *quant=malloc(natoms*nframes*3*sizeof *quant);
1407 fix_t prec_hi, prec_lo;
1408 Ptngc_d_to_i32x2((double)desired_precision,&prec_hi,&prec_lo);
1409 if (quantize_float(vel,natoms,nframes,(float)PRECISION(prec_hi,prec_lo),quant))
1410 data=NULL; /* Error occured. Too large input values. */
1412 data=tng_compress_vel_int(quant,natoms,nframes,prec_hi,prec_lo,speed,algo,nitems);
1417 char DECLSPECDLLEXPORT *tng_compress_vel_find_algo(double *vel, int natoms, int nframes,
1418 double desired_precision,
1427 return tng_compress_vel(vel,natoms,nframes,desired_precision,speed,algo,nitems);
1430 char DECLSPECDLLEXPORT *tng_compress_vel_float_find_algo(float *vel, int natoms, int nframes,
1431 float desired_precision,
1440 return tng_compress_vel_float(vel,natoms,nframes,desired_precision,speed,algo,nitems);
1443 char DECLSPECDLLEXPORT *tng_compress_vel_int_find_algo(int *vel, int natoms, int nframes,
1444 unsigned long prec_hi, unsigned long prec_lo,
1453 return tng_compress_vel_int(vel,natoms,nframes,prec_hi,prec_lo,speed,algo,nitems);
1456 int DECLSPECDLLEXPORT tng_compress_inquire(char *data,int *vel, int *natoms,
1457 int *nframes, double *precision,
1461 fix_t prec_hi, prec_lo;
1462 int initial_coding, initial_coding_parameter;
1463 int coding, coding_parameter;
1465 magic_int=(int)readbufferfix((unsigned char *)data+bufloc,4);
1467 if (magic_int==MAGIC_INT_POS)
1469 else if (magic_int==MAGIC_INT_VEL)
1473 /* Number of atoms. */
1474 *natoms=(int)readbufferfix((unsigned char *)data+bufloc,4);
1476 /* Number of frames. */
1477 *nframes=(int)readbufferfix((unsigned char *)data+bufloc,4);
1479 /* Initial coding. */
1480 initial_coding=(int)readbufferfix((unsigned char *)data+bufloc,4);
1482 /* Initial coding parameter. */
1483 initial_coding_parameter=(int)readbufferfix((unsigned char *)data+bufloc,4);
1486 coding=(int)readbufferfix((unsigned char *)data+bufloc,4);
1488 /* Coding parameter. */
1489 coding_parameter=(int)readbufferfix((unsigned char *)data+bufloc,4);
1492 prec_lo=readbufferfix((unsigned char *)data+bufloc,4);
1494 prec_hi=readbufferfix((unsigned char *)data+bufloc,4);
1495 *precision=PRECISION(prec_hi, prec_lo);
1496 algo[0]=initial_coding;
1497 algo[1]=initial_coding_parameter;
1499 algo[3]=coding_parameter;
1503 static int tng_compress_uncompress_pos_gen(char *data,double *posd,float *posf,int *posi,unsigned long *prec_hi, unsigned long *prec_lo)
1507 int natoms, nframes;
1508 int initial_coding, initial_coding_parameter;
1509 int coding, coding_parameter;
1511 struct coder *coder=NULL;
1514 /* Magic integer for positions */
1515 magic_int=(int)readbufferfix((unsigned char *)data+bufloc,4);
1517 if (magic_int!=MAGIC_INT_POS)
1522 /* Number of atoms. */
1523 natoms=(int)readbufferfix((unsigned char *)data+bufloc,4);
1525 /* Number of frames. */
1526 nframes=(int)readbufferfix((unsigned char *)data+bufloc,4);
1528 /* Initial coding. */
1529 initial_coding=(int)readbufferfix((unsigned char *)data+bufloc,4);
1531 /* Initial coding parameter. */
1532 initial_coding_parameter=(int)readbufferfix((unsigned char *)data+bufloc,4);
1535 coding=(int)readbufferfix((unsigned char *)data+bufloc,4);
1537 /* Coding parameter. */
1538 coding_parameter=(int)readbufferfix((unsigned char *)data+bufloc,4);
1541 *prec_lo=readbufferfix((unsigned char *)data+bufloc,4);
1543 *prec_hi=readbufferfix((unsigned char *)data+bufloc,4);
1545 /* Allocate the memory for the quantized positions */
1546 quant=malloc(natoms*nframes*3*sizeof *quant);
1547 /* The data block length. */
1548 length=(int)readbufferfix((unsigned char *)data+bufloc,4);
1550 /* The initial frame */
1551 coder=Ptngc_coder_init();
1552 rval=Ptngc_unpack_array(coder,(unsigned char*)data+bufloc,quant,natoms*3,
1553 initial_coding,initial_coding_parameter,natoms);
1554 Ptngc_coder_deinit(coder);
1557 /* Skip past the actual data block. */
1559 /* Obtain the actual positions for the initial block. */
1560 if ((initial_coding==TNG_COMPRESS_ALGO_POS_XTC2) ||
1561 (initial_coding==TNG_COMPRESS_ALGO_POS_TRIPLET_ONETOONE) ||
1562 (initial_coding==TNG_COMPRESS_ALGO_POS_XTC3))
1565 unquantize(posd,natoms,1,PRECISION(*prec_hi,*prec_lo),quant);
1567 unquantize_float(posf,natoms,1,(float)PRECISION(*prec_hi,*prec_lo),quant);
1569 memcpy(posi,quant,natoms*3*sizeof *posi);
1571 else if ((initial_coding==TNG_COMPRESS_ALGO_POS_TRIPLET_INTRA) ||
1572 (initial_coding==TNG_COMPRESS_ALGO_POS_BWLZH_INTRA))
1575 unquantize_intra_differences(posd,natoms,1,PRECISION(*prec_hi,*prec_lo),quant);
1577 unquantize_intra_differences_float(posf,natoms,1,(float)PRECISION(*prec_hi,*prec_lo),quant);
1579 unquantize_intra_differences_int(posi,natoms,1,quant);
1580 unquant_intra_differences_first_frame(quant,natoms);
1582 /* The remaining frames. */
1586 coder=Ptngc_coder_init();
1587 rval=Ptngc_unpack_array(coder,(unsigned char *)data+bufloc,quant+natoms*3,(nframes-1)*natoms*3,
1588 coding,coding_parameter,natoms);
1589 Ptngc_coder_deinit(coder);
1592 if ((coding==TNG_COMPRESS_ALGO_POS_STOPBIT_INTER) ||
1593 (coding==TNG_COMPRESS_ALGO_POS_TRIPLET_INTER) ||
1594 (coding==TNG_COMPRESS_ALGO_POS_BWLZH_INTER))
1596 /* This requires that the first frame is already in one-to-one format, even if intra-frame
1597 compression was done there. Therefore the unquant_intra_differences_first_frame should be called
1598 before to convert it correctly. */
1600 unquantize_inter_differences(posd,natoms,nframes,PRECISION(*prec_hi,*prec_lo),quant);
1602 unquantize_inter_differences_float(posf,natoms,nframes,(float)PRECISION(*prec_hi,*prec_lo),quant);
1604 unquantize_inter_differences_int(posi,natoms,nframes,quant);
1606 else if ((coding==TNG_COMPRESS_ALGO_POS_XTC2) ||
1607 (coding==TNG_COMPRESS_ALGO_POS_XTC3) ||
1608 (coding==TNG_COMPRESS_ALGO_POS_TRIPLET_ONETOONE))
1611 unquantize(posd+natoms*3,natoms,nframes-1,PRECISION(*prec_hi,*prec_lo),quant+natoms*3);
1613 unquantize_float(posf+natoms*3,natoms,nframes-1,(float)PRECISION(*prec_hi,*prec_lo),quant+natoms*3);
1615 memcpy(posi+natoms*3,quant+natoms*3,natoms*3*(nframes-1)*sizeof *posi);
1617 else if ((coding==TNG_COMPRESS_ALGO_POS_TRIPLET_INTRA) ||
1618 (coding==TNG_COMPRESS_ALGO_POS_BWLZH_INTRA))
1621 unquantize_intra_differences(posd+natoms*3,natoms,nframes-1,PRECISION(*prec_hi,*prec_lo),quant+natoms*3);
1623 unquantize_intra_differences_float(posf+natoms*3,natoms,nframes-1,(float)PRECISION(*prec_hi,*prec_lo),quant+natoms*3);
1625 unquantize_intra_differences_int(posi+natoms*3,natoms,nframes-1,quant+natoms*3);
1633 static int tng_compress_uncompress_pos(char *data,double *pos)
1635 unsigned long prec_hi, prec_lo;
1636 return tng_compress_uncompress_pos_gen(data,pos,NULL,NULL,&prec_hi,&prec_lo);
1639 static int tng_compress_uncompress_pos_float(char *data,float *pos)
1641 unsigned long prec_hi, prec_lo;
1642 return tng_compress_uncompress_pos_gen(data,NULL,pos,NULL,&prec_hi,&prec_lo);
1645 static int tng_compress_uncompress_pos_int(char *data,int *pos, unsigned long *prec_hi, unsigned long *prec_lo)
1647 return tng_compress_uncompress_pos_gen(data,NULL,NULL,pos,prec_hi,prec_lo);
1650 static int tng_compress_uncompress_vel_gen(char *data,double *veld,float *velf,int *veli,unsigned long *prec_hi, unsigned long *prec_lo)
1654 int natoms, nframes;
1655 int initial_coding, initial_coding_parameter;
1656 int coding, coding_parameter;
1658 struct coder *coder=NULL;
1661 /* Magic integer for velocities */
1662 magic_int=(int)readbufferfix((unsigned char *)data+bufloc,4);
1664 if (magic_int!=MAGIC_INT_VEL)
1669 /* Number of atoms. */
1670 natoms=(int)readbufferfix((unsigned char *)data+bufloc,4);
1672 /* Number of frames. */
1673 nframes=(int)readbufferfix((unsigned char *)data+bufloc,4);
1675 /* Initial coding. */
1676 initial_coding=(int)readbufferfix((unsigned char *)data+bufloc,4);
1678 /* Initial coding parameter. */
1679 initial_coding_parameter=(int)readbufferfix((unsigned char *)data+bufloc,4);
1682 coding=(int)readbufferfix((unsigned char *)data+bufloc,4);
1684 /* Coding parameter. */
1685 coding_parameter=(int)readbufferfix((unsigned char *)data+bufloc,4);
1688 *prec_lo=readbufferfix((unsigned char *)data+bufloc,4);
1690 *prec_hi=readbufferfix((unsigned char *)data+bufloc,4);
1692 /* Allocate the memory for the quantized positions */
1693 quant=malloc(natoms*nframes*3*sizeof *quant);
1694 /* The data block length. */
1695 length=(int)readbufferfix((unsigned char *)data+bufloc,4);
1697 /* The initial frame */
1698 coder=Ptngc_coder_init();
1699 rval=Ptngc_unpack_array(coder,(unsigned char*)data+bufloc,quant,natoms*3,
1700 initial_coding,initial_coding_parameter,natoms);
1701 Ptngc_coder_deinit(coder);
1704 /* Skip past the actual data block. */
1706 /* Obtain the actual positions for the initial block. */
1707 if ((initial_coding==TNG_COMPRESS_ALGO_VEL_STOPBIT_ONETOONE) ||
1708 (initial_coding==TNG_COMPRESS_ALGO_VEL_TRIPLET_ONETOONE) ||
1709 (initial_coding==TNG_COMPRESS_ALGO_VEL_BWLZH_ONETOONE))
1712 unquantize(veld,natoms,1,PRECISION(*prec_hi,*prec_lo),quant);
1714 unquantize_float(velf,natoms,1,(float)PRECISION(*prec_hi,*prec_lo),quant);
1716 memcpy(veli,quant,natoms*3*sizeof *veli);
1718 /* The remaining frames. */
1722 coder=Ptngc_coder_init();
1723 rval=Ptngc_unpack_array(coder,(unsigned char *)data+bufloc,quant+natoms*3,(nframes-1)*natoms*3,
1724 coding,coding_parameter,natoms);
1725 Ptngc_coder_deinit(coder);
1728 /* Inter-frame compression? */
1729 if ((coding==TNG_COMPRESS_ALGO_VEL_TRIPLET_INTER) ||
1730 (coding==TNG_COMPRESS_ALGO_VEL_STOPBIT_INTER) ||
1731 (coding==TNG_COMPRESS_ALGO_VEL_BWLZH_INTER))
1733 /* This requires that the first frame is already in one-to-one format. */
1735 unquantize_inter_differences(veld,natoms,nframes,PRECISION(*prec_hi,*prec_lo),quant);
1737 unquantize_inter_differences_float(velf,natoms,nframes,(float)PRECISION(*prec_hi,*prec_lo),quant);
1739 unquantize_inter_differences_int(veli,natoms,nframes,quant);
1741 /* One-to-one compression? */
1742 else if ((coding==TNG_COMPRESS_ALGO_VEL_STOPBIT_ONETOONE) ||
1743 (coding==TNG_COMPRESS_ALGO_VEL_TRIPLET_ONETOONE) ||
1744 (coding==TNG_COMPRESS_ALGO_VEL_BWLZH_ONETOONE))
1747 unquantize(veld+natoms*3,natoms,nframes-1,PRECISION(*prec_hi,*prec_lo),quant+natoms*3);
1749 unquantize_float(velf+natoms*3,natoms,nframes-1,(float)PRECISION(*prec_hi,*prec_lo),quant+natoms*3);
1751 memcpy(veli+natoms*3,quant+natoms*3,natoms*3*(nframes-1)*sizeof *veli);
1759 static int tng_compress_uncompress_vel(char *data,double *vel)
1761 unsigned long prec_hi, prec_lo;
1762 return tng_compress_uncompress_vel_gen(data,vel,NULL,NULL,&prec_hi,&prec_lo);
1765 static int tng_compress_uncompress_vel_float(char *data,float *vel)
1767 unsigned long prec_hi, prec_lo;
1768 return tng_compress_uncompress_vel_gen(data,NULL,vel,NULL,&prec_hi,&prec_lo);
1771 static int tng_compress_uncompress_vel_int(char *data,int *vel, unsigned long *prec_hi, unsigned long *prec_lo)
1773 return tng_compress_uncompress_vel_gen(data,NULL,NULL,vel,prec_hi,prec_lo);
1776 /* Uncompresses any tng compress block, positions or velocities. It determines whether it is positions or velocities from the data buffer. The return value is 0 if ok, and 1 if not.
1778 int DECLSPECDLLEXPORT tng_compress_uncompress(char *data,double *posvel)
1781 magic_int=(int)readbufferfix((unsigned char *)data,4);
1782 if (magic_int==MAGIC_INT_POS)
1783 return tng_compress_uncompress_pos(data,posvel);
1784 else if (magic_int==MAGIC_INT_VEL)
1785 return tng_compress_uncompress_vel(data,posvel);
1790 int DECLSPECDLLEXPORT tng_compress_uncompress_float(char *data,float *posvel)
1793 magic_int=(int)readbufferfix((unsigned char *)data,4);
1794 if (magic_int==MAGIC_INT_POS)
1795 return tng_compress_uncompress_pos_float(data,posvel);
1796 else if (magic_int==MAGIC_INT_VEL)
1797 return tng_compress_uncompress_vel_float(data,posvel);
1802 int DECLSPECDLLEXPORT tng_compress_uncompress_int(char *data,int *posvel, unsigned long *prec_hi, unsigned long *prec_lo)
1805 magic_int=(int)readbufferfix((unsigned char *)data,4);
1806 if (magic_int==MAGIC_INT_POS)
1807 return tng_compress_uncompress_pos_int(data,posvel,prec_hi,prec_lo);
1808 else if (magic_int==MAGIC_INT_VEL)
1809 return tng_compress_uncompress_vel_int(data,posvel,prec_hi,prec_lo);
1814 void DECLSPECDLLEXPORT tng_compress_int_to_double(int *posvel_int,unsigned long prec_hi, unsigned long prec_lo,
1815 int natoms,int nframes,
1816 double *posvel_double)
1818 unquantize(posvel_double,natoms,nframes,PRECISION(prec_hi,prec_lo),posvel_int);
1821 void DECLSPECDLLEXPORT tng_compress_int_to_float(int *posvel_int,unsigned long prec_hi, unsigned long prec_lo,
1822 int natoms,int nframes,
1823 float *posvel_float)
1825 unquantize_float(posvel_float,natoms,nframes,(float)PRECISION(prec_hi,prec_lo),posvel_int);
1828 static char *compress_algo_pos[TNG_COMPRESS_ALGO_MAX]={
1829 "Positions invalid algorithm",
1830 "Positions stopbits interframe",
1831 "Positions triplet interframe",
1832 "Positions triplet intraframe",
1833 "Positions invalid algorithm",
1835 "Positions invalid algorithm",
1836 "Positions triplet one to one",
1837 "Positions BWLZH interframe",
1838 "Positions BWLZH intraframe",
1842 static char *compress_algo_vel[TNG_COMPRESS_ALGO_MAX]={
1843 "Velocities invalid algorithm",
1844 "Velocities stopbits one to one",
1845 "Velocities triplet interframe",
1846 "Velocities triplet one to one",
1847 "Velocities invalid algorithm",
1848 "Velocities invalid algorithm",
1849 "Velocities stopbits interframe",
1850 "Velocities invalid algorithm",
1851 "Velocities BWLZH interframe",
1852 "Velocities BWLZH one to one",
1853 "Velocities invalid algorithm"
1856 char DECLSPECDLLEXPORT *tng_compress_initial_pos_algo(int *algo)
1861 if (i>=TNG_COMPRESS_ALGO_MAX)
1863 return compress_algo_pos[i];
1866 char DECLSPECDLLEXPORT *tng_compress_pos_algo(int *algo)
1871 if (i>=TNG_COMPRESS_ALGO_MAX)
1873 return compress_algo_pos[i];
1876 char DECLSPECDLLEXPORT *tng_compress_initial_vel_algo(int *algo)
1881 if (i>=TNG_COMPRESS_ALGO_MAX)
1883 return compress_algo_vel[i];
1886 char DECLSPECDLLEXPORT *tng_compress_vel_algo(int *algo)
1891 if (i>=TNG_COMPRESS_ALGO_MAX)
1893 return compress_algo_vel[i];