diff options
| -rw-r--r-- | src/compression/coder.c | 36 | ||||
| -rw-r--r-- | src/lib/tng_io.c | 253 | ||||
| -rw-r--r-- | src/tests/tng_io_read_pos_util.c | 18 |
3 files changed, 193 insertions, 114 deletions
diff --git a/src/compression/coder.c b/src/compression/coder.c index a64e741..a9a9117 100644 --- a/src/compression/coder.c +++ b/src/compression/coder.c @@ -18,6 +18,20 @@ #include "coder.h" #include "warnmalloc.h" +#ifndef USE_WINDOWS +#if defined(WIN32) || defined(_WIN32) || defined(WIN64) || defined(_WIN64) +#define USE_WINDOWS +#endif /* win32... */ +#endif /* not defined USE_WINDOWS */ + +#ifdef USE_WINDOWS +#define TNG_INLINE __inline +#define TNG_SNPRINTF _snprintf +#else +#define TNG_INLINE inline +#define TNG_SNPRINTF snprintf +#endif + struct coder *Ptngc_coder_init(void) { struct coder *coder=warnmalloc(sizeof *coder); @@ -30,17 +44,21 @@ void Ptngc_coder_deinit(struct coder *coder) free(coder); } -void Ptngc_out8bits(struct coder *coder, unsigned char **output) +TNG_INLINE void Ptngc_out8bits(struct coder *coder, unsigned char **output) { - while (coder->pack_temporary_bits>=8) + int pack_temporary_bits=coder->pack_temporary_bits; + unsigned int pack_temporary=coder->pack_temporary; + while (pack_temporary_bits>=8) { - unsigned int mask=~(0xFFU<<(coder->pack_temporary_bits-8)); - unsigned char out=(unsigned char)(coder->pack_temporary>>(coder->pack_temporary_bits-8)); + unsigned int mask=~(0xFFU<<(pack_temporary_bits-8)); + unsigned char out=(unsigned char)(pack_temporary>>(pack_temporary_bits-8)); **output=out; (*output)++; - coder->pack_temporary_bits-=8; - coder->pack_temporary&=mask; - } + pack_temporary_bits-=8; + pack_temporary&=mask; + } + coder->pack_temporary_bits=pack_temporary_bits; + coder->pack_temporary=pack_temporary; } void Ptngc_write_pattern(struct coder *coder,unsigned int pattern, int nbits, unsigned char **output) @@ -62,7 +80,7 @@ void Ptngc_write_pattern(struct coder *coder,unsigned int pattern, int nbits, un } /* Write up to 24 bits */ -void Ptngc_writebits(struct coder *coder,unsigned int value,int nbits, unsigned char **output_ptr) +TNG_INLINE void Ptngc_writebits(struct coder *coder,unsigned int value,int nbits, unsigned char **output_ptr) { /* Make room for the bits. */ coder->pack_temporary<<=nbits; @@ -360,7 +378,7 @@ static int unpack_array_stop_bits(struct coder *coder,unsigned char *packed,int if ((pattern%2)==0) s=-s; output[i]=s; - } + } return 0; } diff --git a/src/lib/tng_io.c b/src/lib/tng_io.c index a931ed6..28bb9e6 100644 --- a/src/lib/tng_io.c +++ b/src/lib/tng_io.c @@ -3766,6 +3766,7 @@ static tng_function_status tng_compress(tng_trajectory_t tng_data, int nalgo; int new_len; char *dest, *temp; + int algo_find_n_frames; if(block->id != TNG_TRAJ_POSITIONS && block->id != TNG_TRAJ_VELOCITIES) @@ -3780,45 +3781,54 @@ static tng_function_status tng_compress(tng_trajectory_t tng_data, return(TNG_FAILURE); } -// #if 1 -// { -// int iframe,i,j; -// for (iframe=0; iframe<n_frames; iframe++) -// for (i=0; i<n_particles; i++) -// { -// printf("tng_io: %d %d: ",iframe,i); -// for (j=0; j<3; j++) -// printf(" %g",((float*)start_pos)[iframe*n_particles*3+i*3+j]); -// printf("\n"); -// } -// } -// #endif - if(block->id == TNG_TRAJ_POSITIONS) { if(!tng_data->compress_algo_pos) { + if(n_frames > 10) + { + algo_find_n_frames = 5; + } + else + { + algo_find_n_frames = n_frames; + } + nalgo = tng_compress_nalgo(); tng_data->compress_algo_pos=malloc(nalgo * sizeof *tng_data->compress_algo_pos); - if(type == TNG_FLOAT_DATA) { dest = tng_compress_pos_float_find_algo(start_pos, n_particles, - n_frames, + algo_find_n_frames, 0.001, 0, tng_data-> compress_algo_pos, &new_len); + + if(algo_find_n_frames < n_frames) + { + dest = tng_compress_pos_float(start_pos, n_particles, + n_frames, 0.001, 0, + tng_data->compress_algo_pos, + &new_len); + } } else { dest = tng_compress_pos_find_algo(start_pos, n_particles, - n_frames, + algo_find_n_frames, 0.001, 0, tng_data-> compress_algo_pos, &new_len); + if(algo_find_n_frames < n_frames) + { + dest = tng_compress_pos(start_pos, n_particles, + n_frames, 0.001, 0, + tng_data->compress_algo_pos, + &new_len); + } } } else @@ -3841,6 +3851,15 @@ static tng_function_status tng_compress(tng_trajectory_t tng_data, { if(!tng_data->compress_algo_vel) { + if(n_frames > 10) + { + algo_find_n_frames = 5; + } + else + { + algo_find_n_frames = n_frames; + } + nalgo = tng_compress_nalgo(); tng_data->compress_algo_vel=malloc(nalgo * sizeof *tng_data->compress_algo_vel); @@ -3848,20 +3867,34 @@ static tng_function_status tng_compress(tng_trajectory_t tng_data, if(type == TNG_FLOAT_DATA) { dest = tng_compress_vel_float_find_algo(start_pos, n_particles, - n_frames, + algo_find_n_frames, 0.001, 0, tng_data-> compress_algo_vel, &new_len); + if(algo_find_n_frames < n_frames) + { + dest = tng_compress_vel_float(start_pos, n_particles, + n_frames, 0.001, 0, + tng_data->compress_algo_pos, + &new_len); + } } else { dest = tng_compress_vel_find_algo(start_pos, n_particles, - n_frames, + algo_find_n_frames, 0.001, 0, tng_data-> compress_algo_vel, &new_len); + if(algo_find_n_frames < n_frames) + { + dest = tng_compress_vel(start_pos, n_particles, + n_frames, 0.001, 0, + tng_data->compress_algo_pos, + &new_len); + } } } else @@ -4165,7 +4198,7 @@ static tng_function_status tng_allocate_particle_data_mem n_frames = tng_max(1, n_frames); data->stride_length = tng_max(1, stride_length); data->n_values_per_frame = n_values_per_frame; - frame_alloc = tng_max(1, (n_frames % stride_length) ? n_frames / stride_length + 1 : n_frames / stride_length); + frame_alloc = (n_frames % stride_length) ? n_frames / stride_length + 1 : n_frames / stride_length; if(data->datatype == TNG_CHAR_DATA) { @@ -4474,7 +4507,7 @@ static tng_function_status tng_particle_data_read tot_n_particles = tng_data->n_particles; } - n_frames_div = tng_max(1, (n_frames % stride_length) ? n_frames / stride_length + 1 : n_frames / stride_length); + n_frames_div = (n_frames % stride_length) ? n_frames / stride_length + 1 : n_frames / stride_length; if(codec_id != TNG_UNCOMPRESSED) { @@ -4706,11 +4739,8 @@ static tng_function_status tng_particle_data_block_write n_frames = tng_min(n_frames, frame_set->n_frames); } - frame_step = n_frames / stride_length; - if(n_frames % stride_length == 1) - { - frame_step += 1; - } + frame_step = (n_frames % stride_length) ? n_frames / stride_length + 1: + n_frames / stride_length; if(mapping && mapping->n_particles != 0) { @@ -5206,7 +5236,7 @@ static tng_function_status tng_allocate_data_mem data->stride_length = tng_max(1, stride_length); n_frames = tng_max(1, n_frames); data->n_values_per_frame = n_values_per_frame; - frame_alloc = tng_max(1, (n_frames % stride_length) ? n_frames / stride_length + 1 : n_frames / stride_length); + frame_alloc = (n_frames % stride_length) ? n_frames / stride_length + 1 : n_frames / stride_length; if(data->datatype == TNG_CHAR_DATA) { @@ -5366,7 +5396,7 @@ static tng_function_status tng_data_read(tng_trajectory_t tng_data, data->compression_multiplier = multiplier; } - n_frames_div = tng_max(1, (n_frames % stride_length) ? n_frames / stride_length + 1 : n_frames / stride_length); + n_frames_div = (n_frames % stride_length) ? n_frames / stride_length + 1 : n_frames / stride_length; if(codec_id != TNG_UNCOMPRESSED) { @@ -5572,11 +5602,8 @@ static tng_function_status tng_data_block_write(tng_trajectory_t tng_data, n_frames = tng_min(n_frames, frame_set->n_frames); } - frame_step = n_frames / stride_length; - if(n_frames % stride_length == 1) - { - frame_step += 1; - } + frame_step = (n_frames % stride_length) ? n_frames / stride_length + 1: + n_frames / stride_length; block->block_contents_size = sizeof(char) * 2 + sizeof(data->n_values_per_frame) + @@ -10666,6 +10693,7 @@ tng_function_status DECLSPECDLLEXPORT tng_data_block_add tng_non_particle_data_t data; char **first_dim_values; char *new_data_c=new_data; + int64_t n_frames_div; frame_set = &tng_data->current_trajectory_frame_set; @@ -10746,9 +10774,13 @@ tng_function_status DECLSPECDLLEXPORT tng_data_block_add frame_set->n_written_frames = n_frames; } + n_frames_div = (n_frames % stride_length) ? + n_frames / stride_length + 1: + n_frames / stride_length; + if(datatype == TNG_CHAR_DATA) { - for(i = 0; i < n_frames / stride_length; i++) + for(i = 0; i < n_frames_div; i++) { first_dim_values = data->strings[i]; for(j = 0; j < n_values_per_frame; j++) @@ -10774,7 +10806,7 @@ tng_function_status DECLSPECDLLEXPORT tng_data_block_add } else { - memcpy(data->values, new_data, size * n_frames / stride_length * + memcpy(data->values, new_data, size * n_frames_div * n_values_per_frame); } } @@ -10797,7 +10829,7 @@ tng_function_status DECLSPECDLLEXPORT tng_particle_data_block_add void *new_data) { int i, j, k, size, len; - int64_t tot_n_particles; + int64_t tot_n_particles, n_frames_div; char ***first_dim_values, **second_dim_values; tng_trajectory_frame_set_t frame_set; tng_particle_data_t data; @@ -10882,9 +10914,13 @@ tng_function_status DECLSPECDLLEXPORT tng_particle_data_block_add frame_set->n_written_frames = n_frames; } + n_frames_div = (n_frames % stride_length) ? + n_frames / stride_length + 1: + n_frames / stride_length; + if(datatype == TNG_CHAR_DATA) { - for(i = 0; i < n_frames / stride_length; i++) + for(i = 0; i < n_frames_div; i++) { first_dim_values = data->strings[i]; for(j = num_first_particle; j < num_first_particle + n_particles; @@ -10928,7 +10964,7 @@ tng_function_status DECLSPECDLLEXPORT tng_particle_data_block_add size = sizeof(double); } - memcpy(data->values, new_data, size * n_frames / stride_length * + memcpy(data->values, new_data, size * n_frames_div * n_particles * n_values_per_frame); } } @@ -11855,7 +11891,7 @@ tng_function_status DECLSPECDLLEXPORT tng_frame_particle_data_write else { fwrite(values, val_n_particles * n_values_per_frame, size, - tng_data->output_file); + tng_data->output_file); } fflush(tng_data->output_file); @@ -12212,7 +12248,7 @@ tng_function_status tng_data_vector_get(tng_trajectory_t tng_data, int64_t *n_values_per_frame, char *type) { - int64_t file_pos, data_size, frames_div; + int64_t file_pos, data_size, n_frames_div; int i, block_index, size; tng_non_particle_data_t data; tng_trajectory_frame_set_t frame_set = @@ -12286,12 +12322,13 @@ tng_function_status tng_data_vector_get(tng_trajectory_t tng_data, size = sizeof(double); } - *n_frames = tng_max(1, data->n_frames); + *n_frames = data->n_frames; *n_values_per_frame = data->n_values_per_frame; *stride_length = data->stride_length; - frames_div = (*n_frames % *stride_length) ? *n_frames / *stride_length + 1 : *n_frames / *stride_length; + n_frames_div = (*n_frames % *stride_length) ? *n_frames / *stride_length + 1: + *n_frames / *stride_length; - data_size = frames_div * size * + data_size = n_frames_div * size * *n_values_per_frame; temp = realloc(*values, data_size); @@ -12507,7 +12544,7 @@ tng_function_status DECLSPECDLLEXPORT tng_data_vector_interval_get int64_t *n_values_per_frame, char *type) { - int64_t n_frames, tot_n_frames, frames_div, frames_div_2; + int64_t n_frames, tot_n_frames, n_frames_div, n_frames_div_2; int64_t current_frame_pos, data_size, frame_size; int64_t last_frame_pos; int size; @@ -12561,10 +12598,10 @@ tng_function_status DECLSPECDLLEXPORT tng_data_vector_interval_get size = sizeof(double); } - frames_div = (tot_n_frames % *stride_length) ? - tot_n_frames / *stride_length + 1 : tot_n_frames / *stride_length; - data_size = frames_div * size * - (*n_values_per_frame); + n_frames_div = (tot_n_frames % *stride_length) ? + tot_n_frames / *stride_length + 1: + tot_n_frames / *stride_length; + data_size = n_frames_div * size * (*n_values_per_frame); temp = realloc(*values, data_size); if(!temp) @@ -12587,18 +12624,17 @@ tng_function_status DECLSPECDLLEXPORT tng_data_vector_interval_get frame_size = size * (*n_values_per_frame); - last_frame_pos = tng_min(frame_set->first_frame + n_frames - - *stride_length, - end_frame_nr - frame_set->first_frame); + last_frame_pos = tng_min(n_frames, + end_frame_nr - current_frame_pos); + + n_frames_div = current_frame_pos / *stride_length; + n_frames_div_2 = (last_frame_pos % *stride_length) ? + last_frame_pos / *stride_length + 1: + last_frame_pos / *stride_length; + n_frames_div_2 = tng_max(1, n_frames_div_2); - frames_div = (current_frame_pos % *stride_length) ? - current_frame_pos / *stride_length + 1 : - current_frame_pos / *stride_length; - frames_div_2 = ((last_frame_pos - current_frame_pos) % *stride_length) ? - (last_frame_pos - current_frame_pos) / *stride_length + 1 : - (last_frame_pos - current_frame_pos) / *stride_length; - memcpy(*values, current_values + frames_div * frame_size, - (frames_div_2 + 1) * frame_size); + memcpy(*values, current_values + n_frames_div * frame_size, + n_frames_div_2 * frame_size); current_frame_pos += n_frames - frame_set->first_frame - current_frame_pos; @@ -12626,20 +12662,18 @@ tng_function_status DECLSPECDLLEXPORT tng_data_vector_interval_get return(stat); } - last_frame_pos = tng_min(frame_set->first_frame + n_frames - - *stride_length, - end_frame_nr - frame_set->first_frame); + last_frame_pos = tng_min(n_frames, + end_frame_nr - current_frame_pos); - frames_div = ((current_frame_pos - start_frame_nr) % *stride_length) ? - (current_frame_pos - start_frame_nr) / *stride_length +1 : - (current_frame_pos - start_frame_nr) / *stride_length; - frames_div_2 = ((last_frame_pos - current_frame_pos) % *stride_length) ? - (last_frame_pos - current_frame_pos) / *stride_length + 1 : - (last_frame_pos - current_frame_pos) / *stride_length; + n_frames_div = current_frame_pos / *stride_length; + n_frames_div_2 = (last_frame_pos % *stride_length) ? + last_frame_pos / *stride_length + 1: + last_frame_pos / *stride_length; + n_frames_div_2 = tng_max(1, n_frames_div_2); - memcpy(*values + frames_div * frame_size, + memcpy(*values + n_frames_div * frame_size, current_values, - (frames_div_2 + 1) * frame_size); + n_frames_div_2 * frame_size); current_frame_pos += frame_set->n_frames; } @@ -12853,7 +12887,7 @@ tng_function_status DECLSPECDLLEXPORT tng_particle_data_vector_get int64_t *n_values_per_frame, char *type) { - int64_t i, j, mapping, file_pos, i_step, data_size, frames_div; + int64_t i, j, mapping, file_pos, i_step, data_size, n_frames_div; int block_index, size; tng_particle_data_t data; tng_trajectory_frame_set_t frame_set = @@ -12952,10 +12986,11 @@ tng_function_status DECLSPECDLLEXPORT tng_particle_data_vector_get *n_values_per_frame = data->n_values_per_frame; *stride_length = data->stride_length; - frames_div = (*n_frames % *stride_length) ? *n_frames / *stride_length + 1 : - *n_frames / *stride_length; + n_frames_div = (*n_frames % *stride_length) ? + *n_frames / *stride_length + 1: + *n_frames / *stride_length; - data_size = frames_div * size * (*n_particles) * + data_size = n_frames_div * size * (*n_particles) * (*n_values_per_frame); temp = realloc(*values, data_size); @@ -13020,7 +13055,6 @@ tng_function_status DECLSPECDLLEXPORT tng_particle_data_interval_get return(stat); } - tng_block_init(&block); file_pos = ftell(tng_data->input_file); /* Read all blocks until next frame set block */ @@ -13230,10 +13264,11 @@ tng_function_status DECLSPECDLLEXPORT tng_particle_data_vector_interval_get int64_t *n_values_per_frame, char *type) { - int64_t n_frames, tot_n_frames, frames_div, frames_div_2; - int64_t current_frame_pos, last_frame_pos, data_size, frame_size; + int64_t n_frames, tot_n_frames, n_frames_div, n_frames_div_2; + int64_t file_pos, current_frame_pos, last_frame_pos, data_size, frame_size; int size; tng_trajectory_frame_set_t frame_set; + tng_gen_block_t block; tng_function_status stat; void *current_values = 0, *temp; @@ -13243,6 +13278,33 @@ tng_function_status DECLSPECDLLEXPORT tng_particle_data_vector_interval_get return(stat); } + tng_block_init(&block); + file_pos = ftell(tng_data->input_file); + /* Read all blocks until next frame set block */ + stat = tng_block_header_read(tng_data, block); + while(file_pos < tng_data->input_file_len && + stat != TNG_CRITICAL && + block->id != TNG_TRAJECTORY_FRAME_SET) + { + stat = tng_block_read_next(tng_data, block, + hash_mode); + if(stat != TNG_CRITICAL) + { + file_pos = ftell(tng_data->input_file); + if(file_pos < tng_data->input_file_len) + { + stat = tng_block_header_read(tng_data, block); + } + } + } + tng_block_destroy(&block); + if(stat == TNG_CRITICAL) + { + printf("Cannot read block header at pos %"PRId64". %s: %d\n", + file_pos, __FILE__, __LINE__); + return(stat); + } + frame_set = &tng_data->current_trajectory_frame_set; stat = tng_particle_data_vector_get(tng_data, block_id, ¤t_values, @@ -13264,7 +13326,7 @@ tng_function_status DECLSPECDLLEXPORT tng_particle_data_vector_interval_get } else { - tot_n_frames = end_frame_nr - start_frame_nr + 1 *(*stride_length); + tot_n_frames = end_frame_nr - start_frame_nr + 1; } switch(*type) @@ -13282,10 +13344,11 @@ tng_function_status DECLSPECDLLEXPORT tng_particle_data_vector_interval_get size = sizeof(double); } - frames_div = (tot_n_frames % *stride_length) ? tot_n_frames / *stride_length + 1: + n_frames_div = (tot_n_frames % *stride_length) ? + tot_n_frames / *stride_length + 1: tot_n_frames / *stride_length; - data_size = frames_div * size * (*n_particles) * + data_size = n_frames_div * size * (*n_particles) * (*n_values_per_frame); temp = realloc(*values, data_size); @@ -13311,17 +13374,16 @@ tng_function_status DECLSPECDLLEXPORT tng_particle_data_vector_interval_get frame_size = size * (*n_particles) * (*n_values_per_frame); last_frame_pos = tng_min(n_frames, - end_frame_nr - current_frame_pos + 1); + end_frame_nr - current_frame_pos); - frames_div = (current_frame_pos % *stride_length) ? - current_frame_pos / *stride_length + 1 : - current_frame_pos / *stride_length; - frames_div_2 = (last_frame_pos % *stride_length) ? - last_frame_pos / *stride_length + 1 : + n_frames_div = current_frame_pos / *stride_length; + n_frames_div_2 = (last_frame_pos % *stride_length) ? + last_frame_pos / *stride_length + 1: last_frame_pos / *stride_length; + n_frames_div_2 = tng_max(1, n_frames_div_2); - memcpy(*values, current_values + frames_div * frame_size, - frames_div_2 * frame_size); + memcpy(*values, current_values + n_frames_div * frame_size, + n_frames_div_2 * frame_size); current_frame_pos += n_frames - frame_set->first_frame - current_frame_pos; @@ -13350,18 +13412,17 @@ tng_function_status DECLSPECDLLEXPORT tng_particle_data_vector_interval_get } last_frame_pos = tng_min(n_frames, - end_frame_nr - current_frame_pos + 1); + end_frame_nr - current_frame_pos); - frames_div = ((current_frame_pos - start_frame_nr) % *stride_length) ? - (current_frame_pos - start_frame_nr) / *stride_length + 1 : - (current_frame_pos - start_frame_nr) / *stride_length; - frames_div_2 = (last_frame_pos % *stride_length) ? - last_frame_pos / *stride_length + 1 : + n_frames_div = current_frame_pos / *stride_length; + n_frames_div_2 = (last_frame_pos % *stride_length) ? + last_frame_pos / *stride_length + 1: last_frame_pos / *stride_length; + n_frames_div_2 = tng_max(1, n_frames_div_2); - memcpy(*values + frames_div * frame_size, + memcpy(*values + n_frames_div * frame_size, current_values, - frames_div_2 * frame_size); + n_frames_div_2 * frame_size); current_frame_pos += n_frames; } diff --git a/src/tests/tng_io_read_pos_util.c b/src/tests/tng_io_read_pos_util.c index 22e75d5..a80ae53 100644 --- a/src/tests/tng_io_read_pos_util.c +++ b/src/tests/tng_io_read_pos_util.c @@ -27,8 +27,7 @@ int main(int argc, char **argv) float *positions = 0, *box_shape = 0; int64_t n_particles, n_frames, tot_n_frames, stride_length, i, j; // Set a default frame range - int64_t first_frame = 0, last_frame = 5000; - int k; + int64_t first_frame = 0, last_frame = 5000, n_strides; if(argc <= 1) { @@ -91,6 +90,9 @@ int main(int argc, char **argv) printf("Simulation box shape not set in the file (or could not be read)\n"); } + n_strides = (n_frames % stride_length) ? n_frames / stride_length + 1: + n_frames / stride_length; + // Get the positions of all particles in the requested frame range. // The positions are stored in the positions array. // N.B. No proper error checks. @@ -98,17 +100,15 @@ int main(int argc, char **argv) == TNG_SUCCESS) { // Print the positions of the wanted particle (zero based) - for(i=0; i < n_frames; i += stride_length) + for(i=0; i < n_strides; i++) { - printf("\nFrame %"PRId64":\n", first_frame + i); + printf("\nFrame %"PRId64":\n", first_frame + i*stride_length); for(j=0; j < n_particles; j++) { printf("Atom nr: %"PRId64"", j); - for(k=0; k < 3; k++) - { - printf("\t%f", positions[i/stride_length*n_particles* - 3+j*3+k]); - } + printf("\t%f", positions[i*n_particles*3+j*3]); + printf("\t%f", positions[i*n_particles*3+j*3+1]); + printf("\t%f", positions[i*n_particles*3+j*3+2]); printf("\n"); } } |