Last active
February 10, 2021 17:54
-
-
Save data-panda/93d3ad48e098ade08f4154b016527dd3 to your computer and use it in GitHub Desktop.
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
/* minimal code example showing how to call the zfp (de)compressor */ | |
#ifdef DOUBLE_DT | |
typedef double md; | |
#else | |
typedef float md; | |
#endif | |
#include <math.h> | |
#include <stdio.h> | |
#include <stdlib.h> | |
#include <string.h> | |
#include <zfp.h> | |
/* fill grid with a floating point value */ | |
static inline void | |
fillGridWithVal(md *array, int nx, int ny, int nz, md val) | |
{ | |
#pragma omp parallel for | |
for (int k = 0; k < nz; k++) | |
for (int j = 0; j < ny; j++) | |
for (int i = 0; i < nx; i++) { | |
array[ i + nx * (j + ny * k) ] = val; | |
} | |
} | |
/* fill the domain with cosine wave */ | |
/* Fill the 3d domain with cosine wave */ | |
static inline md | |
input_func_cos(int i, int j, int k, int nx, int ny, int nz) | |
{ | |
md delta = 1.; | |
md x = i * delta; | |
md y = j * delta; | |
md z = k * delta; | |
// source of seismic wave at top | |
md xs = (nx * delta)/2.; | |
md ys = (ny * delta)/2.; | |
md zs = 0; | |
md wl = 348.23; | |
md amp = 25000; | |
x -= xs; | |
y -= ys; | |
z -= zs; | |
md r = sqrtf( x * x + y * y + z * z); | |
#ifdef DOUBLE_DT | |
md val = amp * fmin(1., 1./r) * cos(r / wl * 2.* M_PI); | |
#else | |
md val = amp * fminf(1., 1./r) * cos(r / wl * 2.* M_PI); | |
#endif | |
// printf("%e, %e, %e\n",r, fminf(1., 1./r), cos(r / wl * 2.* M_PI)); | |
return val; | |
} | |
/* array visualization functions */ | |
static void | |
force_big_endian(unsigned char *bytes) | |
{ | |
static int doneTest = 0; | |
static int shouldSwap = 0; | |
if(!doneTest) | |
{ | |
int tmp1 = 1; | |
unsigned char *tmp2 = (unsigned char *) &tmp1; | |
if (*tmp2 != 0) | |
shouldSwap = 1; | |
doneTest = 1; | |
} | |
if (shouldSwap) | |
{ | |
unsigned char tmp = bytes[0]; | |
bytes[0] = bytes[3]; | |
bytes[3] = tmp; | |
tmp = bytes[1]; | |
bytes[1] = bytes[2]; | |
bytes[2] = tmp; | |
} | |
} | |
void | |
writeStructuredPoints(int nx, int ny, int nz, float *data, int binary, const char *filename) | |
{ | |
long int nt = nx * ny * nz; | |
char full_filename[1024]; | |
if (strstr(filename, ".vtk") != NULL) | |
{ | |
strcpy(full_filename, filename); | |
} | |
else | |
{ | |
sprintf(full_filename, "%s.vtk", filename); | |
} | |
FILE *fp = fopen(full_filename, "w+"); | |
/* write header */ | |
fprintf(fp, "# vtk DataFile Version 2.0\n"); | |
fprintf(fp, "Volume data from array\n"); | |
if (binary) fprintf(fp, "BINARY\n"); | |
else fprintf(fp, "ASCII\n"); | |
fprintf(fp, "DATASET STRUCTURED_POINTS\n"); | |
fprintf(fp, "DIMENSIONS %d %d %d\n", nx, ny, nz); | |
fprintf(fp, "SPACING %f %f %f\n", 1., 1., 1.); | |
fprintf(fp, "ORIGIN %f %f %f\n", 0., 0., 0.); | |
fprintf(fp, "POINT_DATA %ld\n", nt); | |
fprintf(fp, "SCALARS value float 1\n"); //hardcode variable name as value | |
fprintf(fp, "LOOKUP_TABLE default\n"); | |
/* write the values */ | |
if(binary) // ignore big endian and little endian check | |
{ | |
for(int i = 0; i < nt; i++) | |
{ | |
float val = data[i]; | |
force_big_endian((unsigned char *) &val); | |
fwrite(&val, sizeof(float), 1, fp); | |
} | |
fprintf(fp, "\n"); | |
} | |
else | |
{ | |
int numInColumn = 0; | |
int values = 0; | |
for(long int i = 0; i < nt; i++) | |
{ | |
fprintf(fp, "%f\n", data[i]); values++; | |
} | |
printf("array dim is %d and value written is %d\n", nt, values); | |
} | |
fclose(fp); | |
} | |
/* compress or decompress array */ | |
static int | |
compress(md* array, int nx, int ny, int nz, double bpv, int decompress, size_t *cmp_size) | |
{ | |
int status = 0; /* return value: 0 = success */ | |
zfp_type type; /* array scalar type */ | |
zfp_field* field; /* array meta data */ | |
zfp_stream* zfp; /* compressed stream */ | |
void* buffer; /* storage for compressed stream */ | |
size_t bufsize; /* byte size of compressed buffer */ | |
bitstream* stream; /* bit stream to write to or read from */ | |
size_t zfpsize; /* byte size of compressed stream */ | |
/* allocate meta data for the 3D array a[nz][ny][nx] */ | |
#ifdef DOUBLE_DT | |
type = zfp_type_double; | |
#else | |
type = zfp_type_float; | |
#endif | |
field = zfp_field_3d(array, type, nx, ny, nz); | |
/* allocate meta data for a compressed stream */ | |
zfp = zfp_stream_open(NULL); | |
/* set compression mode and parameters via one of four functions */ | |
zfp_stream_set_rate(zfp, bpv, type, zfp_field_dimensionality(field), 0); | |
//zfp_stream_set_accuracy(zfp, 1e-3); | |
zfp_stream_set_execution(zfp, zfp_exec_serial); | |
/* allocate buffer for compressed data */ | |
bufsize = zfp_stream_maximum_size(zfp, field); | |
buffer = malloc(bufsize); | |
/* associate bit stream with allocated buffer */ | |
stream = stream_open(buffer, bufsize); | |
zfp_stream_set_bit_stream(zfp, stream); | |
zfp_stream_rewind(zfp); | |
/* compress or decompress entire array */ | |
if (decompress) { | |
zfp_stream_set_execution(zfp, zfp_exec_serial); | |
zfp_field_set_pointer(field, (void *)array); | |
/* read compressed stream and decompress and output array */ | |
FILE *fp = fopen("compressed_file", "rb+"); | |
zfpsize = fread(buffer, 1, *cmp_size, fp); | |
fclose(fp); | |
zfpsize = zfp_decompress(zfp, field); | |
if (!zfpsize) { | |
fprintf(stderr, "decompression failed\n"); | |
status = EXIT_FAILURE; | |
} | |
} | |
else { | |
/* compress array and output compressed stream */ | |
zfpsize = zfp_compress(zfp, field); | |
if (!zfpsize) { | |
fprintf(stderr, "compression failed\n"); | |
status = EXIT_FAILURE; | |
} | |
else { | |
FILE *fp = fopen("compressed_file", "wb"); | |
fwrite(buffer, 1, zfpsize, fp); | |
fclose(fp); | |
*cmp_size = zfpsize; | |
} | |
} | |
/* clean up */ | |
zfp_field_free(field); | |
zfp_stream_close(zfp); | |
stream_close(stream); | |
free(buffer); | |
return status; | |
} | |
/* Comprares values of decompressed array with raw array */ | |
double | |
calculate_rmse(md *c_array, md *d_array, int nx, int ny, int nz) | |
{ | |
long int nt = nx * ny * nz; /* total num values */ | |
double rmse, error2 = 0.; | |
//long int i, j, k; /* array index vars */ | |
// #pragma omp parallel for reduction(+:error2) | |
for (int k = 0; k < nz; k++) | |
for (int j = 0; j < ny; j++) | |
for (int i = 0; i < nx; i++) { | |
error2 += (d_array[i + nx * (j + ny * k)] - c_array[i + nx * (j + ny * k)]) * | |
(d_array[i + nx * (j + ny * k)] - c_array[i + nx * (j + ny * k)]); | |
} | |
printf("\n \n"); | |
rmse = sqrt(error2/nt); | |
return rmse; | |
} | |
int main(int argc, char* argv[]) | |
{ | |
int status = 0; | |
int decompress = 0; | |
double rmse = 0.; | |
size_t cmp_size = 0; | |
double bpv = 16.0; // bits per value 16, compression rate 2 | |
/* allocate 100x100x100 array of floats */ | |
int nx = 768; | |
int ny = 768; | |
int nz = 768; | |
md* c_array = malloc(nx * ny * nz * sizeof(md)); | |
md* d_array = malloc(nx * ny * nz * sizeof(md)); | |
fillGridWithVal(c_array, nx, ny, nz, 0.); | |
fillGridWithVal(d_array, nx, ny, nz, 0.); | |
/* initialize array to be compressed */ | |
int i, j, k; | |
for (k = 0; k < nz; k++) | |
for (j = 0; j < ny; j++) | |
for (i = 0; i < nx; i++) { | |
md x = 2.0 * i / nx; | |
md y = 2.0 * j / ny; | |
md z = 2.0 * k / nz; | |
//c_array[i + nx * (j + ny * k)] = exp(-(x * x + y * y + z * z)); | |
c_array[i + nx * (j + ny * k)] = input_func_cos(i, j, k, nx, ny, nz); | |
} | |
/* compress */ | |
status = compress(c_array, nx, ny, nz, bpv, 0, &cmp_size); | |
/* decomrpess */ | |
if (!status) { | |
decompress = 1; | |
status = compress(d_array, nx, ny, nz, bpv, decompress, &cmp_size); | |
//writeStructuredPoints(nx, ny, nz, d_array, 1, "uncomp_cos_serial"); | |
rmse = calculate_rmse(c_array, d_array, nx, ny, nz); | |
printf("RMSE = %e\n", rmse); | |
} | |
else { | |
printf("! Error in Compression \n"); | |
} | |
free(c_array); | |
free(d_array); | |
return status; | |
} |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment