diff --git a/Makefile.am b/Makefile.am index c71cc8d00c797f0e2afc034cb1abfff7eba14c88..40ba64dcdd1c7270712288bce938ab56e918694d 100644 --- a/Makefile.am +++ b/Makefile.am @@ -23,6 +23,9 @@ SUBDIRS = src argparse examples doc tests tools if HAVEEAGLECOOLING SUBDIRS += examples/Cooling/CoolingRates endif +if HAVELOGGER +SUBDIRS += logger +endif # Non-standard files that should be part of the distribution. EXTRA_DIST = INSTALL.swift .clang-format format.sh diff --git a/configure.ac b/configure.ac index 338edec60f956c37f666f7592a931b2c20a9f6e8..ba119835680f3679809f3d4d3f0133f74ea04454 100644 --- a/configure.ac +++ b/configure.ac @@ -86,6 +86,7 @@ AC_ARG_ENABLE([logger], if test "$with_logger" = "yes"; then AC_DEFINE([WITH_LOGGER], 1, [logger enabled]) fi +AM_CONDITIONAL([HAVELOGGER],[test $with_logger = "yes"]) # Interprocedural optimization support. Needs special handling for linking and # archiving as well as compilation with Intels, needs to be done before @@ -996,6 +997,50 @@ fi AC_SUBST([TBBMALLOC_LIBS]) AM_CONDITIONAL([HAVETBBMALLOC],[test -n "$TBBMALLOC_LIBS"]) +# Check for python. +have_python="no" +AC_ARG_WITH([python], + [AS_HELP_STRING([--with-python=PATH], + [root directory where python is installed @<:@yes/no@:>@] + )], + [with_python="$withval"], + [with_python="no"] +) +if test "x$with_python" != "xno"; then + AC_CACHE_CHECK( + "python version", + ac_cv_ver_python, + [ac_cv_ver_python=`python -c 'import sys;print(sys.version[[:3]])' 2> /dev/null`]) + + if test "x$with_python" != "xyes" -a "x$with_python" != "x"; then + PYTHON_LIBS="" + PYTHON_INCS="-I$with_python/include/python$ac_cv_ver_python" + else + PYTHON_LIBS="" + PYTHON_INCS="" + fi + have_python="yes" + + AC_CHECK_PROGS( + [PYTHON_BIN], + [python$ac_cv_ver_python], + [AC_MSG_ERROR(Cannot find python binary!)], + [$with_python/bin] + ) + + AC_CHECK_LIB( + [python${ac_cv_ver_python}m], + [PyArg_ParseTuple], + [AC_DEFINE([HAVE_PYTHON],1,[The python library appears to be present.]) + ], + [AC_MSG_ERROR(Cannot find python library!)] + ) +fi +AC_SUBST([PYTHON_LIBS]) +AC_SUBST([PYTHON_INCS]) +AM_CONDITIONAL([HAVEPYTHON],[test -n "$PYTHON_LIBS"]) + + # Check for HDF5. This is required. AX_LIB_HDF5 if test "$with_hdf5" != "yes"; then @@ -1991,7 +2036,7 @@ AM_CONDITIONAL([HAVEEAGLEFEEDBACK], [test $with_feedback = "EAGLE"]) # Handle .in files. AC_CONFIG_FILES([Makefile src/Makefile examples/Makefile examples/Cooling/CoolingRates/Makefile doc/Makefile doc/Doxyfile tests/Makefile]) -AC_CONFIG_FILES([argparse/Makefile tools/Makefile]) +AC_CONFIG_FILES([argparse/Makefile tools/Makefile logger/Makefile]) AC_CONFIG_FILES([tests/testReading.sh], [chmod +x tests/testReading.sh]) AC_CONFIG_FILES([tests/testActivePair.sh], [chmod +x tests/testActivePair.sh]) AC_CONFIG_FILES([tests/test27cells.sh], [chmod +x tests/test27cells.sh]) @@ -2046,6 +2091,7 @@ AC_MSG_RESULT([ VELOCIraptor enabled : $have_velociraptor Particle Logger : $with_logger FoF activated: : $enable_fof + Python enabled : $have_python Hydro scheme : $with_hydro Dimensionality : $with_dimension diff --git a/logger/Makefile.am b/logger/Makefile.am new file mode 100644 index 0000000000000000000000000000000000000000..d2cc4a8a72f4aa27cdda7b9809c90f727a461ba8 --- /dev/null +++ b/logger/Makefile.am @@ -0,0 +1,92 @@ +# This file is part of SWIFT. +# Copyright (c) 2012 Pedro Gonnet (pedro.gonnet@durham.ac.uk), +# Matthieu Schaller (matthieu.schaller@durham.ac.uk). +# +# This program is free software: you can redistribute it and/or modify +# it under the terms of the GNU General Public License as published by +# the Free Software Foundation, either version 3 of the License, or +# (at your option) any later version. +# +# This program is distributed in the hope that it will be useful, +# but WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +# GNU General Public License for more details. +# +# You should have received a copy of the GNU General Public License +# along with this program. If not, see <http://www.gnu.org/licenses/>. + +# Add the non-standard paths to the included library headers +AM_CFLAGS = $(PYTHON_INCS) + +# Assign a "safe" version number +AM_LDFLAGS = -version-info 0:0:0 + +# The git command, if available. +GIT_CMD = @GIT_CMD@ + +# Additional dependencies for shared libraries. +EXTRA_LIBS = $(PROFILER_LIBS) $(TCMALLOC_LIBS) $(JEMALLOC_LIBS) $(TBBMALLOC_LIBS) $(PYTHON_LIBS) + +# MPI libraries. +# MPI_LIBS = $(MPI_THREAD_LIBS) +# MPI_FLAGS = -DWITH_MPI + +# Build the libswiftsim library +lib_LTLIBRARIES = libswiftlogger.la +# Build a MPI-enabled version too? +# if HAVEMPI +# lib_LTLIBRARIES += libswiftlogger_mpi.la +# endif + +# List required headers +include_HEADERS = header.h io.h particle.h timeline.h tools.h + +# Common source files +AM_SOURCES = header.c io.c logger_loader.c particle.c \ + timeline.c tools.c + +# Include files for distribution, not installation. +nobase_noinst_HEADERS = + +# Sources and flags for regular library +libswiftlogger_la_SOURCES = $(AM_SOURCES) +libswiftlogger_la_CFLAGS = $(AM_CFLAGS) +libswiftlogger_la_LDFLAGS = $(AM_LDFLAGS) $(EXTRA_LIBS) + +# Sources and flags for MPI library +# libswiftlogger_mpi_la_SOURCES = $(AM_SOURCES) +# libswiftlogger_mpi_la_CFLAGS = $(AM_CFLAGS) $(MPI_FLAGS) +# libswiftlogger_mpi_la_LDFLAGS = $(AM_LDFLAGS) $(MPI_LIBS) $(EXTRA_LIBS) +# libswiftlogger_mpi_la_SHORTNAME = mpi +# libswiftlogger_mpi_la_LIBADD = + + +# Versioning. If any sources change then update the version_string.h file with +# the current git revision and package version. +# May have a checkout without a version_string.h file and no git command (tar/zip +# download), allow that, but make sure we know it. +version_string.h: ../src/version_string.h.in $(AM_SOURCES) $(include_HEADERS) $(noinst_HEADERS) + if test "X$(GIT_CMD)" != "X"; then \ + GIT_REVISION=`$(GIT_CMD) describe --abbrev=8 --always --tags --dirty`; \ + GIT_BRANCH=`$(GIT_CMD) branch | sed -n 's/^\* \(.*\)/\1/p'`; \ + GIT_DATE=`$(GIT_CMD) log -1 --format=%ci`; \ + sed -e "s,@PACKAGE_VERSION\@,$(PACKAGE_VERSION)," \ + -e "s,@GIT_REVISION\@,$${GIT_REVISION}," \ + -e "s|@GIT_BRANCH\@|$${GIT_BRANCH}|" \ + -e "s|@GIT_DATE\@|$${GIT_DATE}|" \ + -e "s|@SWIFT_CFLAGS\@|$(CFLAGS)|" $< > version_string.h; \ + else \ + if test ! -f version_string.h; then \ + sed -e "s,@PACKAGE_VERSION\@,$(PACKAGE_VERSION)," \ + -e "s,@GIT_REVISION\@,unknown," \ + -e "s,@GIT_BRANCH\@,unknown," \ + -e "s,@GIT_DATE\@,unknown," \ + -e "s|@SWIFT_CFLAGS\@|$(CFLAGS)|" $< > version_string.h; \ + fi; \ + fi + +# Make sure version_string.h is built first. +BUILT_SOURCES = version_string.h + +# And distribute the built files. +EXTRA_DIST = version_string.h ../src/version_string.h.in diff --git a/logger/header.c b/logger/header.c new file mode 100644 index 0000000000000000000000000000000000000000..6edfd6d5e0e97a9429f83969552a4aa35c891099 --- /dev/null +++ b/logger/header.c @@ -0,0 +1,138 @@ +#include "header.h" + +#include "io.h" +#include "tools.h" + +#include <Python.h> +#include <stdio.h> +#include <stdlib.h> + +void header_print(const struct header *h) { +#ifdef SWIFT_DEBUG_CHECKS + printf("Debug checks enabled\n"); +#endif + printf("Version: %s\n", h->version); + printf("First Offset: %lu\n", h->offset_first); + char direction[20]; + if (h->forward_offset) + strcpy(direction, "Forward"); + else + strcpy(direction, "Backward"); + printf("Offset direction: %s\n", direction); + printf("Number masks: %lu\n", h->nber_mask); + + for (size_t i = 0; i < h->nber_mask; i++) { + printf("\tMask: %s\n", h->masks_name[i]); + printf("\tValue: %lu\n", h->masks[i]); + printf("\tSize: %lu\n", h->masks_size[i]); + printf("\n"); + } + + /* mask contains... TODO */ + +}; + +void header_free(struct header *h) { + for (size_t i = 0; i < h->nber_mask; i++) { + free(h->masks_name[i]); + } + free(h->masks_name); + free(h->masks); + free(h->masks_size); +}; + +int header_is_present_and_get_index(const struct header *h, const char *field, + size_t *ind) { + for (size_t i = 0; i < h->nber_mask; i++) { + if (strcmp(h->masks_name[i], field) == 0) { + if (ind != NULL) { + *ind = i; + } + return 1; + } + } + + return 0; +}; + +int header_is_present(const struct header *h, const char *field) { + return header_is_present_and_get_index(h, field, NULL); +}; + +int header_change_offset_direction(struct header *h, void *map) { + h->forward_offset = !h->forward_offset; + size_t offset = LOGGER_VERSION_SIZE; + + return io_write_data(map, LOGGER_NBER_SIZE, &h->forward_offset, + &offset); +} + +int header_read(struct header *h, void *map) { + size_t offset = 0; + + /* read version */ + io_read_data(map, LOGGER_VERSION_SIZE, &h->version, &offset); + + /* read offset direction */ + h->forward_offset = 0; + io_read_data(map, LOGGER_NBER_SIZE, &h->forward_offset, &offset); + + if (h->forward_offset != 0 && h->forward_offset != 1) + error(EIO, "Non boolean value for the offset direction (%i)", h->forward_offset); + + /* read offset to first data */ + h->offset_first = 0; + io_read_data(map, LOGGER_OFFSET_SIZE, &h->offset_first, &offset); + + /* read name size */ + h->name = 0; + io_read_data(map, LOGGER_NBER_SIZE, &h->name, &offset); + + /* check if value defined in this file is large enough */ + if (STRING_SIZE < h->name) { + error(EOVERFLOW, "Name too large in dump file"); + } + + /* read number of masks */ + h->nber_mask = 0; + io_read_data(map, LOGGER_NBER_SIZE, &h->nber_mask, &offset); + + /* allocate memory */ + h->masks = malloc(sizeof(size_t) * h->nber_mask); + h->masks_name = malloc(sizeof(void *) * h->nber_mask); + h->masks_size = malloc(sizeof(size_t) * h->nber_mask); + + /* loop over all masks */ + for (size_t i = 0; i < h->nber_mask; i++) { + /* read mask name */ + h->masks_name[i] = malloc(h->name); + io_read_data(map, h->name, h->masks_name[i], &offset); + + /* get mask value */ + h->masks[i] = 1 << i; + + /* read mask data size */ + h->masks_size[i] = 0; + io_read_data(map, LOGGER_NBER_SIZE, &h->masks_size[i], &offset); + } + + if (offset != h->offset_first) { +#ifdef SWIFT_DEBUG_CHECKS + header_print(h); +#endif + error(EIO, "Wrong header size (in header %li, current %li)", + h->offset_first, offset); + } + + return 0; +}; + +size_t header_get_mask_size(const struct header *h, const size_t mask) { + size_t count = 0; + for (size_t i = 0; i < h->nber_mask; i++) { + if (mask & h->masks[i]) { + count += h->masks_size[i]; + } + } + return count; +} diff --git a/logger/header.h b/logger/header.h new file mode 100644 index 0000000000000000000000000000000000000000..c2e89850418b94eb4d05381e698534a3719652c9 --- /dev/null +++ b/logger/header.h @@ -0,0 +1,95 @@ +#ifndef __LOGGER_HEADER_H__ +#define __LOGGER_HEADER_H__ + +#include "tools.h" + +#include <Python.h> +#include <stdio.h> +#include <stdlib.h> + +#define LOGGER_VERSION_SIZE 20 +#define LOGGER_NBER_SIZE 4 +#define LOGGER_OFFSET_SIZE 7 +#define LOGGER_MASK_SIZE 1 + +struct header { + /* Logger version */ + char version[STRING_SIZE]; + + /* offset of the first header */ + size_t offset_first; + + /* Number of bytes for names */ + size_t name; + + /* number of masks */ + size_t nber_mask; + + /* list of masks */ + size_t *masks; + + /* list of mask name */ + char **masks_name; + + /* size of data in a mask */ + size_t *masks_size; + + /* offset direction */ + int forward_offset; +}; + +/** + * @brief print a header struct + */ +void header_print(const struct header *h); + +/** + * @brief free allocated memory + */ +void header_free(struct header *h); + +/** + * @brief check if field is present in header + * + * @param field name of the requested field + * @param ind (return value) indice of the requested field + */ +int header_is_present_and_get_index(const struct header *h, const char *field, + size_t *ind); + +/** + * @brief check if field is present in header + * + * @param field name of the requested field + */ +int header_is_present(const struct header *h, const char *field); + +/** + * @brief read the logger header + * + * @param h out: header + * @param map file mapping + */ +int header_read(struct header *h, void *map); + +/** + * @brief count number of bits in a given mask + * + * @param h #header file structure + * @param mask mask to compute + * + * @return number of bits in mask + */ +size_t header_get_mask_size(const struct header *h, const size_t mask); + +/** + * @brief Inverse the offset direction + * + * @param h #header file structure + * @param map file mapping + * + * @return error code + */ +int header_change_offset_direction(struct header *h, void *map); + +#endif // __LOGGER_HEADER_H__ diff --git a/logger/io.c b/logger/io.c new file mode 100644 index 0000000000000000000000000000000000000000..9ed62e7fc905e816199d17122aa993fa4e27785d --- /dev/null +++ b/logger/io.c @@ -0,0 +1,84 @@ +#include "io.h" +#include "header.h" +#include "tools.h" + +/* file size */ +#include <sys/stat.h> +/* mapping */ +#include <sys/mman.h> +/* open */ +#include <fcntl.h> + +int io_get_file_size(int fd, size_t *size) { + struct stat s; + int status = fstat(fd, &s); + if (status != 0) error(errno, "Unable to get file size"); + *size = s.st_size; + + return 0; +} + +int io_open_file(char *filename, int *fd, void **map) { + /* open file */ + *fd = open(filename, O_RDWR); + if (*fd == -1) error(errno, "Unable to open file %s", filename); + + /* get file size */ + size_t size; + int status = io_get_file_size(*fd, &size); + if (status != 0) return status; + + /* map memory */ + *map = mmap(NULL, size, PROT_WRITE | PROT_READ, MAP_SHARED, *fd, 0); + if (map == MAP_FAILED) + error(errno, "Failed to allocate map of size %zi bytes.", size); + + return 0; +} + +int io_close_file(int *fd, void **map) { + /* get file size */ + size_t size; + int status = io_get_file_size(*fd, &size); + if (status != 0) return status; + + /* unmap */ + if (munmap(*map, size) != 0) error(errno, "Unable to unmap the file"); + + close(*fd); + + return 0; +} + +int io_read_mask(const struct header *h, void *map, size_t *offset, + size_t *mask, size_t *diff_offset) { + /* read mask */ + if (mask) { + *mask = 0; + memcpy(mask, map + *offset, LOGGER_MASK_SIZE); + } + *offset += LOGGER_MASK_SIZE; + + /* read offset */ + if (diff_offset) { + *diff_offset = 0; + memcpy(diff_offset, map + *offset, LOGGER_OFFSET_SIZE); + } + *offset += LOGGER_OFFSET_SIZE; + + return 0; +} + +int io_read_data(void *map, const size_t size, void *p, size_t *offset) { + memcpy(p, map + *offset, size); + *offset += size; + + return 0; +}; + +int io_write_data(void *map, const size_t size, const void *p, size_t *offset) { + memcpy(map + *offset, p, size); + *offset += size; + + return 0; +}; diff --git a/logger/io.h b/logger/io.h new file mode 100644 index 0000000000000000000000000000000000000000..da4cf34e5367372e3c641f1081a280972faecd66 --- /dev/null +++ b/logger/io.h @@ -0,0 +1,80 @@ +#ifndef __SWIFT_LOGGER_IO_H__ +#define __SWIFT_LOGGER_IO_H__ + +#include "header.h" +#include "tools.h" + +#include <stdio.h> +#include <stdlib.h> + +/** + * @brief get the size of a file + * + * @param fd file id + * @param size out: file size + * + * @return error code + */ +int io_get_file_size(int fd, size_t *size); + +/** + * @brief Open a file and map it + * + * @param filename file to read + * @param fd out: file id + * @param map out: file mapping + * + * @return error code + */ +int io_open_file(char *filename, int *fd, void **map); + +/** + * @brief Close a file and unmap it + * + * @param fd file id + * @param map file mapping + * + * @return error code + */ +int io_close_file(int *fd, void **map); + +/** + * @brief read a single value in a file + * + * @param map file mapping + * @param size size of the chunk to read + * @param p pointer where to store the data + * @param offset In: position to read, Out: shifted by size + * + * @return error code + */ +int io_read_data(void *map, const size_t size, void *p, size_t *offset); + +/** + * @brief write a single value in a file + * + * @param map file mapping + * @param size size of the chunk to write + * @param p pointer to the data + * @param offset In: position to write, Out: shifted by size + * + * @return error code + */ +int io_write_data(void *map, const size_t size, const void *p, size_t *offset); + +/** + * @brief read a maks with its offset + * + * @param h #header file structure + * @param map file mapping + * @param offset In: position in the file, Out: shifted by the mask + offset + * size + * @param mask mask read + * @param diff_offset offset difference to previous/next corresponding chunk + * + * @return error code + */ +int io_read_mask(const struct header *h, void *map, size_t *offset, + size_t *mask, size_t *diff_offset); + +#endif // __SWIFT_LOGGER_IO_H__ diff --git a/logger/logger_loader.c b/logger/logger_loader.c new file mode 100644 index 0000000000000000000000000000000000000000..108f94a89c660147551744d22cffe01c1cf59e2d --- /dev/null +++ b/logger/logger_loader.c @@ -0,0 +1,367 @@ +#include "header.h" +#include "io.h" +#include "particle.h" +#include "timeline.h" + +#include <Python.h> +#include <errno.h> +#include <numpy/arrayobject.h> +#include <stdio.h> +#include <stdlib.h> + +/** + * @brief Reverse offset in dump file + * + * @param filename string filename of the dump file + */ +int reverseOffset(char *filename) { + struct header h; + + /* open file */ + int fd; + void *map; + if (io_open_file(filename, &fd, &map) != 0) return 1; + + /* read header */ + if (header_read(&h, map) != 0) return 1; + + header_print(&h); + + /* check offset direction */ + if (h.forward_offset) { + error_no_return(EIO, "Offset are already reversed"); + return 1; + } + + /* compute file size */ + size_t sz; + int status = io_get_file_size(fd, &sz); + if (status != 0) return 1; + + size_t offset; + int error_code; + +#ifdef SWIFT_DEBUG_CHECKS + /* check offset */ + printf("Check offsets...\n"); + offset = h.offset_first; + while (offset < sz) { + error_code = tools_check_offset(&h, map, &offset); + if (error_code != 0) return 1; + } + printf("Check done\n"); +#endif + + /* reverse header offset */ + header_change_offset_direction(&h, map); + + offset = h.offset_first; + + /* reverse chunks */ + printf("Reversing offsets...\n"); + while (offset < sz) { + error_code = tools_reverse_offset(&h, map, &offset); + if (error_code != 0) return 1; + } + printf("Reversing done\n"); + +#ifdef SWIFT_DEBUG_CHECKS + /* check offset */ + printf("Check offsets...\n"); + offset = h.offset_first; + while (offset < sz) { + error_code = tools_check_offset(&h, map, &offset); + + if (error_code != 0) return 1; + } + printf("Check done\n"); +#endif + + /* free internal variables */ + header_free(&h); + + if (io_close_file(&fd, &map) != 0) return 1; + + return 0; +} + +/** + * @brief load data from the offset without any interpolation + * + * @param offset PyArrayObject list of offset for each particle + * @param filename string filename of the dump file + * @return dictionnary containing the data read + */ +static PyObject *loadFromIndex(__attribute__((unused)) PyObject *self, + PyObject *args) { + + struct header h; + + /* input */ + PyArrayObject *offset = NULL; + char *filename = NULL; + + /* output */ + PyArrayObject *pos = NULL; + PyArrayObject *vel = NULL; + PyArrayObject *acc = NULL; + PyArrayObject *entropy = NULL; + PyArrayObject *h_sph = NULL; + PyArrayObject *rho = NULL; + PyArrayObject *mass = NULL; + PyArrayObject *id = NULL; + + size_t time_offset; + + /* parse arguments */ + + if (!PyArg_ParseTuple(args, "OsL", &offset, &filename, &time_offset)) + return NULL; + + if (!PyArray_Check(offset)) { + error_no_return(ENOTSUP, "Offset is not a numpy array"); + return NULL; + } + if (PyArray_NDIM(offset) != 1) { + error_no_return(ENOTSUP, "Offset is not a 1 dimensional array"); + return NULL; + } + if (PyArray_TYPE(offset) != NPY_UINT64) { + error_no_return(ENOTSUP, "Offset does not contain unsigned int"); + return NULL; + } + + /* open file */ + int fd; + void *map; + if (io_open_file(filename, &fd, &map) != 0) return NULL; + + /* read header */ + if (header_read(&h, map) != 0) return NULL; + + /* reverse offset if needed */ + if (!h.forward_offset) { + if (io_close_file(&fd, &map) != 0) return NULL; + + if (reverseOffset(filename) != 0) return NULL; + + if (io_open_file(filename, &fd, &map) != 0) return NULL; + + /* Reset header */ + header_free(&h); + + if (header_read(&h, map) != 0) return NULL; + } + + /* read timestamps */ + struct time_array times; + + if (time_array_init(×, &h, map, fd) != 0) return NULL; + + time_array_print(×); + /* get required time */ + double time = time_array_get_time(×, time_offset); + + /* init array */ + npy_intp dim[2]; + dim[0] = PyArray_DIMS(offset)[0]; + dim[1] = DIM; + + /* init output */ + if (header_is_present(&h, "position")) { + pos = (PyArrayObject *)PyArray_SimpleNew(2, dim, NPY_DOUBLE); + } + + if (header_is_present(&h, "velocity")) { + vel = (PyArrayObject *)PyArray_SimpleNew(2, dim, NPY_FLOAT); + } + + if (header_is_present(&h, "acceleration")) { + acc = (PyArrayObject *)PyArray_SimpleNew(2, dim, NPY_FLOAT); + } + + if (header_is_present(&h, "entropy")) { + entropy = + (PyArrayObject *)PyArray_SimpleNew(1, PyArray_DIMS(offset), NPY_FLOAT); + } + + if (header_is_present(&h, "cutoff radius")) { + h_sph = + (PyArrayObject *)PyArray_SimpleNew(1, PyArray_DIMS(offset), NPY_FLOAT); + } + + if (header_is_present(&h, "density")) { + rho = + (PyArrayObject *)PyArray_SimpleNew(1, PyArray_DIMS(offset), NPY_FLOAT); + } + + if (header_is_present(&h, "consts")) { + mass = + (PyArrayObject *)PyArray_SimpleNew(1, PyArray_DIMS(offset), NPY_FLOAT); + id = (PyArrayObject *)PyArray_SimpleNew(1, PyArray_DIMS(offset), NPY_ULONG); + } + + int error_code = 0; + + /* check data type in particles */ + if (!particle_check_data_type(&h)) { + error_no_return(ENOTSUP, "Particle data type are not compatible"); + return NULL; + } + + /* loop over all particles */ + for (npy_intp i = 0; i < PyArray_DIMS(offset)[0]; i++) { + struct particle part; + + size_t *offset_particle = (size_t *)PyArray_GETPTR1(offset, i); + + error_code = particle_read(&part, &h, map, offset_particle, time, + reader_lin, ×); + if (error_code != 0) return NULL; + + double *dtmp; + float *ftmp; + size_t *stmp; + + /* copy data */ + for (size_t k = 0; k < DIM; k++) { + if (pos) { + dtmp = PyArray_GETPTR2(pos, i, k); + *dtmp = part.pos[k]; + } + + if (vel) { + ftmp = PyArray_GETPTR2(vel, i, k); + *ftmp = part.vel[k]; + } + + if (acc) { + ftmp = PyArray_GETPTR2(acc, i, k); + *ftmp = part.acc[k]; + } + } + + if (entropy) { + ftmp = PyArray_GETPTR1(entropy, i); + *ftmp = part.entropy; + } + + if (rho) { + ftmp = PyArray_GETPTR1(rho, i); + *ftmp = part.density; + } + + if (h_sph) { + ftmp = PyArray_GETPTR1(h_sph, i); + *ftmp = part.h; + } + + if (mass) { + ftmp = PyArray_GETPTR1(mass, i); + *ftmp = part.mass; + } + + if (id) { + stmp = PyArray_GETPTR1(id, i); + *stmp = part.id; + } + } + + header_free(&h); + + /* construct return */ + PyObject *dict = PyDict_New(); + PyObject *key = PyUnicode_FromString("position"); + PyDict_SetItem(dict, key, PyArray_Return(pos)); + + if (vel) { + key = PyUnicode_FromString("velocity"); + PyDict_SetItem(dict, key, PyArray_Return(vel)); + } + + if (acc) { + key = PyUnicode_FromString("acceleration"); + PyDict_SetItem(dict, key, PyArray_Return(acc)); + } + + if (entropy) { + key = PyUnicode_FromString("entropy"); + PyDict_SetItem(dict, key, PyArray_Return(entropy)); + } + + if (rho) { + key = PyUnicode_FromString("rho"); + PyDict_SetItem(dict, key, PyArray_Return(rho)); + } + + if (h_sph) { + key = PyUnicode_FromString("h_sph"); + PyDict_SetItem(dict, key, PyArray_Return(h_sph)); + } + + if (mass) { + key = PyUnicode_FromString("mass"); + PyDict_SetItem(dict, key, PyArray_Return(mass)); + } + + if (id) { + key = PyUnicode_FromString("id"); + PyDict_SetItem(dict, key, PyArray_Return(id)); + } + + io_close_file(&fd, &map); + + return dict; +} + +/** + * @brief Reverse offset in dump file + * + * @param filename string filename of the dump file + */ +static PyObject *pyReverseOffset(__attribute__((unused)) PyObject *self, + PyObject *args) { + /* input */ + char *filename = NULL; + + /* parse arguments */ + + if (!PyArg_ParseTuple(args, "s", &filename)) return NULL; + + if (reverseOffset(filename) != 0) return NULL; + + return Py_BuildValue(""); +} + +/* definition of the method table */ + +static PyMethodDef logger_loaderMethods[] = { + {"loadFromIndex", loadFromIndex, METH_VARARGS, + "Load snapshot directly from the offset in an index file."}, + {"reverseOffset", pyReverseOffset, METH_VARARGS, + "Reverse the offset (from pointing backward to forward)."}, + + {NULL, NULL, 0, NULL} /* Sentinel */ +}; + +static struct PyModuleDef logger_loadermodule = { + PyModuleDef_HEAD_INIT, + "logger_loader", + "Module reading a SWIFTsim logger snapshot", + -1, + logger_loaderMethods, + NULL, /* m_slots */ + NULL, /* m_traverse */ + NULL, /* m_clear */ + NULL /* m_free */ +}; + +PyMODINIT_FUNC PyInit_logger_loader(void) { + PyObject *m; + m = PyModule_Create(&logger_loadermodule); + if (m == NULL) return NULL; + + import_array(); + + return m; +} diff --git a/logger/particle.c b/logger/particle.c new file mode 100644 index 0000000000000000000000000000000000000000..26b20a8fd5bdb7a38f744d32a7730e88996b5cd8 --- /dev/null +++ b/logger/particle.c @@ -0,0 +1,175 @@ +#include "particle.h" +#include "header.h" +#include "io.h" +#include "timeline.h" +#include "tools.h" + +#include <stdio.h> +#include <stdlib.h> + +void particle_print(const struct particle *p) { + printf("ID: %lu\n", p->id); + printf("Mass: %g\n", p->mass); + printf("Time: %g\n", p->time); + printf("Cutoff Radius: %g\n", p->h); + printf("Position: (%g, %g, %g)\n", p->pos[0], p->pos[1], p->pos[2]); + printf("Velocity: (%g, %g, %g)\n", p->vel[0], p->vel[1], p->vel[2]); + printf("Acceleration: (%g, %g, %g)\n", p->acc[0], p->acc[1], p->acc[2]); + printf("Entropy: %g\n", p->entropy); + printf("Density: %g\n", p->density); +} + +int particle_check_data_type(__attribute__((unused)) const struct header *h) { + printf("TODO check_data_type\n"); + return 1; +} + +void particle_init(struct particle *part) { + for (size_t k = 0; k < DIM; k++) { + part->pos[k] = 0; + part->vel[k] = 0; + part->acc[k] = 0; + } + + part->entropy = -1; + part->density = -1; + part->h = -1; + part->mass = -1; + part->id = SIZE_MAX; +} + +int particle_read_field(struct particle *part, void *map, size_t *offset, + const char *field, const size_t size) { + void *p = NULL; + + if (strcmp("position", field) == 0) { + p = &part->pos; + } else if (strcmp("velocity", field) == 0) { + p = &part->vel; + } else if (strcmp("acceleration", field) == 0) { + p = &part->acc; + } else if (strcmp("entropy", field) == 0) { + p = &part->entropy; + } else if (strcmp("cutoff radius", field) == 0) { + p = &part->h; + } else if (strcmp("density", field) == 0) { + p = &part->density; + } else if (strcmp("consts", field) == 0) { + p = malloc(size); + } else { + error(ENOTSUP, "Type %s not defined", field); + } + + int error_code = io_read_data(map, size, p, offset); + + if (error_code != 0) return error_code; + + if (strcmp("consts", field) == 0) { + part->mass = 0; + part->id = 0; + memcpy(&part->mass, p, sizeof(float)); + p += sizeof(float); + memcpy(&part->id, p, sizeof(size_t)); + p -= sizeof(float); + free(p); + } + return error_code; +} + +int particle_read(struct particle *part, const struct header *h, void *map, + size_t *offset, const double time, const int reader, + struct time_array *times) { + int error_code = 0; + size_t mask = 0; + size_t h_offset = 0; + + particle_init(part); + + error_code = io_read_mask(h, map, offset, &mask, &h_offset); + if (error_code != 0) return error_code; + + if (mask != 127) error(EIO, "Unexpected mask: %lu", mask); + + for (size_t i = 0; i < h->nber_mask; i++) { + if (mask & h->masks[i]) { + error_code = particle_read_field(part, map, offset, h->masks_name[i], + h->masks_size[i]); + if (error_code != 0) return error_code; + } + } + + if (times) /* move offset by 1 in order to be in the required chunk */ + part->time = time_array_get_time(times, *offset - 1); + else + part->time = -1; + + /* end of const case */ + if (reader == reader_const) return 0; + + /* read next particle */ + struct particle part_next; + + if (!h->forward_offset) error(ENOSYS, "TODO"); + + if (h_offset == 0) return 0; + /* get absolute offset of next particle */ + h_offset += *offset - header_get_mask_size(h, mask) - LOGGER_MASK_SIZE - LOGGER_OFFSET_SIZE; + + part_next.time = time_array_get_time(times, h_offset); + + /* previous part exists */ + error_code = particle_read(&part_next, h, map, &h_offset, part_next.time, + reader_const, times); + if (error_code != 0) return error_code; + + error_code = particle_interpolate(part, &part_next, time); + if (error_code != 0) return error_code; + + return 0; +} + +int particle_interpolate(struct particle *part_curr, + const struct particle *part_next, + const double time) { + + if (!part_curr) error(EFAULT, "part_curr is NULL"); + if (!part_next) error(EFAULT, "part_next is NULL"); + +#ifdef SWIFT_DEBUG_CHECKS + if (part_next->time <= part_curr->time) + error(EIO, "Wrong particle order (next before current)"); + if ((time < part_curr->time) || (part_next->time < time)) + error(EIO, + "Interpolating, not extrapolating (particle time: %f, " + "interpolating time: %f, next particle time: %f)", + part_curr->time, time, part_next->time); +#endif + + double scaling = part_next->time - part_curr->time; + + scaling = (time - part_curr->time) / scaling; + + double tmp; + float ftmp; + + /* interpolate vectors */ + for (size_t i = 0; i < DIM; i++) { + tmp = (part_next->pos[i] - part_curr->pos[i]); + part_curr->pos[i] += tmp * scaling; + + ftmp = (part_next->vel[i] - part_curr->vel[i]); + part_curr->vel[i] += ftmp * scaling; + + ftmp = (part_next->acc[i] - part_curr->acc[i]); + part_curr->acc[i] += ftmp * scaling; + } + + /* interpolate scalars */ + ftmp = (part_next->entropy - part_curr->entropy); + part_curr->entropy += ftmp * scaling; + + /* set time */ + part_curr->time = time; + + return 0; +} diff --git a/logger/particle.h b/logger/particle.h new file mode 100644 index 0000000000000000000000000000000000000000..07712cd369f6e284ec9d39338ac2773194ff665f --- /dev/null +++ b/logger/particle.h @@ -0,0 +1,116 @@ +#ifndef __PARTICLE_H__ +#define __PARTICLE_H__ + +#include "header.h" +#include "timeline.h" +#include "tools.h" + +#include <stdio.h> +#include <stdlib.h> + +#define DIM 3 + +struct particle { + /* position */ + double pos[DIM]; + + /* velocity */ + float vel[DIM]; + + /* acceleration */ + float acc[DIM]; + + /* entropy */ + float entropy; + + /* cutoff radius */ + float h; + + /* density */ + float density; + + /* mass */ + float mass; + + /* id */ + size_t id; + + /* time */ + double time; +}; + +enum reader_type { + reader_const, + reader_lin, +}; + +/** + * @brief print a particle + * + * @param p #particle particle to print + */ +void particle_print(const struct particle *p); + +/** + * @brief read a particle in the dump file + * + * @param part #particle particle to update + * @param h #header structure of the file + * @param map file mapping + * @param offset offset of the chunk to read (update it to the end of the chunk) + * @param time time to interpolate (if #reader_type is an interpolating one) + * @param reader #reader_type + * @param times #time_array times in the dump + * + * @return error code + */ +int particle_read(struct particle *part, const struct header *h, void *map, + size_t *offset, const double time, const int reader, + struct time_array *times); + +/** + * @brief Check if dump data type are compatible with the particle type + * + * @param h #header structure of the file + * + * @return error code + */ +int particle_check_data_type(const struct header *h); + +/** + * @brief initialize a particle + * + * @param part #particle particle to initialize + */ +void particle_init(struct particle *part); + +/** + * @brief read a single field for a particle + * + * @param part @particle particle to update + * @param map file mapping + * @param offset In: read position, Out: input shifted by the required amount of + * data + * @param field field to read + * @param size number of bits to read + * + * @return error code + */ +int particle_read_field(struct particle *part, void *map, size_t *offset, + const char *field, const size_t size); + +/** + * @brief interpolate two particles at a given time + * + * @param part_curr #particle In: current particle (before time), Out: + * interpolated particle + * @param part_next #particle next particle (after time) + * @param time interpolation time + * + * @return error code + */ +int particle_interpolate(struct particle *part_curr, + const struct particle *part_next, + const double time); + +#endif //__PARTICLE_H__ diff --git a/logger/timeline.c b/logger/timeline.c new file mode 100644 index 0000000000000000000000000000000000000000..7eec898d54c0152231bc28ccb7b1ac211ec08a15 --- /dev/null +++ b/logger/timeline.c @@ -0,0 +1,197 @@ +#include "timeline.h" +#include "io.h" + +double time_convert_to_double(const integertime_t ti, const double timeBase) { + return ti * timeBase; // should add timebegin +} + +integertime_t time_convert_to_integer(const double_t ti, + const double timeBase) { + return ti / timeBase; // should add timebegin +} + +int time_read(integertime_t *timestamp, double *time, const struct header *h, void *map, + size_t *offset) { + int error_code = 0; + size_t mask = 0; + size_t ind = 0; + size_t prev_offset = 0; + *timestamp = 0; + *time = 0; + + /* read chunck header */ + error_code = io_read_mask(h, map, offset, &mask, &prev_offset); + if (error_code != 0) return error_code; + +#ifdef SWIFT_DEBUG_CHECKS + /* check if timestamp is present */ + if (!header_is_present_and_get_index(h, "timestamp", &ind)) + error(EIO, "Header does not contain a timestamp"); + + /* check if timestamp */ + if (h->masks[ind] != mask) error(EIO, "Not a timestamp"); +#endif + + /* read data */ + // TODO + error_code = io_read_data(map, sizeof(unsigned long long int), timestamp, offset); + error_code = io_read_data(map, sizeof(double), time, offset); + + return error_code; +} + +int time_first_timestamp(const struct header *h, void *map, size_t *offset) { + *offset = h->offset_first; + int test; + + size_t i; + + if (!header_is_present_and_get_index(h, "timestamp", &i)) + error(EIO, "Time stamp not present in header"); + + size_t tmp = *offset; + size_t mask = 0; + test = io_read_mask(h, map, offset, &mask, NULL); + if (test != 0) return test; + + if (mask != h->masks[i]) error(EIO, "Dump should begin by timestep"); + + *offset = tmp; + return 0; +} + +int time_array_init(struct time_array *t, const struct header *h, void *map, + int fd) { + + t->next = NULL; + t->prev = NULL; + + /* get first time stamp */ + size_t offset = 0; + time_first_timestamp(h, map, &offset); + + integertime_t timestamp = 0; + double time = 0; + + /* get file size */ + size_t file_size; + int test = io_get_file_size(fd, &file_size); + + while (offset < file_size) { + + /* read time */ + t->offset = offset; + size_t tmp_offset = offset; + time_read(×tamp, &time, h, map, &tmp_offset); + t->timestamp = timestamp; + t->time = time; + + /* get next chunk */ + test = tools_get_next_chunk(h, map, &offset, fd); + if (test == -1) + break; + else if (test != 0) + return test; + + /* allocate next time_array */ + struct time_array *tmp = malloc(sizeof(struct time_array)); + tmp->prev = t; + tmp->next = NULL; + t->next = tmp; + t = tmp; + } + + /* free unused time_array */ + struct time_array *tmp = t->prev; + tmp->next = NULL; + free(t); + + return 0; +} + +integertime_t time_array_get_integertime(struct time_array *t, const size_t offset) { + const struct time_array *tmp = time_array_get_time_array(t, offset); + return tmp->timestamp; +} + +double time_array_get_time(struct time_array *t, const size_t offset) { + const struct time_array *tmp = time_array_get_time_array(t, offset); + return tmp->time; +} + +struct time_array *time_array_get_time_array(struct time_array *t, + const size_t offset) { + +#ifdef SWIFT_DEBUG_CHECKS + if (!t) return 0; +#endif + while (t->next) { + if (t->offset > offset) break; + + t = t->next; + } + + if (t->prev == NULL) return t; + + struct time_array *tmp = t->prev; + return tmp; +} + +void time_array_free(struct time_array *t) { + struct time_array *tmp; + while (t->next) { + tmp = t->next; + free(t); + t = tmp; + } +} + +void time_array_print(const struct time_array *t) { + size_t threshold = 4; + + size_t n = time_array_count(t); + size_t i = 0; + size_t up_threshold = n - threshold; + + printf("Times (size %lu): [%lli (%g)", n, t->timestamp, t->time); + + while (t->next) { + i += 1; + t = t->next; + if (i < threshold || i > up_threshold) printf(", %lli (%g)", t->timestamp, t->time); + + if (i == threshold) printf(", ..."); + } + + printf("]\n"); +} + +void time_array_print_offset(const struct time_array *t) { + size_t threshold = 4; + + size_t n = time_array_count(t); + size_t i = 0; + size_t up_threshold = n - threshold; + + printf("Times (size %lu): [%lu", n, t->offset); + + while (t->next) { + i += 1; + t = t->next; + if (i < threshold || i > up_threshold) printf(", %lu", t->offset); + + if (i == threshold) printf(", ..."); + } + + printf("]\n"); +} + +size_t time_array_count(const struct time_array *t) { + size_t count = 1; + while (t->next) { + t = t->next; + count += 1; + } + + return count; +} diff --git a/logger/timeline.h b/logger/timeline.h new file mode 100644 index 0000000000000000000000000000000000000000..c5231dcaf78fe4816adbf1bb87d7dd799013142e --- /dev/null +++ b/logger/timeline.h @@ -0,0 +1,133 @@ +#ifndef __TIMELINE_H__ +#define __TIMELINE_H__ + +#include "header.h" +#include "tools.h" + +typedef char timebin_t; +typedef long long integertime_t; + +struct time_array { + void *next; + void *prev; + integertime_t timestamp; + double time; + size_t offset; +}; + +/** + * @brief convert an integer time to a real time + * + * @param ti integer time to convert + * @param timeBase time base of the integer time + * + * @return converted time + */ +double time_convert_to_double(const integertime_t ti, const double timeBase); + +/** + * @brief convert an integer time to a real time + * + * @param ti integer time to convert + * @param timeBase time base of the integer time + * + * @return converted time + */ +integertime_t time_convert_to_integer(const double ti, const double timeBase); + +/** + * @brief read a time stamp + * + * @param timestamp timestamp read + * @param time time read + * @param h #header file structure + * @param map file mapping + * @param offset In: position in the file, Out: shifted by the timestamp + */ +int time_read(integertime_t *timestamp, double *time, const struct header *h, void *map, + size_t *offset); + +/** + * @brief Initialize a time array + * + * @param t #time_array to initialize + * @param h #header file structure + * @param map file mapping + * @param fd file id + */ +int time_array_init(struct time_array *t, const struct header *h, void *map, + int fd); + +/** + * @brief access the time of a given chunk (by its offset) + * + * @param t #time_array to access + * @param offset offset of the chunk + * + * @return integer time of the chunk + */ +integertime_t time_array_get_integertime(struct time_array *t, const size_t offset); + +/** + * @brief access the time of a given chunk (by its offset) + * + * @param t #time_array to access + * @param offset offset of the chunk + * + * @return time of the chunk + */ +double time_array_get_time(struct time_array *t, const size_t offset); + +/** + * @brief access the #time_array of a given chunk (by its offset) + * + * @param t #time_array to access + * @param offset offset of the chunk + * + * @return pointer to the requested #time_array + */ +struct time_array *time_array_get_time_array(struct time_array *t, + const size_t offset); + +/** + * @brief free memory of a #time_array + * + * @param t #time_array to free + */ +void time_array_free(struct time_array *t); + +/** + * @brief print a #time_array + * + * @param t #time_array to print + */ +void time_array_print(const struct time_array *t); + +/** + * @brief print a #time_array (offset) + * + * @param t #time_array to print + */ +void time_array_print_offset(const struct time_array *t); + +/** + * @brief count number of element in #time_array + * + * @param t #time_array to count + * + * @return number of element + */ +size_t time_array_count(const struct time_array *t); + +/** + * @brief get offset of first timestamp + * + * @param h file #header + * @param map file mapping + * @param offset out: offset of first timestamp + * + * @return error code + */ +int time_first_timestamp(const struct header *h, void *map, size_t *offset); + +#endif // __TIMELINE_H__ diff --git a/logger/tools.c b/logger/tools.c new file mode 100644 index 0000000000000000000000000000000000000000..ce2fde8f4c5f12fc0491a2a867e88c2420d8c80b --- /dev/null +++ b/logger/tools.c @@ -0,0 +1,161 @@ +#include "tools.h" +#include "header.h" +#include "io.h" + +#include "particle.h" + +#include <stdio.h> + +int tools_get_next_chunk(const struct header *h, void *map, size_t *offset, + int fd) { + if (h->forward_offset) + return _tools_get_next_chunk_forward(h, map, offset); + else + return _tools_get_next_chunk_backward(h, map, offset, fd); +} + +int _tools_get_next_chunk_forward(const struct header *h, void *map, + size_t *offset) { + size_t diff_offset = 0; + + /* read offset */ + int test = io_read_mask(h, map, offset, NULL, &diff_offset); + if (test != 0) return test; + + if (diff_offset == 0) return -1; + + /* set absolute offset */ + *offset += diff_offset - LOGGER_MASK_SIZE - LOGGER_OFFSET_SIZE; + return 0; +} + +int _tools_get_next_chunk_backward(const struct header *h, void *map, + size_t *offset, int fd) { +#ifndef SWIFT_DEBUG_CHECKS + error(ENOTSUP, "Should not be used, method too slow"); +#endif + size_t current_offset = *offset; + size_t chunk_header = LOGGER_MASK_SIZE + LOGGER_OFFSET_SIZE; + + size_t file_size = 0; + int test = io_get_file_size(fd, &file_size); + + if (test != 0) return test; + + while (current_offset < file_size) { + size_t mask = 0; + size_t prev_offset; + int test = io_read_mask(h, map, ¤t_offset, &mask, &prev_offset); + if (test != 0) return test; + + prev_offset = current_offset - prev_offset - chunk_header; + if (*offset == prev_offset) { + *offset = current_offset - chunk_header; + return 0; + } + + current_offset += header_get_mask_size(h, mask); + } + + return -1; +} + +int tools_reverse_offset(const struct header *h, void *map, size_t *offset) { + size_t mask = 0; + size_t prev_offset = 0; + const size_t cur_offset = *offset; + + /* read mask + offset */ + int error_code = io_read_mask(h, map, offset, &mask, &prev_offset); + if (error_code != 0) return error_code; + + /* write offset of zero (in case it is the last chunk) */ + const size_t zero = 0; + *offset -= LOGGER_OFFSET_SIZE; + error_code = io_write_data(map, LOGGER_OFFSET_SIZE, &zero, offset); + if (error_code != 0) error(errno, "Unable to write data while reversing"); + + /* set offset after current chunk */ + *offset += header_get_mask_size(h, mask); + + /* first chunks do not have a previous partner */ + if (prev_offset == cur_offset) return 0; + + if (prev_offset > cur_offset) + error(EIO, "Unexpected offset, header %lu, current %lu", prev_offset, + cur_offset); + + /* modify previous offset */ + size_t abs_prev_offset = cur_offset - prev_offset + LOGGER_MASK_SIZE; + error_code = io_write_data(map, LOGGER_OFFSET_SIZE, &prev_offset, &abs_prev_offset); + if (error_code != 0) error(errno, "Unable to write offset"); + +#ifdef SWIFT_DEBUG_CHECKS + size_t prev_mask = 0; + abs_prev_offset -= LOGGER_MASK_SIZE + LOGGER_OFFSET_SIZE; + error_code = io_read_mask(h, map, &abs_prev_offset, &prev_mask, NULL); + + if (prev_mask != mask) + error(EIO, "Unexpected mask: %lu (at %lu), got %lu (at %lu)", mask, *offset, + prev_mask, prev_offset); + +#endif // SWIFT_DEBUG_CHECKS + + return 0; +} + +int tools_check_offset(const struct header *h, void *map, size_t *offset) { +#ifndef SWIFT_DEBUG_CHECKS + error(ENOTSUP, "Should not check in non debug mode"); +#endif + + size_t tmp = *offset; + + size_t mask; + size_t pointed_offset; + + /* read mask + offset */ + int error_code = io_read_mask(h, map, offset, &mask, &pointed_offset); + if (error_code != 0) return error_code; + + /* get absolute offset */ + if (h->forward_offset) + pointed_offset += tmp; + else { + if (tmp < pointed_offset) + error(EIO, "Offset too large (%lu) at %lu with mask %lu", pointed_offset, + tmp, mask); + pointed_offset = tmp - pointed_offset; + } + + /* set offset after current chunk */ + *offset += header_get_mask_size(h, mask); + + if (pointed_offset == tmp || pointed_offset == 0) return 0; + + /* read mask of the pointed chunk */ + size_t pointed_mask = 0; + error_code = io_read_mask(h, map, &pointed_offset, &pointed_mask, NULL); + + /* check masks */ + if (pointed_mask != mask) + error(EIO, "Error in the offset (mask %lu != %lu) at %lu and %lu", mask, + pointed_mask, + *offset - header_get_mask_size(h, mask) - LOGGER_MASK_SIZE - LOGGER_OFFSET_SIZE, + pointed_offset - LOGGER_MASK_SIZE - LOGGER_OFFSET_SIZE); + + if (pointed_mask == 128) return 0; + + struct particle part; + error_code = particle_read(&part, h, map, &tmp, 0, reader_const, NULL); + + size_t id = part.id; + tmp = pointed_offset - LOGGER_MASK_SIZE - LOGGER_OFFSET_SIZE; + error_code = particle_read(&part, h, map, &tmp, 0, reader_const, NULL); + + if (id != part.id) + error(EIO, "Offset wrong, id incorrect (%lu != %lu) at %lu", id, part.id, + tmp); + + return 0; +} diff --git a/logger/tools.h b/logger/tools.h new file mode 100644 index 0000000000000000000000000000000000000000..f3b666485fa000be0724d51e60127aefcbeaa973 --- /dev/null +++ b/logger/tools.h @@ -0,0 +1,103 @@ +#ifndef __TOOLS_H__ +#define __TOOLS_H__ + +#define NPY_NO_DEPRECATED_API NPY_1_7_API_VERSION + +#include <Python.h> +#include <stdio.h> +#include <stdlib.h> + +#define STRING_SIZE 200 + +struct header; + +/* Set the error message for python and return the error code + * WARNING for python, you need to return NULL and not the error code + */ +#define error(err, s, ...) \ + ({ \ + error_no_return(err, s, ##__VA_ARGS__); \ + return err; \ + }) + +/* Same as error but does not return the error code. + * Should be used for python functions */ +#define error_no_return(err, s, ...) \ + ({ \ + char error_msg[STRING_SIZE]; \ + sprintf(error_msg, "%s:%s():%i: " s ": %s\n", __FILE__, __FUNCTION__, \ + __LINE__, ##__VA_ARGS__, strerror(err)); \ + PyErr_SetString(PyExc_RuntimeError, error_msg); \ + }) + +#define message(s, ...) \ + ({ \ + printf("%s:%s():%i: " s "\n", __FILE__, __FUNCTION__, __LINE__, \ + ##__VA_ARGS__); \ + }) + +/** + * @brief get the offset of the next corresponding chunk + * + * @param h #header structure of the file + * @param map file mapping + * @param offset In: initial offset, Out: offset of the next chunk + * @param fd file id + * + * @return error code, -1 if no next chunk + */ +int tools_get_next_chunk(const struct header *h, void *map, size_t *offset, + int fd); + +/** + * @brief internal function of #tools_get_next_chunk. Should not be used (very + * slow) + * + * @param h #header structure of the file + * @param map file mapping + * @param offset In: initial offset, Out: offset of the next chunk + * @param fd file id + * + * @return error code, -1 if no next chunk + */ +int _tools_get_next_chunk_backward(const struct header *h, void *map, + size_t *offset, int fd); + +/** + * @brief internal function of #tools_get_next_chunk. Should not be used outside + * + * @param h #header structure of the file + * @param map file mapping + * @param offset In: initial offset, Out: offset of the next chunk + * + * @return error code, -1 if no next chunk + */ +int _tools_get_next_chunk_forward(const struct header *h, void *map, + size_t *offset); + +/** + * @brief switch side offset + * + * From current chunk, switch side of the offset of the previous one. + * @param h #header structure of the file + * @param map file mapping + * @param offset In: initial offset, Out: offset of next chunk + * + * @return error code + */ +int tools_reverse_offset(const struct header *h, void *map, size_t *offset); + +/** + * @brief debugging function checking the offset of a chunk + * + * Compare the mask with the one pointed by the header. + * if the chunk is a particle, check the id too. + * @param h #header structure of the file + * @param map file mapping + * @param offset In: chunk offset, Out: offset after the chunk + * + * @return error code + */ +int tools_check_offset(const struct header *h, void *map, size_t *offset); + +#endif //__TOOLS_H__