diff --git a/.github/workflows/build.yml b/.github/workflows/build.yml new file mode 100644 index 0000000..d3bc4d4 --- /dev/null +++ b/.github/workflows/build.yml @@ -0,0 +1,36 @@ +name: build +on: [push, pull_request] + +jobs: + build_debug: + runs-on: ${{ matrix.os }} + strategy: + fail-fast: false + matrix: + os: [ubuntu-latest,macos-latest] + toolchain: + - {compiler: gcc, cxx: g++-12} + - {compiler: clang, cxx: clang++} + env: + CXX: ${{ matrix.toolchain.cxx }} + + steps: + - name: Checkout source code + uses: actions/checkout@v3 + + - name: Install argp-standalone (MacOS) + if: contains(matrix.os, 'macos') + run: | + brew install argp-standalone + + - name: Install libtirpc (Ubuntu) + if: contains(matrix.os,'ubuntu') + run: | + sudo apt-get install libtirpc-dev + export TIRPC_DIR=/usr + + - name: Build hemeXtract + run: | + mkdir build && cd build + cmake .. -DCMAKE_BUILD_TYPE=Debug + make hemeXtract \ No newline at end of file diff --git a/CMakeLists.txt b/CMakeLists.txt new file mode 100644 index 0000000..db41814 --- /dev/null +++ b/CMakeLists.txt @@ -0,0 +1,31 @@ +cmake_minimum_required(VERSION 3.15 FATAL_ERROR) + +project(hemeXtract + VERSION 0.1.0 + LANGUAGES CXX) + +set(CMAKE_CXX_STANDARD 14) + +set(CMAKE_MODULE_PATH "${CMAKE_SOURCE_DIR}/cmake") + +find_package(argp REQUIRED) + +if(CMAKE_SYSTEM_NAME STREQUAL "Linux") + find_package(TIRPC REQUIRED) +endif() + +add_executable(hemeXtract + src/hemeXtract.cc + src/Mapping.cc + src/Snapshot.cc + src/Compare.cc + src/HemeLBExtractionFile.cc) +target_include_directories(hemeXtract PRIVATE include/) + +target_include_directories(hemeXtract PRIVATE ${argp_INCLUDE_DIRS}) +target_link_libraries(hemeXtract PRIVATE ${argp_LIBRARIES}) + +if(TIRPC_FOUND) + target_include_directories(hemeXtract PRIVATE ${TIRPC_INCLUDE_DIRS}) + target_link_libraries(hemeXtract PRIVATE ${TIRPC_LIBRARY}) +endif() diff --git a/HemeLBExtractionFile.h b/HemeLBExtractionFile.h deleted file mode 100644 index 5131eb0..0000000 --- a/HemeLBExtractionFile.h +++ /dev/null @@ -1,606 +0,0 @@ -#ifndef INCLUDED_HEMELBEXTRACTIONFILE_H -#define INCLUDED_HEMELBEXTRACTIONFILE_H - -#include -#include -#include -#include -#include -#include -#include -#include -//#include -#include -#include -#include - -#include "HemeLBExtractionFileTypes.h" -#include "Snapshot.h" -#include "Vector3.h" - -class HemeLBExtractionFile -{ - public: - HemeLBExtractionFile(char *fname, double step_length, double scaling, bool verbose, Vector3 *translate) - { - this->step_length = step_length; - this->scaling = scaling; - this->bool_verbose = verbose; - bool_correctly_initialised = false; - snapshot = NULL; - in = fopen(fname, "rb"); - if (in == NULL) { - fprintf(stderr, "Could not open file '%s'. Does it exist?\n", fname); - return; - } - xdrstdio_create(&xdrs, in, XDR_DECODE); - - has_velocity = false; - has_pressure = false; - has_shearstress = false; - - // Read the file header, populating the 'header' class variable. Determines whether file is "normal" or colloids - int r; - r = read_header(); - if(r != 0) { - fprintf(stderr, "Problem reading header for file: '%s'.\n", fname); - return; - } - - // If file is a colloids file, return now (no more info to be read in) - if(header->is_colloids_file == true) { - bool_correctly_initialised = true; - return; - } - - /// Otherwise, if file is "normal", read the field header - r = read_field_header(); - if(r != 0) { - fprintf(stderr, "Problem reading field header for file: '%s'.\n", fname); - return; - } - - // If specified, convert lattice translation vector (real units) into grid units (via rounding based on voxel size) - if(translate != NULL) { - this->translate_grid_x = (int32_t)(translate->get_x()/header->voxelsz); - this->translate_grid_y = (int32_t)(translate->get_y()/header->voxelsz); - this->translate_grid_z = (int32_t)(translate->get_z()/header->voxelsz); - if(bool_verbose == true) { - fprintf(stderr, "Lattice translation vector (%e,%e,%e) in real units, converted to (%d,%d,%d) in grid units (voxel size = %e).\n", translate->get_x(), translate->get_y(), translate->get_z(), translate_grid_x, translate_grid_y, translate_grid_z, header->voxelsz); - } - } else { - this->translate_grid_x = 0; - this->translate_grid_y = 0; - this->translate_grid_z = 0; - } - - // Allocate a snapshot of the right dimensions - snapshot = new Snapshot(header, field_header); - - // Get the time of the next snapshot - bool next = read_time_next(); - if(next == false) { - bool_no_more_snapshots = true; - } else { - bool_no_more_snapshots = false; - load_next_snapshot(); // Load in the first snapshot - bool_correctly_initialised = true; - } - } - - ~HemeLBExtractionFile() - { - xdr_destroy(&xdrs); - fclose(in); - } - - bool read_time_next() - { - // Check if we have reached the end of the file - if (feof(in)) { - if(bool_verbose == true) { - fprintf(stderr, "# Reached end of file.\n"); - } - return false; - } - - // Work around for problem with xdr_u_long() not working - uint32_t dt1, dt2; - uint64_t step; - xdr_u_int(&xdrs, &dt1); - xdr_u_int(&xdrs, &dt2); - step = ((uint64_t)dt1)<<32 | dt2; - - time_next = step * step_length; - return true; - } - - int load_next_snapshot() - { - if(bool_no_more_snapshots == true) { - if(bool_verbose == true) { - fprintf(stderr, "Warning: No more snapshots left in file.\n"); - } - return 1; - } - - uint32_t gridposx, gridposy, gridposz; - float value; - - snapshot->set_timestep(time_next); - for(unsigned int s = 0; s < header->num_sites; s++) { - // Read in grid position of site - xdr_u_int(&xdrs, &gridposx); - xdr_u_int(&xdrs, &gridposy); - xdr_u_int(&xdrs, &gridposz); - - // Apply any requested translation of this xtr file directly to the site positions - gridposx += translate_grid_x; - gridposy += translate_grid_y; - gridposz += translate_grid_z; - - snapshot->site_set(s, gridposx, gridposy, gridposz); - for(unsigned int i = 0; i < header->num_columns; i++) { - // Read value and add it to the appropriate data column (scaling it by the specified scaling) - xdr_float(&xdrs, &value); - snapshot->column_set_plus_offset(i, s, value * this->scaling); - } - } - bool next = read_time_next(); - if(next == false) { - bool_no_more_snapshots = true; - if(bool_verbose == true) {fprintf(stderr, "No more snapshots. (Note that final snapshot may have been incomplete)\n");} - } - return 0; - } - - // Workaround for XDR's broken long reading ability - void xdr_long(XDR *xdrs, uint64_t *ret) - { - uint32_t half1, half2; - xdr_u_int(xdrs, &half1); - xdr_u_int(xdrs, &half2); - *ret = ((uint64_t)half1)<<32 | half2; // Assume big-endianness - } - - - void read_and_print_colloids(FILE *outfile) - { - uint32_t headerLen, recordLen; - uint64_t dsetLen, timeStep; - uint64_t id, rank; - double X, Y, Z, VX, VY, VZ; - - while(!feof(in)) { - - // Reinitialise reading vars to catch partially written header - headerLen = 0; - recordLen = 0; - dsetLen = 0; - timeStep = 0; - - // Read timestep header - xdr_u_int(&xdrs, &headerLen); - xdr_u_int(&xdrs, &recordLen); - xdr_long(&xdrs, &dsetLen); - xdr_long(&xdrs, &timeStep); - - // Check for end of file here (to avoid spurious EOF output) - if(feof(in)) { - break; - } - - if(recordLen <= 0) { - break; - } - - if(bool_verbose == true) { - fprintf(outfile, "# headerLen: %u recordLen %u dsetLen %ld timeStep: %ld\n", headerLen, recordLen, dsetLen, timeStep); - } - - // Calc. num particles from the record length data - uint32_t num_particles = dsetLen/recordLen; - - for (uint32_t i = 0; i < num_particles; i++) { - // Read particle data - xdr_long(&xdrs, &id); - xdr_long(&xdrs, &rank); - xdr_double(&xdrs, &X); - xdr_double(&xdrs, &Y); - xdr_double(&xdrs, &Z); - xdr_double(&xdrs, &VX); - xdr_double(&xdrs, &VY); - xdr_double(&xdrs, &VZ); - - // Rescale the positions by the scaling factor provided to hemeXtract - X *= this->scaling; - Y *= this->scaling; - Z *= this->scaling; - -// fprintf(outfile, "TIME: %ld ID: %ld RANK: %ld A0: %e Ah: %e X: %e Y: %e Z: %e\n", timeStep, id, rank, A0, Ah, X, Y, Z); - fprintf(outfile, "%ld %ld %ld %.13e %.13e %.13e %.13e %.13e %.13e\n", timeStep, id, rank, X, Y, Z, VX, VY, VZ); - } - } - if(bool_verbose == true) { - fprintf(stderr, "# Reached end of file.\n"); - } - } - - bool correctly_initialised() - { - return bool_correctly_initialised; - } - - bool is_colloids_file() - { - return header->is_colloids_file; - } - - bool no_more_snapshots() { - return bool_no_more_snapshots; - } - - uint64_t get_num_sites() - { - return header->num_sites; - } - - double get_voxelsz() - { - return header->voxelsz; - } - - bool hasVelocity() - { - return has_velocity; - } - - bool hasShearStress() - { - return has_shearstress; - } - - bool hasPressure() - { - return has_pressure; - } - - double get_scaling() - { - return scaling; - } - - double get_scalar_quantity(uint32_t column_index, uint64_t site_index) - { - return snapshot->get(column_index, site_index); - } - - double get_interpolated_scalar_quantity(uint32_t column_index, lattice_map *map) - { - double average_existing = 0; - uint32_t num_existing = 0; - for(uint32_t i = 0; i < 8; i++) { - if(map->index[i].exists == true) { - average_existing += get_scalar_quantity(column_index, map->index[i].index); - num_existing++; - } - } - if(num_existing == 0) { - return 0; - } - average_existing /= num_existing; - - // Get the scalars at the 8 lattice sites forming the cube - double s[8]; - for(uint32_t i = 0; i < 8; i++) { - // If there was no corresponding site_index, the site was not measured, so set it to the average - // value of the existing sites in the map - if(map->index[i].exists == false) { - s[i] = average_existing; - } else { - s[i] = get_scalar_quantity(column_index, map->index[i].index); - } - } - - // Trilinear interpolation within the cube - return s[0] * (1 - map->a) * (1 - map->b) * (1 - map->c) - + s[1] * (map->a * (1 - map->b) * (1 - map->c)) - + s[2] * ((1 - map->a) * map->b * (1 - map->c)) - + s[3] * ((1 - map->a) * (1 - map->b) * map->c) - + s[4] * (map->a * (1 - map->b) * map->c) - + s[5] * ((1 - map->a) * map->b * map->c) - + s[6] * (map->a * map->b * (1 - map->c)) - + s[7] * (map->a * map->b * map->c); - } - void get_vector_quantity(uint32_t column_index, uint64_t site_index, Vector3 *returned_val) - { - double vx, vy, vz; - vx = get_scalar_quantity(column_index, site_index); - vy = get_scalar_quantity(column_index + 1, site_index); - vz = get_scalar_quantity(column_index + 2, site_index); - returned_val->set(vx, vy, vz); - } - - void get_interpolated_vector_quantity(uint32_t column_index, lattice_map *map, Vector3 *returned_val) - { - double vx, vy, vz; - vx = get_interpolated_scalar_quantity(column_index, map); - vy = get_interpolated_scalar_quantity(column_index + 1, map); - vz = get_interpolated_scalar_quantity(column_index + 2, map); - returned_val->set(vx, vy, vz); - } - - void get_pressure(uint64_t site_index, double *returned_val) - { - *returned_val = get_scalar_quantity(column_pressure, site_index); - } - - void get_velocity(uint64_t site_index, Vector3 *returned_val) - { - get_vector_quantity(column_velocity, site_index, returned_val); - } - - void get_shearstress(uint64_t site_index, double *returned_val) - { - *returned_val = get_scalar_quantity(column_shearstress, site_index); - } - - void get_interpolated_pressure(lattice_map *map, double *returned_val) - { - *returned_val = get_interpolated_scalar_quantity(column_pressure, map); - } - - void get_interpolated_velocity(lattice_map *map, Vector3 *returned_val) - { - get_interpolated_vector_quantity(column_velocity, map, returned_val); - } - - void get_interpolated_shearstress(lattice_map *map, double *returned_val) - { - *returned_val = get_interpolated_scalar_quantity(column_shearstress, map); - } - - Site * get_sites() - { - return snapshot->get_sites(); - } - - SiteIndex * get_site_indices(Site *list, uint64_t list_size) - { - return snapshot->get_site_indices(list, list_size); - } - - SiteIndex * get_site_indices_hashed_lookup(Site *list, uint64_t list_size) - { - return snapshot->get_site_indices_hashed_lookup(list, list_size, bool_verbose); - } - - double get_time() - { - return snapshot->get_timestep(); - } - - double get_time_next() - { - return time_next; - } - - void print_header(FILE *outfile) - { - fprintf(outfile, "\n# HEADER:\n# -------\n"); - fprintf(outfile, "# version=%u\n", header->version); - fprintf(outfile, "# voxelsz=%f\n", header->voxelsz); - fprintf(outfile, "# originx=%f\n", header->originx); - fprintf(outfile, "# originy=%f\n", header->originy); - fprintf(outfile, "# originz=%f\n", header->originz); - fprintf(outfile, "# num_sites=%lu\n", header->num_sites); - fprintf(outfile, "# field_count=%u\n", header->field_count); - fprintf(outfile, "# field_header_length=%u\n", header->field_header_length); - } - - void print_field_header(FILE *outfile) - { - fprintf(outfile, "\n# FIELD HEADERS:\n# -------\n"); - for(unsigned int i = 0; i < header->field_count; i++) { - fprintf(outfile, "# fieldname=%s\n# num_floats=%u\n# offset=%lf\n# ---\n", field_header[i].name, field_header[i].num_floats, field_header[i].offset); - } - fprintf(outfile, "# Number of 'columns' per snapshot = %d\n", header->num_columns); - } - - void print_column_headings(FILE *outfile) - { - fprintf(outfile, "# step | grid_x | grid_y | grid_z"); - for(unsigned int i = 0; i < header->field_count; i++) { - if(field_header[i].num_floats == 1) { - fprintf(outfile, " | %s", field_header[i].name); - } else { - for(unsigned int j = 0; j < field_header[i].num_floats; j++) { - fprintf(outfile, " | %s(%d)", field_header[i].name, j); - } - } - } - fprintf(outfile, "\n"); - } - void print_stats_column_headings(FILE *outfile) - { - fprintf(outfile, "# step"); - for(unsigned int i = 0; i < header->field_count; i++) { - if(field_header[i].num_floats == 1) { - fprintf(outfile, " | %s [average | stdev | max | min]", field_header[i].name); - } else { - for(unsigned int j = 0; j < field_header[i].num_floats; j++) { - fprintf(outfile, " | %s(%d) [average | stdev | max | min]", field_header[i].name, j); - } - } - } - fprintf(outfile, "\n"); - } - - void print_all(FILE *outfile) - { - snapshot->print(outfile); - } - - void print_stats(FILE *outfile) - { - snapshot->print_stats(outfile); - } - - private: - /* True if file has been opened and the headers read without problem */ - bool bool_correctly_initialised; - - /* True if user wants verbose output. If false, means quiet. */ - bool bool_verbose; - - /* All info contained in the file header */ - HEADER *header; - FIELD_HEADER *field_header; - - /** The step length for this file (has to be provided by the user since Heme extraction files don't store this value...) */ - double step_length; - - /** The velocity scaling factor */ - double scaling; - - /** The currently loaded snapshot */ - Snapshot *snapshot; - - /** The open extraction file's stream */ - FILE *in; - - /** The XDR reader for this file */ - XDR xdrs; - - /** The time of the next snapshot (next to be loaded) */ - double time_next; - - /** True when no more snapshots can be loaded */ - bool bool_no_more_snapshots; - - /** True if the file contains these field types */ - bool has_velocity, has_shearstress, has_pressure; - - /** The column indices in which to find each field type */ - uint32_t column_velocity; - uint32_t column_pressure; - uint32_t column_shearstress; - - /** The translation vector (in grid units) to be applied to this lattice (if requested by --translate option) */ - int32_t translate_grid_x; - int32_t translate_grid_y; - int32_t translate_grid_z; - - /** Magic numbers designating a heme file, a heme extraction file and a colloid file */ - static const uint32_t heme_magic = 0x686C6221; // hlb! - static const uint32_t extract_magic = 0x78747204; // xtr - static const uint32_t colloid_magic = 0x636F6C04; // colloid - - /** Reads a HemeLB extraction file header. Checks that the magic numbers are correct. */ - int read_header() - { - uint32_t magic1 = 0, magic2 = 0; - uint32_t num_sites1, num_sites2; - - // Read heme magic number - xdr_u_int(&xdrs, &magic1); - if(magic1 != heme_magic) { - fprintf(stderr, "First uint32_t does not match heme magic number.\n"); - return 1; - } - - header = new HEADER(); - - // Read extraction magic number - xdr_u_int(&xdrs, &magic2); - if(magic2 == extract_magic) { - if(bool_verbose == true) { - fprintf(stderr, "Reading a normal extraction file.\n"); - } - - header->is_colloids_file = false; - - xdr_u_int(&xdrs, &(header->version)); - xdr_double(&xdrs, &(header->voxelsz)); - xdr_double(&xdrs, &(header->originx)); - xdr_double(&xdrs, &(header->originy)); - xdr_double(&xdrs, &(header->originz)); - - // There is something wrong with either Heme or XDR which means that xdr_u_long is not working (possibly even reading more/less than 8 bytes...) - xdr_u_int(&xdrs, &num_sites1); - xdr_u_int(&xdrs, &num_sites2); - - // Combine the two ints into a long (assuming bigendian-ness) - header->num_sites = ((uint64_t)num_sites1)<<32 | num_sites2; - - xdr_u_int(&xdrs, &(header->field_count)); - xdr_u_int(&xdrs, &(header->field_header_length)); - - } else if (magic2 == colloid_magic) { - if(bool_verbose == true) { - fprintf(stderr, "Reading a colloids extraction file.\n"); - } - header->is_colloids_file = true; - - // Check version - uint32_t version=0; - xdr_u_int(&xdrs, &version); - if(bool_verbose == true) { - printf("Colloids version %u\n", version); - } - } else { - fprintf(stderr, "Unknown format: Second uint32_t does not match extraction file or colloids file magic number.\n"); - return 1; - } - - // Everything went well - return 0; - } - - /** Reads the field headers and remembers which columns correspond to velocity, pressure and shearstress. */ - int read_field_header() - { - if(bool_verbose == true) { - printf("Reading field header...\n"); - } - int const maxlength = 100; // surely it can't be bigger than this...? - field_header = new FIELD_HEADER[header->field_count]; - for(unsigned int i = 0; i < header->field_count; i++) { - field_header[i].name = new char[maxlength]; - xdr_string(&xdrs, &(field_header[i].name), maxlength); - xdr_u_int(&xdrs, &(field_header[i].num_floats)); - xdr_double(&xdrs, &(field_header[i].offset)); - } - if(bool_verbose == true) { - printf("Finished reading field header...\n"); - } - - // Calc. number of 'columns' needed to represent a snapshot - uint32_t num_columns = 0; - for(unsigned int i = 0; i < header->field_count; i++) { - // If file contains velocity, pressure or shearstress, remember the corresponding column index for that field type. - if(strcmp(field_header[i].name, "velocity") == 0) { - has_velocity = true; - column_velocity = num_columns; - } else if(strcmp(field_header[i].name, "pressure") == 0) { - has_pressure = true; - column_pressure = num_columns; - } else if(strcmp(field_header[i].name, "shearstress") == 0) { - has_shearstress = true; - column_shearstress = num_columns; - } else if(strcmp(field_header[i].name, "d_shearstress") == 0) { // For "legacy" reasons only. Some old MCA runs have this field. - has_shearstress = true; - column_shearstress = num_columns; - } - - // This field has 'num_floats' columns to it - num_columns += field_header[i].num_floats; - } - header->num_columns = num_columns; - - // Everything went well - return 0; - } -}; - -#endif diff --git a/Makefile b/Makefile deleted file mode 100644 index 22ff3f5..0000000 --- a/Makefile +++ /dev/null @@ -1,4 +0,0 @@ -hemeXtract: - g++ -Wall -O3 hemeXtract.cc -o hemeXtract -debug: - g++ -Wall -g hemeXtract.cc -o hemeXtract diff --git a/Mapping.h b/Mapping.h deleted file mode 100644 index 6a951fd..0000000 --- a/Mapping.h +++ /dev/null @@ -1,69 +0,0 @@ -#ifndef INCLUDED_MAPPING_H -#define INCLUDED_MAPPING_H - -#include -#include -#include -#include -#include -#include -#include - -#include "HemeLBExtractionFileTypes.h" -#include "HemeLBExtractionFile.h" - -void map(double old_voxelsz, double new_voxelsz, uint32_t in_coord, uint32_t *out_root, double *out_interpolate) -{ - double real_pos = in_coord * old_voxelsz; - *out_root = floor(real_pos/new_voxelsz); - *out_interpolate = (real_pos - *out_root * new_voxelsz) / new_voxelsz; -} - -/** Returns a mapping for each site in A to a corresponding trilinearly interpolated position in B. */ -lattice_map * get_mapping(HemeLBExtractionFile *A, HemeLBExtractionFile *B) -{ - // Work out mapping between the two lattices - uint64_t num_sites_A = A->get_num_sites(); - Site *sitesA = A->get_sites(); - double voxelA = A->get_voxelsz(); - double voxelB = B->get_voxelsz(); - uint32_t mapped_x, mapped_y, mapped_z; - lattice_map *mapping_A_B = new lattice_map[num_sites_A]; - if(mapping_A_B == NULL) { - fprintf(stderr, "Could not allocate memory for lattice_map\n"); - return NULL; - } - Site *mapped_sites = new Site[num_sites_A * 8]; - if(mapped_sites == NULL) { - fprintf(stderr, "Could not allocate memory for mapped_sites\n"); - return NULL; - } - - for(uint64_t i = 0; i < num_sites_A; i++) { - // Get the interpolation constants for this site in terms of the other lattice - map(voxelA, voxelB, sitesA[i].x, &mapped_x, &mapping_A_B[i].a); - map(voxelA, voxelB, sitesA[i].y, &mapped_y, &mapping_A_B[i].b); - map(voxelA, voxelB, sitesA[i].z, &mapped_z, &mapping_A_B[i].c); - - // Add the 8 sites forming the cube this point lies within - mapped_sites[i * 8].set(mapped_x, mapped_y, mapped_z); // root (000) - mapped_sites[i * 8 + 1].set(mapped_x + 1, mapped_y, mapped_z); // (100) - mapped_sites[i * 8 + 2].set(mapped_x, mapped_y + 1, mapped_z); // (010) - mapped_sites[i * 8 + 3].set(mapped_x, mapped_y, mapped_z + 1); // (001) - mapped_sites[i * 8 + 4].set(mapped_x + 1, mapped_y, mapped_z + 1); // (101) - mapped_sites[i * 8 + 5].set(mapped_x, mapped_y + 1, mapped_z + 1); // (011) - mapped_sites[i * 8 + 6].set(mapped_x + 1, mapped_y + 1, mapped_z); // (110) - mapped_sites[i * 8 + 7].set(mapped_x + 1, mapped_y + 1, mapped_z + 1); // (111) - } - - SiteIndex *site_indices = B->get_site_indices_hashed_lookup(mapped_sites, num_sites_A * 8); - - for(uint64_t i = 0; i < num_sites_A; i++) { - for(uint32_t j = 0; j < 8; j++) { - mapping_A_B[i].index[j] = site_indices[i * 8 + j]; - } - } - - return mapping_A_B; -} -#endif diff --git a/Site.h b/Site.h deleted file mode 100644 index 9afd5c9..0000000 --- a/Site.h +++ /dev/null @@ -1,62 +0,0 @@ -#ifndef INCLUDED_SITE_H -#define INCLUDED_SITE_H - -#include -#include -#include -#include -#include -#include -#include -#include -//#include -#include -#include -#include - -class Site { - public: - Site() - { - this->x = 0; - this->y = 0; - this->z = 0; - } - - void set(uint32_t x, uint32_t y, uint32_t z) { - this->x = x; - this->y = y; - this->z = z; - } - - uint32_t x, y, z; - - void print(FILE *outfile) { - fprintf(outfile, "%u %u %u", x, y, z); - } - - void print(FILE *outfile, double voxelsz) { - fprintf(outfile, "%f %f %f", x * voxelsz, y * voxelsz, z * voxelsz); - } - - bool equals(uint32_t a, uint32_t b, uint32_t c) - { - if(a != x) { - return false; - } - if(b != y) { - return false; - } - if(c != z) { - return false; - } - return true; - } - - bool equals(Site *s) - { - return equals(s->x, s->y, s->z); - } -}; - -#endif diff --git a/Snapshot.h b/Snapshot.h deleted file mode 100644 index 13e14c5..0000000 --- a/Snapshot.h +++ /dev/null @@ -1,238 +0,0 @@ -#ifndef INCLUDED_SNAPSHOT_H -#define INCLUDED_SNAPSHOT_H - -#include -#include -#include -#include -#include -#include -#include -#include -//#include -#include -#include -#include - -#include "HemeLBExtractionFileTypes.h" -#include "Site.h" -#include "Mapping.h" - -class Column { - public: - Column(unsigned int num_records, double offset) - { - this->num_records = num_records; - this->offset = offset; - stats_calcd = false; - records = new double[num_records]; - } - ~Column() - { - delete[] records; - } - - double get(uint64_t index) - { - if(index >= num_records) { - fprintf(stderr, "Error: index (%lu) exceeds number of records in column (%lu)\n", index, num_records); - exit(1); - } - return records[index]; - } - void set_plus_offset(uint64_t index, float value) - { - if(index >= num_records) { - fprintf(stderr, "Error: index (%lu) exceeds number of records in column (%lu)\n", index, num_records); - exit(1); - } - records[index] = value + offset; - stats_calcd = false; - } - void print(FILE *outfile, uint64_t index) { - if(index >= num_records) { - fprintf(stderr, "Error: index (%lu) exceeds number of records in column (%lu)\n", index, num_records); - exit(1); - } - fprintf(outfile, "%e", records[index]); - } - void print_stats(FILE *outfile) { - calc_stats(); - fprintf(outfile, "%e %e %e %e", average, stdev, max, min); - } - - private: - uint64_t num_records; - double *records; - double offset; - bool stats_calcd; - double average; - double stdev; - double max, min; - void calc_stats() - { - if(!stats_calcd) { - double sum = 0, sum2 = 0; - max = -INFINITY; - min = INFINITY; - for(uint64_t i = 0; i < num_records; i++) { - sum += records[i]; - sum2 += records[i] * records[i]; - if(max < records[i]) { - max = records[i]; - } - if(min > records[i]) { - min = records[i]; - } - } - average = sum / num_records; - stdev = sqrt(sum2/num_records - average*average); - stats_calcd = true; - } - } -}; - -class Snapshot { - public: - Snapshot(HEADER *header, FIELD_HEADER *field_header) - { - this->num_sites = header->num_sites; - this->num_columns = header->num_columns; - this->voxelsz = header->voxelsz; - sites = new Site[num_sites]; - columns = new Column*[num_columns]; - unsigned int c = 0; - for(unsigned int i = 0; i < header->field_count; i++) { - for(unsigned int bro=0; bro < field_header[i].num_floats; bro++) { - columns[c] = new Column(num_sites, field_header[i].offset); - c++; - } - } - } - - void set_timestep(double timestep) - { - this->timestep = timestep; - } - double get_timestep() - { - return timestep; - } - - void site_set(unsigned int site_index, uint32_t x, uint32_t y, uint32_t z) - { - sites[site_index].set(x, y, z); - } - - void column_set_plus_offset(uint32_t column_index, uint64_t site_index, double value) - { - if(column_index >= num_columns) { - fprintf(stderr, "Error: column index (%u) exceeds number of columns (%u)\n", column_index, num_columns); - exit(1); - } - columns[column_index]->set_plus_offset(site_index, value); - } - - double get(uint32_t column_index, uint64_t site_index) - { - return columns[column_index]->get(site_index); - } - - bool get_site_index(uint32_t a, uint32_t b, uint32_t c, uint64_t *site_index) - { - for(unsigned int i = 0; i < num_sites; i++) { - if(sites[i].equals(a, b, c)) { - *site_index = i; - return true; - } - } - - // If no such site can be found, return false - return false; - } - - int get_site_index(Site *s, uint64_t *site_index) - { - return get_site_index(s->x, s->y, s->z, site_index); - } - - Site * get_sites() - { - return sites; - } - - /* Very inefficient calculation of site_index corresponding to the given site. Only suitable for small data sets. */ - SiteIndex * get_site_indices(Site *list, uint64_t list_size) - { - SiteIndex *indices = new SiteIndex[list_size]; - for(uint64_t i = 0; i < list_size; i++) { - indices[i].exists = get_site_index(&list[i], &(indices[i].index)); - } - return indices; - } - - /** Builds hashtable to get indices corresponding to the given site. Necessary for very large data sets. */ - SiteIndex * get_site_indices_hashed_lookup(Site *list, uint64_t list_size, bool bool_verbose) - { - // Build the hash table - if(bool_verbose == true) { fprintf(stderr, "# Building hashtable...\n"); } - //std::tr1::unordered_map hashtable; //JM for C++11 capability - std::unordered_map hashtable; - for(uint64_t i = 0; i < num_sites; i++) { - std::ostringstream oss; - oss << sites[i].x << "," << sites[i].y << "," << sites[i].z; - std::string hashstr = oss.str(); - hashtable[hashstr] = i; - } - if(bool_verbose == true) { fprintf(stderr, "# ...done building hashtable\n"); } - - if(bool_verbose == true) { fprintf(stderr, "# Calculating mapping...\n"); } - SiteIndex *indices = new SiteIndex[list_size]; - for(uint64_t i = 0; i < list_size; i++) { - std::ostringstream oss; - oss << list[i].x << "," << list[i].y << "," << list[i].z; - std::string hashstr = oss.str(); - if(hashtable.find(hashstr) != hashtable.end()) { - indices[i].exists = true; - indices[i].index = hashtable[hashstr]; - } else { - indices[i].exists = false; - } - } - if(bool_verbose == true) { fprintf(stderr, "# ... done calculating mapping\n"); } - - return indices; - } - - void print(FILE *outfile) - { - for(uint64_t s = 0; s < num_sites; s++) { - fprintf(outfile, "%f ", timestep); - sites[s].print(outfile, this->voxelsz); - for(uint32_t c = 0; c < num_columns; c++) { - fprintf(outfile, " "); - columns[c]->print(outfile, s); - } - fprintf(outfile, "\n"); - } - fprintf(outfile, "\n"); - } - void print_stats(FILE *outfile) - { - fprintf(outfile, "%f", timestep); - for(uint32_t c = 0; c < num_columns; c++) { - fprintf(outfile, " "); - columns[c]->print_stats(outfile); - } - fprintf(outfile, "\n"); - } - private: - double timestep; - uint64_t num_sites; - uint32_t num_columns; - double voxelsz; - Column **columns; - Site *sites; -}; - -#endif diff --git a/cmake/FindTIRPC.cmake b/cmake/FindTIRPC.cmake new file mode 100644 index 0000000..c5b0037 --- /dev/null +++ b/cmake/FindTIRPC.cmake @@ -0,0 +1,24 @@ +find_package(PkgConfig QUIET) +pkg_check_modules(PC_TIRPC libtirpc) + +# check for header file +find_path(TIRPC_INCLUDE_DIRS + NAMES netconfig.h + HINTS ${TIRPC_DIR}/include $ENV{TIRPC_DIR}/include ${PC_TIRPC_INCLUDE_DIRS} + PATH_SUFFIXES tirpc + DOC "directory where the TIRPC header is located") + +# check for tirpc +find_library(TIRPC_LIBRARY + NAMES tirpc + HINTS ${TIRPC_DIR}/lib $ENV{TIRPC_DIR}/lib ${PC_TIRPC_LIBRARY_DIRS} + NO_DEFAULT_PATH + DOC "the TIRPC library") +find_library(TIRPC_LIBRARY + NAMES tirpc + DOC "the TIRPC library") + +include(FindPackageHandleStandardArgs) +find_package_handle_standard_args(TIRPC + REQUIRED_VARS TIRPC_INCLUDE_DIRS TIRPC_LIBRARY + VERSION_VAR TIRPC_VERSION) \ No newline at end of file diff --git a/cmake/Findargp.cmake b/cmake/Findargp.cmake new file mode 100644 index 0000000..f9a3cf6 --- /dev/null +++ b/cmake/Findargp.cmake @@ -0,0 +1,78 @@ +#[=======================================================================[.rst: +Findargp +------- + +Finds the argp library. + +Result Variables +^^^^^^^^^^^^^^^^ + +This will define the following variables: + +``argp_FOUND`` + True if the system has the argp library. +``argp_INCLUDE_DIRS`` + Include directories needed to use argp. +``argp_LIBRARIES`` + Libraries needed to link to argp. + +Cache Variables +^^^^^^^^^^^^^^^ + +The following cache variables may also be set: + +``argp_INCLUDE_DIR`` + The directory containing ``argp.h``. +``argp_LIBRARY`` + The path to the argp library. + +#]=======================================================================] + +include(FindPackageHandleStandardArgs) +include(CheckCXXSymbolExists) + +set(CMAKE_REQUIRED_QUITE TRUE) + +# Check if header file exists +find_path(argp_INCLUDE_DIR + NAMES argp.h + HINTS /usr/include /usr/local/include) + +set(CMAKE_REQUIRED_INCLUDES ${argp_INCLUDE_DIR}) + +if(argp_INCLUDE_DIR) + # Check for function argp_parse in default libraries (libc) + check_cxx_symbol_exists(argp_parse "argp.h" ARGP_LIBC) + if(ARGP_LIBC) + set(argp_LIBRARY c CACHE STRING + "argp library file") + else() + find_library(argp_LIBRARY "argp" + HINTS /usr/local/lib + DOC "argp library file") + if (argp_LIBRARY) + set(CMAKE_REQUIRED_LIBRARIES ${argp_LIBRARY}) + + # Check for function argp_parse + check_cxx_symbol_exists(argp_parse "argp.h" ARGP_EXTERNAL) + if(NOT ARGP_EXTERNAL) + message(FATAL_ERROR + "The argp library was found in ${argp_LIBRARY}," + "but it doesn't contain a symbol named argp_parse.") + endif() + endif() + endif() +endif() + +find_package_handle_standard_args(argp + DEFAULT_MSG argp_LIBRARY argp_INCLUDE_DIR) + +if(argp_FOUND) + set(argp_LIBRARIES ${argp_LIBRARY}) + set(argp_INCLUDE_DIRS ${argp_INCLUDE_DIR}) +endif() + +mark_as_advanced( + argp_LIBRARY + argp_INCLUDE_DIR +) \ No newline at end of file diff --git a/include/Compare.h b/include/Compare.h new file mode 100644 index 0000000..c587863 --- /dev/null +++ b/include/Compare.h @@ -0,0 +1,28 @@ +#ifndef INCLUDED_COMPARE_H +#define INCLUDED_COMPARE_H + +#include + +#include "Vector3.h" +#include "HemeLBExtractionFileTypes.h" // lattice_map +#include "HemeLBExtractionFile.h" // HemeLBExtractionFile + +/** Compares two lattices, A and B, with mapping mapA_to_B, by calculating their crosscorrelation, and their L2 norm difference. */ +void compare(FILE *outfile, + const lattice_map *mapA_to_B, + const HemeLBExtractionFile *A, + const HemeLBExtractionFile *B, + int minexistent, + bool normalize_correl); + +/** Calculates the difference between each site in A, and the corresponding (trilinearly interpolated) point in B */ +void diff(FILE *outfile, + const lattice_map *mapA_to_B, + const HemeLBExtractionFile *A, + const HemeLBExtractionFile *B, + int minexistent, + Vector3 *project, + bool relativeErr, + bool verbose); + +#endif diff --git a/include/HemeLBExtractionFile.h b/include/HemeLBExtractionFile.h new file mode 100644 index 0000000..5abe5ef --- /dev/null +++ b/include/HemeLBExtractionFile.h @@ -0,0 +1,147 @@ +#ifndef INCLUDED_HEMELBEXTRACTIONFILE_H +#define INCLUDED_HEMELBEXTRACTIONFILE_H + +#include +#include +#include // unique_ptr + +#include // must precede xdr.h +#include // XDR + +#include "Vector3.h" // Vector3 +#include "Site.h" // Site +#include "HemeLBExtractionFileTypes.h" // lattice_map, SiteIndex +#include "Snapshot.h" // Snapshot + +class HemeLBExtractionFile +{ + public: + // Constructor + HemeLBExtractionFile(char *fname, + double step_length, + double scaling, + bool verbose, + Vector3 *translate); + // Destructor + ~HemeLBExtractionFile(); + + bool read_time_next(); + + int load_next_snapshot(); + + void read_and_print_colloids(FILE *outfile); + + bool correctly_initialised() const; + + bool is_colloids_file() const; + + bool no_more_snapshots() const; + + uint64_t get_num_sites() const; + + double get_voxelsz() const; + + bool hasVelocity() const; + bool hasShearStress() const ; + bool hasPressure() const; + + double get_scaling() const; + + double get_scalar_quantity(uint32_t column_index, uint64_t site_index) const; + + double get_interpolated_scalar_quantity( + uint32_t column_index, + const lattice_map &map) const; + + void get_vector_quantity( + uint32_t column_index, + uint64_t site_index, + Vector3 *returned_val) const; + + void get_interpolated_vector_quantity( + uint32_t column_index, + const lattice_map &map, + Vector3 *returned_val) const; + + void get_pressure(uint64_t site_index, double *returned_val) const; + + void get_velocity(uint64_t site_index, Vector3 *returned_val) const; + + void get_shearstress(uint64_t site_index, double *returned_val) const; + + void get_interpolated_pressure(const lattice_map &map, double *returned_val) const; + + void get_interpolated_velocity(const lattice_map &map, Vector3 *returned_val) const; + + void get_interpolated_shearstress(const lattice_map &map, double *returned_val) const; + + const Site * get_sites() const; + + SiteIndex * get_site_indices(Site *list, uint64_t list_size); + SiteIndex * get_site_indices_hashed_lookup(Site *list, uint64_t list_size); + + double get_time() const; + double get_time_next() const; + + void print_header(FILE *outfile) const; + void print_field_header(FILE *outfile) const; + void print_column_headings(FILE *outfile) const; + void print_stats_column_headings(FILE *outfile) const; + void print_all(FILE *outfile) const; + void print_stats(FILE *outfile) const; + + private: + /* True if file has been opened and the headers read without problem */ + bool bool_correctly_initialised{false}; + + /* True if user wants verbose output. If false, means quiet. */ + bool bool_verbose; + + /* All info contained in the file header */ + std::unique_ptr
header{nullptr}; + std::vector field_header; + + /** The step length for this file (has to be provided by the user since Heme extraction files don't store this value...) */ + double step_length; + + /** The velocity scaling factor */ + double scaling; + + /** The currently loaded snapshot */ + std::unique_ptr snapshot{nullptr}; + + /** The open extraction file's stream */ + FILE *in; + + /** The XDR reader for this file */ + XDR xdrs; + + /** The time of the next snapshot (next to be loaded) */ + double time_next; + + /** True when no more snapshots can be loaded */ + bool bool_no_more_snapshots; + + /** True if the file contains these field types */ + bool has_velocity{false}; + bool has_shearstress{false}; + bool has_pressure{false}; + + /** The column indices in which to find each field type */ + uint32_t column_velocity; + uint32_t column_pressure; + uint32_t column_shearstress; + + /** The translation vector (in grid units) to be applied to this lattice (if requested by --translate option) */ + int32_t translate_grid_x{0}; + int32_t translate_grid_y{0}; + int32_t translate_grid_z{0}; + + /** Reads a HemeLB extraction file header. Checks that the magic numbers are correct. */ + int read_header(); + + /** Reads the field headers and remembers which columns correspond to velocity, pressure and shearstress. */ + int read_field_header(); +}; + +#endif diff --git a/HemeLBExtractionFileTypes.h b/include/HemeLBExtractionFileTypes.h similarity index 78% rename from HemeLBExtractionFileTypes.h rename to include/HemeLBExtractionFileTypes.h index d9ab48c..11a980c 100644 --- a/HemeLBExtractionFileTypes.h +++ b/include/HemeLBExtractionFileTypes.h @@ -1,6 +1,8 @@ #ifndef INCLUDED_HEMELBEXTRACTIONFILETYPES_H #define INCLUDED_HEMELBEXTRACTIONFILETYPES_H +#include + typedef struct { uint32_t version; double voxelsz; @@ -25,19 +27,15 @@ typedef struct { class SiteIndex { public: - SiteIndex() - { - index = 0; - exists = false; - } - uint64_t index; - bool exists; + SiteIndex() {} + uint64_t index{0}; + bool exists{false}; }; typedef struct { - SiteIndex index[8]; - double a, b, c; + SiteIndex index[8]; + double a, b, c; } lattice_map; #endif diff --git a/include/Mapping.h b/include/Mapping.h new file mode 100644 index 0000000..4f7c41b --- /dev/null +++ b/include/Mapping.h @@ -0,0 +1,12 @@ +#ifndef INCLUDED_MAPPING_H +#define INCLUDED_MAPPING_H + +#include "HemeLBExtractionFileTypes.h" // lattice_map +#include "HemeLBExtractionFile.h" // HemeLBExtractionFile + +/* Returns a mapping for each site in A to a corresponding + * trilinearly interpolated position in B. + */ +lattice_map *get_mapping(HemeLBExtractionFile *A, HemeLBExtractionFile *B); + +#endif diff --git a/include/Site.h b/include/Site.h new file mode 100644 index 0000000..a3c8e49 --- /dev/null +++ b/include/Site.h @@ -0,0 +1,51 @@ +#ifndef INCLUDED_SITE_H +#define INCLUDED_SITE_H + +#include +#include + +class Site +{ + public: + + Site() {} + Site(uint32_t x_, uint32_t y_, uint32_t z_) + : x(x_), y(y_), z(z_) {} + + void set(uint32_t x, uint32_t y, uint32_t z) { + this->x = x; + this->y = y; + this->z = z; + } + + void print(FILE *outfile) const { + fprintf(outfile, "%u %u %u", x, y, z); + } + + void print(FILE *outfile, double voxelsz) const { + fprintf(outfile, "%f %f %f", x * voxelsz, y * voxelsz, z * voxelsz); + } + + inline bool equals(uint32_t a, uint32_t b, uint32_t c) const + { + if(a != x) { + return false; + } + if(b != y) { + return false; + } + if(c != z) { + return false; + } + return true; + } + + bool equals(const Site &s) const + { + return equals(s.x, s.y, s.z); + } + + uint32_t x{0}, y{0}, z{0}; +}; + +#endif diff --git a/include/Snapshot.h b/include/Snapshot.h new file mode 100644 index 0000000..a72882d --- /dev/null +++ b/include/Snapshot.h @@ -0,0 +1,51 @@ +#ifndef INCLUDED_SNAPSHOT_H +#define INCLUDED_SNAPSHOT_H + +#include +#include +#include + +#include "Site.h" // Site +#include "HemeLBExtractionFileTypes.h" // HEADER, FIELD_HEADER, SiteIndex + +class Column; // Forward declaration (could be private potentially) + +class Snapshot { + public: + Snapshot(const HEADER *header, const FIELD_HEADER *field_header); + + void set_timestep(double timestep); + double get_timestep() const; + + void site_set(unsigned int site_index, uint32_t x, uint32_t y, uint32_t z); + + void column_set_plus_offset(uint32_t column_index, uint64_t site_index, double value); + + double get(uint32_t column_index, uint64_t site_index) const; + + bool get_site_index(const Site &s, uint64_t &site_index) const; + + const Site *get_sites() const; + + /* Very inefficient calculation of site_index corresponding to the + * given site. Only suitable for small data sets. */ + SiteIndex *get_site_indices(const Site *list, uint64_t list_size); + + /* Builds hashtable to get indices corresponding to the given + * site. Necessary for very large data sets. */ + SiteIndex *get_site_indices_hashed_lookup(const Site *list, + uint64_t list_size, bool bool_verbose); + + void print(FILE *outfile); + void print_stats(FILE *outfile); + + private: + double timestep; + uint64_t num_sites; + uint32_t num_columns; + double voxelsz; + Column **columns; // Memory Leak!! + std::vector sites; +}; + +#endif diff --git a/Vector3.h b/include/Vector3.h similarity index 59% rename from Vector3.h rename to include/Vector3.h index 5ad8cb3..400560a 100644 --- a/Vector3.h +++ b/include/Vector3.h @@ -1,35 +1,18 @@ #ifndef INCLUDED_VECTOR3_H #define INCLUDED_VECTOR3_H -#include -#include -#include -#include -#include -#include -#include -#include -//#include -#include -#include -#include +#include +#include /** Reinventing the wheel... */ class Vector3 { public: - Vector3() - { - x = 0; - y = 0; - z = 0; - } - Vector3(double x, double y, double z) - { - this->x = x; - this->y = y; - this->z = z; - } + + Vector3() {} + + Vector3(double x_, double y_, double z_) + : x(x_), y(y_), z(z_) {} void set(double x, double y, double z) { @@ -45,12 +28,12 @@ class Vector3 this->z *= s; } - double dot(Vector3 *v) + double dot(Vector3 *v) const { return this->x*v->x + this->y*v->y + this->z*v->z; } - double length() + double length() const { return sqrt(x*x + y*y + z*z); } @@ -71,7 +54,7 @@ class Vector3 z /= l; } - double abs_diff(Vector3 *v) + double abs_diff(Vector3 *v) const { double dx = x - v->x; double dy = y - v->y; @@ -86,30 +69,29 @@ class Vector3 this->z += v->z * scale; } - double get_x() + double get_x() const { return this->x; } - double get_y() + double get_y() const { return this->y; } - double get_z() + double get_z() const { return this->z; } - void print(FILE *outfile) + void print(FILE *outfile) const { fprintf(outfile, "(%f %f %f)\n", x, y, z); } private: - //static const double tol = 1e-10; - static constexpr double tol = 1e-10; //JM for C++ 11 - double x, y, z; + static constexpr double tol{1e-10}; //JM for C++ 11 + double x{0}, y{0}, z{0}; }; #endif diff --git a/Compare.h b/src/Compare.cc similarity index 82% rename from Compare.h rename to src/Compare.cc index 05896d6..7a7f9b1 100644 --- a/Compare.h +++ b/src/Compare.cc @@ -1,16 +1,15 @@ -#ifndef INCLUDED_COMPARE_H -#define INCLUDED_COMPARE_H +#include "Compare.h" -#include -#include -#include -#include -#include -#include -#include +#include /** Compares two lattices, A and B, with mapping mapA_to_B, by calculating their crosscorrelation, and their L2 norm difference. */ -void compare(FILE *outfile, lattice_map *mapA_to_B, HemeLBExtractionFile *A, HemeLBExtractionFile *B, int minexistent, bool normalize_correl) +void compare( + FILE *outfile, + const lattice_map *mapA_to_B, + const HemeLBExtractionFile *A, + const HemeLBExtractionFile *B, + int minexistent, + bool normalize_correl) { uint64_t num_sites_A = A->get_num_sites(); double timeA = A->get_time(); @@ -48,7 +47,7 @@ void compare(FILE *outfile, lattice_map *mapA_to_B, HemeLBExtractionFile *A, Hem // Velocity A->get_velocity(i, &velA); - B->get_interpolated_velocity(&mapA_to_B[i], &velB); + B->get_interpolated_velocity(mapA_to_B[i], &velB); // Max amplitude if(velA.length() > max_vel_A) { @@ -84,7 +83,7 @@ void compare(FILE *outfile, lattice_map *mapA_to_B, HemeLBExtractionFile *A, Hem for(uint64_t i = 0; i < num_sites_A; i++) { // Cross correlation shearstress A->get_shearstress(i, &shearstressA); - B->get_interpolated_shearstress(&mapA_to_B[i], &shearstressB); + B->get_interpolated_shearstress(mapA_to_B[i], &shearstressB); crosscorrelAB += shearstressA * shearstressB; // Autocorrelations @@ -94,8 +93,8 @@ void compare(FILE *outfile, lattice_map *mapA_to_B, HemeLBExtractionFile *A, Hem // L2 shearstressL2 += (shearstressA - shearstressB) * (shearstressA - shearstressB); } - correl_shearstress = crosscorrelAB / (sqrt(autocorrelA) * sqrt(autocorrelB)); - shearstressL2 = sqrt(shearstressL2); + correl_shearstress = crosscorrelAB / (std::sqrt(autocorrelA) * std::sqrt(autocorrelB)); + shearstressL2 = std::sqrt(shearstressL2); } if(do_pressure) { @@ -105,22 +104,30 @@ void compare(FILE *outfile, lattice_map *mapA_to_B, HemeLBExtractionFile *A, Hem for(uint64_t i = 0; i < num_sites_A; i++) { // Pressure A->get_pressure(i, &presA); - B->get_interpolated_pressure(&mapA_to_B[i], &presB); + B->get_interpolated_pressure(mapA_to_B[i], &presB); autocorrelA += presA * presA; autocorrelB += presB * presB; crosscorrelAB += presA * presB; presL2 += (presA - presB) * (presA - presB); } - correl_pres = crosscorrelAB / (sqrt(autocorrelA) * sqrt(autocorrelB)); - presL2 = sqrt(presL2); + correl_pres = crosscorrelAB / (std::sqrt(autocorrelA) * std::sqrt(autocorrelB)); + presL2 = std::sqrt(presL2); } fprintf(outfile, "%f %f %f %f %f %f %f %f %f\n", timeA, timeB, correl_vel, max_vel_A, max_vel_B, correl_shearstress, shearstressL2, correl_pres, presL2); } /** Calculates the difference between each site in A, and the corresponding (trilinearly interpolated) point in B */ -void diff(FILE *outfile, lattice_map *mapA_to_B, HemeLBExtractionFile *A, HemeLBExtractionFile *B, int minexistent, Vector3 *project, bool relativeErr, bool verbose) +void diff( + FILE *outfile, + const lattice_map *mapA_to_B, + const HemeLBExtractionFile *A, + const HemeLBExtractionFile *B, + int minexistent, + Vector3 *project, + bool relativeErr, + bool verbose) { uint64_t num_sites_A = A->get_num_sites(); double timeA = A->get_time(); @@ -132,7 +139,7 @@ void diff(FILE *outfile, lattice_map *mapA_to_B, HemeLBExtractionFile *A, HemeLB if(do_velocity) { - if(verbose == true) { + if(verbose) { fprintf(stderr, "Velocity difference calc\n"); fprintf(outfile, "# timeA=%f timeB=%f\n", timeA, timeB); } @@ -154,7 +161,7 @@ void diff(FILE *outfile, lattice_map *mapA_to_B, HemeLBExtractionFile *A, HemeLB // Get velocity from A and corresponding interpolated velocity from B A->get_velocity(i, &velA); - B->get_interpolated_velocity(&mapA_to_B[i], &velB); + B->get_interpolated_velocity(mapA_to_B[i], &velB); // fprintf(stderr, "velA "); // velA.print(stderr); @@ -187,7 +194,7 @@ void diff(FILE *outfile, lattice_map *mapA_to_B, HemeLBExtractionFile *A, HemeLB } if(do_shearstress) { - if(verbose == true) { + if(verbose) { fprintf(stderr, "WSS difference calc\n"); fprintf(outfile, "# timeA=%f timeB=%f\n", timeA, timeB); } @@ -209,15 +216,16 @@ void diff(FILE *outfile, lattice_map *mapA_to_B, HemeLBExtractionFile *A, HemeLB // Get shearstress from A and corresponding interpolated shearstress from B A->get_shearstress(i, &shearA); - B->get_interpolated_shearstress(&mapA_to_B[i], &shearB); + B->get_interpolated_shearstress(mapA_to_B[i], &shearB); // Print out site coords - A->get_sites()[i].print(outfile, voxelA); + auto sites = A->get_sites(); + sites[i].print(outfile, voxelA); double sheardiff = shearA - shearB; // If requested, give difference relative to local WSS from input file A - if(relativeErr == true) { + if (relativeErr) { sheardiff /= shearA; } @@ -225,5 +233,3 @@ void diff(FILE *outfile, lattice_map *mapA_to_B, HemeLBExtractionFile *A, HemeLB } } } - -#endif diff --git a/src/HemeLBExtractionFile.cc b/src/HemeLBExtractionFile.cc new file mode 100644 index 0000000..51b571a --- /dev/null +++ b/src/HemeLBExtractionFile.cc @@ -0,0 +1,541 @@ +#include "HemeLBExtractionFile.h" + +#include + +/* Magic numbers designating: + * - a heme file, + * - a heme extraction file, + * - and a colloid file + */ +static const uint32_t heme_magic = 0x686C6221; // hlb! +static const uint32_t extract_magic = 0x78747204; // xtr +static const uint32_t colloid_magic = 0x636F6C04; // colloid + + +// Workaround for XDR's broken long reading ability +void xdr_long(XDR *xdrs, uint64_t *ret) +{ + uint32_t half1, half2; + xdr_u_int(xdrs, &half1); + xdr_u_int(xdrs, &half2); + *ret = ((uint64_t)half1)<<32 | half2; // Assume big-endianness +} + +HemeLBExtractionFile::HemeLBExtractionFile( + char *fname, + double step_length_, + double scaling_, + bool verbose, + Vector3 *translate) + : step_length(step_length_), scaling(scaling_), bool_verbose(verbose) +{ + in = fopen(fname, "rb"); + if (in == NULL) { + fprintf(stderr, "Could not open file '%s'. Does it exist?\n", fname); + return; + } + xdrstdio_create(&xdrs, in, XDR_DECODE); + + // Read the file header, populating the 'header' class variable. Determines whether file is "normal" or colloids + int r; + r = read_header(); + if(r != 0) { + fprintf(stderr, "Problem reading header for file: '%s'.\n", fname); + return; + } + + // If file is a colloids file, return now (no more info to be read in) + if (header->is_colloids_file) { + bool_correctly_initialised = true; + return; + } + + /// Otherwise, if file is "normal", read the field header + r = read_field_header(); + if(r != 0) { + fprintf(stderr, "Problem reading field header for file: '%s'.\n", fname); + return; + } + + // If specified, convert lattice translation vector (real units) into grid units (via rounding based on voxel size) + if(translate != nullptr) { + translate_grid_x = (int32_t)(translate->get_x()/header->voxelsz); + translate_grid_y = (int32_t)(translate->get_y()/header->voxelsz); + translate_grid_z = (int32_t)(translate->get_z()/header->voxelsz); + if(bool_verbose == true) { + fprintf(stderr, "Lattice translation vector (%e,%e,%e) in real units, converted to (%d,%d,%d) in grid units (voxel size = %e).\n", translate->get_x(), translate->get_y(), translate->get_z(), translate_grid_x, translate_grid_y, translate_grid_z, header->voxelsz); + } + } + + // Allocate a snapshot of the right dimensions + snapshot = std::make_unique(header.get(), field_header.data()); + + // Get the time of the next snapshot + bool next = read_time_next(); + if(next == false) { + bool_no_more_snapshots = true; + } else { + bool_no_more_snapshots = false; + load_next_snapshot(); // Load in the first snapshot + bool_correctly_initialised = true; + } +} + +HemeLBExtractionFile::~HemeLBExtractionFile() +{ + xdr_destroy(&xdrs); + fclose(in); +} + +bool HemeLBExtractionFile::read_time_next() +{ + // Check if we have reached the end of the file + if (feof(in)) { + if(bool_verbose) { + fprintf(stderr, "# Reached end of file.\n"); + } + return false; + } + + uint64_t step; + xdr_long(&xdrs,&step); + + time_next = step * step_length; + return true; +} + +int HemeLBExtractionFile::load_next_snapshot() +{ + if(bool_no_more_snapshots) { + if(bool_verbose) { + fprintf(stderr, "Warning: No more snapshots left in file.\n"); + } + return 1; + } + + uint32_t gridposx, gridposy, gridposz; + float value; + + snapshot->set_timestep(time_next); + for(unsigned int s = 0; s < header->num_sites; s++) { + // Read in grid position of site + xdr_u_int(&xdrs, &gridposx); + xdr_u_int(&xdrs, &gridposy); + xdr_u_int(&xdrs, &gridposz); + + // Apply any requested translation of this xtr file directly to the site positions + gridposx += translate_grid_x; + gridposy += translate_grid_y; + gridposz += translate_grid_z; + + snapshot->site_set(s, gridposx, gridposy, gridposz); + for(unsigned int i = 0; i < header->num_columns; i++) { + // Read value and add it to the appropriate data column (scaling it by the specified scaling) + xdr_float(&xdrs, &value); + snapshot->column_set_plus_offset(i, s, value * this->scaling); + } + } + + if(!read_time_next()) { + bool_no_more_snapshots = true; + if(bool_verbose) { + fprintf(stderr, "No more snapshots. (Note that final snapshot may have been incomplete)\n"); + } + } + return 0; +} + + + + +void HemeLBExtractionFile::read_and_print_colloids(FILE *outfile) +{ + uint32_t headerLen, recordLen; + uint64_t dsetLen, timeStep; + uint64_t id, rank; + double X, Y, Z, VX, VY, VZ; + + while(!feof(in)) { + + // Reinitialise reading vars to catch partially written header + headerLen = 0; + recordLen = 0; + dsetLen = 0; + timeStep = 0; + + // Read timestep header + xdr_u_int(&xdrs, &headerLen); + xdr_u_int(&xdrs, &recordLen); + xdr_long(&xdrs, &dsetLen); + xdr_long(&xdrs, &timeStep); + + // Check for end of file here (to avoid spurious EOF output) + if (feof(in)) { + break; + } + + if (recordLen <= 0) { + break; + } + + if (bool_verbose) { + fprintf(outfile, "# headerLen: %u recordLen %u dsetLen %ld timeStep: %ld\n", headerLen, recordLen, dsetLen, timeStep); + } + + // Calc. num particles from the record length data + uint32_t num_particles = dsetLen/recordLen; + + for (uint32_t i = 0; i < num_particles; i++) { + // Read particle data + xdr_long(&xdrs, &id); + xdr_long(&xdrs, &rank); + xdr_double(&xdrs, &X); + xdr_double(&xdrs, &Y); + xdr_double(&xdrs, &Z); + xdr_double(&xdrs, &VX); + xdr_double(&xdrs, &VY); + xdr_double(&xdrs, &VZ); + + // Rescale the positions by the scaling factor provided to hemeXtract + X *= this->scaling; + Y *= this->scaling; + Z *= this->scaling; + +// fprintf(outfile, "TIME: %ld ID: %ld RANK: %ld A0: %e Ah: %e X: %e Y: %e Z: %e\n", timeStep, id, rank, A0, Ah, X, Y, Z); + fprintf(outfile, "%ld %ld %ld %.13e %.13e %.13e %.13e %.13e %.13e\n", timeStep, id, rank, X, Y, Z, VX, VY, VZ); + } + } + if (bool_verbose) { + fprintf(stderr, "# Reached end of file.\n"); + } +} + +bool HemeLBExtractionFile::correctly_initialised() const +{ + return bool_correctly_initialised; +} + +bool HemeLBExtractionFile::is_colloids_file() const +{ + return header->is_colloids_file; +} + +bool HemeLBExtractionFile::no_more_snapshots() const +{ + return bool_no_more_snapshots; +} + +uint64_t HemeLBExtractionFile::get_num_sites() const +{ + return header->num_sites; +} + +double HemeLBExtractionFile::get_voxelsz() const +{ + return header->voxelsz; +} + +bool HemeLBExtractionFile::hasVelocity() const +{ + return has_velocity; +} + +bool HemeLBExtractionFile::hasShearStress() const +{ + return has_shearstress; +} + +bool HemeLBExtractionFile::hasPressure() const +{ + return has_pressure; +} + +double HemeLBExtractionFile::get_scaling() const +{ + return scaling; +} + +double HemeLBExtractionFile::get_scalar_quantity(uint32_t column_index, uint64_t site_index) const +{ + return snapshot->get(column_index, site_index); +} + +double HemeLBExtractionFile::get_interpolated_scalar_quantity(uint32_t column_index, const lattice_map &map) const +{ + double average_existing = 0; + int num_existing = 0; + for (int i = 0; i < 8; i++) { + if(map.index[i].exists) { + average_existing += get_scalar_quantity(column_index, map.index[i].index); + num_existing++; + } + } + if (num_existing == 0) { + return 0; + } + average_existing /= num_existing; + + // Get the scalars at the 8 lattice sites forming the cube + double s[8]; + for (int i = 0; i < 8; i++) { + // If there was no corresponding site_index, the site was not measured, so set it to the average + // value of the existing sites in the map + if(map.index[i].exists == false) { + s[i] = average_existing; + } else { + s[i] = get_scalar_quantity(column_index, map.index[i].index); + } + } + + // Trilinear interpolation within the cube + return s[0] * (1 - map.a) * (1 - map.b) * (1 - map.c) + + s[1] * (map.a * (1 - map.b) * (1 - map.c)) + + s[2] * ((1 - map.a) * map.b * (1 - map.c)) + + s[3] * ((1 - map.a) * (1 - map.b) * map.c) + + s[4] * (map.a * (1 - map.b) * map.c) + + s[5] * ((1 - map.a) * map.b * map.c) + + s[6] * (map.a * map.b * (1 - map.c)) + + s[7] * (map.a * map.b * map.c); +} + +void HemeLBExtractionFile::get_vector_quantity(uint32_t column_index, uint64_t site_index, Vector3 *returned_val) const +{ + double vx, vy, vz; + vx = get_scalar_quantity(column_index, site_index); + vy = get_scalar_quantity(column_index + 1, site_index); + vz = get_scalar_quantity(column_index + 2, site_index); + returned_val->set(vx, vy, vz); +} + +void HemeLBExtractionFile::get_interpolated_vector_quantity(uint32_t column_index, const lattice_map &map, Vector3 *returned_val) const +{ + double vx, vy, vz; + vx = get_interpolated_scalar_quantity(column_index, map); + vy = get_interpolated_scalar_quantity(column_index + 1, map); + vz = get_interpolated_scalar_quantity(column_index + 2, map); + returned_val->set(vx, vy, vz); +} + +void HemeLBExtractionFile::get_pressure(uint64_t site_index, double *returned_val) const +{ + *returned_val = get_scalar_quantity(column_pressure, site_index); +} + +void HemeLBExtractionFile::get_velocity(uint64_t site_index, Vector3 *returned_val) const +{ + get_vector_quantity(column_velocity, site_index, returned_val); +} + +void HemeLBExtractionFile::get_shearstress(uint64_t site_index, double *returned_val) const +{ + *returned_val = get_scalar_quantity(column_shearstress, site_index); +} + +void HemeLBExtractionFile::get_interpolated_pressure(const lattice_map &map, double *returned_val) const +{ + *returned_val = get_interpolated_scalar_quantity(column_pressure, map); +} + +void HemeLBExtractionFile::get_interpolated_velocity(const lattice_map &map, Vector3 *returned_val) const +{ + get_interpolated_vector_quantity(column_velocity, map, returned_val); +} + +void HemeLBExtractionFile::get_interpolated_shearstress(const lattice_map &map, double *returned_val) const +{ + *returned_val = get_interpolated_scalar_quantity(column_shearstress, map); +} + +const Site * HemeLBExtractionFile::get_sites() const +{ + return snapshot->get_sites(); +} + +SiteIndex * HemeLBExtractionFile::get_site_indices(Site *list, uint64_t list_size) +{ + return snapshot->get_site_indices(list, list_size); +} + +SiteIndex * HemeLBExtractionFile::get_site_indices_hashed_lookup(Site *list, uint64_t list_size) +{ + return snapshot->get_site_indices_hashed_lookup(list, list_size, bool_verbose); +} + +double HemeLBExtractionFile::get_time() const +{ + return snapshot->get_timestep(); +} + +double HemeLBExtractionFile::get_time_next() const +{ + return time_next; +} + +void HemeLBExtractionFile::print_header(FILE *outfile) const +{ + fprintf(outfile, "\n# HEADER:\n# -------\n"); + fprintf(outfile, "# version=%u\n", header->version); + fprintf(outfile, "# voxelsz=%f\n", header->voxelsz); + fprintf(outfile, "# originx=%f\n", header->originx); + fprintf(outfile, "# originy=%f\n", header->originy); + fprintf(outfile, "# originz=%f\n", header->originz); + fprintf(outfile, "# num_sites=%lu\n", header->num_sites); + fprintf(outfile, "# field_count=%u\n", header->field_count); + fprintf(outfile, "# field_header_length=%u\n", header->field_header_length); +} + +void HemeLBExtractionFile::print_field_header(FILE *outfile) const +{ + fprintf(outfile, "\n# FIELD HEADERS:\n# -------\n"); + for(unsigned int i = 0; i < header->field_count; i++) { + fprintf(outfile, "# fieldname=%s\n# num_floats=%u\n# offset=%lf\n# ---\n", field_header[i].name, field_header[i].num_floats, field_header[i].offset); + } + fprintf(outfile, "# Number of 'columns' per snapshot = %d\n", header->num_columns); +} + +void HemeLBExtractionFile::print_column_headings(FILE *outfile) const +{ + fprintf(outfile, "# step | grid_x | grid_y | grid_z"); + for(unsigned int i = 0; i < header->field_count; i++) { + if(field_header[i].num_floats == 1) { + fprintf(outfile, " | %s", field_header[i].name); + } else { + for(unsigned int j = 0; j < field_header[i].num_floats; j++) { + fprintf(outfile, " | %s(%d)", field_header[i].name, j); + } + } + } + fprintf(outfile, "\n"); +} +void HemeLBExtractionFile::print_stats_column_headings(FILE *outfile) const +{ + fprintf(outfile, "# step"); + for(unsigned int i = 0; i < header->field_count; i++) { + if(field_header[i].num_floats == 1) { + fprintf(outfile, " | %s [average | stdev | max | min]", field_header[i].name); + } else { + for(unsigned int j = 0; j < field_header[i].num_floats; j++) { + fprintf(outfile, " | %s(%d) [average | stdev | max | min]", field_header[i].name, j); + } + } + } + fprintf(outfile, "\n"); +} + +void HemeLBExtractionFile::print_all(FILE *outfile) const +{ + snapshot->print(outfile); +} + +void HemeLBExtractionFile::print_stats(FILE *outfile) const +{ + snapshot->print_stats(outfile); +} + +/** Reads a HemeLB extraction file header. Checks that the magic numbers are correct. */ +int HemeLBExtractionFile::read_header() +{ + uint32_t magic1 = 0, magic2 = 0; + uint32_t num_sites1, num_sites2; + + // Read heme magic number + xdr_u_int(&xdrs, &magic1); + if(magic1 != heme_magic) { + fprintf(stderr, "First uint32_t does not match heme magic number.\n"); + return 1; + } + + header = std::make_unique
(); + + // Read extraction magic number + xdr_u_int(&xdrs, &magic2); + if(magic2 == extract_magic) { + if(bool_verbose) { + fprintf(stderr, "Reading a normal extraction file.\n"); + } + + header->is_colloids_file = false; + + xdr_u_int(&xdrs, &(header->version)); + xdr_double(&xdrs, &(header->voxelsz)); + xdr_double(&xdrs, &(header->originx)); + xdr_double(&xdrs, &(header->originy)); + xdr_double(&xdrs, &(header->originz)); + + // There is something wrong with either Heme or XDR which means that xdr_u_long is not working (possibly even reading more/less than 8 bytes...) + xdr_u_int(&xdrs, &num_sites1); + xdr_u_int(&xdrs, &num_sites2); + + // Combine the two ints into a long (assuming bigendian-ness) + header->num_sites = ((uint64_t)num_sites1)<<32 | num_sites2; + + xdr_u_int(&xdrs, &(header->field_count)); + xdr_u_int(&xdrs, &(header->field_header_length)); + + } else if (magic2 == colloid_magic) { + if(bool_verbose) { + fprintf(stderr, "Reading a colloids extraction file.\n"); + } + header->is_colloids_file = true; + + // Check version + uint32_t version=0; + xdr_u_int(&xdrs, &version); + if(bool_verbose) { + printf("Colloids version %u\n", version); + } + } else { + fprintf(stderr, "Unknown format: Second uint32_t does not match extraction file or colloids file magic number.\n"); + return 1; + } + + // Everything went well + return 0; +} + +/** Reads the field headers and remembers which columns correspond to velocity, pressure and shearstress. */ +int HemeLBExtractionFile::read_field_header() +{ + if(bool_verbose) { + printf("Reading field header...\n"); + } + int const maxlength = 100; // surely it can't be bigger than this...? + + field_header.resize(header->field_count); + + for(unsigned int i = 0; i < header->field_count; i++) { + field_header[i].name = new char[maxlength]; + xdr_string(&xdrs, &(field_header[i].name), maxlength); + xdr_u_int(&xdrs, &(field_header[i].num_floats)); + xdr_double(&xdrs, &(field_header[i].offset)); + } + if(bool_verbose) { + printf("Finished reading field header...\n"); + } + + // Calculate number of 'columns' needed to represent a snapshot + uint32_t num_columns = 0; + for(unsigned int i = 0; i < header->field_count; i++) { + // If file contains velocity, pressure or shearstress, remember the corresponding column index for that field type. + if(strcmp(field_header[i].name, "velocity") == 0) { + has_velocity = true; + column_velocity = num_columns; + } else if(strcmp(field_header[i].name, "pressure") == 0) { + has_pressure = true; + column_pressure = num_columns; + } else if(strcmp(field_header[i].name, "shearstress") == 0) { + has_shearstress = true; + column_shearstress = num_columns; + } else if(strcmp(field_header[i].name, "d_shearstress") == 0) { // For "legacy" reasons only. Some old MCA runs have this field. + has_shearstress = true; + column_shearstress = num_columns; + } + + // This field has 'num_floats' columns to it + num_columns += field_header[i].num_floats; + } + header->num_columns = num_columns; + + // Everything went well + return 0; +} diff --git a/src/Mapping.cc b/src/Mapping.cc new file mode 100644 index 0000000..78c5fc0 --- /dev/null +++ b/src/Mapping.cc @@ -0,0 +1,66 @@ +#include "Mapping.h" + +#include +#include +#include + +inline void map( + double old_voxelsz, + double new_voxelsz, + uint32_t in_coord, + uint32_t &out_root, + double &out_interpolate) +{ + double real_pos = in_coord * old_voxelsz; + out_root = std::floor(real_pos/new_voxelsz); + out_interpolate = (real_pos - out_root * new_voxelsz) / new_voxelsz; +} + +lattice_map * get_mapping(HemeLBExtractionFile *A, HemeLBExtractionFile *B) +{ + // Work out mapping between the two lattices + uint64_t num_sites_A = A->get_num_sites(); + const Site *sitesA = A->get_sites(); + double voxelA = A->get_voxelsz(); + double voxelB = B->get_voxelsz(); + uint32_t mapped_x, mapped_y, mapped_z; + lattice_map *mapping_A_B = new lattice_map[num_sites_A]; + if(mapping_A_B == NULL) { + fprintf(stderr, "Could not allocate memory for lattice_map\n"); + return NULL; + } + //Site *mapped_sites = new Site[num_sites_A * 8]; + //if(mapped_sites == NULL) { + // fprintf(stderr, "Could not allocate memory for mapped_sites\n"); + // return NULL; + //} + std::vector mapped_sites(num_sites_A * 8); + + for(uint64_t i = 0; i < num_sites_A; i++) { + // Get the interpolation constants for this site in terms of the other lattice + map(voxelA, voxelB, sitesA[i].x, mapped_x, mapping_A_B[i].a); + map(voxelA, voxelB, sitesA[i].y, mapped_y, mapping_A_B[i].b); + map(voxelA, voxelB, sitesA[i].z, mapped_z, mapping_A_B[i].c); + + // Add the 8 sites forming the cube this point lies within + mapped_sites[i * 8 ].set(mapped_x , mapped_y , mapped_z ); // (000) root + mapped_sites[i * 8 + 1].set(mapped_x + 1, mapped_y , mapped_z ); // (100) + mapped_sites[i * 8 + 2].set(mapped_x , mapped_y + 1, mapped_z ); // (010) + mapped_sites[i * 8 + 3].set(mapped_x , mapped_y , mapped_z + 1); // (001) + mapped_sites[i * 8 + 4].set(mapped_x + 1, mapped_y , mapped_z + 1); // (101) + mapped_sites[i * 8 + 5].set(mapped_x , mapped_y + 1, mapped_z + 1); // (011) + mapped_sites[i * 8 + 6].set(mapped_x + 1, mapped_y + 1, mapped_z ); // (110) + mapped_sites[i * 8 + 7].set(mapped_x + 1, mapped_y + 1, mapped_z + 1); // (111) + } + + SiteIndex *site_indices = B->get_site_indices_hashed_lookup( + mapped_sites.data(), num_sites_A * 8); + + for(uint64_t i = 0; i < num_sites_A; i++) { + for(int j = 0; j < 8; j++) { + mapping_A_B[i].index[j] = site_indices[i * 8 + j]; + } + } + + return mapping_A_B; +} diff --git a/src/Snapshot.cc b/src/Snapshot.cc new file mode 100644 index 0000000..521794c --- /dev/null +++ b/src/Snapshot.cc @@ -0,0 +1,209 @@ +#include "Snapshot.h" + +#include +#include +#include +#include +#include +#include +#include +#include + +// Calculate statistics (helper function) +inline void calc_stats( + const std::vector &records, + double &average, double &stdev, double &max, double &min) +{ + double sum{0}, sum2{0}; + max = -INFINITY; + min = INFINITY; + for(const double& rec : records) { + sum += rec; + sum2 += rec * rec; + if(max < rec) max = rec; + if(min > rec) min = rec; + } + average = sum/records.size(); + stdev = std::sqrt(sum2/records.size() - average*average); +} + +// +// Column (helper class) +// +class Column +{ + public: + Column(unsigned int num_records_, double offset_) + : records(num_records_), offset(offset_) {} + + double get(uint64_t index) const + { + if(index >= records.size()) { + fprintf(stderr, "Error: index (%lu) exceeds number of records in column (%lu)\n", index, records.size()); + exit(1); + } + return records[index]; + } + void set_plus_offset(uint64_t index, float value) + { + if(index >= records.size()) { + fprintf(stderr, "Error: index (%lu) exceeds number of records in column (%lu)\n", index, records.size()); + exit(1); + } + records[index] = value + offset; + stats_calcd = false; + } + void print(FILE *outfile, uint64_t index) { + if(index >= records.size()) { + fprintf(stderr, "Error: index (%lu) exceeds number of records in column (%lu)\n", index, records.size()); + exit(1); + } + fprintf(outfile, "%e", records[index]); + } + void print_stats(FILE *outfile) { + if (!stats_calcd) { + calc_stats(records,average,stdev,max,min); + stats_calcd = true; + } + fprintf(outfile, "%e %e %e %e", average, stdev, max, min); + } + + private: + std::vector records; + double offset; + bool stats_calcd{false}; + double average; + double stdev; + double max, min; +}; + +Snapshot::Snapshot(const HEADER *header, const FIELD_HEADER *field_header) + : num_sites(header->num_sites), + num_columns(header->num_columns), + voxelsz(header->voxelsz), + sites(num_sites) +{ + columns = new Column*[num_columns]; + unsigned int c = 0; + for(unsigned int i = 0; i < header->field_count; i++) { + for(unsigned int bro=0; bro < field_header[i].num_floats; bro++) { + columns[c] = new Column(num_sites, field_header[i].offset); + c++; + } + } +} + +void Snapshot::set_timestep(double timestep) +{ + this->timestep = timestep; +} +double Snapshot::get_timestep() const +{ + return timestep; +} + +void Snapshot::site_set(unsigned int site_index, uint32_t x, uint32_t y, uint32_t z) +{ + sites[site_index].set(x, y, z); +} + +void Snapshot::column_set_plus_offset(uint32_t column_index, uint64_t site_index, double value) +{ + if(column_index >= num_columns) { + fprintf(stderr, "Error: column index (%u) exceeds number of columns (%u)\n", column_index, num_columns); + exit(1); + } + columns[column_index]->set_plus_offset(site_index, value); +} + +double Snapshot::get(uint32_t column_index, uint64_t site_index) const +{ + return columns[column_index]->get(site_index); +} + +bool Snapshot::get_site_index(const Site &s, uint64_t &site_index) const +{ + auto res = std::find_if(sites.begin(),sites.end(), + [=](Site st){ return st.equals(s); }); + + if (res != sites.end()) { + site_index = std::distance(sites.begin(),res); + return true; + } else { + // If no such site can be found, return false + return false; + } +} + +const Site *Snapshot::get_sites() const +{ + return sites.data(); +} + +/* Very inefficient calculation of site_index corresponding to the given site. + * Only suitable for small data sets. */ +SiteIndex * Snapshot::get_site_indices(const Site *list, uint64_t list_size) +{ + SiteIndex *indices = new SiteIndex[list_size]; + for(uint64_t i = 0; i < list_size; i++) { + indices[i].exists = get_site_index(list[i], indices[i].index); + } + return indices; +} + +/* Builds hashtable to get indices corresponding to the given site. + * Necessary for very large data sets. */ +SiteIndex * Snapshot::get_site_indices_hashed_lookup(const Site *list, uint64_t list_size, bool bool_verbose) +{ + // Build the hash table + if(bool_verbose == true) { fprintf(stderr, "# Building hashtable...\n"); } + std::unordered_map hashtable; + std::ostringstream oss; + for(uint64_t i = 0; i < num_sites; i++) { + oss << sites[i].x << "," << sites[i].y << "," << sites[i].z; + std::string hashstr = oss.str(); + hashtable[hashstr] = i; + } + if(bool_verbose == true) { fprintf(stderr, "# ...done building hashtable\n"); } + + if(bool_verbose == true) { fprintf(stderr, "# Calculating mapping...\n"); } + SiteIndex *indices = new SiteIndex[list_size]; + for(uint64_t i = 0; i < list_size; i++) { + oss << list[i].x << "," << list[i].y << "," << list[i].z; + std::string hashstr = oss.str(); + if(hashtable.find(hashstr) != hashtable.end()) { + indices[i].exists = true; + indices[i].index = hashtable[hashstr]; + } else { + indices[i].exists = false; + } + } + if(bool_verbose == true) { fprintf(stderr, "# ... done calculating mapping\n"); } + + return indices; +} + +void Snapshot::print(FILE *outfile) +{ + for(uint64_t s = 0; s < num_sites; s++) { + fprintf(outfile, "%f ", timestep); + sites[s].print(outfile, this->voxelsz); + for(uint32_t c = 0; c < num_columns; c++) { + fprintf(outfile, " "); + columns[c]->print(outfile, s); + } + fprintf(outfile, "\n"); + } + fprintf(outfile, "\n"); +} + +void Snapshot::print_stats(FILE *outfile) +{ + fprintf(outfile, "%f", timestep); + for(uint32_t c = 0; c < num_columns; c++) { + fprintf(outfile, " "); + columns[c]->print_stats(outfile); + } + fprintf(outfile, "\n"); +} + diff --git a/hemeXtract.cc b/src/hemeXtract.cc similarity index 97% rename from hemeXtract.cc rename to src/hemeXtract.cc index 1fb2999..2e29617 100644 --- a/hemeXtract.cc +++ b/src/hemeXtract.cc @@ -1,16 +1,14 @@ -#include -#include -#include -#include -#include +#include #include -#include -#include #include +#include + +#include -#include "Mapping.h" -#include "Compare.h" -#include "HemeLBExtractionFile.h" +#include "Vector3.h" +#include "HemeLBExtractionFile.h" // HemeLBExtractionFile +#include "Compare.h" // compare, diff +#include "Mapping.h" // get_mapping /* Command line parsing stuff */ const char *argp_program_version = "hemeXtract 0.0.1"; @@ -54,6 +52,7 @@ struct arguments { bool relativeErr; bool normalize_correl; }; + static Vector3 * parse_vector3(char* arg) { std::stringstream ss(arg); std::vector result; @@ -68,6 +67,7 @@ static Vector3 * parse_vector3(char* arg) { } return new Vector3(atof(result[0].c_str()), atof(result[1].c_str()), atof(result[2].c_str())); } + static error_t parse_opt(int key, char *arg, struct argp_state *state) { struct arguments *arggs = reinterpret_cast(state->input); @@ -133,10 +133,11 @@ static error_t parse_opt(int key, char *arg, struct argp_state *state) { return 0; } + static struct argp argp = { options, parse_opt, args_doc, doc, 0, 0, 0 }; /** Simple pointer swapper */ -void swap(HemeLBExtractionFile **i, HemeLBExtractionFile **j) { +inline void swap(HemeLBExtractionFile **i, HemeLBExtractionFile **j) { HemeLBExtractionFile *temp = *i; *i = *j; *j = temp;