src/external/tng_io/src/tests/compression/testsuite.c

   1 /* tng compression routines */
   2
   3 /* Only modify testsuite.c
   4  *Then* run testsuite.sh to perform the test.
   5  */
   6
   7 #include <stdio.h>
   8 #include <stdlib.h>
   9 #include <string.h>
  10 #include <math.h>
  11 #include "../../../include/compression/tng_compress.h"
  12 #include TESTPARAM
  13
  14 #ifdef TEST_FLOAT
  15 #define REAL float
  16 #else
  17 #define REAL double
  18 #endif
  19
  20 #ifndef TNG_COMPRESS_FILES_DIR
  21 #define TNG_COMPRESS_FILES_DIR ""
  22 #endif
  23
  24 #define FUDGE 1.1 /* 10% off target precision is acceptable */
  25
  26 static void keepinbox(int *val)
  27 {
  28   while (val[0]>INTMAX1)
  29     val[0]-=(INTMAX1-INTMIN1+1);
  30   while (val[0]<INTMIN1)
  31     val[0]+=(INTMAX1-INTMIN1+1);
  32   while (val[1]>INTMAX2)
  33     val[1]-=(INTMAX2-INTMIN2+1);
  34   while (val[1]<INTMIN2)
  35     val[1]+=(INTMAX2-INTMIN2+1);
  36   while (val[2]>INTMAX3)
  37     val[2]-=(INTMAX3-INTMIN3+1);
  38   while (val[2]<INTMIN3)
  39     val[2]+=(INTMAX3-INTMIN3+1);
  40 }
  41
  42 static int intsintable[128]={
  43 0 , 3215 , 6423 , 9615 , 12785 , 15923 , 19023 , 22078 ,
  44 25079 , 28019 , 30892 , 33691 , 36409 , 39039 , 41574 , 44010 ,
  45 46340 , 48558 , 50659 , 52638 , 54490 , 56211 , 57796 , 59242 ,
  46 60546 , 61704 , 62713 , 63570 , 64275 , 64825 , 65219 , 65456 ,
  47 65535 , 65456 , 65219 , 64825 , 64275 , 63570 , 62713 , 61704 ,
  48 60546 , 59242 , 57796 , 56211 , 54490 , 52638 , 50659 , 48558 ,
  49 46340 , 44010 , 41574 , 39039 , 36409 , 33691 , 30892 , 28019 ,
  50 25079 , 22078 , 19023 , 15923 , 12785 , 9615 , 6423 , 3215 ,
  51 0 , -3215 , -6423 , -9615 , -12785 , -15923 , -19023 , -22078 ,
  52 -25079 , -28019 , -30892 , -33691 , -36409 , -39039 , -41574 , -44010 ,
  53 -46340 , -48558 , -50659 , -52638 , -54490 , -56211 , -57796 , -59242 ,
  54 -60546 , -61704 , -62713 , -63570 , -64275 , -64825 , -65219 , -65456 ,
  55 -65535 , -65456 , -65219 , -64825 , -64275 , -63570 , -62713 , -61704 ,
  56 -60546 , -59242 , -57796 , -56211 , -54490 , -52638 , -50659 , -48558 ,
  57 -46340 , -44010 , -41574 , -39039 , -36409 , -33691 , -30892 , -28019 ,
  58 -25079 , -22078 , -19023 , -15923 , -12785 , -9615 , -6423 , -3215 ,
  59 };
  60
  61 static int intsin(int i)
  62 {
  63   int sign=1;
  64   if (i<0)
  65     {
  66       i=0;
  67       sign=-1;
  68     }
  69   return sign*intsintable[i%128];
  70 }
  71
  72 static int intcos(int i)
  73 {
  74   if (i<0)
  75     i=0;
  76   return intsin(i+32);
  77 }
  78
  79 static void molecule(int *target,
  80                      int *base,
  81                      int length,
  82                      int scale, int *direction,
  83                      int flip,
  84                      int iframe)
  85 {
  86   int i;
  87   for (i=0; i<length; i++)
  88     {
  89       int ifl=i;
  90       if ((i==0) && (flip) && (length>1))
  91         ifl=1;
  92       else if ((i==1) && (flip) && (length>1))
  93         ifl=0;
  94       target[ifl*3]=base[0]+(intsin((i+iframe)*direction[0])*scale)/256;
  95       target[ifl*3+1]=base[1]+(intcos((i+iframe)*direction[1])*scale)/256;
  96       target[ifl*3+2]=base[2]+(intcos((i+iframe)*direction[2])*scale)/256;
  97       keepinbox(target+ifl*3);
  98     }
  99 }
 100
 101 #ifndef FRAMESCALE
 102 #define FRAMESCALE 1
 103 #endif
 104
 105 static void genibox(int *intbox, int iframe)
 106 {
 107   int molecule_length=1;
 108   int molpos[3];
 109   int direction[3]={1,1,1};
 110   int scale=1;
 111   int flip=0;
 112   int i=0;
 113   molpos[0]=intsin(iframe*FRAMESCALE)/32;
 114   molpos[1]=1+intcos(iframe*FRAMESCALE)/32;
 115   molpos[2]=2+intsin(iframe*FRAMESCALE)/16;
 116   keepinbox(molpos);
 117   while (i<NATOMS)
 118     {
 119       int this_mol_length=molecule_length;
 120       int dir;
 121 #ifdef REGULAR
 122       this_mol_length=4;
 123       flip=0;
 124       scale=1;
 125 #endif
 126       if (i+this_mol_length>NATOMS)
 127         this_mol_length=NATOMS-i;
 128       /* We must test the large rle as well. This requires special
 129          sequencies to get triggered. So insert these from time to
 130          time */
 131 #ifndef REGULAR
 132       if ((i%10)==0)
 133         {
 134           int j;
 135           intbox[i*3]=molpos[0];
 136           intbox[i*3+1]=molpos[1];
 137           intbox[i*3+2]=molpos[2];
 138           for (j=1; j<this_mol_length; j++)
 139             {
 140               intbox[(i+j)*3]=intbox[(i+j-1)*3]+(INTMAX1-INTMIN1+1)/5;
 141               intbox[(i+j)*3+1]=intbox[(i+j-1)*3+1]+(INTMAX2-INTMIN2+1)/5;
 142               intbox[(i+j)*3+2]=intbox[(i+j-1)*3+2]+(INTMAX3-INTMIN3+1)/5;
 143               keepinbox(intbox+(i+j)*3);
 144             }
 145         }
 146       else
 147 #endif
 148         molecule(intbox+i*3,molpos,this_mol_length,scale,direction,flip,iframe*FRAMESCALE);
 149       i+=this_mol_length;
 150       dir=1;
 151       if (intsin(i*3)<0)
 152         dir=-1;
 153       molpos[0]+=dir*(INTMAX1-INTMIN1+1)/20;
 154       dir=1;
 155       if (intsin(i*5)<0)
 156         dir=-1;
 157       molpos[1]+=dir*(INTMAX2-INTMIN2+1)/20;
 158       dir=1;
 159       if (intsin(i*7)<0)
 160         dir=-1;
 161       molpos[2]+=dir*(INTMAX3-INTMIN3+1)/20;
 162       keepinbox(molpos);
 163
 164       direction[0]=((direction[0]+1)%7)+1;
 165       direction[1]=((direction[1]+1)%3)+1;
 166       direction[2]=((direction[2]+1)%6)+1;
 167
 168       scale++;
 169       if (scale>5)
 170         scale=1;
 171
 172       molecule_length++;
 173       if (molecule_length>30)
 174         molecule_length=1;
 175       if (i%9)
 176         flip=1-flip;
 177     }
 178 }
 179
 180 static void genivelbox(int *intvelbox, int iframe)
 181 {
 182   int i;
 183   for (i=0; i<NATOMS; i++)
 184     {
 185 #ifdef VELINTMUL
 186       intvelbox[i*3]=((intsin((i+iframe*FRAMESCALE)*3))/10)*VELINTMUL+i;
 187       intvelbox[i*3+1]=1+((intcos((i+iframe*FRAMESCALE)*5))/10)*VELINTMUL+i;
 188       intvelbox[i*3+2]=2+((intsin((i+iframe*FRAMESCALE)*7)+intcos((i+iframe*FRAMESCALE)*9))/20)*VELINTMUL+i;
 189 #else
 190       intvelbox[i*3]=((intsin((i+iframe*FRAMESCALE)*3))/10);
 191       intvelbox[i*3+1]=1+((intcos((i+iframe*FRAMESCALE)*5))/10);
 192       intvelbox[i*3+2]=2+((intsin((i+iframe*FRAMESCALE)*7)+intcos((i+iframe*FRAMESCALE)*9))/20);
 193 #endif
 194     }
 195 }
 196
 197 #ifndef STRIDE1
 198 #define STRIDE1 3
 199 #endif
 200
 201 #ifndef STRIDE2
 202 #define STRIDE2 3
 203 #endif
 204
 205 #ifndef GENPRECISION
 206 #define GENPRECISION PRECISION
 207 #endif
 208
 209 #ifndef GENVELPRECISION
 210 #define GENVELPRECISION VELPRECISION
 211 #endif
 212
 213 static void realbox(int *intbox, REAL *realbox, int stride)
 214 {
 215   int i,j;
 216   for (i=0; i<NATOMS; i++)
 217     {
 218       for (j=0; j<3; j++)
 219         realbox[i*stride+j]=(REAL)(intbox[i*3+j]*GENPRECISION*SCALE);
 220       for (j=3; j<stride; j++)
 221         realbox[i*stride+j]=0.;
 222     }
 223 }
 224
 225 static void realvelbox(int *intbox, REAL *realbox, int stride)
 226 {
 227   int i,j;
 228   for (i=0; i<NATOMS; i++)
 229     {
 230       for (j=0; j<3; j++)
 231         realbox[i*stride+j]=(REAL)(intbox[i*3+j]*GENVELPRECISION*SCALE);
 232       for (j=3; j<stride; j++)
 233         realbox[i*stride+j]=0.;
 234     }
 235 }
 236
 237 static int equalarr(REAL *arr1, REAL *arr2, REAL prec, int len, int itemlen, int stride1, int stride2)
 238 {
 239   REAL maxdiff=0.;
 240   int i,j;
 241   for (i=0; i<len; i++)
 242     {
 243       for (j=0; j<itemlen; j++)
 244         if (fabs(arr1[i*stride1+j]-arr2[i*stride2+j])>maxdiff)
 245           maxdiff=(REAL)fabs(arr1[i*stride1+j]-arr2[i*stride2+j]);
 246     }
 247 #if 0
 248   for (i=0; i<len; i++)
 249     {
 250       for (j=0; j<itemlen; j++)
 251         printf("%d %d: %g %g\n",i,j,arr1[i*stride1+j],arr2[i*stride2+j]);
 252     }
 253 #endif
 254 #if 0
 255   fprintf(stderr,"Error is %g. Acceptable error is %g.\n",maxdiff,prec*0.5*FUDGE);
 256 #endif
 257   if (maxdiff>prec*0.5*FUDGE)
 258     {
 259       return 0;
 260     }
 261   else
 262     return 1;
 263 }
 264
 265 struct tng_file
 266 {
 267   FILE *f;
 268   int natoms;
 269   int chunky;
 270   REAL precision;
 271   REAL velprecision;
 272   int initial_coding;
 273   int initial_coding_parameter;
 274   int coding;
 275   int coding_parameter;
 276   int initial_velcoding;
 277   int initial_velcoding_parameter;
 278   int velcoding;
 279   int velcoding_parameter;
 280   int speed;
 281   int nframes;
 282   int nframes_delivered;
 283   int writevel;
 284   REAL *pos;
 285   REAL *vel;
 286   int *ipos;
 287   int *ivel;
 288   unsigned int prec_hi, prec_lo;
 289   unsigned int velprec_hi, velprec_lo;
 290 };
 291
 292 static size_t fwrite_int_le(int *x,FILE *f)
 293 {
 294   unsigned char c[4];
 295   unsigned int i=(unsigned int)*x;
 296   c[0]=(unsigned char)(i&0xFFU);
 297   c[1]=(unsigned char)((i>>8)&0xFFU);
 298   c[2]=(unsigned char)((i>>16)&0xFFU);
 299   c[3]=(unsigned char)((i>>24)&0xFFU);
 300   return fwrite(c,1,4,f);
 301 }
 302
 303 static size_t fread_int_le(int *x,FILE *f)
 304 {
 305   unsigned char c[4];
 306   unsigned int i;
 307   size_t n=fread(c,1,4,f);
 308   if (n)
 309     {
 310       i=(((unsigned int)c[3])<<24)|(((unsigned int)c[2])<<16)|(((unsigned int)c[1])<<8)|((unsigned int)c[0]);
 311       *x=(int)i;
 312     }
 313   return n;
 314 }
 315
 316 static struct tng_file *open_tng_file_write(char *filename,
 317                                             int natoms,int chunky,
 318                                             REAL precision,
 319                                             int writevel,
 320                                             REAL velprecision,
 321                                             int initial_coding,
 322                                             int initial_coding_parameter,
 323                                             int coding,
 324                                             int coding_parameter,
 325                                             int initial_velcoding,
 326                                             int initial_velcoding_parameter,
 327                                             int velcoding,
 328                                             int velcoding_parameter,
 329                                             int speed)
 330 {
 331   struct tng_file *tng_file=malloc(sizeof *tng_file);
 332   tng_file->pos=NULL;
 333   tng_file->vel=NULL;
 334   tng_file->ipos=NULL;
 335   tng_file->ivel=NULL;
 336   tng_file->nframes=0;
 337   tng_file->chunky=chunky;
 338   tng_file->precision=precision;
 339   tng_file->natoms=natoms;
 340   tng_file->writevel=writevel;
 341   tng_file->velprecision=velprecision;
 342   tng_file->initial_coding=initial_coding;
 343   tng_file->initial_coding_parameter=initial_coding_parameter;
 344   tng_file->coding=coding;
 345   tng_file->coding_parameter=coding_parameter;
 346   tng_file->initial_velcoding=initial_velcoding;
 347   tng_file->initial_velcoding_parameter=initial_velcoding_parameter;
 348   tng_file->velcoding=velcoding;
 349   tng_file->velcoding_parameter=velcoding_parameter;
 350   tng_file->speed=speed;
 351   tng_file->pos=malloc(natoms*chunky*3*sizeof *tng_file->pos);
 352   tng_file->f=fopen(filename,"wb");
 353   if (writevel)
 354     tng_file->vel=malloc(natoms*chunky*3*sizeof *tng_file->vel);
 355   fwrite_int_le(&natoms,tng_file->f);
 356   return tng_file;
 357 }
 358
 359 static struct tng_file *open_tng_file_write_int(char *filename,
 360                                                 int natoms,int chunky,
 361                                                 int writevel,
 362                                                 int initial_coding,
 363                                                 int initial_coding_parameter,
 364                                                 int coding,
 365                                                 int coding_parameter,
 366                                                 int initial_velcoding,
 367                                                 int initial_velcoding_parameter,
 368                                                 int velcoding,
 369                                                 int velcoding_parameter,
 370                                                 int speed)
 371 {
 372   struct tng_file *tng_file=malloc(sizeof *tng_file);
 373   tng_file->pos=NULL;
 374   tng_file->vel=NULL;
 375   tng_file->ipos=NULL;
 376   tng_file->ivel=NULL;
 377   tng_file->nframes=0;
 378   tng_file->chunky=chunky;
 379   tng_file->natoms=natoms;
 380   tng_file->writevel=writevel;
 381   tng_file->initial_coding=initial_coding;
 382   tng_file->initial_coding_parameter=initial_coding_parameter;
 383   tng_file->coding=coding;
 384   tng_file->coding_parameter=coding_parameter;
 385   tng_file->initial_velcoding=initial_velcoding;
 386   tng_file->initial_velcoding_parameter=initial_velcoding_parameter;
 387   tng_file->velcoding=velcoding;
 388   tng_file->velcoding_parameter=velcoding_parameter;
 389   tng_file->speed=speed;
 390   tng_file->ipos=malloc(natoms*chunky*3*sizeof *tng_file->ipos);
 391   tng_file->f=fopen(filename,"wb");
 392   if (writevel)
 393     tng_file->ivel=malloc(natoms*chunky*3*sizeof *tng_file->ivel);
 394   fwrite_int_le(&natoms,tng_file->f);
 395   return tng_file;
 396 }
 397
 398 static void flush_tng_frames(struct tng_file *tng_file,
 399                              unsigned long prec_hi, unsigned long prec_lo,
 400                              unsigned long velprec_hi, unsigned long velprec_lo)
 401 {
 402   int algo[4];
 403   char *buf;
 404   int nitems;
 405
 406   /* Make sure these variables are used to avoid compilation warnings */
 407   (void)prec_hi;
 408   (void)prec_lo;
 409   (void)velprec_hi;
 410   (void)velprec_lo;
 411
 412   fwrite_int_le(&tng_file->nframes,tng_file->f);
 413   algo[0]=tng_file->initial_coding;
 414   algo[1]=tng_file->initial_coding_parameter;
 415   algo[2]=tng_file->coding;
 416   algo[3]=tng_file->coding_parameter;
 417 #ifdef RECOMPRESS
 418   buf=tng_compress_pos_int(tng_file->ipos,
 419                              tng_file->natoms,
 420                              tng_file->nframes,
 421                              prec_hi,prec_lo,
 422                              tng_file->speed,algo,&nitems);
 423 #else /* RECOMPRESS */
 424 #ifdef TEST_FLOAT
 425   buf=tng_compress_pos_float(tng_file->pos,
 426                              tng_file->natoms,
 427                              tng_file->nframes,
 428                              tng_file->precision,
 429                              tng_file->speed,algo,&nitems);
 430 #else  /* TEST_FLOAT */
 431   buf=tng_compress_pos(tng_file->pos,
 432                        tng_file->natoms,
 433                        tng_file->nframes,
 434                        tng_file->precision,
 435                        tng_file->speed,algo,&nitems);
 436 #endif  /* TEST_FLOAT */
 437 #endif /* RECOMPRESS */
 438   tng_file->initial_coding=algo[0];
 439   tng_file->initial_coding_parameter=algo[1];
 440   tng_file->coding=algo[2];
 441   tng_file->coding_parameter=algo[3];
 442   fwrite_int_le(&nitems,tng_file->f);
 443   fwrite(buf,1,nitems,tng_file->f);
 444   free(buf);
 445   if (tng_file->writevel)
 446     {
 447       algo[0]=tng_file->initial_velcoding;
 448       algo[1]=tng_file->initial_velcoding_parameter;
 449       algo[2]=tng_file->velcoding;
 450       algo[3]=tng_file->velcoding_parameter;
 451 #ifdef RECOMPRESS
 452       buf=tng_compress_vel_int(tng_file->ivel,
 453                                tng_file->natoms,
 454                                tng_file->nframes,
 455                                velprec_hi,velprec_lo,
 456                                tng_file->speed,algo,&nitems);
 457 #else /* RECOMPRESS */
 458 #ifdef TEST_FLOAT
 459       buf=tng_compress_vel_float(tng_file->vel,
 460                            tng_file->natoms,
 461                            tng_file->nframes,
 462                            tng_file->velprecision,
 463                            tng_file->speed,algo,&nitems);
 464 #else /* TEST_FLOAT */
 465       buf=tng_compress_vel(tng_file->vel,
 466                            tng_file->natoms,
 467                            tng_file->nframes,
 468                            tng_file->velprecision,
 469                            tng_file->speed,algo,&nitems);
 470 #endif /* TEST_FLOAT */
 471 #endif /* RECOMPRESS */
 472       tng_file->initial_velcoding=algo[0];
 473       tng_file->initial_velcoding_parameter=algo[1];
 474       tng_file->velcoding=algo[2];
 475       tng_file->velcoding_parameter=algo[3];
 476       fwrite_int_le(&nitems,tng_file->f);
 477       fwrite(buf,1,nitems,tng_file->f);
 478       free(buf);
 479     }
 480   tng_file->nframes=0;
 481 }
 482
 483 static void write_tng_file(struct tng_file *tng_file,
 484                            REAL *pos,REAL *vel)
 485 {
 486   memcpy(tng_file->pos+tng_file->nframes*tng_file->natoms*3,pos,tng_file->natoms*3*sizeof *tng_file->pos);
 487   if (tng_file->writevel)
 488     memcpy(tng_file->vel+tng_file->nframes*tng_file->natoms*3,vel,tng_file->natoms*3*sizeof *tng_file->vel);
 489   tng_file->nframes++;
 490   if (tng_file->nframes==tng_file->chunky)
 491     flush_tng_frames(tng_file,0,0,0,0);
 492 }
 493
 494 static void write_tng_file_int(struct tng_file *tng_file,
 495                                int *ipos,int *ivel,
 496                                unsigned long prec_hi, unsigned long prec_lo,
 497                                unsigned long velprec_hi, unsigned long velprec_lo)
 498 {
 499   memcpy(tng_file->ipos+tng_file->nframes*tng_file->natoms*3,ipos,tng_file->natoms*3*sizeof *tng_file->ipos);
 500   if (tng_file->writevel)
 501     memcpy(tng_file->ivel+tng_file->nframes*tng_file->natoms*3,ivel,tng_file->natoms*3*sizeof *tng_file->ivel);
 502   tng_file->nframes++;
 503   if (tng_file->nframes==tng_file->chunky)
 504     flush_tng_frames(tng_file,prec_hi,prec_lo,velprec_hi,velprec_lo);
 505   tng_file->prec_hi=prec_hi;
 506   tng_file->prec_lo=prec_lo;
 507   tng_file->velprec_hi=velprec_hi;
 508   tng_file->velprec_lo=velprec_lo;
 509 }
 510
 511 static void close_tng_file_write(struct tng_file *tng_file)
 512 {
 513   if (tng_file->nframes)
 514     flush_tng_frames(tng_file,tng_file->prec_hi,tng_file->prec_lo,tng_file->velprec_hi,tng_file->velprec_lo);
 515   fclose(tng_file->f);
 516   free(tng_file->pos);
 517   free(tng_file->vel);
 518   free(tng_file->ipos);
 519   free(tng_file->ivel);
 520   free(tng_file);
 521 }
 522
 523 static struct tng_file *open_tng_file_read(char *filename, int writevel)
 524 {
 525   struct tng_file *tng_file=malloc(sizeof *tng_file);
 526   tng_file->pos=NULL;
 527   tng_file->vel=NULL;
 528   tng_file->ipos=NULL;
 529   tng_file->ivel=NULL;
 530   tng_file->f=fopen(filename,"rb");
 531   tng_file->nframes=0;
 532   tng_file->nframes_delivered=0;
 533   tng_file->writevel=writevel;
 534   if (tng_file->f)
 535     fread_int_le(&tng_file->natoms,tng_file->f);
 536   else
 537     {
 538       free(tng_file);
 539       tng_file=NULL;
 540     }
 541   return tng_file;
 542 }
 543
 544 static struct tng_file *open_tng_file_read_int(char *filename, int writevel)
 545 {
 546   struct tng_file *tng_file=malloc(sizeof *tng_file);
 547   tng_file->pos=NULL;
 548   tng_file->vel=NULL;
 549   tng_file->ipos=NULL;
 550   tng_file->ivel=NULL;
 551   tng_file->f=fopen(filename,"rb");
 552   tng_file->nframes=0;
 553   tng_file->nframes_delivered=0;
 554   tng_file->writevel=writevel;
 555   if (tng_file->f)
 556     fread_int_le(&tng_file->natoms,tng_file->f);
 557   else
 558     {
 559       free(tng_file);
 560       tng_file=NULL;
 561     }
 562   return tng_file;
 563 }
 564
 565 static int read_tng_file(struct tng_file *tng_file,
 566                          REAL *pos,
 567                          REAL *vel)
 568 {
 569   if (tng_file->nframes==tng_file->nframes_delivered)
 570     {
 571       int nitems;
 572       char *buf;
 573       free(tng_file->pos);
 574       free(tng_file->vel);
 575       if (!fread_int_le(&tng_file->nframes,tng_file->f))
 576         return 1;
 577       if (!fread_int_le(&nitems,tng_file->f))
 578         return 1;
 579       buf=malloc(nitems);
 580       if (!fread(buf,1,nitems,tng_file->f))
 581       {
 582           free(buf);
 583           return 1;
 584       }
 585       tng_file->pos=malloc(tng_file->natoms*tng_file->nframes*3*sizeof *tng_file->pos);
 586       if (tng_file->writevel)
 587         tng_file->vel=malloc(tng_file->natoms*tng_file->nframes*3*sizeof *tng_file->vel);
 588 #if 0
 589       {
 590         int natoms, nframes, algo[4];
 591         double precision;
 592         int ivel;
 593         char *initial_coding, *coding;
 594         tng_compress_inquire(buf,&ivel,&natoms,&nframes,&precision,algo);
 595         initial_coding=tng_compress_initial_pos_algo(algo);
 596         coding=tng_compress_pos_algo(algo);
 597         printf("ivel=%d natoms=%d nframes=%d precision=%g initial pos=%s pos=%s\n",ivel,natoms,nframes,precision,initial_coding,coding);
 598       }
 599 #endif
 600 #ifdef TEST_FLOAT
 601       tng_compress_uncompress_float(buf,tng_file->pos);
 602 #else
 603       tng_compress_uncompress(buf,tng_file->pos);
 604 #endif
 605       free(buf);
 606       if (tng_file->writevel)
 607         {
 608           if (!fread_int_le(&nitems,tng_file->f))
 609             return 1;
 610           buf=malloc(nitems);
 611           if (!fread(buf,1,nitems,tng_file->f))
 612           {
 613               free(buf);
 614               return 1;
 615           }
 616 #if 0
 617           {
 618             int natoms, nframes, algo[4];
 619             double precision;
 620             int ivel;
 621             char *initial_coding, *coding;
 622             tng_compress_inquire(buf,&ivel,&natoms,&nframes,&precision,algo);
 623             initial_coding=tng_compress_initial_vel_algo(algo);
 624             coding=tng_compress_vel_algo(algo);
 625             printf("ivel=%d natoms=%d nframes=%d precision=%g initial vel=%s vel=%s\n",ivel,natoms,nframes,precision,initial_coding,coding);
 626           }
 627 #endif
 628 #ifdef TEST_FLOAT
 629           tng_compress_uncompress_float(buf,tng_file->vel);
 630 #else
 631           tng_compress_uncompress(buf,tng_file->vel);
 632 #endif
 633           free(buf);
 634         }
 635       tng_file->nframes_delivered=0;
 636     }
 637   memcpy(pos,tng_file->pos+tng_file->nframes_delivered*tng_file->natoms*3,tng_file->natoms*3*sizeof *pos);
 638   if (tng_file->writevel)
 639     memcpy(vel,tng_file->vel+tng_file->nframes_delivered*tng_file->natoms*3,tng_file->natoms*3*sizeof *vel);
 640   tng_file->nframes_delivered++;
 641   return 0;
 642 }
 643
 644 static int read_tng_file_int(struct tng_file *tng_file,
 645                              int *ipos,
 646                              int *ivel,
 647                              unsigned long *prec_hi, unsigned long *prec_lo,
 648                              unsigned long *velprec_hi, unsigned long *velprec_lo)
 649 {
 650   if (tng_file->nframes==tng_file->nframes_delivered)
 651     {
 652       int nitems;
 653       char *buf;
 654       free(tng_file->ipos);
 655       free(tng_file->ivel);
 656       if (!fread_int_le(&tng_file->nframes,tng_file->f))
 657         return 1;
 658       if (!fread_int_le(&nitems,tng_file->f))
 659         return 1;
 660       buf=malloc(nitems);
 661       if (!fread(buf,1,nitems,tng_file->f))
 662       {
 663           free(buf);
 664           return 1;
 665       }
 666       tng_file->ipos=malloc(tng_file->natoms*tng_file->nframes*3*sizeof *tng_file->ipos);
 667       if (tng_file->writevel)
 668         tng_file->ivel=malloc(tng_file->natoms*tng_file->nframes*3*sizeof *tng_file->ivel);
 669       tng_compress_uncompress_int(buf,tng_file->ipos,prec_hi,prec_lo);
 670       free(buf);
 671       if (tng_file->writevel)
 672         {
 673           if (!fread_int_le(&nitems,tng_file->f))
 674             return 1;
 675           buf=malloc(nitems);
 676           if (!fread(buf,1,nitems,tng_file->f))
 677           {
 678               free(buf);
 679               return 1;
 680           }
 681           tng_compress_uncompress_int(buf,tng_file->ivel,velprec_hi,velprec_lo);
 682           free(buf);
 683         }
 684       tng_file->nframes_delivered=0;
 685     }
 686   memcpy(ipos,tng_file->ipos+tng_file->nframes_delivered*tng_file->natoms*3,tng_file->natoms*3*sizeof *ipos);
 687   if (tng_file->writevel)
 688     memcpy(ivel,tng_file->ivel+tng_file->nframes_delivered*tng_file->natoms*3,tng_file->natoms*3*sizeof *ivel);
 689   tng_file->nframes_delivered++;
 690   return 0;
 691 }
 692
 693 static void close_tng_file_read(struct tng_file *tng_file)
 694 {
 695   free(tng_file->vel);
 696   free(tng_file->pos);
 697   free(tng_file->ivel);
 698   free(tng_file->ipos);
 699   fclose(tng_file->f);
 700   free(tng_file);
 701 }
 702
 703
 704
 705 #ifndef EXPECTED_FILESIZE
 706 #define EXPECTED_FILESIZE 1
 707 #endif
 708
 709 #ifndef INITIALVELCODING
 710 #define INITIALVELCODING -1
 711 #endif
 712 #ifndef INITIALVELCODINGPARAMETER
 713 #define INITIALVELCODINGPARAMETER -1
 714 #endif
 715
 716 #ifndef SPEED
 717 #define SPEED 5
 718 #endif
 719
 720 /* Return value 1 means file error.
 721    Return value 4 means coding error in coordinates.
 722    Return value 5 means coding error in velocities.
 723    Return value 9 means filesize seems too off.
 724
 725    Return value 100+ means test specific error.
 726  */
 727 static int algotest()
 728 {
 729   int i;
 730   int *intbox=malloc(NATOMS*3*sizeof *intbox);
 731   int *intvelbox=malloc(NATOMS*3*sizeof *intvelbox);
 732 #ifdef RECOMPRESS
 733   unsigned long pos_prec_hi,pos_prec_lo;
 734   unsigned long vel_prec_hi,vel_prec_lo;
 735 #endif
 736   REAL *box1=malloc(NATOMS*STRIDE1*sizeof *box1);
 737   REAL *velbox1=malloc(NATOMS*STRIDE1*sizeof *velbox1);
 738   int startframe=0;
 739   int endframe=NFRAMES;
 740 #ifdef GEN
 741   FILE *file;
 742   REAL filesize;
 743   REAL *box2=0;
 744   REAL *velbox2=0;
 745 #else
 746   int readreturn;
 747   REAL *box2=malloc(NATOMS*STRIDE2*sizeof *box2);
 748   REAL *velbox2=malloc(NATOMS*STRIDE2*sizeof *velbox2);
 749 #endif
 750 #ifdef RECOMPRESS
 751   void *dumpfile=open_tng_file_write_int(TNG_COMPRESS_FILES_DIR FILENAME,NATOMS,CHUNKY,
 752                                      WRITEVEL,
 753                                      INITIALCODING,
 754                                      INITIALCODINGPARAMETER,CODING,CODINGPARAMETER,
 755                                      INITIALVELCODING,INITIALVELCODINGPARAMETER,
 756                                      VELCODING,VELCODINGPARAMETER,SPEED);
 757   void *dumpfile_recompress=open_tng_file_read_int(TNG_COMPRESS_FILES_DIR RECOMPRESS,WRITEVEL);
 758   if (!dumpfile_recompress)
 759     {
 760       free(intbox);
 761       free(intvelbox);
 762       free(box1);
 763       free(velbox1);
 764       if(box2)
 765         free(box2);
 766       if(velbox2)
 767         free(velbox2);
 768       return 1;
 769     }
 770 #else /* RECOMPRESS */
 771 #ifdef GEN
 772   void *dumpfile=open_tng_file_write(TNG_COMPRESS_FILES_DIR FILENAME,NATOMS,CHUNKY,
 773                                      PRECISION,WRITEVEL,VELPRECISION,
 774                                      INITIALCODING,
 775                                      INITIALCODINGPARAMETER,CODING,CODINGPARAMETER,
 776                                      INITIALVELCODING,INITIALVELCODINGPARAMETER,
 777                                      VELCODING,VELCODINGPARAMETER,SPEED);
 778 #else
 779   void *dumpfile=open_tng_file_read(TNG_COMPRESS_FILES_DIR FILENAME,WRITEVEL);
 780 #endif
 781 #endif /* RECOMPRESS */
 782   if (!dumpfile)
 783     {
 784       free(intbox);
 785       free(intvelbox);
 786       free(box1);
 787       free(velbox1);
 788       if(box2)
 789         free(box2);
 790       if(velbox2)
 791         free(velbox2);
 792       return 1;
 793     }
 794   for (i=startframe; i<endframe; i++)
 795     {
 796 #ifdef RECOMPRESS
 797       unsigned long prec_hi, prec_lo;
 798       unsigned long velprec_hi, velprec_lo;
 799       if (read_tng_file_int(dumpfile_recompress,intbox,intvelbox,&prec_hi,&prec_lo,&velprec_hi,&velprec_lo))
 800         return 1;
 801       write_tng_file_int(dumpfile,intbox,intvelbox,prec_hi,prec_lo,velprec_hi,velprec_lo);
 802 #else /* RECOMPRESS */
 803       genibox(intbox,i);
 804       realbox(intbox,box1,STRIDE1);
 805 #if WRITEVEL
 806       genivelbox(intvelbox,i);
 807       realvelbox(intvelbox,velbox1,STRIDE1);
 808 #endif
 809 #ifdef GEN
 810       write_tng_file(dumpfile,box1,velbox1);
 811 #else /* GEN */
 812 #ifdef INTTOFLOAT
 813       {
 814         unsigned long prec_hi, prec_lo;
 815         unsigned long velprec_hi, velprec_lo;
 816         readreturn=read_tng_file_int(dumpfile,intbox,intvelbox,&prec_hi,&prec_lo,&velprec_hi,&velprec_lo);
 817         if (!readreturn)
 818           {
 819             tng_compress_int_to_float(intbox,prec_hi,prec_lo,NATOMS,1,box2);
 820 #if WRITEVEL
 821             tng_compress_int_to_float(intvelbox,velprec_hi,velprec_lo,NATOMS,1,velbox2);
 822 #endif
 823           }
 824       }
 825 #else /* INTTOFLOAT */
 826 #ifdef INTTODOUBLE
 827       {
 828         unsigned long prec_hi, prec_lo;
 829         unsigned long velprec_hi, velprec_lo;
 830         readreturn=read_tng_file_int(dumpfile,intbox,intvelbox,&prec_hi,&prec_lo,&velprec_hi,&velprec_lo);
 831         if (!readreturn)
 832           {
 833             tng_compress_int_to_double(intbox,prec_hi,prec_lo,NATOMS,1,box2);
 834 #if WRITEVEL
 835             tng_compress_int_to_double(intvelbox,velprec_hi,velprec_lo,NATOMS,1,velbox2);
 836 #endif
 837           }
 838       }
 839 #else /* INTTODOUBLE */
 840       readreturn=read_tng_file(dumpfile,box2,velbox2);
 841 #endif /* INTTODOUBLE */
 842 #endif /* INTTOFLOAT */
 843       if (readreturn==1) /* general read error  */
 844         {
 845           free(intbox);
 846           free(intvelbox);
 847           free(box1);
 848           free(velbox1);
 849           if(box2)
 850             free(box2);
 851           if(velbox2)
 852             free(velbox2);
 853           return 1;
 854         }
 855 #endif /* GEN */
 856 #ifndef GEN
 857       /* Check for equality of boxes. */
 858       if (!equalarr(box1,box2,(REAL)PRECISION,NATOMS,3,STRIDE1,STRIDE2))
 859         {
 860           free(intbox);
 861           free(intvelbox);
 862           free(box1);
 863           free(velbox1);
 864           if(box2)
 865             free(box2);
 866           if(velbox2)
 867             free(velbox2);
 868           return 4;
 869         }
 870 #if WRITEVEL
 871       if (!equalarr(velbox1,velbox2,(REAL)VELPRECISION,NATOMS,3,STRIDE1,STRIDE2))
 872         {
 873           free(intbox);
 874           free(intvelbox);
 875           free(box1);
 876           free(velbox1);
 877           if(box2)
 878             free(box2);
 879           if(velbox2)
 880             free(velbox2);
 881           return 5;
 882         }
 883 #endif
 884 #endif /* GEN */
 885 #endif /* RECOMPRESS */
 886     }
 887 #ifdef GEN
 888   close_tng_file_write(dumpfile);
 889 #else
 890   close_tng_file_read(dumpfile);
 891 #endif
 892 #ifdef RECOMPRESS
 893   close_tng_file_read(dumpfile_recompress);
 894 #endif
 895 #ifdef GEN
 896   /* Check against expected filesize for this test. */
 897   if (!(file=fopen(TNG_COMPRESS_FILES_DIR FILENAME,"rb")))
 898     {
 899       fprintf(stderr,"ERROR: Cannot open file "TNG_COMPRESS_FILES_DIR FILENAME"\n");
 900       exit(EXIT_FAILURE);
 901     }
 902   filesize=0;
 903   while(1)
 904     {
 905       char b;
 906       if (!fread(&b,1,1,file))
 907         break;
 908       filesize++;
 909     }
 910   fclose(file);
 911   if (filesize>0)
 912     {
 913       if ((fabs(filesize-EXPECTED_FILESIZE)/EXPECTED_FILESIZE)>0.05)
 914       {
 915           free(intvelbox);
 916           free(box1);
 917           free(velbox1);
 918           if(box2)
 919             free(box2);
 920           if(velbox2)
 921             free(velbox2);
 922           return 9;
 923       }
 924     }
 925 #endif
 926   free(intvelbox);
 927   free(box1);
 928   free(velbox1);
 929   if(box2)
 930     free(box2);
 931   if(velbox2)
 932     free(velbox2);
 933   return 0;
 934 }
 935
 936 int main()
 937 {
 938   int testval;
 939   if (sizeof(int)<4)
 940     {
 941       fprintf(stderr,"ERROR: sizeof(int) is too small: %d<4\n",(int)sizeof(int));
 942       exit(EXIT_FAILURE);
 943     }
 944 #ifdef GEN
 945   printf("Tng compress testsuite generating (writing) test: %s\n",TESTNAME);
 946 #else
 947   printf("Tng compress testsuite running (reading) test: %s\n",TESTNAME);
 948 #endif
 949   testval=algotest();
 950   if (testval==0)
 951     printf("Passed.\n");
 952   else if (testval==1)
 953     {
 954       printf("ERROR: File error.\n");
 955       exit(EXIT_FAILURE);
 956     }
 957   else if (testval==4)
 958     {
 959       printf("ERROR: Read coding error in coordinates.\n");
 960       exit(EXIT_FAILURE);
 961     }
 962   else if (testval==5)
 963     {
 964       printf("ERROR: Read coding error in velocities.\n");
 965       exit(EXIT_FAILURE);
 966     }
 967   else if (testval==9)
 968     {
 969       printf("ERROR: Generated filesize differs too much.\n");
 970       exit(EXIT_FAILURE);
 971     }
 972   else
 973     {
 974       printf("ERROR: Unknown error.\n");
 975       exit(EXIT_FAILURE);
 976     }
 977   return 0;
 978 }