/* tng compression routines */ /* Only modify testsuite.c *Then* run testsuite.sh to perform the test. */ #include #include #include #include #include "compression/tng_compress.h" #include TESTPARAM #ifdef TEST_FLOAT #define REAL float #else #define REAL double #endif #ifndef TNG_COMPRESS_FILES_DIR #define TNG_COMPRESS_FILES_DIR "" #endif #define FUDGE 1.1 /* 10% off target precision is acceptable */ static void keepinbox(int *val) { while (val[0]>INTMAX1) val[0]-=(INTMAX1-INTMIN1+1); while (val[0]INTMAX2) val[1]-=(INTMAX2-INTMIN2+1); while (val[1]INTMAX3) val[2]-=(INTMAX3-INTMIN3+1); while (val[2]1)) ifl=1; else if ((i==1) && (flip) && (length>1)) ifl=0; target[ifl*3]=base[0]+(intsin((i+iframe)*direction[0])*scale)/256; target[ifl*3+1]=base[1]+(intcos((i+iframe)*direction[1])*scale)/256; target[ifl*3+2]=base[2]+(intcos((i+iframe)*direction[2])*scale)/256; keepinbox(target+ifl*3); } } #ifndef FRAMESCALE #define FRAMESCALE 1 #endif static void genibox(int *intbox, int iframe) { int molecule_length=1; int molpos[3]; int direction[3]={1,1,1}; int scale=1; int flip=0; int i=0; molpos[0]=intsin(iframe*FRAMESCALE)/32; molpos[1]=1+intcos(iframe*FRAMESCALE)/32; molpos[2]=2+intsin(iframe*FRAMESCALE)/16; keepinbox(molpos); while (iNATOMS) this_mol_length=NATOMS-i; /* We must test the large rle as well. This requires special sequencies to get triggered. So insert these from time to time */ #ifndef REGULAR if ((i%10)==0) { int j; intbox[i*3]=molpos[0]; intbox[i*3+1]=molpos[1]; intbox[i*3+2]=molpos[2]; for (j=1; j5) scale=1; molecule_length++; if (molecule_length>30) molecule_length=1; if (i%9) flip=1-flip; } } static void genivelbox(int *intvelbox, int iframe) { int i; for (i=0; imaxdiff) maxdiff=(REAL)fabs(arr1[i*stride1+j]-arr2[i*stride2+j]); } #if 0 for (i=0; iprec*0.5*FUDGE) { return 0; } else return 1; } struct tng_file { FILE *f; int natoms; int chunky; REAL precision; REAL velprecision; int initial_coding; int initial_coding_parameter; int coding; int coding_parameter; int initial_velcoding; int initial_velcoding_parameter; int velcoding; int velcoding_parameter; int speed; int nframes; int nframes_delivered; int writevel; REAL *pos; REAL *vel; int *ipos; int *ivel; unsigned int prec_hi, prec_lo; unsigned int velprec_hi, velprec_lo; }; static size_t fwrite_int_le(int *x,FILE *f) { unsigned char c[4]; unsigned int i=(unsigned int)*x; c[0]=(unsigned char)(i&0xFFU); c[1]=(unsigned char)((i>>8)&0xFFU); c[2]=(unsigned char)((i>>16)&0xFFU); c[3]=(unsigned char)((i>>24)&0xFFU); return fwrite(c,1,4,f); } static size_t fread_int_le(int *x,FILE *f) { unsigned char c[4]; unsigned int i; size_t n=fread(c,1,4,f); if (n) { i=(((unsigned int)c[3])<<24)|(((unsigned int)c[2])<<16)|(((unsigned int)c[1])<<8)|((unsigned int)c[0]); *x=(int)i; } return n; } static struct tng_file *open_tng_file_write(char *filename, int natoms,int chunky, REAL precision, int writevel, REAL velprecision, int initial_coding, int initial_coding_parameter, int coding, int coding_parameter, int initial_velcoding, int initial_velcoding_parameter, int velcoding, int velcoding_parameter, int speed) { struct tng_file *tng_file=malloc(sizeof *tng_file); tng_file->pos=NULL; tng_file->vel=NULL; tng_file->ipos=NULL; tng_file->ivel=NULL; tng_file->nframes=0; tng_file->chunky=chunky; tng_file->precision=precision; tng_file->natoms=natoms; tng_file->writevel=writevel; tng_file->velprecision=velprecision; tng_file->initial_coding=initial_coding; tng_file->initial_coding_parameter=initial_coding_parameter; tng_file->coding=coding; tng_file->coding_parameter=coding_parameter; tng_file->initial_velcoding=initial_velcoding; tng_file->initial_velcoding_parameter=initial_velcoding_parameter; tng_file->velcoding=velcoding; tng_file->velcoding_parameter=velcoding_parameter; tng_file->speed=speed; tng_file->pos=malloc(natoms*chunky*3*sizeof *tng_file->pos); tng_file->f=fopen(filename,"wb"); if (writevel) tng_file->vel=malloc(natoms*chunky*3*sizeof *tng_file->vel); fwrite_int_le(&natoms,tng_file->f); return tng_file; } static struct tng_file *open_tng_file_write_int(char *filename, int natoms,int chunky, int writevel, int initial_coding, int initial_coding_parameter, int coding, int coding_parameter, int initial_velcoding, int initial_velcoding_parameter, int velcoding, int velcoding_parameter, int speed) { struct tng_file *tng_file=malloc(sizeof *tng_file); tng_file->pos=NULL; tng_file->vel=NULL; tng_file->ipos=NULL; tng_file->ivel=NULL; tng_file->nframes=0; tng_file->chunky=chunky; tng_file->natoms=natoms; tng_file->writevel=writevel; tng_file->initial_coding=initial_coding; tng_file->initial_coding_parameter=initial_coding_parameter; tng_file->coding=coding; tng_file->coding_parameter=coding_parameter; tng_file->initial_velcoding=initial_velcoding; tng_file->initial_velcoding_parameter=initial_velcoding_parameter; tng_file->velcoding=velcoding; tng_file->velcoding_parameter=velcoding_parameter; tng_file->speed=speed; tng_file->ipos=malloc(natoms*chunky*3*sizeof *tng_file->ipos); tng_file->f=fopen(filename,"wb"); if (writevel) tng_file->ivel=malloc(natoms*chunky*3*sizeof *tng_file->ivel); fwrite_int_le(&natoms,tng_file->f); return tng_file; } static void flush_tng_frames(struct tng_file *tng_file, unsigned long prec_hi, unsigned long prec_lo, unsigned long velprec_hi, unsigned long velprec_lo) { int algo[4]; char *buf; int nitems; /* Make sure these variables are used to avoid compilation warnings */ (void)prec_hi; (void)prec_lo; (void)velprec_hi; (void)velprec_lo; fwrite_int_le(&tng_file->nframes,tng_file->f); algo[0]=tng_file->initial_coding; algo[1]=tng_file->initial_coding_parameter; algo[2]=tng_file->coding; algo[3]=tng_file->coding_parameter; #ifdef RECOMPRESS buf=tng_compress_pos_int(tng_file->ipos, tng_file->natoms, tng_file->nframes, prec_hi,prec_lo, tng_file->speed,algo,&nitems); #else /* RECOMPRESS */ #ifdef TEST_FLOAT buf=tng_compress_pos_float(tng_file->pos, tng_file->natoms, tng_file->nframes, tng_file->precision, tng_file->speed,algo,&nitems); #else /* TEST_FLOAT */ buf=tng_compress_pos(tng_file->pos, tng_file->natoms, tng_file->nframes, tng_file->precision, tng_file->speed,algo,&nitems); #endif /* TEST_FLOAT */ #endif /* RECOMPRESS */ tng_file->initial_coding=algo[0]; tng_file->initial_coding_parameter=algo[1]; tng_file->coding=algo[2]; tng_file->coding_parameter=algo[3]; fwrite_int_le(&nitems,tng_file->f); fwrite(buf,1,nitems,tng_file->f); free(buf); if (tng_file->writevel) { algo[0]=tng_file->initial_velcoding; algo[1]=tng_file->initial_velcoding_parameter; algo[2]=tng_file->velcoding; algo[3]=tng_file->velcoding_parameter; #ifdef RECOMPRESS buf=tng_compress_vel_int(tng_file->ivel, tng_file->natoms, tng_file->nframes, velprec_hi,velprec_lo, tng_file->speed,algo,&nitems); #else /* RECOMPRESS */ #ifdef TEST_FLOAT buf=tng_compress_vel_float(tng_file->vel, tng_file->natoms, tng_file->nframes, tng_file->velprecision, tng_file->speed,algo,&nitems); #else /* TEST_FLOAT */ buf=tng_compress_vel(tng_file->vel, tng_file->natoms, tng_file->nframes, tng_file->velprecision, tng_file->speed,algo,&nitems); #endif /* TEST_FLOAT */ #endif /* RECOMPRESS */ tng_file->initial_velcoding=algo[0]; tng_file->initial_velcoding_parameter=algo[1]; tng_file->velcoding=algo[2]; tng_file->velcoding_parameter=algo[3]; fwrite_int_le(&nitems,tng_file->f); fwrite(buf,1,nitems,tng_file->f); free(buf); } tng_file->nframes=0; } static void write_tng_file(struct tng_file *tng_file, REAL *pos,REAL *vel) { memcpy(tng_file->pos+tng_file->nframes*tng_file->natoms*3,pos,tng_file->natoms*3*sizeof *tng_file->pos); if (tng_file->writevel) memcpy(tng_file->vel+tng_file->nframes*tng_file->natoms*3,vel,tng_file->natoms*3*sizeof *tng_file->vel); tng_file->nframes++; if (tng_file->nframes==tng_file->chunky) flush_tng_frames(tng_file,0,0,0,0); } static void write_tng_file_int(struct tng_file *tng_file, int *ipos,int *ivel, unsigned long prec_hi, unsigned long prec_lo, unsigned long velprec_hi, unsigned long velprec_lo) { memcpy(tng_file->ipos+tng_file->nframes*tng_file->natoms*3,ipos,tng_file->natoms*3*sizeof *tng_file->ipos); if (tng_file->writevel) memcpy(tng_file->ivel+tng_file->nframes*tng_file->natoms*3,ivel,tng_file->natoms*3*sizeof *tng_file->ivel); tng_file->nframes++; if (tng_file->nframes==tng_file->chunky) flush_tng_frames(tng_file,prec_hi,prec_lo,velprec_hi,velprec_lo); tng_file->prec_hi=prec_hi; tng_file->prec_lo=prec_lo; tng_file->velprec_hi=velprec_hi; tng_file->velprec_lo=velprec_lo; } static void close_tng_file_write(struct tng_file *tng_file) { if (tng_file->nframes) flush_tng_frames(tng_file,tng_file->prec_hi,tng_file->prec_lo,tng_file->velprec_hi,tng_file->velprec_lo); fclose(tng_file->f); free(tng_file->pos); free(tng_file->vel); free(tng_file->ipos); free(tng_file->ivel); free(tng_file); } static struct tng_file *open_tng_file_read(char *filename, int writevel) { struct tng_file *tng_file=malloc(sizeof *tng_file); tng_file->pos=NULL; tng_file->vel=NULL; tng_file->ipos=NULL; tng_file->ivel=NULL; tng_file->f=fopen(filename,"rb"); tng_file->nframes=0; tng_file->nframes_delivered=0; tng_file->writevel=writevel; if (tng_file->f) fread_int_le(&tng_file->natoms,tng_file->f); else { free(tng_file); tng_file=NULL; } return tng_file; } static struct tng_file *open_tng_file_read_int(char *filename, int writevel) { struct tng_file *tng_file=malloc(sizeof *tng_file); tng_file->pos=NULL; tng_file->vel=NULL; tng_file->ipos=NULL; tng_file->ivel=NULL; tng_file->f=fopen(filename,"rb"); tng_file->nframes=0; tng_file->nframes_delivered=0; tng_file->writevel=writevel; if (tng_file->f) fread_int_le(&tng_file->natoms,tng_file->f); else { free(tng_file); tng_file=NULL; } return tng_file; } static int read_tng_file(struct tng_file *tng_file, REAL *pos, REAL *vel) { if (tng_file->nframes==tng_file->nframes_delivered) { int nitems; char *buf; free(tng_file->pos); free(tng_file->vel); if (!fread_int_le(&tng_file->nframes,tng_file->f)) return 1; if (!fread_int_le(&nitems,tng_file->f)) return 1; buf=malloc(nitems); if (!fread(buf,1,nitems,tng_file->f)) { free(buf); return 1; } tng_file->pos=malloc(tng_file->natoms*tng_file->nframes*3*sizeof *tng_file->pos); if (tng_file->writevel) tng_file->vel=malloc(tng_file->natoms*tng_file->nframes*3*sizeof *tng_file->vel); #if 0 { int natoms, nframes, algo[4]; double precision; int ivel; char *initial_coding, *coding; tng_compress_inquire(buf,&ivel,&natoms,&nframes,&precision,algo); initial_coding=tng_compress_initial_pos_algo(algo); coding=tng_compress_pos_algo(algo); printf("ivel=%d natoms=%d nframes=%d precision=%g initial pos=%s pos=%s\n",ivel,natoms,nframes,precision,initial_coding,coding); } #endif #ifdef TEST_FLOAT tng_compress_uncompress_float(buf,tng_file->pos); #else tng_compress_uncompress(buf,tng_file->pos); #endif free(buf); if (tng_file->writevel) { if (!fread_int_le(&nitems,tng_file->f)) return 1; buf=malloc(nitems); if (!fread(buf,1,nitems,tng_file->f)) { free(buf); return 1; } #if 0 { int natoms, nframes, algo[4]; double precision; int ivel; char *initial_coding, *coding; tng_compress_inquire(buf,&ivel,&natoms,&nframes,&precision,algo); initial_coding=tng_compress_initial_vel_algo(algo); coding=tng_compress_vel_algo(algo); printf("ivel=%d natoms=%d nframes=%d precision=%g initial vel=%s vel=%s\n",ivel,natoms,nframes,precision,initial_coding,coding); } #endif #ifdef TEST_FLOAT tng_compress_uncompress_float(buf,tng_file->vel); #else tng_compress_uncompress(buf,tng_file->vel); #endif free(buf); } tng_file->nframes_delivered=0; } memcpy(pos,tng_file->pos+tng_file->nframes_delivered*tng_file->natoms*3,tng_file->natoms*3*sizeof *pos); if (tng_file->writevel) memcpy(vel,tng_file->vel+tng_file->nframes_delivered*tng_file->natoms*3,tng_file->natoms*3*sizeof *vel); tng_file->nframes_delivered++; return 0; } static int read_tng_file_int(struct tng_file *tng_file, int *ipos, int *ivel, unsigned long *prec_hi, unsigned long *prec_lo, unsigned long *velprec_hi, unsigned long *velprec_lo) { if (tng_file->nframes==tng_file->nframes_delivered) { int nitems; char *buf; free(tng_file->ipos); free(tng_file->ivel); if (!fread_int_le(&tng_file->nframes,tng_file->f)) return 1; if (!fread_int_le(&nitems,tng_file->f)) return 1; buf=malloc(nitems); if (!fread(buf,1,nitems,tng_file->f)) { free(buf); return 1; } tng_file->ipos=malloc(tng_file->natoms*tng_file->nframes*3*sizeof *tng_file->ipos); if (tng_file->writevel) tng_file->ivel=malloc(tng_file->natoms*tng_file->nframes*3*sizeof *tng_file->ivel); tng_compress_uncompress_int(buf,tng_file->ipos,prec_hi,prec_lo); free(buf); if (tng_file->writevel) { if (!fread_int_le(&nitems,tng_file->f)) return 1; buf=malloc(nitems); if (!fread(buf,1,nitems,tng_file->f)) { free(buf); return 1; } tng_compress_uncompress_int(buf,tng_file->ivel,velprec_hi,velprec_lo); free(buf); } tng_file->nframes_delivered=0; } memcpy(ipos,tng_file->ipos+tng_file->nframes_delivered*tng_file->natoms*3,tng_file->natoms*3*sizeof *ipos); if (tng_file->writevel) memcpy(ivel,tng_file->ivel+tng_file->nframes_delivered*tng_file->natoms*3,tng_file->natoms*3*sizeof *ivel); tng_file->nframes_delivered++; return 0; } static void close_tng_file_read(struct tng_file *tng_file) { free(tng_file->vel); free(tng_file->pos); free(tng_file->ivel); free(tng_file->ipos); fclose(tng_file->f); free(tng_file); } #ifndef EXPECTED_FILESIZE #define EXPECTED_FILESIZE 1 #endif #ifndef INITIALVELCODING #define INITIALVELCODING -1 #endif #ifndef INITIALVELCODINGPARAMETER #define INITIALVELCODINGPARAMETER -1 #endif #ifndef SPEED #define SPEED 5 #endif /* Return value 1 means file error. Return value 4 means coding error in coordinates. Return value 5 means coding error in velocities. Return value 9 means filesize seems too off. Return value 100+ means test specific error. */ static int algotest() { int i; int *intbox=malloc(NATOMS*3*sizeof *intbox); int *intvelbox=malloc(NATOMS*3*sizeof *intvelbox); #ifdef RECOMPRESS unsigned long pos_prec_hi,pos_prec_lo; unsigned long vel_prec_hi,vel_prec_lo; #endif REAL *box1=malloc(NATOMS*STRIDE1*sizeof *box1); REAL *velbox1=malloc(NATOMS*STRIDE1*sizeof *velbox1); int startframe=0; int endframe=NFRAMES; #ifdef GEN FILE *file; REAL filesize; REAL *box2=0; REAL *velbox2=0; #else int readreturn; REAL *box2=malloc(NATOMS*STRIDE2*sizeof *box2); REAL *velbox2=malloc(NATOMS*STRIDE2*sizeof *velbox2); #endif #ifdef RECOMPRESS void *dumpfile=open_tng_file_write_int(TNG_COMPRESS_FILES_DIR FILENAME,NATOMS,CHUNKY, WRITEVEL, INITIALCODING, INITIALCODINGPARAMETER,CODING,CODINGPARAMETER, INITIALVELCODING,INITIALVELCODINGPARAMETER, VELCODING,VELCODINGPARAMETER,SPEED); void *dumpfile_recompress=open_tng_file_read_int(TNG_COMPRESS_FILES_DIR RECOMPRESS,WRITEVEL); if (!dumpfile_recompress) { free(intbox); free(intvelbox); free(box1); free(velbox1); if(box2) free(box2); if(velbox2) free(velbox2); return 1; } #else /* RECOMPRESS */ #ifdef GEN void *dumpfile=open_tng_file_write(TNG_COMPRESS_FILES_DIR FILENAME,NATOMS,CHUNKY, PRECISION,WRITEVEL,VELPRECISION, INITIALCODING, INITIALCODINGPARAMETER,CODING,CODINGPARAMETER, INITIALVELCODING,INITIALVELCODINGPARAMETER, VELCODING,VELCODINGPARAMETER,SPEED); #else void *dumpfile=open_tng_file_read(TNG_COMPRESS_FILES_DIR FILENAME,WRITEVEL); #endif #endif /* RECOMPRESS */ if (!dumpfile) { free(intbox); free(intvelbox); free(box1); free(velbox1); if(box2) free(box2); if(velbox2) free(velbox2); return 1; } for (i=startframe; i0) { if ((fabs(filesize-EXPECTED_FILESIZE)/EXPECTED_FILESIZE)>0.05) { free(intvelbox); free(box1); free(velbox1); if(box2) free(box2); if(velbox2) free(velbox2); return 9; } } #endif free(intvelbox); free(box1); free(velbox1); if(box2) free(box2); if(velbox2) free(velbox2); return 0; } int main() { int testval; if (sizeof(int)<4) { fprintf(stderr,"ERROR: sizeof(int) is too small: %d<4\n",(int)sizeof(int)); exit(EXIT_FAILURE); } #ifdef GEN printf("Tng compress testsuite generating (writing) test: %s\n",TESTNAME); #else printf("Tng compress testsuite running (reading) test: %s\n",TESTNAME); #endif testval=algotest(); if (testval==0) printf("Passed.\n"); else if (testval==1) { printf("ERROR: File error.\n"); exit(EXIT_FAILURE); } else if (testval==4) { printf("ERROR: Read coding error in coordinates.\n"); exit(EXIT_FAILURE); } else if (testval==5) { printf("ERROR: Read coding error in velocities.\n"); exit(EXIT_FAILURE); } else if (testval==9) { printf("ERROR: Generated filesize differs too much.\n"); exit(EXIT_FAILURE); } else { printf("ERROR: Unknown error.\n"); exit(EXIT_FAILURE); } return 0; }