star-ck

Describe the radiative properties of gas mixtures
git clone git://git.meso-star.fr/star-ck.git
Log | Files | Refs | README | LICENSE

commit d5e520d214c340bb9173a01b2c732d694e9a4f7e
parent 6987976c813f315dc45dd3339d27d0f1b5b27de7
Author: Vincent Forest <vincent.forest@meso-star.com>
Date:   Thu, 24 Mar 2022 09:27:21 +0100

Upd the fileformat and rename library in Star-CK

Diffstat:
MREADME.md | 23++++++++++++++---------
Mcmake/CMakeLists.txt | 82+++++++++++++++++++++++++++++++++++++++----------------------------------------
Ddoc/atrck.5.txt | 112-------------------------------------------------------------------------------
Adoc/sck.5.scd | 110+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
Dsrc/atrck.c | 544-------------------------------------------------------------------------------
Dsrc/atrck.h | 126-------------------------------------------------------------------------------
Dsrc/atrck_c.h | 139-------------------------------------------------------------------------------
Dsrc/atrck_log.c | 125-------------------------------------------------------------------------------
Dsrc/atrck_log.h | 68--------------------------------------------------------------------
Asrc/sck.c | 581+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
Asrc/sck.h | 136+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
Asrc/sck_c.h | 136+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
Asrc/sck_log.c | 125+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
Asrc/sck_log.h | 68++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
Dsrc/test_atrck.c | 70----------------------------------------------------------------------
Dsrc/test_atrck_load.c | 405-------------------------------------------------------------------------------
Asrc/test_sck.c | 71+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
Asrc/test_sck_load.c | 411+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
18 files changed, 1692 insertions(+), 1640 deletions(-)

diff --git a/README.md b/README.md @@ -1,15 +1,16 @@ -# AsToRia: Correlated K +# Star-CorrelatedK This C library loads the radiative properties of a gas mixture saved wrt the -AtrCK fileformat. +Star-CK fileformat. ## How to build -This library is compatible GNU/Linux 64-bits. It relies the +This library is compatible with 64-bits POSIX systems. It relies the [CMake](http://www.cmake.org) and the -[RCMake](https://gitlab.com/vaplv/rcmake/) packages to build. -It also depends on the -[RSys](https://gitlab.com/vaplv/rsys/) library. +[RCMake](https://gitlab.com/vaplv/rcmake/) packages to build. It also depends +on the [RSys](https://gitlab.com/vaplv/rsys/) library. It optionally depends on +[scdoc](https://sr.ht/~sircmpwn/scdoc/) which, if available, is used to +generate the man page of the Star-CK file format. First ensure that CMake is installed on your system. Then install the RCMake package as well as the aforementioned prerequisites. Finally generate the @@ -19,9 +20,13 @@ resulting project can be edited, built, tested and installed as any CMake project. Refer to the [CMake documentation](https://cmake.org/documentation) for further informations on CMake. -## License +## Copyright notice +Copyright (C) 2022 [|Meso|Star>](https://www.meso-star.com) (<contact@meso-star.com>). Copyright (C) 2020, 2021 Centre National de la Recherche Scientifique (CNRS). -AtrCK is free software released under the GPL v3+ license: GNU GPL version 3 or -later. You are welcome to redistribute it under certain conditions; refer to + +## License + +Star-CK is free software released under the GPL v3+ license: GNU GPL version 3 +or later. You are welcome to redistribute it under certain conditions; refer to the COPYING file for details. diff --git a/cmake/CMakeLists.txt b/cmake/CMakeLists.txt @@ -13,11 +13,11 @@ # You should have received a copy of the GNU General Public License # along with this program. If not, see <http://www.gnu.org/licenses/>. -cmake_minimum_required(VERSION 2.8) -project(atrck C) +cmake_minimum_required(VERSION 3.1) +project(sck C) enable_testing() -set(ATRCK_SOURCE_DIR ${PROJECT_SOURCE_DIR}/../src) +set(SCK_SOURCE_DIR ${PROJECT_SOURCE_DIR}/../src) option(NO_TEST "Do not build tests" OFF) ################################################################################ @@ -40,81 +40,79 @@ set(VERSION_MINOR 0) set(VERSION_PATCH 0) set(VERSION ${VERSION_MAJOR}.${VERSION_MINOR}.${VERSION_PATCH}) -set(ATRCK_FILES_SRC - atrck.c - atrck_log.c) -set(ATRCK_FILES_INC - atrck_c.h - atrck_log.h) -set(ATRCK_FILES_INC_API - atrck.h) +set(SCK_FILES_SRC + sck.c + sck_log.c) +set(SCK_FILES_INC + sck_c.h + sck_log.h) +set(SCK_FILES_INC_API + sck.h) -set(ATRCK_FILES_DOC COPYING README.md) +set(SCK_FILES_DOC COPYING README.md) -# Prepend each file in the `ATRCK_FILES_<SRC|INC>' list by `ATRCK_SOURCE_DIR' -rcmake_prepend_path(ATRCK_FILES_SRC ${ATRCK_SOURCE_DIR}) -rcmake_prepend_path(ATRCK_FILES_INC ${ATRCK_SOURCE_DIR}) -rcmake_prepend_path(ATRCK_FILES_INC_API ${ATRCK_SOURCE_DIR}) -rcmake_prepend_path(ATRCK_FILES_DOC ${PROJECT_SOURCE_DIR}/../) +# Prepend each file in the `SCK_FILES_<SRC|INC>' list by `SCK_SOURCE_DIR' +rcmake_prepend_path(SCK_FILES_SRC ${SCK_SOURCE_DIR}) +rcmake_prepend_path(SCK_FILES_INC ${SCK_SOURCE_DIR}) +rcmake_prepend_path(SCK_FILES_INC_API ${SCK_SOURCE_DIR}) +rcmake_prepend_path(SCK_FILES_DOC ${PROJECT_SOURCE_DIR}/../) -add_library(atrck SHARED ${ATRCK_FILES_SRC} ${ATRCK_FILES_INC} ${ATRCK_FILES_INC_API}) -target_link_libraries(atrck RSys) +add_library(sck SHARED ${SCK_FILES_SRC} ${SCK_FILES_INC} ${SCK_FILES_INC_API}) +target_link_libraries(sck RSys) if(CMAKE_COMPILER_IS_GNUCC) - target_link_libraries(atrck m) + target_link_libraries(sck m) endif() -set_target_properties(atrck PROPERTIES - DEFINE_SYMBOL ATRCK_SHARED_BUILD +set_target_properties(sck PROPERTIES + DEFINE_SYMBOL SCK_SHARED_BUILD VERSION ${VERSION} SOVERSION ${VERSION_MAJOR}) -rcmake_setup_devel(atrck AtrCK ${VERSION} astoria/atrck_version.h) +rcmake_setup_devel(sck AtrCK ${VERSION} astoria/sck_version.h) ################################################################################ # Add tests ################################################################################ if(NOT NO_TEST) function(new_test _name) - add_executable(${_name} ${ATRCK_SOURCE_DIR}/${_name}.c) - target_link_libraries(${_name} atrck RSys ${ARGN}) + add_executable(${_name} ${SCK_SOURCE_DIR}/${_name}.c) + target_link_libraries(${_name} sck RSys ${ARGN}) add_test(${_name} ${_name}) endfunction() - new_test(test_atrck) - new_test(test_atrck_load) + new_test(test_sck) + new_test(test_sck_load) endif() ################################################################################ # Man page ############################################################################### -find_program(A2X NAMES a2x a2x.py) -if(NOT A2X) +find_program(SCDOC NAMES scdoc) +if(NOT SCDOC) message(WARNING - "The `a2x' program is missing. " - "The atrck man page cannot be generated.") + "The `scdoc' program is missing. " + "The Star-CK man page cannot be generated.") else() - set(_src ${PROJECT_SOURCE_DIR}/../doc/atrck.5.txt) - set(_txt ${CMAKE_CURRENT_BINARY_DIR}/atrck.5.txt) + set(_src ${PROJECT_SOURCE_DIR}/../doc/sck.5.scd) add_custom_command( - OUTPUT atrck.5 - COMMAND ${CMAKE_COMMAND} -E copy ${_src} ${_txt} - COMMAND ${A2X} -dmanpage -fmanpage ${_txt} + OUTPUT sck.5 + COMMAND ${SCDOC} < ${_src} > sck.5 DEPENDS ${_src} WORKING_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR} - COMMENT "Buid ROFF man page atrck.5" + COMMENT "Buid ROFF man page sck.5" VERBATIM) - add_custom_target(man-roff ALL DEPENDS atrck.5) - install(FILES ${CMAKE_CURRENT_BINARY_DIR}/atrck.5 DESTINATION share/man/man5) + add_custom_target(man-roff ALL DEPENDS sck.5) + install(FILES ${CMAKE_CURRENT_BINARY_DIR}/sck.5 DESTINATION share/man/man5) endif() ################################################################################ # Define output & install directories ################################################################################ -install(TARGETS atrck +install(TARGETS sck ARCHIVE DESTINATION bin LIBRARY DESTINATION lib RUNTIME DESTINATION bin) -install(FILES ${ATRCK_FILES_INC_API} DESTINATION include/astoria) -install(FILES ${ATRCK_FILES_DOC} DESTINATION share/doc/atrck) +install(FILES ${SCK_FILES_INC_API} DESTINATION include/star) +install(FILES ${SCK_FILES_DOC} DESTINATION share/doc/star-ck) diff --git a/doc/atrck.5.txt b/doc/atrck.5.txt @@ -1,112 +0,0 @@ -// Copyright (C) 2020, 2021 CNRS -// -// 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/>. -:toc: - -atrck(5) -======= - -NAME ----- -atrck - AsToRia: Correlated K - -DESCRIPTION ------------ - -*atrck* is a binary file format for storing the correlated K of a gas mixture. -The volumetric mesh to which the CKs are attached is _not_ described there but -must be defined in a separated file via, for example, the *sth*(5) file -format. - -An *atrck* file begins by a header that describes its data layout. Next the -spectral bands on which the CKs were evaluated are listed. Finally, the entire -list of CKs are stored. - -The header is made up of 3 64-bit integers. The first integer is a power of -two (usually 4096) that defines the _<pagesize>_ in bytes to which the per -quadrature point CKs are aligned. The two remaining integers store -respectively the number of bands and the number of nodes, i.e. the number of -CKs per quadrature point which is in fact the number of vertices of the -volumetric mesh to which the K are mapped. - -After the header comes the list of spectral bands sorted in ascending order -with respect to their spectral interval. Each spectral band begins with a -unique identifier stored on a 64-bit signed integer. Then 2 double precision -floating point numbers store the lower and upper limit of the band in -nanometers. Another 64-bit signed integer then gives the number of quadrature -points in the band. Finally the band's quadrature points are listed where each -point is in fact 2 double precision floating numbers: the first gives its -abscissa and the second its weight. - -Padding bytes follow the list of spectral bands to ensure alignment of the CK -to _<pagesize>_. - -Then CKs are listed as a set of _<ka-list>_ each consisting of _<#nodes>_ -single precision floating point numbers per list. Each _<ka-list>_ is followed -by pad bytes to ensure the alignment of the next list to _<pagesize>_. A -_<ka-list>_ is in fact the list of the absorption coefficients of a given -quadrature point in a specific spectral band. These lists are sorted according -to the order of the aforementioned spectral bands and their associated -quadrature points. - -BINARY FILE FORMAT ------------------- - -Data are encoded with respect to the little endian bytes ordering, i.e. least -significant bytes are stored first. - -[verse] -------- -<atrck> ::= <pagesize> <#bands> <#nodes> - <bands-list> - <padding> - <meshes-list> - -<pagesize> ::= INT64 -<#bands> ::= INT64 -<#nodes> ::= INT64 -<padding> ::= [ BYTE ... ] #ensure alignment on <pagesize> - ---- - -<bands-list> ::= <band> [ <band> ... ] -<band> ::= <band-id> <band-low> <band-upp> <#quad-pts> - <quad-pts-list> - -<band-id> ::= INT64 -<band-low> ::= DOUBLE # In nm -<band-upp> ::= DOUBLE # In nĂ¹ -<#quad-pts> ::= INT64 - ---- - -<quad-pts-list> ::= <quad-pt> [ <quad-pt> ... ] -<quad-pt> ::= <quad-abscissa> <quad-weight> - -<quad-abscissa> ::= DOUBLE -<quad-weight> ::= DOUBLE - ---- - -<meshes-list> ::= <ka-list> - <padding> - [ <ka-list> ... ] -<ka-list> ::= <ka> [ <ka> ... ] -<ka> ::= FLOAT - -------- - -SEE ALSO --------- -*mmap*(2), *sth*(5) diff --git a/doc/sck.5.scd b/doc/sck.5.scd @@ -0,0 +1,110 @@ +sck(5) + +; Copyright (C) 2022 |Meso|Star> (contact@meso-star.com) +; Copyright (C) 2020, 2021 CNRS +; +; 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/>. + +# NAME + +sck - Star Correlated K file format + +# DESCRIPTION + +*sck* is a binary file format for storing the correlated K of a gas mixture. +The volumetric mesh to which the CKs are attached is _not_ described there but +must be defined in a separated file via, for example, the *smsh*(5) file +format. + +An *sck* file begins by a header that describes the layout of the data. Next the +spectral bands on which the CKs were evaluated are listed. Finally, the entire +list of CKs are stored. + +A *sck* file begins with a header that describes the layout of the data. Then, +the spectral bands on which the CKs have been evaluated are listed. Finally, +the entire list of radiative coefficients is stored. + +The header consists of 3 64-bit integers. The first integer is a power of two +(usually 4096) that defines the _pagesize_ in bytes on which the radiative +coefficients are aligned. The remaining two integers store the number of bands +and the number of nodes, that is, the number of radiative coefficients per band +or per quadrature point, respectively, which is actually the number of nodes in +the volumetric mesh to which these coefficients are attached. + +After the header comes the list of spectral bands sorted in ascending order +relative to their interval. Each spectral band begins with 2 double-precision +floating-point numbers that represent the lower and upper limits of the band in +nanometers. Another 64-bit integer then gives the number of quadrature points in +the band. Finally, the weights of the quadrature points are stored in a list of +floating numbers with double precision. + +Fill bytes follow the list of spectral bands to ensure alignment of the first +radiative coefficients on _pagesize_. Then, for each band, the diffusion +coefficients are listed before the list of absorption coefficients for each +quadrature point. Each list of radiative coefficients is followed by a +_padding_, which is a list of bytes that provides memory alignment of the +following data to _pagesize_. Bands and quadrature points are sorted according +to the order in which they were previously declared. + +# BINARY FILE FORMAT + +Data are encoded with respect to the little endian bytes ordering, i.e. least +significant bytes are stored first. + +``` +<sck> ::= <pagesize> <#bands> <#nodes> + <bands> + <padding> + <rad-coefs> + +<pagesize> ::= UINT64 +<#bands> ::= UINT64 +<#nodes> ::= UINT64 +<padding> ::= [ BYTE ... ] #ensure alignment on <pagesize> + +--- + +<band-list> ::= <band> [ <band> ... ] +<band> ::= <band-low> <band-upp> <#quad-pts> <quad-pt-list> + +<band-id> ::= UINT64 +<band-low> ::= DOUBLE # In nm +<band-upp> ::= DOUBLE # In nm +<#quad-pts> ::= UINT64 + +--- + +<quad-pt-list> ::= <quad-weight> [ <quad-weight> ... ] +<quad-weight> ::= DOUBLE + +--- + +<rad-coefs> ::= <per-band-k> + [ <per-band-k> ... ] + +<per-band-k> ::= <ks-list> + <per-quad-pt-ka> + +<per-quad-pt-ka> ::= <ka-list> + [ <ka-list>... ] + +<ka-list> ::= <ka> [ <ka> ... ] <padding> +<ks-list> ::= <ks> [ <ks> ... ] <padding> +<ka> ::= FLOAT # in m^-1 +<ks> ::= FLOAT # in m^-1 +``` + +# SEE ALSO + +*mmap*(2), *smsh*(5) diff --git a/src/atrck.c b/src/atrck.c @@ -1,544 +0,0 @@ -/* Copyright (C) 2020, 2021 CNRS - * - * 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/>. */ - -#define _POSIX_C_SOURCE 200809L /* mmap support */ -#define _DEFAULT_SOURCE 1 /* MAP_POPULATE support */ -#define _BSD_SOURCE 1 /* MAP_POPULATE for glibc < 2.19 */ - -#include "atrck.h" -#include "atrck_c.h" -#include "atrck_log.h" - -#include <rsys/algorithm.h> - -#include <unistd.h> /* sysconf support */ - -#include <errno.h> -#include <string.h> -#include <sys/mman.h> /* mmap */ - -/******************************************************************************* - * Helper functions - ******************************************************************************/ -static void -reset_atrck(struct atrck* atrck) -{ - ASSERT(atrck); - atrck->pagesize = 0; - atrck->nnodes = 0; - darray_band_purge(&atrck->bands); -} - -static INLINE int -cmp_dbl(const void* a, const void* b) -{ - const double d0 = *((const double*)a); - const double d1 = *((const double*)b); - return d0 < d1 ? -1 : (d0 > d1 ? 1 : 0); -} - -static res_T -read_quad_pt - (struct atrck* atrck, - struct quad_pt* quad_pt, - FILE* stream, - const char* stream_name, - size_t iband, - size_t iquad_pt) -{ - res_T res = RES_OK; - ASSERT(atrck && quad_pt && stream_name); - - #define READ(Var, Name) { \ - if(fread((Var), sizeof(*(Var)), 1, stream) != 1) { \ - log_err(atrck, \ - "%s: band %lu: quadrature point %lu: could not read the %s.\n", \ - stream_name, iband, iquad_pt, Name); \ - res = RES_IO_ERR; \ - goto error; \ - } \ - } (void)0 - READ(&quad_pt->abscissa, "quadrature point abscissa"); - READ(&quad_pt->weight, "quadrature point weight"); - #undef READ - - quad_pt->ka_list = NULL; - -exit: - return res; -error: - goto exit; -} - -static res_T -read_band - (struct atrck* atrck, - struct band* band, - FILE* stream, - const char* stream_name) -{ - double cumul; - size_t iquad_pt; - uint64_t nquad_pts; - res_T res = RES_OK; - ASSERT(atrck && band && stream_name); - - /* Read band definition */ - #define READ(Var, Name) { \ - if(fread((Var), sizeof(*(Var)), 1, stream) != 1) { \ - log_err(atrck, "%s: band %lu: could not read the %s.\n", \ - stream_name, (unsigned long)band->id, (Name)); \ - res = RES_IO_ERR; \ - goto error; \ - } \ - } (void)0 - READ(&band->id, "band identifier"); - READ(&band->low, "band lower bound"); - READ(&band->upp, "band upper bound"); - READ(&nquad_pts, "#quadrature points"); - #undef READ - - /* Check band description */ - if(band->low < 0 || band->low > band->upp) { - log_err(atrck, - "%s: band %lu: invalid band range [%g, %g].\n", - stream_name, (unsigned long)band->id, band->low, band->upp); - res = RES_BAD_ARG; - goto error; - } - if(nquad_pts == 0) { - log_err(atrck, - "%s: band %lu: invalid number fo quadrature points (#points=%lu).\n", - stream_name, (unsigned long)band->id, (unsigned long)nquad_pts); - res = RES_BAD_ARG; - goto error; - } - - /* Allocate the list of quadrature points */ - res = darray_quad_pt_resize(&band->quad_pts, nquad_pts); - if(res != RES_OK) { - log_err(atrck, - "%s: band %lu: could not allocate the list of quadrature points " - "(#points=%lu).\n", - stream_name, (unsigned long)band->id, (unsigned long)nquad_pts); - goto error; - } - - /* Allocate the cumulative used to sample the quadrature points */ - res = darray_double_resize(&band->quad_pts_cumul, nquad_pts); - if(res != RES_OK) { - log_err(atrck, - "%s: band %lu: could not allocate the cumulative of quadrature points " - "(#points=%lu).\n", - stream_name, (unsigned long)band->id, (unsigned long)nquad_pts); - goto error; - } - - /* Read the data of the quadrature points of the band */ - cumul = 0; - FOR_EACH(iquad_pt, 0, nquad_pts) { - struct quad_pt* quad_pt = darray_quad_pt_data_get(&band->quad_pts)+iquad_pt; - - /* Read current quadrature point */ - res = read_quad_pt(atrck, quad_pt, stream, stream_name, band->id, iquad_pt); - if(res != RES_OK) goto error; - - /* Compute the quadrature point cumulative */ - cumul += quad_pt->weight; - darray_double_data_get(&band->quad_pts_cumul)[iquad_pt] = cumul; - } - - /* The weights are not normalized */ - if(!eq_eps(cumul, 1.0, 1.e-6)) { - log_warn(atrck, - "%s: band %lu: the weights of the quadrature points are not normalised.\n", - stream_name, (unsigned long)band->id); - - /* Renormalize the cumulative */ - FOR_EACH(iquad_pt, 0, nquad_pts) { - darray_double_data_get(&band->quad_pts_cumul)[iquad_pt] /= cumul; - } - } - - /* Handle imprecision issue */ - darray_double_data_get(&band->quad_pts_cumul)[nquad_pts-1] = 1.0; - -exit: - return res; -error: - goto exit; -} - -static res_T -load_stream(struct atrck* atrck, FILE* stream, const char* stream_name) -{ - size_t map_len; - size_t iband; - uint64_t nbands; - off_t offset = 0; - res_T res = RES_OK; - ASSERT(atrck && stream && stream_name); - - reset_atrck(atrck); - - /* Read file header */ - #define READ(Var, Name) { \ - if(fread((Var), sizeof(*(Var)), 1, stream) != 1) { \ - log_err(atrck, "%s: could not read the %s.\n", stream_name, (Name)); \ - res = RES_IO_ERR; \ - goto error; \ - } \ - } (void)0 - READ(&atrck->pagesize, "page size"); - READ(&nbands, "number of bands"); - READ(&atrck->nnodes, "number of nodes"); - #undef READ - - /* Check band description */ - if(!IS_ALIGNED(atrck->pagesize, atrck->pagesize_os)) { - log_err(atrck, - "%s: invalid page size %lu. The page size attribute must be aligned on " - "the page size of the operating system (%lu).\n", - stream_name, - (unsigned long)atrck->pagesize, - (unsigned long)atrck->pagesize_os); - res = RES_BAD_ARG; - goto error; - } - if(!nbands) { - log_err(atrck, "%s: invalid number of bands %lu.\n", - stream_name, (unsigned long)nbands); - res = RES_BAD_ARG; - goto error; - } - if(!atrck->nnodes) { - log_err(atrck, "%s: invalid number of nodes %lu.\n", - stream_name, (unsigned long)atrck->nnodes); - res = RES_BAD_ARG; - goto error; - } - - /* Allocate the bands */ - res = darray_band_resize(&atrck->bands, nbands); - if(res != RES_OK) { - log_err(atrck, "%s: could not allocate the list of bands (#bands=%lu).\n", - stream_name, (unsigned long)nbands); - goto error; - } - - /* Read the band description */ - FOR_EACH(iband, 0, nbands) { - struct band* band = darray_band_data_get(&atrck->bands) + iband; - res = read_band(atrck, band, stream, stream_name); - if(res != RES_OK) goto error; - } - - /* Compute the length in bytes of the ka to map for each quadrature point */ - map_len = ALIGN_SIZE(atrck->nnodes * sizeof(float), atrck->pagesize); - - /* Compute the offset toward the 1st list of ka */ - offset = ftell(stream); - offset = (off_t)ALIGN_SIZE((uint64_t)offset, atrck->pagesize); - - /* Load the per band, per quadrature point and per node correlated K */ - FOR_EACH(iband, 0, nbands) { - struct band* band = darray_band_data_get(&atrck->bands) + iband; - size_t iquad_pt; - FOR_EACH(iquad_pt, 0, darray_quad_pt_size_get(&band->quad_pts)) { - struct quad_pt* quad_pt = NULL; - - quad_pt = darray_quad_pt_data_get(&band->quad_pts)+iquad_pt; - quad_pt->map_len = map_len; - quad_pt->ka_list = mmap(NULL, quad_pt->map_len, PROT_READ, - MAP_PRIVATE|MAP_POPULATE, fileno(stream), offset); - if(quad_pt->ka_list == MAP_FAILED) { - log_err(atrck, - "%s: band %lu: quadrature point %lu: coult not map the correlated K " - "-- %s.\n", - stream_name, - (unsigned long)band->id, - (unsigned long)iquad_pt, - strerror(errno)); - res = RES_IO_ERR; - goto error; - } - offset = (off_t)((size_t)offset + map_len); - } - } - -exit: - return res; -error: - reset_atrck(atrck); - goto exit; -} - -static void -release_atrck(ref_T* ref) -{ - struct atrck* atrck; - ASSERT(ref); - atrck = CONTAINER_OF(ref, struct atrck, ref); - if(atrck->logger == &atrck->logger__) logger_release(&atrck->logger__); - darray_band_release(&atrck->bands); - MEM_RM(atrck->allocator, atrck); -} - -/******************************************************************************* - * Exported functions - ******************************************************************************/ -res_T -atrck_create - (struct logger* logger, /* NULL <=> use default logger */ - struct mem_allocator* mem_allocator, /* NULL <=> use default allocator */ - const int verbose, /* Verbosity level */ - struct atrck** out_atrck) -{ - struct atrck* atrck = NULL; - struct mem_allocator* allocator = NULL; - res_T res = RES_OK; - - if(!out_atrck) { - res = RES_BAD_ARG; - goto error; - } - - allocator = mem_allocator ? mem_allocator : &mem_default_allocator; - atrck = MEM_CALLOC(allocator, 1, sizeof(*atrck)); - if(!atrck) { - if(verbose) { - #define ERR_STR "Could not allocate the AtrKC device.\n" - if(logger) { - logger_print(logger, LOG_ERROR, ERR_STR); - } else { - fprintf(stderr, MSG_ERROR_PREFIX ERR_STR); - } - #undef ERR_STR - } - res = RES_MEM_ERR; - goto error; - } - ref_init(&atrck->ref); - atrck->allocator = allocator; - atrck->verbose = verbose; - atrck->pagesize_os = (size_t)sysconf(_SC_PAGESIZE); - darray_band_init(allocator, &atrck->bands); - if(logger) { - atrck->logger = logger; - } else { - setup_log_default(atrck); - } - -exit: - if(out_atrck) *out_atrck = atrck; - return res; -error: - if(atrck) { - ATRCK(ref_put(atrck)); - atrck = NULL; - } - goto exit; -} - -res_T -atrck_ref_get(struct atrck* atrck) -{ - if(!atrck) return RES_BAD_ARG; - ref_get(&atrck->ref); - return RES_OK; -} - -res_T -atrck_ref_put(struct atrck* atrck) -{ - if(!atrck) return RES_BAD_ARG; - ref_put(&atrck->ref, release_atrck); - return RES_OK; -} - -res_T -atrck_load(struct atrck* atrck, const char* path) -{ - FILE* file = NULL; - res_T res = RES_OK; - - if(!atrck || !path) { - res = RES_BAD_ARG; - goto error; - } - - file = fopen(path, "r"); - if(!file) { - log_err(atrck, "%s: error opening file `%s'.\n", FUNC_NAME, path); - res = RES_IO_ERR; - goto error; - } - - res = load_stream(atrck, file, path); - if(res != RES_OK) goto error; - -exit: - if(file) fclose(file); - return res; -error: - goto exit; -} - -res_T -atrck_load_stream - (struct atrck* atrck, - FILE* stream, - const char* stream_name) -{ - if(!atrck || !stream) return RES_BAD_ARG; - return load_stream(atrck, stream, stream_name ? stream_name : "<stream>"); -} - -size_t -atrck_get_bands_count(const struct atrck* atrck) -{ - ASSERT(atrck); - return darray_band_size_get(&atrck->bands); -} - -size_t -atrck_get_nodes_count(const struct atrck* atrck) -{ - ASSERT(atrck); - return atrck->nnodes; -} - -res_T -atrck_get_band - (const struct atrck* atrck, - const size_t iband, - struct atrck_band* atrck_band) -{ - const struct band* band = NULL; - res_T res = RES_OK; - - if(!atrck || !atrck_band) { - res = RES_BAD_ARG; - goto error; - } - - if(iband >= atrck_get_bands_count(atrck)) { - log_err(atrck, "%s: invalid band index %lu.\n", - FUNC_NAME, (unsigned long)iband); - res = RES_BAD_ARG; - goto error; - } - - band = darray_band_cdata_get(&atrck->bands) + iband; - atrck_band->lower = band->low; - atrck_band->upper = band->upp; - atrck_band->quad_pts_count = darray_quad_pt_size_get(&band->quad_pts); - atrck_band->id = (size_t)band->id; - atrck_band->atrck__ = atrck; - atrck_band->band__ = band; - -exit: - return res; -error: - goto exit; -} - -res_T -atrck_band_get_quad_pt - (const struct atrck_band* atrck_band, - const size_t iquad_pt, - struct atrck_quad_pt* atrck_quad_pt) -{ - const struct band* band = NULL; - const struct quad_pt* quad_pt = NULL; - res_T res = RES_OK; - - if(!atrck_band || !atrck_quad_pt) { - res = RES_BAD_ARG; - goto error; - } - - band = atrck_band->band__; - if(iquad_pt >= atrck_band->quad_pts_count) { - log_err(atrck_band->atrck__, - "%s: band %lu: invalid quadrature point index %lu.\n", - FUNC_NAME, (unsigned long)atrck_band->id, (unsigned long)iquad_pt); - res = RES_BAD_ARG; - goto error; - } - - quad_pt = darray_quad_pt_cdata_get(&band->quad_pts) + iquad_pt; - atrck_quad_pt->ka_list = quad_pt->ka_list; - atrck_quad_pt->abscissa = quad_pt->abscissa; - atrck_quad_pt->weight = quad_pt->weight; - atrck_quad_pt->id = iquad_pt; - -exit: - return res; -error: - goto exit; -} - -res_T -atrck_band_sample_quad_pt - (const struct atrck_band* atrck_band, - const double r, /* Canonical random number in [0, 1[ */ - struct atrck_quad_pt* quad_pt) -{ - const struct band* band = NULL; - const double* cumul = NULL; - const double* find = NULL; - double r_next = nextafter(r, DBL_MAX); - size_t i; - - if(!atrck_band || !quad_pt) return RES_BAD_ARG; - - if(r < 0 || r >= 1) { - log_err(atrck_band->atrck__, - "%s: invalid canonical random number `%g'.\n", FUNC_NAME, r); - return RES_BAD_ARG; - } - - band = atrck_band->band__; - cumul = darray_double_cdata_get(&band->quad_pts_cumul); - - /* Use r_next rather than r in order to find the first entry that is not less - * than *or equal* to r */ - find = search_lower_bound(&r_next, cumul, atrck_band->quad_pts_count, - sizeof(double), cmp_dbl); - ASSERT(find); - - /* Compute the index of the found quadrature point */ - i = (size_t)(find - cumul); - ASSERT(i < atrck_band->quad_pts_count); - ASSERT(cumul[i] > r); - ASSERT(!i || cumul[i-1] <= r); - - /* Fetch the sampled quadrature point */ - atrck_band_get_quad_pt(atrck_band, i, quad_pt); - - return RES_OK; -} - -/******************************************************************************* - * Local functions - ******************************************************************************/ -void -quad_pt_release(struct quad_pt* quad_pt) -{ - ASSERT(quad_pt); - if(quad_pt->ka_list && quad_pt->ka_list != MAP_FAILED) - munmap(quad_pt->ka_list, quad_pt->map_len); -} diff --git a/src/atrck.h b/src/atrck.h @@ -1,126 +0,0 @@ -/* Copyright (C) 2020, 2021 CNRS - * - * 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/>. */ - -#ifndef ATRCK_H -#define ATRCK_H - -#include <rsys/rsys.h> - -/* Library symbol management */ -#if defined(ATRCK_SHARED_BUILD) /* Build shared library */ - #define ATRCK_API extern EXPORT_SYM -#elif defined(ATRCK_STATIC) /* Use/build static library */ - #define ATRCK_API extern LOCAL_SYM -#else /* Use shared library */ - #define ATRCK_API extern IMPORT_SYM -#endif - -/* Helper macro that asserts if the invocation of the atrck function `Func' - * returns an error. One should use this macro on sth function calls for - * which no explicit error checking is performed */ -#ifndef NDEBUG - #define ATRCK(Func) ASSERT(atrck_ ## Func == RES_OK) -#else - #define ATRCK(Func) atrck_ ## Func -#endif - -struct atrck_band { - double lower; /* Lower band wavenumber in cm^-1 */ - double upper; /* Upper band wavenumber in cm^-1 */ - size_t quad_pts_count; /* #quadrature points */ - size_t id; - - /* Internal data */ - const struct atrck* atrck__; - const void* band__; -}; -static const struct atrck_band ATRCK_BAND_NULL; - -struct atrck_quad_pt { - float* ka_list; /* Per node ka */ - double abscissa; /* m^-1 */ - double weight; - size_t id; -}; -static const struct atrck_quad_pt ATRCK_QUAD_PT_NULL; - -/* Forward declaration of external data types */ -struct logger; -struct mem_allocator; - -/* Forward declaration of opaque data types */ -struct atrck; - -BEGIN_DECLS - -/******************************************************************************* - * AtrKC API - ******************************************************************************/ -ATRCK_API res_T -atrck_create - (struct logger* logger, /* NULL <=> use default logger */ - struct mem_allocator* allocator, /* NULL <=> use default allocator */ - const int verbose, /* Verbosity level */ - struct atrck** atrck); - -ATRCK_API res_T -atrck_ref_get - (struct atrck* atrck); - -ATRCK_API res_T -atrck_ref_put - (struct atrck* atrck); - -ATRCK_API res_T -atrck_load - (struct atrck* atrck, - const char* path); - -ATRCK_API res_T -atrck_load_stream - (struct atrck* atrck, - FILE* stream, - const char* stream_name); /* Can be NULL */ - -ATRCK_API size_t -atrck_get_bands_count - (const struct atrck* atrck); - -ATRCK_API size_t -atrck_get_nodes_count - (const struct atrck* atrck); - -ATRCK_API res_T -atrck_get_band - (const struct atrck* atrck, - const size_t iband, - struct atrck_band* band); - -ATRCK_API res_T -atrck_band_get_quad_pt - (const struct atrck_band* band, - const size_t iquad_pt, - struct atrck_quad_pt* quad_pt); - -ATRCK_API res_T -atrck_band_sample_quad_pt - (const struct atrck_band* band, - const double r, /* Canonical random number in [0, 1[ */ - struct atrck_quad_pt* quad_pt); - -END_DECLS - -#endif /* ATRCK_H */ - diff --git a/src/atrck_c.h b/src/atrck_c.h @@ -1,139 +0,0 @@ -/* Copyright (C) 2020, 2021 CNRS - * - * 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/>. */ - -#ifndef ATRCK_C_H -#define ATRCK_C_H - -#include <rsys/dynamic_array_double.h> -#include <rsys/logger.h> -#include <rsys/ref_count.h> - -struct mem_allocator; - -struct quad_pt { - float* ka_list; /* Per node ka */ - size_t map_len; - double abscissa; - double weight; -}; - -static INLINE void -quad_pt_init(struct mem_allocator* allocator, struct quad_pt* quad) -{ - ASSERT(quad); - (void)allocator; - quad->ka_list = NULL; - quad->map_len = 0;; - quad->abscissa = 0; - quad->weight = 0; -} - -extern LOCAL_SYM INLINE void -quad_pt_release - (struct quad_pt* quad); - -/* Define the dynamic array of quadrature points */ -#define DARRAY_NAME quad_pt -#define DARRAY_DATA struct quad_pt -#include <rsys/dynamic_array.h> - -struct band { - uint64_t id; - double low; /* Lower bound in cm^-1 */ - double upp; /* Upper bound in cm^-1 */ - struct darray_quad_pt quad_pts; - struct darray_double quad_pts_cumul; -}; - -static INLINE void -band_init(struct mem_allocator* allocator, struct band* band) -{ - ASSERT(band); - band->id = UINT64_MAX; - band->low = DBL_MAX; - band->upp =-DBL_MAX; - darray_quad_pt_init(allocator, &band->quad_pts); - darray_double_init(allocator, &band->quad_pts_cumul); -} - -static INLINE void -band_release(struct band* band) -{ - ASSERT(band); - darray_quad_pt_release(&band->quad_pts); - darray_double_release(&band->quad_pts_cumul); -} - -static INLINE res_T -band_copy(struct band* dst, const struct band* src) -{ - res_T res = RES_OK; - ASSERT(dst && src); - - dst->id = src->id; - dst->low = src->low; - dst->upp = src->upp; - res = darray_quad_pt_copy(&dst->quad_pts, &src->quad_pts); - if(res != RES_OK) return res; - res = darray_double_copy(&dst->quad_pts_cumul, &src->quad_pts_cumul); - if(res != RES_OK) return res; - return RES_OK; -} - -static INLINE res_T -band_copy_and_release(struct band* dst, struct band* src) -{ - res_T res = RES_OK; - ASSERT(dst && src); - - dst->id = src->id; - dst->low = src->low; - dst->upp = src->upp; - res = darray_quad_pt_copy_and_release(&dst->quad_pts, &src->quad_pts); - if(res != RES_OK) return res; - res = darray_double_copy_and_release - (&dst->quad_pts_cumul, &src->quad_pts_cumul); - if(res != RES_OK) return res; - return RES_OK; - -} - -/* Define the dynamic array of bands */ -#define DARRAY_NAME band -#define DARRAY_DATA struct band -#define DARRAY_FUNCTOR_INIT band_init -#define DARRAY_FUNCTOR_RELEASE band_release -#define DARRAY_FUNCTOR_COPY band_copy -#define DARRAY_FUNCTOR_COPY_AND_RELEASE band_copy_and_release -#include <rsys/dynamic_array.h> - -struct mem_allocator; - -struct atrck { - /* Loaded data */ - uint64_t pagesize; - uint64_t nnodes; - struct darray_band bands; - - size_t pagesize_os; - - struct mem_allocator* allocator; - struct logger* logger; - struct logger logger__; /* Default logger */ - int verbose; - ref_T ref; -}; - -#endif /* ATRCK_C_H */ diff --git a/src/atrck_log.c b/src/atrck_log.c @@ -1,125 +0,0 @@ -/* Copyright (C) 2020, 2021 CNRS - * - * 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/>. */ - -#include "atrck_c.h" -#include "atrck_log.h" - -#include <rsys/cstr.h> -#include <rsys/logger.h> - -#include <stdarg.h> - -/******************************************************************************* - * Helper functions - ******************************************************************************/ -static INLINE void -log_msg - (const struct atrck* atrck, - const enum log_type stream, - const char* msg, - va_list vargs) -{ - ASSERT(atrck && msg); - if(atrck->verbose) { - res_T res; (void)res; - res = logger_vprint(atrck->logger, stream, msg, vargs); - ASSERT(res == RES_OK); - } -} - -static void -print_info(const char* msg, void* ctx) -{ - (void)ctx; - fprintf(stderr, MSG_INFO_PREFIX"%s", msg); -} - -static void -print_err(const char* msg, void* ctx) -{ - (void)ctx; - fprintf(stderr, MSG_ERROR_PREFIX"%s", msg); -} - -static void -print_warn(const char* msg, void* ctx) -{ - (void)ctx; - fprintf(stderr, MSG_WARNING_PREFIX"%s", msg); -} - -/******************************************************************************* - * Local functions - ******************************************************************************/ -res_T -setup_log_default(struct atrck* atrck) -{ - res_T res = RES_OK; - ASSERT(atrck); - - res = logger_init(atrck->allocator, &atrck->logger__); - if(res != RES_OK) { - if(atrck->verbose) { - fprintf(stderr, - MSG_ERROR_PREFIX - "Could not setup the AtrCK default logger -- %s.\n", - res_to_cstr(res)); - } - goto error; - } - logger_set_stream(&atrck->logger__, LOG_OUTPUT, print_info, NULL); - logger_set_stream(&atrck->logger__, LOG_ERROR, print_err, NULL); - logger_set_stream(&atrck->logger__, LOG_WARNING, print_warn, NULL); - atrck->logger = &atrck->logger__; - -exit: - return res; -error: - goto exit; -} - -void -log_info(const struct atrck* atrck, const char* msg, ...) -{ - va_list vargs_list; - ASSERT(atrck && msg); - - va_start(vargs_list, msg); - log_msg(atrck, LOG_OUTPUT, msg, vargs_list); - va_end(vargs_list); -} - -void -log_err(const struct atrck* atrck, const char* msg, ...) -{ - va_list vargs_list; - ASSERT(atrck && msg); - - va_start(vargs_list, msg); - log_msg(atrck, LOG_ERROR, msg, vargs_list); - va_end(vargs_list); -} - -void -log_warn(const struct atrck* atrck, const char* msg, ...) -{ - va_list vargs_list; - ASSERT(atrck && msg); - - va_start(vargs_list, msg); - log_msg(atrck, LOG_WARNING, msg, vargs_list); - va_end(vargs_list); -} - diff --git a/src/atrck_log.h b/src/atrck_log.h @@ -1,68 +0,0 @@ -/* Copyright (C) 2020, 2021 CNRS - * - * 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/>. */ - -#ifndef ATRCK_LOG_H -#define ATRCK_LOG_H - -#include <rsys/rsys.h> - -#define MSG_INFO_PREFIX "AtrCK:\x1b[1m\x1b[32minfo\x1b[0m: " -#define MSG_ERROR_PREFIX "AtrCK:\x1b[1m\x1b[31merror\x1b[0m: " -#define MSG_WARNING_PREFIX "AtrCK:\x1b[1m\x1b[33mwarning\x1b[0m: " - -struct atrck; -struct logger; - -extern LOCAL_SYM res_T -setup_log_default - (struct atrck* atrck); - -/* Conditionally log a message on the LOG_OUTPUT stream of the atrck logger, - * with respect to its verbose flag */ -extern LOCAL_SYM void -log_info - (const struct atrck* atrck, - const char* msg, - ...) -#ifdef COMPILER_GCC - __attribute((format(printf, 2, 3))) -#endif -; - -/* Conditionally log a message on the LOG_ERROR stream of the atrck logger, - * with respect to its verbose flag */ -extern LOCAL_SYM void -log_err - (const struct atrck* atrck, - const char* msg, - ...) -#ifdef COMPILER_GCC - __attribute((format(printf, 2, 3))) -#endif -; - -/* Conditionally log a message on the LOG_WARNING stream of the atrck logger, - * with respect to its verbose flag */ -extern LOCAL_SYM void -log_warn - (const struct atrck* atrck, - const char* msg, - ...) -#ifdef COMPILER_GCC - __attribute((format(printf, 2, 3))) -#endif -; - -#endif /* ATRCK_LOG_H */ diff --git a/src/sck.c b/src/sck.c @@ -0,0 +1,581 @@ +/* Copyright (C) 2020, 2021 CNRS + * + * 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/>. */ + +#define _POSIX_C_SOURCE 200809L /* mmap support */ +#define _DEFAULT_SOURCE 1 /* MAP_POPULATE support */ +#define _BSD_SOURCE 1 /* MAP_POPULATE for glibc < 2.19 */ + +#include "sck.h" +#include "sck_c.h" +#include "sck_log.h" + +#include <rsys/algorithm.h> + +#include <unistd.h> /* sysconf support */ + +#include <errno.h> +#include <string.h> +#include <sys/mman.h> /* mmap */ + +/******************************************************************************* + * Helper functions + ******************************************************************************/ +static INLINE res_T +check_sck_create_args(const struct sck_create_args* args) +{ + /* Nothing to check. Only return RES_BAD_ARG if args is NULL */ + return args ? RES_OK : RES_BAD_ARG; +} + +static void +reset_sck(struct sck* sck) +{ + ASSERT(sck); + sck->pagesize = 0; + sck->nnodes = 0; + darray_band_purge(&sck->bands); +} + +static INLINE int +cmp_dbl(const void* a, const void* b) +{ + const double d0 = *((const double*)a); + const double d1 = *((const double*)b); + return d0 < d1 ? -1 : (d0 > d1 ? 1 : 0); +} + +static res_T +read_quad_pt + (struct sck* sck, + struct quad_pt* quad_pt, + FILE* stream, + const char* stream_name, + size_t iband, + size_t iquad_pt) +{ + res_T res = RES_OK; + ASSERT(sck && quad_pt && stream_name); + + if(fread(&quad_pt->weight, sizeof(quad_pt->weight), 1, stream) != 1) { + log_err(sck, + "%s: band %lu: quadrature point %lu: could not read the weight.\n", + stream_name, iband, iquad_pt); + res = RES_IO_ERR; + goto error; + } + + if(quad_pt->weight < 0) { + log_err(sck, + "%s: band %lu: quadrature point %lu: invalid weight %g.\n", + stream_name, iband, iquad_pt, quad_pt->weight); + res = RES_BAD_ARG; + goto error; + + } + + quad_pt->ka_list = NULL; + +exit: + return res; +error: + goto exit; +} + +static res_T +read_band + (struct sck* sck, + struct band* band, + FILE* stream, + const char* stream_name) +{ + double cumul; + size_t iquad_pt; + size_t iband; + uint64_t nquad_pts; + res_T res = RES_OK; + ASSERT(sck && band && stream_name); + + iband = (size_t)(band - darray_band_cdata_get(&sck->bands)); + + /* Read band definition */ + #define READ(Var, Name) { \ + if(fread((Var), sizeof(*(Var)), 1, stream) != 1) { \ + log_err(sck, "%s: band %lu: could not read the %s.\n", \ + stream_name, (unsigned long)iband, (Name)); \ + res = RES_IO_ERR; \ + goto error; \ + } \ + } (void)0 + READ(&band->low, "band lower bound"); + READ(&band->upp, "band upper bound"); + READ(&nquad_pts, "#quadrature points"); + #undef READ + + /* Check band description */ + if(band->low < 0 || band->low > band->upp) { + log_err(sck, + "%s: band %lu: invalid band range [%g, %g].\n", + stream_name, (unsigned long)iband, band->low, band->upp); + res = RES_BAD_ARG; + goto error; + } + if(nquad_pts == 0) { + log_err(sck, + "%s: band %lu: invalid number fo quadrature points (#points=%lu).\n", + stream_name, (unsigned long)iband, (unsigned long)nquad_pts); + res = RES_BAD_ARG; + goto error; + } + + /* Allocate the list of quadrature points */ + res = darray_quad_pt_resize(&band->quad_pts, nquad_pts); + if(res != RES_OK) { + log_err(sck, + "%s: band %lu: could not allocate the list of quadrature points " + "(#points=%lu).\n", + stream_name, (unsigned long)iband, (unsigned long)nquad_pts); + goto error; + } + + /* Allocate the cumulative used to sample the quadrature points */ + res = darray_double_resize(&band->quad_pts_cumul, nquad_pts); + if(res != RES_OK) { + log_err(sck, + "%s: band %lu: could not allocate the cumulative of quadrature points " + "(#points=%lu).\n", + stream_name, (unsigned long)iband, (unsigned long)nquad_pts); + goto error; + } + + /* Read the data of the quadrature points of the band */ + cumul = 0; + FOR_EACH(iquad_pt, 0, nquad_pts) { + struct quad_pt* quad_pt = darray_quad_pt_data_get(&band->quad_pts)+iquad_pt; + + /* Read current quadrature point */ + res = read_quad_pt(sck, quad_pt, stream, stream_name, iband, iquad_pt); + if(res != RES_OK) goto error; + + /* Compute the quadrature point cumulative */ + cumul += quad_pt->weight; + darray_double_data_get(&band->quad_pts_cumul)[iquad_pt] = cumul; + } + + /* The weights are not normalized */ + if(!eq_eps(cumul, 1.0, 1.e-6)) { + log_warn(sck, + "%s: band %lu: the weights of the quadrature points are not normalised.\n", + stream_name, (unsigned long)iband); + + /* Renormalize the cumulative */ + FOR_EACH(iquad_pt, 0, nquad_pts) { + darray_double_data_get(&band->quad_pts_cumul)[iquad_pt] /= cumul; + } + } + + /* Handle imprecision issue */ + darray_double_data_get(&band->quad_pts_cumul)[nquad_pts-1] = 1.0; + +exit: + return res; +error: + goto exit; +} + +static res_T +load_stream(struct sck* sck, FILE* stream, const char* stream_name) +{ + size_t map_len; + size_t iband; + uint64_t nbands; + off_t offset = 0; + res_T res = RES_OK; + ASSERT(sck && stream && stream_name); + + reset_sck(sck); + + /* Read file header */ + #define READ(Var, Name) { \ + if(fread((Var), sizeof(*(Var)), 1, stream) != 1) { \ + log_err(sck, "%s: could not read the %s.\n", stream_name, (Name)); \ + res = RES_IO_ERR; \ + goto error; \ + } \ + } (void)0 + READ(&sck->pagesize, "page size"); + READ(&nbands, "number of bands"); + READ(&sck->nnodes, "number of nodes"); + #undef READ + + /* Check band description */ + if(!IS_ALIGNED(sck->pagesize, sck->pagesize_os)) { + log_err(sck, + "%s: invalid page size %lu. The page size attribute must be aligned on " + "the page size of the operating system (%lu).\n", + stream_name, + (unsigned long)sck->pagesize, + (unsigned long)sck->pagesize_os); + res = RES_BAD_ARG; + goto error; + } + if(!nbands) { + log_err(sck, "%s: invalid number of bands %lu.\n", + stream_name, (unsigned long)nbands); + res = RES_BAD_ARG; + goto error; + } + if(!sck->nnodes) { + log_err(sck, "%s: invalid number of nodes %lu.\n", + stream_name, (unsigned long)sck->nnodes); + res = RES_BAD_ARG; + goto error; + } + + /* Allocate the bands */ + res = darray_band_resize(&sck->bands, nbands); + if(res != RES_OK) { + log_err(sck, "%s: could not allocate the list of bands (#bands=%lu).\n", + stream_name, (unsigned long)nbands); + goto error; + } + + /* Read the band description */ + FOR_EACH(iband, 0, nbands) { + struct band* band = darray_band_data_get(&sck->bands) + iband; + res = read_band(sck, band, stream, stream_name); + if(res != RES_OK) goto error; + } + + /* Compute the length in bytes of the k to map for each band/quadrature point */ + map_len = ALIGN_SIZE(sck->nnodes * sizeof(float), sck->pagesize); + + /* Compute the offset toward the 1st list of radiative coefficients */ + offset = ftell(stream); + offset = (off_t)ALIGN_SIZE((uint64_t)offset, sck->pagesize); + + FOR_EACH(iband, 0, nbands) { + struct band* band = NULL; + size_t iquad_pt; + + /* Map the per band scattering coefficient */ + band = darray_band_data_get(&sck->bands) + iband; + band->map_len = map_len; + band->ks_list = mmap(NULL, band->map_len, PROT_READ, + MAP_PRIVATE|MAP_POPULATE, fileno(stream), offset); + if(band->ks_list == MAP_FAILED) { + log_err(sck, + "%s: band %lu: could not map the scattering coefficients -- %s.\n", + stream_name, (unsigned long)iband, strerror(errno)); + res = RES_IO_ERR; + goto error; + } + offset = (off_t)((size_t)offset + map_len); + + FOR_EACH(iquad_pt, 0, darray_quad_pt_size_get(&band->quad_pts)) { + struct quad_pt* quad_pt = NULL; + + /* Map the per band absorption coefficient */ + quad_pt = darray_quad_pt_data_get(&band->quad_pts)+iquad_pt; + quad_pt->map_len = map_len; + quad_pt->ka_list = mmap(NULL, quad_pt->map_len, PROT_READ, + MAP_PRIVATE|MAP_POPULATE, fileno(stream), offset); + if(quad_pt->ka_list == MAP_FAILED) { + log_err(sck, + "%s: band %lu: quadrature point %lu: coult not map the correlated K " + "-- %s.\n", + stream_name, + (unsigned long)iband, + (unsigned long)iquad_pt, + strerror(errno)); + res = RES_IO_ERR; + goto error; + } + offset = (off_t)((size_t)offset + map_len); + } + } + +exit: + return res; +error: + reset_sck(sck); + goto exit; +} + +static void +release_sck(ref_T* ref) +{ + struct sck* sck; + ASSERT(ref); + sck = CONTAINER_OF(ref, struct sck, ref); + if(sck->logger == &sck->logger__) logger_release(&sck->logger__); + darray_band_release(&sck->bands); + MEM_RM(sck->allocator, sck); +} + +/******************************************************************************* + * Exported functions + ******************************************************************************/ +res_T +sck_create + (const struct sck_create_args* args, + struct sck** out_sck) +{ + struct sck* sck = NULL; + struct mem_allocator* allocator = NULL; + res_T res = RES_OK; + + if(!out_sck) { res = RES_BAD_ARG; goto error; } + res = check_sck_create_args(args); + if(res != RES_OK) goto error; + + allocator = args->allocator ? args->allocator : &mem_default_allocator; + sck = MEM_CALLOC(allocator, 1, sizeof(*sck)); + if(!sck) { + if(args->verbose) { + #define ERR_STR "Could not allocate the Star-CK device.\n" + if(args->logger) { + logger_print(args->logger, LOG_ERROR, ERR_STR); + } else { + fprintf(stderr, MSG_ERROR_PREFIX ERR_STR); + } + #undef ERR_STR + } + res = RES_MEM_ERR; + goto error; + } + ref_init(&sck->ref); + sck->allocator = allocator; + sck->verbose = args->verbose; + sck->pagesize_os = (size_t)sysconf(_SC_PAGESIZE); + darray_band_init(allocator, &sck->bands); + if(args->logger) { + sck->logger = args->logger; + } else { + setup_log_default(sck); + } + +exit: + if(out_sck) *out_sck = sck; + return res; +error: + if(sck) { + SCK(ref_put(sck)); + sck = NULL; + } + goto exit; +} + +res_T +sck_ref_get(struct sck* sck) +{ + if(!sck) return RES_BAD_ARG; + ref_get(&sck->ref); + return RES_OK; +} + +res_T +sck_ref_put(struct sck* sck) +{ + if(!sck) return RES_BAD_ARG; + ref_put(&sck->ref, release_sck); + return RES_OK; +} + +res_T +sck_load(struct sck* sck, const char* path) +{ + FILE* file = NULL; + res_T res = RES_OK; + + if(!sck || !path) { + res = RES_BAD_ARG; + goto error; + } + + file = fopen(path, "r"); + if(!file) { + log_err(sck, "%s: error opening file `%s'.\n", FUNC_NAME, path); + res = RES_IO_ERR; + goto error; + } + + res = load_stream(sck, file, path); + if(res != RES_OK) goto error; + +exit: + if(file) fclose(file); + return res; +error: + goto exit; +} + +res_T +sck_load_stream + (struct sck* sck, + FILE* stream, + const char* stream_name) +{ + if(!sck || !stream) return RES_BAD_ARG; + return load_stream(sck, stream, stream_name ? stream_name : "<stream>"); +} + +size_t +sck_get_bands_count(const struct sck* sck) +{ + ASSERT(sck); + return darray_band_size_get(&sck->bands); +} + +size_t +sck_get_nodes_count(const struct sck* sck) +{ + ASSERT(sck); + return sck->nnodes; +} + +res_T +sck_get_band + (const struct sck* sck, + const size_t iband, + struct sck_band* sck_band) +{ + const struct band* band = NULL; + res_T res = RES_OK; + + if(!sck || !sck_band) { + res = RES_BAD_ARG; + goto error; + } + + if(iband >= sck_get_bands_count(sck)) { + log_err(sck, "%s: invalid band index %lu.\n", + FUNC_NAME, (unsigned long)iband); + res = RES_BAD_ARG; + goto error; + } + + band = darray_band_cdata_get(&sck->bands) + iband; + sck_band->lower = band->low; + sck_band->upper = band->upp; + sck_band->ks_list = band->ks_list; + sck_band->quad_pts_count = darray_quad_pt_size_get(&band->quad_pts); + sck_band->id = iband; + sck_band->sck__ = sck; + sck_band->band__ = band; + +exit: + return res; +error: + goto exit; +} + +res_T +sck_band_get_quad_pt + (const struct sck_band* sck_band, + const size_t iquad_pt, + struct sck_quad_pt* sck_quad_pt) +{ + const struct band* band = NULL; + const struct quad_pt* quad_pt = NULL; + res_T res = RES_OK; + + if(!sck_band || !sck_quad_pt) { + res = RES_BAD_ARG; + goto error; + } + + band = sck_band->band__; + if(iquad_pt >= sck_band->quad_pts_count) { + log_err(sck_band->sck__, + "%s: band %lu: invalid quadrature point index %lu.\n", + FUNC_NAME, (unsigned long)sck_band->id, (unsigned long)iquad_pt); + res = RES_BAD_ARG; + goto error; + } + + quad_pt = darray_quad_pt_cdata_get(&band->quad_pts) + iquad_pt; + sck_quad_pt->ka_list = quad_pt->ka_list; + sck_quad_pt->weight = quad_pt->weight; + sck_quad_pt->id = iquad_pt; + +exit: + return res; +error: + goto exit; +} + +res_T +sck_band_sample_quad_pt + (const struct sck_band* sck_band, + const double r, /* Canonical random number in [0, 1[ */ + struct sck_quad_pt* quad_pt) +{ + const struct band* band = NULL; + const double* cumul = NULL; + const double* find = NULL; + double r_next = nextafter(r, DBL_MAX); + size_t i; + + if(!sck_band || !quad_pt) return RES_BAD_ARG; + + if(r < 0 || r >= 1) { + log_err(sck_band->sck__, + "%s: invalid canonical random number `%g'.\n", FUNC_NAME, r); + return RES_BAD_ARG; + } + + band = sck_band->band__; + cumul = darray_double_cdata_get(&band->quad_pts_cumul); + + /* Use r_next rather than r in order to find the first entry that is not less + * than *or equal* to r */ + find = search_lower_bound(&r_next, cumul, sck_band->quad_pts_count, + sizeof(double), cmp_dbl); + ASSERT(find); + + /* Compute the index of the found quadrature point */ + i = (size_t)(find - cumul); + ASSERT(i < sck_band->quad_pts_count); + ASSERT(cumul[i] > r); + ASSERT(!i || cumul[i-1] <= r); + + /* Fetch the sampled quadrature point */ + sck_band_get_quad_pt(sck_band, i, quad_pt); + + return RES_OK; +} + +/******************************************************************************* + * Local functions + ******************************************************************************/ +void +band_release(struct band* band) +{ + ASSERT(band); + darray_quad_pt_release(&band->quad_pts); + darray_double_release(&band->quad_pts_cumul); + if(band->ks_list && band->ks_list != MAP_FAILED) { + munmap(band->ks_list, band->map_len); + } +} + +void +quad_pt_release(struct quad_pt* quad_pt) +{ + ASSERT(quad_pt); + if(quad_pt->ka_list && quad_pt->ka_list != MAP_FAILED) { + munmap(quad_pt->ka_list, quad_pt->map_len); + } +} diff --git a/src/sck.h b/src/sck.h @@ -0,0 +1,136 @@ +/* Copyright (C) 2020, 2021 CNRS + * + * 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/>. */ + +#ifndef SCK_H +#define SCK_H + +#include <rsys/rsys.h> + +/* Library symbol management */ +#if defined(SCK_SHARED_BUILD) /* Build shared library */ + #define SCK_API extern EXPORT_SYM +#elif defined(SCK_STATIC) /* Use/build static library */ + #define SCK_API extern LOCAL_SYM +#else /* Use shared library */ + #define SCK_API extern IMPORT_SYM +#endif + +/* Helper macro that asserts if the invocation of the sck function `Func' + * returns an error. One should use this macro on sck function calls for which + * no explicit error checking is performed */ +#ifndef NDEBUG + #define SCK(Func) ASSERT(sck_ ## Func == RES_OK) +#else + #define SCK(Func) sck_ ## Func +#endif + +/* Forward declaration of external data types */ +struct logger; +struct mem_allocator; + +struct sck_create_args { + struct logger* logger; /* May be NULL <=> default logger */ + struct mem_allocator* allocator; /* NULL <=> use default allocator */ + int verbose; /* Verbosity level */ +}; +#define SCK_CREATE_ARGS_DEFAULT__ {NULL, NULL, 0} +static const struct sck_create_args SCK_CREATE_ARGS_DEFAULT = + SCK_CREATE_ARGS_DEFAULT__; + +struct sck_band { + double lower; /* Lower band wavenumber in cm^-1 */ + double upper; /* Upper band wavenumber in cm^-1 */ + size_t quad_pts_count; /* #quadrature points */ + size_t id; + + float* ks_list; /* Per node ks */ + + /* Internal data */ + const struct sck* sck__; + const void* band__; +}; +#define SCK_BAND_NULL__ {0, 0, 0, 0, NULL, NULL, NULL} +static const struct sck_band SCK_BAND_NULL = SCK_BAND_NULL__; + +struct sck_quad_pt { + float* ka_list; /* Per node ka */ + double weight; + size_t id; +}; +#define SCK_QUAD_PT_NULL__ {NULL, 0, 0} +static const struct sck_quad_pt SCK_QUAD_PT_NULL = SCK_QUAD_PT_NULL__; + +/* Forward declaration of opaque data types */ +struct sck; + +BEGIN_DECLS + +/******************************************************************************* + * AtrKC API + ******************************************************************************/ +SCK_API res_T +sck_create + (const struct sck_create_args* args, + struct sck** sck); + +SCK_API res_T +sck_ref_get + (struct sck* sck); + +SCK_API res_T +sck_ref_put + (struct sck* sck); + +SCK_API res_T +sck_load + (struct sck* sck, + const char* path); + +SCK_API res_T +sck_load_stream + (struct sck* sck, + FILE* stream, + const char* stream_name); /* Can be NULL */ + +SCK_API size_t +sck_get_bands_count + (const struct sck* sck); + +SCK_API size_t +sck_get_nodes_count + (const struct sck* sck); + +SCK_API res_T +sck_get_band + (const struct sck* sck, + const size_t iband, + struct sck_band* band); + +SCK_API res_T +sck_band_get_quad_pt + (const struct sck_band* band, + const size_t iquad_pt, + struct sck_quad_pt* quad_pt); + +SCK_API res_T +sck_band_sample_quad_pt + (const struct sck_band* band, + const double r, /* Canonical random number in [0, 1[ */ + struct sck_quad_pt* quad_pt); + +END_DECLS + +#endif /* SCK_H */ + diff --git a/src/sck_c.h b/src/sck_c.h @@ -0,0 +1,136 @@ +/* Copyright (C) 2020, 2021 CNRS + * + * 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/>. */ + +#ifndef SCK_C_H +#define SCK_C_H + +#include <rsys/dynamic_array_double.h> +#include <rsys/logger.h> +#include <rsys/ref_count.h> + +struct mem_allocator; + +struct quad_pt { + float* ka_list; /* Per node ka */ + size_t map_len; + double weight; +}; + +static INLINE void +quad_pt_init(struct mem_allocator* allocator, struct quad_pt* quad) +{ + ASSERT(quad); + (void)allocator; + quad->ka_list = NULL; + quad->map_len = 0; + quad->weight = 0; +} + +extern LOCAL_SYM void +quad_pt_release + (struct quad_pt* quad); + +/* Define the dynamic array of quadrature points */ +#define DARRAY_NAME quad_pt +#define DARRAY_DATA struct quad_pt +#include <rsys/dynamic_array.h> + +struct band { + double low; /* Lower bound in nm */ + double upp; /* Upper bound in nm */ + size_t map_len; + float* ks_list; /* Per node ks */ + struct darray_quad_pt quad_pts; + struct darray_double quad_pts_cumul; +}; + +static INLINE void +band_init(struct mem_allocator* allocator, struct band* band) +{ + ASSERT(band); + band->low = DBL_MAX; + band->upp =-DBL_MAX; + band->map_len = 0; + band->ks_list = NULL; + darray_quad_pt_init(allocator, &band->quad_pts); + darray_double_init(allocator, &band->quad_pts_cumul); +} + +extern LOCAL_SYM INLINE void +band_release + (struct band* band); + +static INLINE res_T +band_copy(struct band* dst, const struct band* src) +{ + res_T res = RES_OK; + ASSERT(dst && src); + + dst->low = src->low; + dst->upp = src->upp; + dst->map_len = src->map_len; + dst->ks_list = src->ks_list; + res = darray_quad_pt_copy(&dst->quad_pts, &src->quad_pts); + if(res != RES_OK) return res; + res = darray_double_copy(&dst->quad_pts_cumul, &src->quad_pts_cumul); + if(res != RES_OK) return res; + return RES_OK; +} + +static INLINE res_T +band_copy_and_release(struct band* dst, struct band* src) +{ + res_T res = RES_OK; + ASSERT(dst && src); + + dst->low = src->low; + dst->upp = src->upp; + dst->map_len = src->map_len; + dst->ks_list = src->ks_list; + res = darray_quad_pt_copy_and_release(&dst->quad_pts, &src->quad_pts); + if(res != RES_OK) return res; + res = darray_double_copy_and_release + (&dst->quad_pts_cumul, &src->quad_pts_cumul); + if(res != RES_OK) return res; + return RES_OK; +} + +/* Define the dynamic array of bands */ +#define DARRAY_NAME band +#define DARRAY_DATA struct band +#define DARRAY_FUNCTOR_INIT band_init +#define DARRAY_FUNCTOR_RELEASE band_release +#define DARRAY_FUNCTOR_COPY band_copy +#define DARRAY_FUNCTOR_COPY_AND_RELEASE band_copy_and_release +#include <rsys/dynamic_array.h> + +struct mem_allocator; + +struct sck { + /* Loaded data */ + uint64_t pagesize; + uint64_t nnodes; + struct darray_band bands; + + size_t pagesize_os; + + struct mem_allocator* allocator; + struct logger* logger; + struct logger logger__; /* Default logger */ + int verbose; + ref_T ref; +}; + +#endif /* SCK_C_H */ diff --git a/src/sck_log.c b/src/sck_log.c @@ -0,0 +1,125 @@ +/* Copyright (C) 2020, 2021 CNRS + * + * 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/>. */ + +#include "sck_c.h" +#include "sck_log.h" + +#include <rsys/cstr.h> +#include <rsys/logger.h> + +#include <stdarg.h> + +/******************************************************************************* + * Helper functions + ******************************************************************************/ +static INLINE void +log_msg + (const struct sck* sck, + const enum log_type stream, + const char* msg, + va_list vargs) +{ + ASSERT(sck && msg); + if(sck->verbose) { + res_T res; (void)res; + res = logger_vprint(sck->logger, stream, msg, vargs); + ASSERT(res == RES_OK); + } +} + +static void +print_info(const char* msg, void* ctx) +{ + (void)ctx; + fprintf(stderr, MSG_INFO_PREFIX"%s", msg); +} + +static void +print_err(const char* msg, void* ctx) +{ + (void)ctx; + fprintf(stderr, MSG_ERROR_PREFIX"%s", msg); +} + +static void +print_warn(const char* msg, void* ctx) +{ + (void)ctx; + fprintf(stderr, MSG_WARNING_PREFIX"%s", msg); +} + +/******************************************************************************* + * Local functions + ******************************************************************************/ +res_T +setup_log_default(struct sck* sck) +{ + res_T res = RES_OK; + ASSERT(sck); + + res = logger_init(sck->allocator, &sck->logger__); + if(res != RES_OK) { + if(sck->verbose) { + fprintf(stderr, + MSG_ERROR_PREFIX + "Could not setup the AtrCK default logger -- %s.\n", + res_to_cstr(res)); + } + goto error; + } + logger_set_stream(&sck->logger__, LOG_OUTPUT, print_info, NULL); + logger_set_stream(&sck->logger__, LOG_ERROR, print_err, NULL); + logger_set_stream(&sck->logger__, LOG_WARNING, print_warn, NULL); + sck->logger = &sck->logger__; + +exit: + return res; +error: + goto exit; +} + +void +log_info(const struct sck* sck, const char* msg, ...) +{ + va_list vargs_list; + ASSERT(sck && msg); + + va_start(vargs_list, msg); + log_msg(sck, LOG_OUTPUT, msg, vargs_list); + va_end(vargs_list); +} + +void +log_err(const struct sck* sck, const char* msg, ...) +{ + va_list vargs_list; + ASSERT(sck && msg); + + va_start(vargs_list, msg); + log_msg(sck, LOG_ERROR, msg, vargs_list); + va_end(vargs_list); +} + +void +log_warn(const struct sck* sck, const char* msg, ...) +{ + va_list vargs_list; + ASSERT(sck && msg); + + va_start(vargs_list, msg); + log_msg(sck, LOG_WARNING, msg, vargs_list); + va_end(vargs_list); +} + diff --git a/src/sck_log.h b/src/sck_log.h @@ -0,0 +1,68 @@ +/* Copyright (C) 2020, 2021 CNRS + * + * 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/>. */ + +#ifndef SCK_LOG_H +#define SCK_LOG_H + +#include <rsys/rsys.h> + +#define MSG_INFO_PREFIX "StarCK:\x1b[1m\x1b[32minfo\x1b[0m: " +#define MSG_ERROR_PREFIX "StarCK:\x1b[1m\x1b[31merror\x1b[0m: " +#define MSG_WARNING_PREFIX "StarCK:\x1b[1m\x1b[33mwarning\x1b[0m: " + +struct sck; +struct logger; + +extern LOCAL_SYM res_T +setup_log_default + (struct sck* sck); + +/* Conditionally log a message on the LOG_OUTPUT stream of the sck logger, + * with respect to its verbose flag */ +extern LOCAL_SYM void +log_info + (const struct sck* sck, + const char* msg, + ...) +#ifdef COMPILER_GCC + __attribute((format(printf, 2, 3))) +#endif +; + +/* Conditionally log a message on the LOG_ERROR stream of the sck logger, + * with respect to its verbose flag */ +extern LOCAL_SYM void +log_err + (const struct sck* sck, + const char* msg, + ...) +#ifdef COMPILER_GCC + __attribute((format(printf, 2, 3))) +#endif +; + +/* Conditionally log a message on the LOG_WARNING stream of the sck logger, + * with respect to its verbose flag */ +extern LOCAL_SYM void +log_warn + (const struct sck* sck, + const char* msg, + ...) +#ifdef COMPILER_GCC + __attribute((format(printf, 2, 3))) +#endif +; + +#endif /* SCK_LOG_H */ diff --git a/src/test_atrck.c b/src/test_atrck.c @@ -1,70 +0,0 @@ -/* Copyright (C) 2020, 2021 CNRS - * - * 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/>. */ - -#include "atrck.h" - -#include <rsys/logger.h> - -static void -log_stream(const char* msg, void* ctx) -{ - ASSERT(msg); - (void)msg, (void)ctx; - printf("%s\n", msg); -} - -int -main(int argc, char** argv) -{ - struct mem_allocator allocator; - struct logger logger; - struct atrck* atrck; - (void)argc, (void)argv; - - CHK(atrck_create(NULL, NULL, 0, NULL) == RES_BAD_ARG); - CHK(atrck_create(NULL, NULL, 0, &atrck) == RES_OK); - - CHK(atrck_ref_get(NULL) == RES_BAD_ARG); - CHK(atrck_ref_get(atrck) == RES_OK); - CHK(atrck_ref_put(NULL) == RES_BAD_ARG); - CHK(atrck_ref_put(atrck) == RES_OK); - CHK(atrck_ref_put(atrck) == RES_OK); - - CHK(mem_init_proxy_allocator(&allocator, &mem_default_allocator) == RES_OK); - CHK(atrck_create(NULL, &allocator, 1, &atrck) == RES_OK); - CHK(atrck_ref_put(atrck) == RES_OK); - - CHK(logger_init(&allocator, &logger) == RES_OK); - logger_set_stream(&logger, LOG_OUTPUT, log_stream, NULL); - logger_set_stream(&logger, LOG_ERROR, log_stream, NULL); - logger_set_stream(&logger, LOG_WARNING, log_stream, NULL); - - CHK(atrck_create(&logger, &allocator, 0, &atrck) == RES_OK); - CHK(atrck_ref_put(atrck) == RES_OK); - CHK(atrck_create(&logger, NULL, 0, &atrck) == RES_OK); - CHK(atrck_ref_put(atrck) == RES_OK); - - logger_release(&logger); - if(MEM_ALLOCATED_SIZE(&allocator)) { - char dump[512]; - MEM_DUMP(&allocator, dump, sizeof(dump)/sizeof(char)); - fprintf(stderr, "%s\n", dump); - FATAL("Memory leaks\n"); - } - mem_shutdown_proxy_allocator(&allocator); - CHK(mem_allocated_size() == 0); - return 0; -} - diff --git a/src/test_atrck_load.c b/src/test_atrck_load.c @@ -1,405 +0,0 @@ -/* Copyright (C) 2020, 2021 CNRS - * - * 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/>. */ - -#include "atrck.h" - -#include <rsys/math.h> -#include <rsys/mem_allocator.h> -#include <math.h> - -static INLINE float -rand_canonic(void) -{ - int r; - while((r = rand()) == RAND_MAX); - return (float)r / (float)RAND_MAX; -} - -/******************************************************************************* - * Helper functions - ******************************************************************************/ -static void -check_atrck_load - (const struct atrck* atrck, - const size_t nbands, - const size_t nnodes) -{ - #define NMOMENTS 4 - struct atrck_band band = ATRCK_BAND_NULL; - struct atrck_quad_pt qpt = ATRCK_QUAD_PT_NULL; - size_t iband; - const size_t nsamps = 1000000; - size_t isamp; - - CHK(atrck); - CHK(nbands); - CHK(nnodes); - - CHK(atrck_get_bands_count(atrck) == nbands); - CHK(atrck_get_nodes_count(atrck) == nnodes); - - CHK(atrck_get_band(NULL, 0, &band) == RES_BAD_ARG); - CHK(atrck_get_band(atrck, nbands, &band) == RES_BAD_ARG); - CHK(atrck_get_band(atrck, nbands, NULL) == RES_BAD_ARG); - CHK(atrck_get_band(atrck, 0, &band) == RES_OK); - - CHK(band.quad_pts_count == 1); - CHK(atrck_band_get_quad_pt(NULL, 0, &qpt) == RES_BAD_ARG); - CHK(atrck_band_get_quad_pt(&band, 1, &qpt) == RES_BAD_ARG); - CHK(atrck_band_get_quad_pt(&band, 0, NULL) == RES_BAD_ARG); - - CHK(atrck_band_sample_quad_pt(NULL, 0, &qpt) == RES_BAD_ARG); - CHK(atrck_band_sample_quad_pt(&band, 1, &qpt) == RES_BAD_ARG); - CHK(atrck_band_sample_quad_pt(&band, 0, NULL) == RES_BAD_ARG); - - FOR_EACH(iband, 0, nbands) { - const double low = (double)iband; - const double upp = (double)(iband+1); - const size_t nqpts = iband + 1; - size_t iqpt; - double sum[NMOMENTS] = {0}; - double sum2[NMOMENTS] = {0}; - double E[NMOMENTS] = {0}; - double V[NMOMENTS] = {0}; - double SE[NMOMENTS] = {0}; - double ref; - int imoment; - - CHK(atrck_get_band(atrck, iband, &band) == RES_OK); - CHK(band.lower == low); - CHK(band.upper == upp); - CHK(band.id == iband); - CHK(band.quad_pts_count == nqpts); - - FOR_EACH(iqpt, 0, nqpts) { - const double abscissa = (double)iqpt*10.0; - const double weight = 1.0/(double)nqpts; - size_t inode; - - CHK(atrck_band_get_quad_pt(&band, iqpt, &qpt) == RES_OK); - CHK(qpt.abscissa == abscissa); - CHK(qpt.weight == weight); - CHK(qpt.id == iqpt); - - FOR_EACH(inode, 0, nnodes) { - const float ka = (float)(iband*1000 + iqpt*100 + inode); - CHK(qpt.ka_list[inode] == ka); - } - } - - /* Check the sampling of the quadrature points */ - FOR_EACH(isamp, 0, nsamps) { - const double r = rand_canonic(); - double w[NMOMENTS]; - - CHK(atrck_band_sample_quad_pt(&band, r, &qpt) == RES_OK); - FOR_EACH(imoment, 0, NMOMENTS) { - if(imoment == 0) { - w[imoment] = (double)(qpt.id + 1); - } else { - w[imoment] = (double)(qpt.id + 1) * w[imoment-1];; - } - sum[imoment] += w[imoment]; - sum2[imoment] += w[imoment]*w[imoment]; - } - } - - FOR_EACH(imoment, 0, NMOMENTS) { - E[imoment] = sum[imoment] / (double)nsamps; - V[imoment] = sum2[imoment] / (double)nsamps - E[imoment]*E[imoment]; - SE[imoment] = sqrt(V[imoment]/(double)nsamps); - } - - ref = (double)(nqpts+1)/2.0; - CHK(eq_eps(ref, E[0], 3*SE[0])); - - ref = (double)(2*nqpts*nqpts + 3*nqpts + 1) / 6.0; - CHK(eq_eps(ref, E[1], 3*SE[1])); - - ref = ((double)nqpts * pow((double)nqpts + 1.0, 2.0)) / 4.0; - CHK(eq_eps(ref, E[2], 3*SE[2])); - - ref = (double)1.0/30.0 * (double) - (6*nqpts*nqpts*nqpts*nqpts + 15*nqpts*nqpts*nqpts + 10*nqpts*nqpts - 1); - CHK(eq_eps(ref, E[3], 3*SE[3])); - } -} - -static void -test_load(struct atrck* atrck) -{ - FILE* fp = NULL; - const char* filename = "test_file.atrck"; - const uint64_t pagesize = 16384; - const uint64_t nbands = 11; - const uint64_t nnodes = 1000; - uint64_t iband; - const char byte = 0; - - fp = fopen(filename, "w+"); - CHK(fp); - - /* Write the header */ - CHK(fwrite(&pagesize, sizeof(pagesize), 1, fp) == 1); - CHK(fwrite(&nbands, sizeof(nbands), 1, fp) == 1); - CHK(fwrite(&nnodes, sizeof(nnodes), 1, fp) == 1); - - FOR_EACH(iband, 0, nbands) { - const double low = (double)iband; - const double upp = (double)(iband+1); - const uint64_t nqpts = iband + 1; - uint64_t iqpt; - - /* Write band description */ - CHK(fwrite(&iband, sizeof(iband), 1, fp) == 1); - CHK(fwrite(&low, sizeof(low), 1, fp) == 1); - CHK(fwrite(&upp, sizeof(upp), 1, fp) == 1); - CHK(fwrite(&nqpts, sizeof(nqpts), 1, fp) == 1); - - /* Write per band quadrature points */ - FOR_EACH(iqpt, 0, nqpts) { - const double abscissa = (double)iqpt*10.0; - const double weight = 1.0/(double)nqpts; - CHK(fwrite(&abscissa, sizeof(abscissa), 1, fp) == 1); - CHK(fwrite(&weight, sizeof(weight), 1, fp) == 1); - } - } - - /* Write per band and per quadrature point list of ka */ - FOR_EACH(iband, 0, nbands) { - const uint64_t nqpts = iband + 1; - uint64_t iqpt; - FOR_EACH(iqpt, 0, nqpts) { - uint64_t inode; - - /* Padding */ - CHK(fseek(fp, (long)ALIGN_SIZE((size_t)ftell(fp), pagesize), SEEK_SET)==0); - - FOR_EACH(inode, 0, nnodes) { - const float ka = (float)(iband*1000 + iqpt*100 + inode); - CHK(fwrite(&ka, sizeof(ka), 1, fp) == 1); - } - } - } - - /* Padding. Write one char to position the EOF indicator */ - CHK(fseek(fp, (long)ALIGN_SIZE((size_t)ftell(fp), pagesize)-1, SEEK_SET) == 0); - CHK(fwrite(&byte, sizeof(byte), 1, fp) == 1); - rewind(fp); - - CHK(atrck_load_stream(NULL, fp, filename) == RES_BAD_ARG); - CHK(atrck_load_stream(atrck, NULL, filename) == RES_BAD_ARG); - CHK(atrck_load_stream(atrck, fp, NULL) == RES_OK); - check_atrck_load(atrck, nbands, nnodes); - - rewind(fp); - CHK(atrck_load_stream(atrck, fp, filename) == RES_OK); - - CHK(atrck_load(NULL, filename) == RES_BAD_ARG); - CHK(atrck_load(atrck, NULL) == RES_BAD_ARG); - CHK(atrck_load(atrck, "nop") == RES_IO_ERR); - CHK(atrck_load(atrck, filename) == RES_OK); - - fclose(fp); -} - -static void -test_load_fail(struct atrck* atrck) -{ - const char byte = 1; - FILE* fp = NULL; - uint64_t pagesize; - uint64_t nnodes; - uint64_t nbands; - uint64_t iband; - uint64_t nqpts; - double low; - double upp; - double abscissa; - double weight; - double ka; - - /* Wrong pagesize */ - fp = tmpfile(); - CHK(fp); - pagesize = 2048; - nnodes = 1; - nbands = 1; - CHK(fwrite(&pagesize, sizeof(pagesize), 1, fp) == 1); - CHK(fwrite(&nbands, sizeof(nbands), 1, fp) == 1); - CHK(fwrite(&nnodes, sizeof(nnodes), 1, fp) == 1); - iband = 0; - low = 0; - upp = 1; - nqpts = 1; - CHK(fwrite(&iband, sizeof(iband), 1, fp) == 1); - CHK(fwrite(&low, sizeof(low), 1, fp) == 1); - CHK(fwrite(&upp, sizeof(upp), 1, fp) == 1); - CHK(fwrite(&nqpts, sizeof(nqpts), 1, fp) == 1); - abscissa = 1; - weight = 1; - CHK(fwrite(&abscissa, sizeof(abscissa), 1, fp) == 1); - CHK(fwrite(&weight, sizeof(weight), 1, fp) == 1); - CHK(fseek(fp, (long)ALIGN_SIZE((size_t)ftell(fp), pagesize), SEEK_SET) == 0); - ka = 1; - CHK(fwrite(&ka, sizeof(ka), 1, fp) == 1); - CHK(fseek(fp, (long)ALIGN_SIZE((size_t)ftell(fp), pagesize)-1, SEEK_SET) == 0); - CHK(fwrite(&byte, sizeof(byte), 1, fp) == 1); /* Positioned the EOF */ - rewind(fp); - - CHK(atrck_load_stream(atrck, fp, NULL) == RES_BAD_ARG); - fclose(fp); - - /* Wrong #bands */ - fp = tmpfile(); - CHK(fp); - pagesize = 4096; - nnodes = 1; - nbands = 0; - CHK(fwrite(&pagesize, sizeof(pagesize), 1, fp) == 1); - CHK(fwrite(&nbands, sizeof(nbands), 1, fp) == 1); - CHK(fwrite(&nnodes, sizeof(nnodes), 1, fp) == 1); - rewind(fp); - - CHK(atrck_load_stream(atrck, fp, NULL) == RES_BAD_ARG); - fclose(fp); - - /* Wrong #nodes */ - fp = tmpfile(); - CHK(fp); - pagesize = 4096; - nnodes = 0; - nbands = 1; - CHK(fwrite(&pagesize, sizeof(pagesize), 1, fp) == 1); - CHK(fwrite(&nbands, sizeof(nbands), 1, fp) == 1); - CHK(fwrite(&nnodes, sizeof(nnodes), 1, fp) == 1); - iband = 0; - low = 0; - upp = 1; - nqpts = 1; - CHK(fwrite(&iband, sizeof(iband), 1, fp) == 1); - CHK(fwrite(&low, sizeof(low), 1, fp) == 1); - CHK(fwrite(&upp, sizeof(upp), 1, fp) == 1); - CHK(fwrite(&nqpts, sizeof(nqpts), 1, fp) == 1); - abscissa = 1; - weight = 1; - CHK(fwrite(&abscissa, sizeof(abscissa), 1, fp) == 1); - CHK(fwrite(&weight, sizeof(weight), 1, fp) == 1); - CHK(fseek(fp, (long)ALIGN_SIZE((size_t)ftell(fp), pagesize)-1, SEEK_SET) == 0); - CHK(fwrite(&byte, sizeof(byte), 1, fp) == 1); /* Positioned the EOF */ - rewind(fp); - - /* Wrong band boundaries */ - fp = tmpfile(); - CHK(fp); - pagesize = 16384; - nnodes = 1; - nbands = 1; - CHK(fwrite(&pagesize, sizeof(pagesize), 1, fp) == 1); - CHK(fwrite(&nbands, sizeof(nbands), 1, fp) == 1); - CHK(fwrite(&nnodes, sizeof(nnodes), 1, fp) == 1); - iband = 0; - low = 1; - upp = 0; - nqpts = 1; - CHK(fwrite(&iband, sizeof(iband), 1, fp) == 1); - CHK(fwrite(&low, sizeof(low), 1, fp) == 1); - CHK(fwrite(&upp, sizeof(upp), 1, fp) == 1); - CHK(fwrite(&nqpts, sizeof(nqpts), 1, fp) == 1); - abscissa = 1; - weight = 1; - CHK(fwrite(&abscissa, sizeof(abscissa), 1, fp) == 1); - CHK(fwrite(&weight, sizeof(weight), 1, fp) == 1); - CHK(fseek(fp, (long)ALIGN_SIZE((size_t)ftell(fp), pagesize), SEEK_SET) == 0); - ka = 1; - CHK(fwrite(&ka, sizeof(ka), 1, fp) == 1); - CHK(fseek(fp, (long)ALIGN_SIZE((size_t)ftell(fp), pagesize)-1, SEEK_SET) == 0); - CHK(fwrite(&byte, sizeof(byte), 1, fp) == 1); /* Positioned the EOF */ - rewind(fp); - - - CHK(atrck_load_stream(atrck, fp, NULL) == RES_BAD_ARG); - fclose(fp); -} - -static void -test_load_files(struct atrck* atrck, int argc, char** argv) -{ - int i; - CHK(atrck); - FOR_EACH(i, 1, argc) { - size_t nnodes; - size_t nbands; - size_t iband; - - printf("Load %s\n", argv[1]); - - CHK(atrck_load(atrck, argv[i]) == RES_OK); - nbands = atrck_get_bands_count(atrck); - nnodes = atrck_get_nodes_count(atrck); - CHK(nbands); - CHK(nnodes); - - FOR_EACH(iband, 0, nbands) { - struct atrck_band band = ATRCK_BAND_NULL; - size_t iqpt; - - CHK(atrck_get_band(atrck, iband, &band) == RES_OK); - printf("band %lu in [%g, %g] cm^-1 with %lu quadrature points\n", - (unsigned long)band.id, - band.lower, band.upper, - (unsigned long)band.quad_pts_count); - - CHK(band.lower == band.lower); /* !NaN */ - CHK(band.upper == band.upper); /* !NaN */ - CHK(band.lower < band.upper); - CHK(band.quad_pts_count); - - FOR_EACH(iqpt, 0, band.quad_pts_count) { - struct atrck_quad_pt qpt; - size_t inode; - - CHK(atrck_band_get_quad_pt(&band, iqpt, &qpt) == RES_OK); - CHK(qpt.abscissa == qpt.abscissa); /* !NaN */ - CHK(qpt.weight == qpt.weight); /* !NaN */ - CHK(qpt.weight > 0); - - FOR_EACH(inode, 0, nnodes) { - CHK(qpt.ka_list[inode] == qpt.ka_list[inode]); /* !NaN */ - } - } - } - } -} - -/******************************************************************************* - * Main function - ******************************************************************************/ -int -main(int argc, char** argv) -{ - struct atrck* atrck = NULL; - (void)argc, (void)argv; - - CHK(atrck_create(NULL, &mem_default_allocator, 1, &atrck) == RES_OK); - - test_load(atrck); - test_load_fail(atrck); - test_load_files(atrck, argc, argv); - - CHK(atrck_ref_put(atrck) == RES_OK); - CHK(mem_allocated_size() == 0); - return 0; -} diff --git a/src/test_sck.c b/src/test_sck.c @@ -0,0 +1,71 @@ +/* Copyright (C) 2020, 2021 CNRS + * + * 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/>. */ + +#include "sck.h" + +#include <rsys/logger.h> + +static void +log_stream(const char* msg, void* ctx) +{ + ASSERT(msg); + (void)msg, (void)ctx; + printf("%s\n", msg); +} + +int +main(int argc, char** argv) +{ + struct mem_allocator allocator; + struct logger logger; + struct sck_create_args args = SCK_CREATE_ARGS_DEFAULT; + struct sck* sck; + (void)argc, (void)argv; + + CHK(sck_create(NULL, &sck) == RES_BAD_ARG); + CHK(sck_create(&args, NULL) == RES_BAD_ARG); + CHK(sck_create(&args, &sck) == RES_OK); + + CHK(sck_ref_get(NULL) == RES_BAD_ARG); + CHK(sck_ref_get(sck) == RES_OK); + CHK(sck_ref_put(NULL) == RES_BAD_ARG); + CHK(sck_ref_put(sck) == RES_OK); + CHK(sck_ref_put(sck) == RES_OK); + + CHK(mem_init_proxy_allocator(&allocator, &mem_default_allocator) == RES_OK); + args.allocator = &allocator; + args.verbose = 1; + CHK(sck_create(&args, &sck) == RES_OK); + CHK(sck_ref_put(sck) == RES_OK); + + CHK(logger_init(&allocator, &logger) == RES_OK); + logger_set_stream(&logger, LOG_OUTPUT, log_stream, NULL); + logger_set_stream(&logger, LOG_ERROR, log_stream, NULL); + logger_set_stream(&logger, LOG_WARNING, log_stream, NULL); + + args.logger = &logger; + args.verbose = 0; + CHK(sck_create(&args, &sck) == RES_OK); + CHK(sck_ref_put(sck) == RES_OK); + args.allocator = NULL; + CHK(sck_create(&args, &sck) == RES_OK); + CHK(sck_ref_put(sck) == RES_OK); + + logger_release(&logger); + mem_shutdown_proxy_allocator(&allocator); + CHK(mem_allocated_size() == 0); + return 0; +} + diff --git a/src/test_sck_load.c b/src/test_sck_load.c @@ -0,0 +1,411 @@ +/* Copyright (C) 2020, 2021 CNRS + * + * 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/>. */ + +#include "sck.h" + +#include <rsys/math.h> +#include <rsys/mem_allocator.h> +#include <math.h> + +static INLINE float +rand_canonic(void) +{ + int r; + while((r = rand()) == RAND_MAX); + return (float)r / (float)RAND_MAX; +} + +/******************************************************************************* + * Helper functions + ******************************************************************************/ +static void +check_sck_load + (const struct sck* sck, + const size_t nbands, + const size_t nnodes) +{ + #define NMOMENTS 4 + struct sck_band band = SCK_BAND_NULL; + struct sck_quad_pt qpt = SCK_QUAD_PT_NULL; + size_t iband; + const size_t nsamps = 1000000; + size_t isamp; + + CHK(sck); + CHK(nbands); + CHK(nnodes); + + CHK(sck_get_bands_count(sck) == nbands); + CHK(sck_get_nodes_count(sck) == nnodes); + + CHK(sck_get_band(NULL, 0, &band) == RES_BAD_ARG); + CHK(sck_get_band(sck, nbands, &band) == RES_BAD_ARG); + CHK(sck_get_band(sck, nbands, NULL) == RES_BAD_ARG); + CHK(sck_get_band(sck, 0, &band) == RES_OK); + + CHK(band.quad_pts_count == 1); + CHK(sck_band_get_quad_pt(NULL, 0, &qpt) == RES_BAD_ARG); + CHK(sck_band_get_quad_pt(&band, 1, &qpt) == RES_BAD_ARG); + CHK(sck_band_get_quad_pt(&band, 0, NULL) == RES_BAD_ARG); + + CHK(sck_band_sample_quad_pt(NULL, 0, &qpt) == RES_BAD_ARG); + CHK(sck_band_sample_quad_pt(&band, 1, &qpt) == RES_BAD_ARG); + CHK(sck_band_sample_quad_pt(&band, 0, NULL) == RES_BAD_ARG); + + FOR_EACH(iband, 0, nbands) { + const double low = (double)iband; + const double upp = (double)(iband+1); + const size_t nqpts = iband + 1; + size_t iqpt; + double sum[NMOMENTS] = {0}; + double sum2[NMOMENTS] = {0}; + double E[NMOMENTS] = {0}; + double V[NMOMENTS] = {0}; + double SE[NMOMENTS] = {0}; + double ref; + size_t inode; + int imoment; + + CHK(sck_get_band(sck, iband, &band) == RES_OK); + CHK(band.lower == low); + CHK(band.upper == upp); + CHK(band.id == iband); + CHK(band.quad_pts_count == nqpts); + + FOR_EACH(inode, 0, nnodes) { + const float ks = (float)(iband*2000 + inode); + CHK(band.ks_list[inode] == ks); + } + + FOR_EACH(iqpt, 0, nqpts) { + const double weight = 1.0/(double)nqpts; + + CHK(sck_band_get_quad_pt(&band, iqpt, &qpt) == RES_OK); + CHK(qpt.weight == weight); + CHK(qpt.id == iqpt); + + FOR_EACH(inode, 0, nnodes) { + const float ka = (float)(iband*1000 + iqpt*100 + inode); + CHK(qpt.ka_list[inode] == ka); + } + } + + /* Check the sampling of the quadrature points */ + FOR_EACH(isamp, 0, nsamps) { + const double r = rand_canonic(); + double w[NMOMENTS]; + + CHK(sck_band_sample_quad_pt(&band, r, &qpt) == RES_OK); + FOR_EACH(imoment, 0, NMOMENTS) { + if(imoment == 0) { + w[imoment] = (double)(qpt.id + 1); + } else { + w[imoment] = (double)(qpt.id + 1) * w[imoment-1];; + } + sum[imoment] += w[imoment]; + sum2[imoment] += w[imoment]*w[imoment]; + } + } + + FOR_EACH(imoment, 0, NMOMENTS) { + E[imoment] = sum[imoment] / (double)nsamps; + V[imoment] = sum2[imoment] / (double)nsamps - E[imoment]*E[imoment]; + SE[imoment] = sqrt(V[imoment]/(double)nsamps); + } + + ref = (double)(nqpts+1)/2.0; + CHK(eq_eps(ref, E[0], 3*SE[0])); + + ref = (double)(2*nqpts*nqpts + 3*nqpts + 1) / 6.0; + CHK(eq_eps(ref, E[1], 3*SE[1])); + + ref = ((double)nqpts * pow((double)nqpts + 1.0, 2.0)) / 4.0; + CHK(eq_eps(ref, E[2], 3*SE[2])); + + ref = (double)1.0/30.0 * (double) + (6*nqpts*nqpts*nqpts*nqpts + 15*nqpts*nqpts*nqpts + 10*nqpts*nqpts - 1); + CHK(eq_eps(ref, E[3], 3*SE[3])); + } +} + +static void +test_load(struct sck* sck) +{ + FILE* fp = NULL; + const char* filename = "test_file.sck"; + const uint64_t pagesize = 16384; + const uint64_t nbands = 11; + const uint64_t nnodes = 1000; + uint64_t iband; + const char byte = 0; + + fp = fopen(filename, "w+"); + CHK(fp); + + /* Write the header */ + CHK(fwrite(&pagesize, sizeof(pagesize), 1, fp) == 1); + CHK(fwrite(&nbands, sizeof(nbands), 1, fp) == 1); + CHK(fwrite(&nnodes, sizeof(nnodes), 1, fp) == 1); + + FOR_EACH(iband, 0, nbands) { + const double low = (double)iband; + const double upp = (double)(iband+1); + const uint64_t nqpts = iband + 1; + uint64_t iqpt; + + /* Write band description */ + CHK(fwrite(&low, sizeof(low), 1, fp) == 1); + CHK(fwrite(&upp, sizeof(upp), 1, fp) == 1); + CHK(fwrite(&nqpts, sizeof(nqpts), 1, fp) == 1); + + /* Write per band quadrature points */ + FOR_EACH(iqpt, 0, nqpts) { + const double weight = 1.0/(double)nqpts; + CHK(fwrite(&weight, sizeof(weight), 1, fp) == 1); + } + } + + /* Write per band ks and per quadrature point ka */ + FOR_EACH(iband, 0, nbands) { + const uint64_t nqpts = iband + 1; + uint64_t iqpt; + uint64_t inode; + + /* Padding */ + CHK(fseek(fp, (long)ALIGN_SIZE((size_t)ftell(fp), pagesize), SEEK_SET)==0); + + FOR_EACH(inode, 0, nnodes) { + const float ks = (float)(iband*2000 + inode); + CHK(fwrite(&ks, sizeof(ks), 1, fp) == 1); + } + + FOR_EACH(iqpt, 0, nqpts) { + + /* Padding */ + CHK(fseek(fp, (long)ALIGN_SIZE((size_t)ftell(fp), pagesize), SEEK_SET)==0); + + FOR_EACH(inode, 0, nnodes) { + const float ka = (float)(iband*1000 + iqpt*100 + inode); + CHK(fwrite(&ka, sizeof(ka), 1, fp) == 1); + } + } + } + + /* Padding. Write one char to position the EOF indicator */ + CHK(fseek(fp, (long)ALIGN_SIZE((size_t)ftell(fp), pagesize)-1, SEEK_SET) == 0); + CHK(fwrite(&byte, sizeof(byte), 1, fp) == 1); + rewind(fp); + + CHK(sck_load_stream(NULL, fp, filename) == RES_BAD_ARG); + CHK(sck_load_stream(sck, NULL, filename) == RES_BAD_ARG); + CHK(sck_load_stream(sck, fp, NULL) == RES_OK); + check_sck_load(sck, nbands, nnodes); + + rewind(fp); + CHK(sck_load_stream(sck, fp, filename) == RES_OK); + + CHK(sck_load(NULL, filename) == RES_BAD_ARG); + CHK(sck_load(sck, NULL) == RES_BAD_ARG); + CHK(sck_load(sck, "nop") == RES_IO_ERR); + CHK(sck_load(sck, filename) == RES_OK); + + fclose(fp); +} + +static void +test_load_fail(struct sck* sck) +{ + const char byte = 1; + FILE* fp = NULL; + uint64_t pagesize; + uint64_t nnodes; + uint64_t nbands; + uint64_t nqpts; + double low; + double upp; + double weight; + float ks; + float ka; + + /* Wrong pagesize */ + fp = tmpfile(); + CHK(fp); + pagesize = 2048; + nnodes = 1; + nbands = 1; + CHK(fwrite(&pagesize, sizeof(pagesize), 1, fp) == 1); + CHK(fwrite(&nbands, sizeof(nbands), 1, fp) == 1); + CHK(fwrite(&nnodes, sizeof(nnodes), 1, fp) == 1); + low = 0; + upp = 1; + nqpts = 1; + CHK(fwrite(&low, sizeof(low), 1, fp) == 1); + CHK(fwrite(&upp, sizeof(upp), 1, fp) == 1); + CHK(fwrite(&nqpts, sizeof(nqpts), 1, fp) == 1); + weight = 1; + CHK(fwrite(&weight, sizeof(weight), 1, fp) == 1); + CHK(fseek(fp, (long)ALIGN_SIZE((size_t)ftell(fp), pagesize), SEEK_SET) == 0); + ks = 1; + CHK(fwrite(&ks, sizeof(ka), 1, fp) == 1); + CHK(fseek(fp, (long)ALIGN_SIZE((size_t)ftell(fp), pagesize), SEEK_SET) == 0); + ka = 1; + CHK(fwrite(&ka, sizeof(ka), 1, fp) == 1); + CHK(fseek(fp, (long)ALIGN_SIZE((size_t)ftell(fp), pagesize)-1, SEEK_SET) == 0); + CHK(fwrite(&byte, sizeof(byte), 1, fp) == 1); /* Positioned the EOF */ + rewind(fp); + + CHK(sck_load_stream(sck, fp, NULL) == RES_BAD_ARG); + fclose(fp); + + /* Wrong #bands */ + fp = tmpfile(); + CHK(fp); + pagesize = 4096; + nnodes = 1; + nbands = 0; + CHK(fwrite(&pagesize, sizeof(pagesize), 1, fp) == 1); + CHK(fwrite(&nbands, sizeof(nbands), 1, fp) == 1); + CHK(fwrite(&nnodes, sizeof(nnodes), 1, fp) == 1); + rewind(fp); + + CHK(sck_load_stream(sck, fp, NULL) == RES_BAD_ARG); + fclose(fp); + + /* Wrong #nodes */ + fp = tmpfile(); + CHK(fp); + pagesize = 4096; + nnodes = 0; + nbands = 1; + CHK(fwrite(&pagesize, sizeof(pagesize), 1, fp) == 1); + CHK(fwrite(&nbands, sizeof(nbands), 1, fp) == 1); + CHK(fwrite(&nnodes, sizeof(nnodes), 1, fp) == 1); + low = 0; + upp = 1; + nqpts = 1; + CHK(fwrite(&low, sizeof(low), 1, fp) == 1); + CHK(fwrite(&upp, sizeof(upp), 1, fp) == 1); + CHK(fwrite(&nqpts, sizeof(nqpts), 1, fp) == 1); + weight = 1; + CHK(fwrite(&weight, sizeof(weight), 1, fp) == 1); + CHK(fseek(fp, (long)ALIGN_SIZE((size_t)ftell(fp), pagesize)-1, SEEK_SET) == 0); + CHK(fwrite(&byte, sizeof(byte), 1, fp) == 1); /* Positioned the EOF */ + fclose(fp); + + /* Wrong band boundaries */ + fp = tmpfile(); + CHK(fp); + pagesize = 16384; + nnodes = 1; + nbands = 1; + CHK(fwrite(&pagesize, sizeof(pagesize), 1, fp) == 1); + CHK(fwrite(&nbands, sizeof(nbands), 1, fp) == 1); + CHK(fwrite(&nnodes, sizeof(nnodes), 1, fp) == 1); + low = 1; + upp = 0; + nqpts = 1; + CHK(fwrite(&low, sizeof(low), 1, fp) == 1); + CHK(fwrite(&upp, sizeof(upp), 1, fp) == 1); + CHK(fwrite(&nqpts, sizeof(nqpts), 1, fp) == 1); + weight = 1; + CHK(fwrite(&weight, sizeof(weight), 1, fp) == 1); + CHK(fseek(fp, (long)ALIGN_SIZE((size_t)ftell(fp), pagesize), SEEK_SET) == 0); + ks = 1; + CHK(fwrite(&ka, sizeof(ks), 1, fp) == 1); + CHK(fseek(fp, (long)ALIGN_SIZE((size_t)ftell(fp), pagesize), SEEK_SET) == 0); + ka = 1; + CHK(fwrite(&ka, sizeof(ka), 1, fp) == 1); + CHK(fseek(fp, (long)ALIGN_SIZE((size_t)ftell(fp), pagesize)-1, SEEK_SET) == 0); + CHK(fwrite(&byte, sizeof(byte), 1, fp) == 1); /* Positioned the EOF */ + rewind(fp); + + CHK(sck_load_stream(sck, fp, NULL) == RES_BAD_ARG); + fclose(fp); +} + +static void +test_load_files(struct sck* sck, int argc, char** argv) +{ + int i; + CHK(sck); + FOR_EACH(i, 1, argc) { + size_t nnodes; + size_t nbands; + size_t iband; + + printf("Load %s\n", argv[1]); + + CHK(sck_load(sck, argv[i]) == RES_OK); + nbands = sck_get_bands_count(sck); + nnodes = sck_get_nodes_count(sck); + CHK(nbands); + CHK(nnodes); + + FOR_EACH(iband, 0, nbands) { + struct sck_band band = SCK_BAND_NULL; + size_t inode; + size_t iqpt; + + CHK(sck_get_band(sck, iband, &band) == RES_OK); + printf("band %lu in [%g, %g] nm with %lu quadrature points\n", + (unsigned long)band.id, + band.lower, band.upper, + (unsigned long)band.quad_pts_count); + + CHK(band.lower == band.lower); /* !NaN */ + CHK(band.upper == band.upper); /* !NaN */ + CHK(band.lower < band.upper); + CHK(band.quad_pts_count); + + FOR_EACH(inode, 0, nnodes) { + CHK(band.ks_list[inode] == band.ks_list[inode]); /* !NaN */ + } + + FOR_EACH(iqpt, 0, band.quad_pts_count) { + struct sck_quad_pt qpt; + + CHK(sck_band_get_quad_pt(&band, iqpt, &qpt) == RES_OK); + CHK(qpt.weight == qpt.weight); /* !NaN */ + CHK(qpt.weight > 0); + + FOR_EACH(inode, 0, nnodes) { + CHK(qpt.ka_list[inode] == qpt.ka_list[inode]); /* !NaN */ + } + } + } + } +} + +/******************************************************************************* + * Main function + ******************************************************************************/ +int +main(int argc, char** argv) +{ + struct sck* sck = NULL; + struct sck_create_args args = SCK_CREATE_ARGS_DEFAULT; + (void)argc, (void)argv; + + args.verbose = 1; + CHK(sck_create(&args, &sck) == RES_OK); + + test_load(sck); + test_load_fail(sck); + test_load_files(sck, argc, argv); + + CHK(sck_ref_put(sck) == RES_OK); + CHK(mem_allocated_size() == 0); + return 0; +}