diff --git a/.github/workflows/tests.yml b/.github/workflows/tests.yml new file mode 100644 index 0000000..fa51d6f --- /dev/null +++ b/.github/workflows/tests.yml @@ -0,0 +1,29 @@ +name: Tests + +on: + push: + pull_request: + workflow_dispatch: + +jobs: + cmake: + name: Build & Test + runs-on: ${{ matrix.os }} + strategy: + fail-fast: false + matrix: + os: [ubuntu-latest, windows-latest] + + steps: + - uses: actions/checkout@v7 + with: + submodules: recursive + + - name: Configure + run: cmake -S . -B build -DBUILD_TESTING=ON -DCMAKE_BUILD_TYPE=Release + + - name: Build + run: cmake --build build --config Release + + - name: Test + run: ctest --test-dir build --build-config Release --output-on-failure diff --git a/.gitignore b/.gitignore new file mode 100644 index 0000000..f9c5c22 --- /dev/null +++ b/.gitignore @@ -0,0 +1,2 @@ +# default build directory +build/ diff --git a/.gitmodules b/.gitmodules new file mode 100644 index 0000000..800345d --- /dev/null +++ b/.gitmodules @@ -0,0 +1,3 @@ +[submodule "external/Catch2"] + path = external/Catch2 + url = https://github.com/catchorg/Catch2.git diff --git a/CMakeLists.txt b/CMakeLists.txt new file mode 100644 index 0000000..6af4887 --- /dev/null +++ b/CMakeLists.txt @@ -0,0 +1,40 @@ +cmake_minimum_required(VERSION 3.16) + +project(NDTable LANGUAGES C CXX) + +function(enable_project_warnings target) + if(MSVC) + target_compile_options(${target} PRIVATE /W4) + else() + target_compile_options(${target} PRIVATE -Wall -Wextra -Wpedantic) + endif() +endfunction() + +add_library(NDTable + src/Core.c + src/Interpolation.c +) + +enable_project_warnings(NDTable) + +target_include_directories(NDTable + PUBLIC + ${CMAKE_CURRENT_SOURCE_DIR}/include +) + +include(CTest) + +if(BUILD_TESTING) + add_subdirectory(external/Catch2) + + add_executable(NDTable_test + test/CoreTests.cpp + test/InterpolationTests.cpp + ) + + target_compile_features(NDTable_test PRIVATE cxx_std_14) + target_link_libraries(NDTable_test PRIVATE NDTable Catch2::Catch2WithMain) + enable_project_warnings(NDTable_test) + + add_test(NAME NDTable_test COMMAND NDTable_test) +endif() diff --git a/LICENSE b/LICENSE index d895fe4..a6a0200 100644 --- a/LICENSE +++ b/LICENSE @@ -1,21 +1,16 @@ -BSD 3-Clause License +BSD 2-Clause License -Copyright (c) 2017, Scientific Data Format -All rights reserved. +Copyright (c) 2026, Dassault Systemes Redistribution and use in source and binary forms, with or without modification, are permitted provided that the following conditions are met: -* Redistributions of source code must retain the above copyright notice, this - list of conditions and the following disclaimer. +1. Redistributions of source code must retain the above copyright notice, this + list of conditions and the following disclaimer. -* Redistributions in binary form must reproduce the above copyright notice, - this list of conditions and the following disclaimer in the documentation - and/or other materials provided with the distribution. - -* Neither the name of the copyright holder nor the names of its - contributors may be used to endorse or promote products derived from - this software without specific prior written permission. +2. Redistributions in binary form must reproduce the above copyright notice, + this list of conditions and the following disclaimer in the documentation + and/or other materials provided with the distribution. THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE diff --git a/README.md b/README.md index f307203..6204873 100644 --- a/README.md +++ b/README.md @@ -1,2 +1,13 @@ +![Build Status](https://ci.appveyor.com/api/projects/status/github/ScientificDataFormat/NDTable?branch=master&svg=true) + # NDTable -A C library to inter- and extrapolate multi-dimensional data + +NDTable is a C library to inter- and extrapolate multi-dimensional data. + +## License + +The code in this repository is licensed under the [3-clause BSD license](LICENSE) + +---------------------------------------------- + +Copyright © 2017 Dassault Systèmes diff --git a/external/Catch2 b/external/Catch2 new file mode 160000 index 0000000..317ac1e --- /dev/null +++ b/external/Catch2 @@ -0,0 +1 @@ +Subproject commit 317ac1ed4c0bb6e6b91eafc817e05c488feffcb3 diff --git a/include/NDTable.h b/include/NDTable.h new file mode 100644 index 0000000..b0cd4b6 --- /dev/null +++ b/include/NDTable.h @@ -0,0 +1,156 @@ +#ifndef NDTABLE_H_ +#define NDTABLE_H_ + +#ifdef __cplusplus +extern "C" { +#endif + +/*! The maximum number of dimensions */ +#define MAX_NDIMS 32 + +/*! Interpolation methods */ +typedef enum { + NDTABLE_INTERP_HOLD = 1, + NDTABLE_INTERP_NEAREST, + NDTABLE_INTERP_LINEAR, + NDTABLE_INTERP_AKIMA, + NDTABLE_INTERP_FRITSCH_BUTLAND, + NDTABLE_INTERP_STEFFEN +} NDTable_InterpMethod_t; + +/*! Extrapolation methods */ +typedef enum { + NDTABLE_EXTRAP_HOLD = 1, + NDTABLE_EXTRAP_LINEAR, + NDTABLE_EXTRAP_NONE +} NDTable_ExtrapMethod_t; + +/*! The structure that holds the data values */ +typedef struct { + int ndims; //!< the number of dimensions of the table + int dims[MAX_NDIMS]; //!< extents of the dimensions + int numel; //!< the number of data values + int offs[MAX_NDIMS]; //!< the index offsets for the dimensions + double *data; //!< the data values + double *scales[MAX_NDIMS]; //!< array of pointers to the scale values +} NDTable_t; + +typedef NDTable_t * NDTable_h; + +/*! Interpolation status codes */ +typedef enum { + NDTABLE_INTERPSTATUS_UNKNOWN_METHOD = -4, + NDTABLE_INTERPSTATUS_DATASETNOTFOUND = -3, + NDTABLE_INTERPSTATUS_WRONGNPARAMS = -2, + NDTABLE_INTERPSTATUS_OUTOFBOUNS = -1, + NDTABLE_INTERPSTATUS_OK = 0 +} NDTable_InterpolationStatus; + + +/*! Get the last error message + * + * @return the error message + */ +const char * NDTable_get_error_message(); + +/*! Evaluate the value of the table at the given sample point using the specified inter- and extrapolation methods + * + * @param [in] table the table handle + * @param [in] nparams the number of dimensions + * @param [in] params the sample point + * @param [in] interp_method the interpolation method + * @param [in] extrap_method the extrapolation method + * @param [out] value the value at the sample point + * + * @return 0 if the value could be evaluated, -1 otherwise + */ +int NDTable_evaluate(NDTable_h table, int nparams, const double params[], NDTable_InterpMethod_t interp_method, NDTable_ExtrapMethod_t extrap_method, double *value); + +/*! Evalute the total differential of the table at the given sample point and deltas using the specified inter- and extrapolation methods + * + * @param [in] table the table handle + * @param [in] nparams the number of dimensions + * @param [in] params the sample point + * @param [in] delta_params the the deltas + * @param [in] interp_method the interpolation method + * @param [in] extrap_method the extrapolation method + * @param [out] value the total differential at the sample point + * + * @return 0 if the value could be evaluated, -1 otherwise + */ +int NDTable_evaluate_derivative(NDTable_h table, int nparams, const double params[], const double delta_params[], NDTable_InterpMethod_t interp_method, NDTable_ExtrapMethod_t extrap_method, double *value); + +/*! The maximum length of an error message */ +#define MAX_MESSAGE_LENGTH 256 + +/*! Sets the error message */ +void NDTable_set_error_message(const char *msg, ...); + +/*! Allocates a new table + * + * @return a pointer to the new table + */ +NDTable_h NDTable_alloc_table(); + +/*! De-allocates a table + * + * @param [in] pointer to the table to de-allocate + */ +void NDTable_free_table(NDTable_h table); + +/*! Converts index to subscripts + * + * @param [in] index the index to convert + * @param [in] table the table for which to convert the index + * @param [out] subs the subscripts + */ +void NDTable_ind2sub(const int index, const NDTable_h table, int *subs); + +/*! Converts subscripts to index + * + * @param [in] subs the subscripts to convert + * @param [in] table the table for which to convert the subscripts + * @param [out] index the index + */ +void NDTable_sub2ind(const int *subs, const NDTable_h table, int *index); + +double NDTable_get_value_subs(const NDTable_h table, const int subs[]); + +/*! Helper function to the indices for the interpolation + * + * @param [in] value the value to search for + * @param [in] num_values the number of values + * @param [in] values the values + * @param [out] index the smallest index in [0;num_values-2] for which values[index] <= value + * @param [out] t the weight for the linear interpolation s.t. value == (1-t)*values[index] + t*values[index+1] + * + * @return 0 + */ +void NDTable_find_index(double value, int num_values, const double values[], int *index, double *t, NDTable_ExtrapMethod_t extrap_method); + +int NDTable_evaluate_internal(const NDTable_h table, const double *t, const int *subs, int *nsubs, int dim, NDTable_InterpMethod_t interp_method, NDTable_ExtrapMethod_t extrap_method, double *value, double *derivatives); + +NDTable_h NDTable_create_table(int ndims, const int *dims, const double *data, const double **scales); + +/*! Calculate the number of offsets from the dimensions + * + * @param [in] ndims the number of dimensions + * @param [in] dims the extent of the dimensions + * @param [out] offs array to write the offsets + */ +void NDTable_calculate_offsets(int ndims, const int dims[], int offs[]); + +/*! Calculate the number of elements from the dimensions + * + * @param [in] ndims the number of dimensions + * @param [in] dims the extent of the dimensions + * + * @return the number of elements + */ +int NDTable_calculate_numel(int ndims, const int dims[]); + +#ifdef __cplusplus +} +#endif + +#endif /*NDTABLE_H_*/ \ No newline at end of file diff --git a/src/Core.c b/src/Core.c new file mode 100644 index 0000000..81d61bc --- /dev/null +++ b/src/Core.c @@ -0,0 +1,204 @@ +#include +#include +#include +#include +#include +#include + +#include "NDTable.h" + +#ifdef _WIN32 +#define ISFINITE(x) _finite(x) +#else +#define ISFINITE(x) isfinite(x) +#endif + +static char error_message[MAX_MESSAGE_LENGTH] = ""; + +void NDTable_set_error_message(const char *msg, ...) { + va_list vargs; + va_start(vargs, msg); + vsprintf(error_message, msg, vargs); + va_end(vargs); +} + +NDTable_h NDTable_alloc_table() { + return (NDTable_h )calloc(1, sizeof(NDTable_t)); +} + +void NDTable_free_table(NDTable_h table) { + int i; + + if(!table) return; + + free(table->data); + + for(i = 0; i < MAX_NDIMS; i++) { + free(table->scales[i]); + } + + free(table); +} + +const char * NDTable_get_error_message() { + return error_message; +} + +void NDTable_calculate_offsets(int ndims, const int dims[], int *offs) { + int i; + + if(ndims < 1) { + return; + } + + offs[ndims-1] = 1; + + for(i = ndims-2; i >= 0; i--) { + offs[i] = offs[i+1] * dims[i+1]; + } +} + +int NDTable_validate_table(NDTable_h table) { + int i, j, numel, offs[MAX_NDIMS]; + double v; + + // check the rank + if(table->ndims < 0 || table->ndims > 32) { + NDTable_set_error_message("The rank of '%s' in '%s' must be in the range [0;32] but was %d", "table->datasetname", "table->filename", table->ndims); + return -1; + } + + // check the extent of the dimensions + for(i = 0; i < table->ndims; i++) { + if(table->dims[i] < 1) { + NDTable_set_error_message("Extent of dimension %d of '%s' in '%s' must be >=0 but was %d", i, "table->datasetname", "table->filename", table->dims[i]); + return -1; + } + } + + // check the number of values + numel = 1; + for(i = 0; i < table->ndims; i++) { + numel = numel * table->dims[i]; + } + + if(table->numel != numel) { + NDTable_set_error_message("The size of '%s' in '%s' does not match its extent", "table->datasetname", "table->filename"); + return -1; + } + + // check the offsets + NDTable_calculate_offsets(table->ndims, table->dims, offs); + for(i = 0; i < table->ndims; i++) { + if(table->offs[i] != offs[i]) { + NDTable_set_error_message("The offset[%d] of '%s' in '%s' must be %d but was %d", i, "table->datasetname", "table->filename", offs[i], table->offs[i]); + return -1; + } + } + + // check the scales + for(i = 0; i < table->ndims; i++) { + + // make sure a scale is set + if(table->scales[i] == NULL) { + NDTable_set_error_message("Scale for dimension %d of '%s' in '%s' is not set", i, "table->datasetname", "table->filename"); + return -1; + } + + // check strict monotonicity + v = table->scales[i][0]; + for(j = 1; j < table->dims[i]; j++) { + if(v >= table->scales[i][j]) { + NDTable_set_error_message("Scale for dimension %d of '%s' in '%s' is not strictly monotonic increasing at index %d", i, "table->datasetname", "table->filename", j); + return -1; + } + v = table->scales[i][j]; + } + + if(table->offs[i] != offs[i]) { + NDTable_set_error_message("The offset[%d] of '%s' in '%s' must be %d but was %d", i, "table->datasetname", "table->filename", offs[i], table->offs[i]); + return -1; + } + } + + // check the data for non-finite values + for(i = 0; i < table->numel; i++) { + if(!ISFINITE(table->data[i])) { + NDTable_set_error_message("The data value at index %d of '%s' in '%s' is not finite", i, "table->datasetname", "table->filename"); + return -1; + } + } + + return 0; +} + +int NDTable_calculate_numel(int ndims, const int dims[]) { + int i, numel = 1; + + for(i = 0; i < ndims; i++) { + numel = numel * dims[i]; + } + + return numel; +} + +void NDTable_ind2sub(const int index, const NDTable_h table, int *subs) { + int i, n = index; // number of remaining elements + + for(i = 0; i < table->ndims; i++) { + subs[i] = n / table->offs[i]; + n -= subs[i] * table->offs[i]; + } +} + +void NDTable_sub2ind(const int *subs, const NDTable_h table, int *index) { + int i, k = 1; + + (*index) = 0; + + for(i = table->ndims-1; i >= 0; i--) { + (*index) += subs[i] * k; + k *= table->dims[i]; // TODO use pre-calculated offsets + } +} + +double NDTable_get_value_subs(const NDTable_h table, const int subs[]) { + int index; + NDTable_sub2ind(subs, table, &index); + return table->data[index]; +} + +NDTable_h NDTable_create_table(int ndims, const int *dims, const double *data, const double **scales) { + int i, j; + NDTable_h table = NULL; + + // check scales for strict monotonicity + for(i = 0; i < ndims; i++) { + for(j = 0; j < dims[i] - 1; j++) { + if (scales[i][j] >= scales[i][j + 1]) { + NDTable_set_error_message("The scale for dimension %d is not strictly monotonic at index %d", i + 1, j + 1); + goto out; + } + } + } + + table = NDTable_alloc_table(); + + table->ndims = ndims; + + table->numel = NDTable_calculate_numel(ndims, dims); + + NDTable_calculate_offsets(ndims, dims, table->offs); + + table->data = (double *)malloc(table->numel * sizeof(double)); + memcpy(table->data, data, table->numel * sizeof(double)); + + for(i = 0; i < ndims; i++) { + table->dims[i] = dims[i]; + table->scales[i] = (double *)malloc(dims[i] * sizeof(double)); + memcpy(table->scales[i], scales[i], dims[i] * sizeof(double)); + } + +out: + return table; +} diff --git a/src/Interpolation.c b/src/Interpolation.c new file mode 100644 index 0000000..b373f9b --- /dev/null +++ b/src/Interpolation.c @@ -0,0 +1,566 @@ +#include +#include +#include +#include +#include + +#include "NDTable.h" + +#ifndef MAX +#define MAX(a,b) (((a) > (b)) ? (a) : (b)) +#endif + +#ifndef MIN +#define MIN(a,b) (((a) < (b)) ? (a) : (b)) +#endif + +#ifndef NAN +static const unsigned long __nan[2] = { 0xffffffff, 0x7fffffff }; +#define NAN (*(const float *) __nan) +#endif + +#ifdef _WIN32 +#define ISFINITE(x) _finite(x) +#else +#define ISFINITE(x) isfinite(x) +#endif + +/** +Prototype of an interpolation function + +@param table [in] table handle +@param t [in] weights for the interpolation (normalized) +@param subs [in] subscripts of the left sample point +@param nsubs [in,out] subscripts of the right (next) sample point +@param dim [in] index of the current dimension +@param interp_method [in] index of the current dimension +@param extrap_method [in] index of the current dimension +@param value [out] interpoated value +@param derivatives [out] partial derivatives + +@return status code +*/ +typedef int(*interp_fun)(const NDTable_h table, const double *t, const int *subs, int *nsubs, int dim, NDTable_InterpMethod_t interp_method, NDTable_ExtrapMethod_t extrap_method, double *value, double derivatives[]); + +// forward declare inter- and extrapolation functions +static int interp_hold (const NDTable_h table, const double *t, const int *subs, int *nsubs, int dim, NDTable_InterpMethod_t interp_method, NDTable_ExtrapMethod_t extrap_method, double *value, double derivatives[]); +static int interp_nearest (const NDTable_h table, const double *t, const int *subs, int *nsubs, int dim, NDTable_InterpMethod_t interp_method, NDTable_ExtrapMethod_t extrap_method, double *value, double derivatives[]); +static int interp_linear (const NDTable_h table, const double *t, const int *subs, int *nsubs, int dim, NDTable_InterpMethod_t interp_method, NDTable_ExtrapMethod_t extrap_method, double *value, double derivatives[]); +static int interp_akima (const NDTable_h table, const double *t, const int *subs, int *nsubs, int dim, NDTable_InterpMethod_t interp_method, NDTable_ExtrapMethod_t extrap_method, double *value, double derivatives[]); +static int interp_fritsch_butland (const NDTable_h table, const double *t, const int *subs, int *nsubs, int dim, NDTable_InterpMethod_t interp_method, NDTable_ExtrapMethod_t extrap_method, double *value, double derivatives[]); +static int interp_steffen (const NDTable_h table, const double *t, const int *subs, int *nsubs, int dim, NDTable_InterpMethod_t interp_method, NDTable_ExtrapMethod_t extrap_method, double *value, double derivatives[]); +static int extrap_hold (const NDTable_h table, const double *t, const int *subs, int *nsubs, int dim, NDTable_InterpMethod_t interp_method, NDTable_ExtrapMethod_t extrap_method, double *value, double derivatives[]); +static int extrap_linear (const NDTable_h table, const double *t, const int *subs, int *nsubs, int dim, NDTable_InterpMethod_t interp_method, NDTable_ExtrapMethod_t extrap_method, double *value, double derivatives[]); + + +void NDTable_find_index(double value, int nvalues, const double *values, int *index, double *t, NDTable_ExtrapMethod_t extrap_method) { + int i; + double a, b; + double min = values[0]; + double max = values[nvalues - 1]; + double range = max - min; + + if(nvalues < 2) { + *t = 0.0; + *index = 0; + return; + } + + // estimate the index and make sure that i <= 0 and i <= 2nd last + i = MAX(0, MIN((int)(nvalues * (value - min) / range), nvalues - 2)); + + // go up until value < values[i+1] + while (i < nvalues - 2 && value > values[i+1]) { i++; } + + // go down until values[i] < value + while (i > 0 && value < values[i]) { i--; } + + a = values[i]; + b = values[i+1]; + + *t = (value - a) / (b - a); + + *index = i; +} + +int NDTable_evaluate(NDTable_h table, int nparams, const double params[], NDTable_InterpMethod_t interp_method, NDTable_ExtrapMethod_t extrap_method, double *value) { + int i; + double t [MAX_NDIMS]; // the weights for the interpolation + int subs [MAX_NDIMS]; // the subscripts + int nsubs [MAX_NDIMS]; // the neighboring subscripts + double derivatives [MAX_NDIMS]; + + // TODO: add null check + + // if the dataset is scalar return the value + if (table->ndims == 0) { + *value = table->data[0]; + return NDTABLE_INTERPSTATUS_OK; + } + + // find entry point and weights + for (i = 0; i < table->ndims; i++) { + NDTable_find_index(params[i], table->dims[i], table->scales[i], &subs[i], &t[i], extrap_method); + } + + return NDTable_evaluate_internal(table, t, subs, nsubs, 0, interp_method, extrap_method, value, derivatives); +} + +int NDTable_evaluate_derivative(NDTable_h table, int nparams, const double params[], const double delta_params[], NDTable_InterpMethod_t interp_method, NDTable_ExtrapMethod_t extrap_method, double *value) { + int i, err; + double t[MAX_NDIMS]; // the weights for the interpolation + int subs[MAX_NDIMS]; // the subscripts + int nsubs[MAX_NDIMS]; // the neighboring subscripts + double derivatives[MAX_NDIMS]; + + // TODO: add null check + + // if the dataset is scalar return the value + if (table->ndims == 0) { + *value = table->data[0]; + return NDTABLE_INTERPSTATUS_OK; + } + + // find entry point and weights + for (i = 0; i < table->ndims; i++) { + NDTable_find_index(params[i], table->dims[i], table->scales[i], &subs[i], &t[i], extrap_method); + } + + if ((err = NDTable_evaluate_internal(table, t, subs, nsubs, 0, interp_method, extrap_method, value, derivatives)) != 0) { + return err; + } + + *value = 0.0; + + for (i = 0; i < nparams; i++) { + *value += delta_params[i] * derivatives[i]; + } + + return 0; +} + +int NDTable_evaluate_internal(const NDTable_h table, const double *t, const int *subs, int *nsubs, int dim, NDTable_InterpMethod_t interp_method, NDTable_ExtrapMethod_t extrap_method, double *value, double derivatives[]) { + + interp_fun func; + + // check arguments + if (table == NULL || t == NULL || subs == NULL || nsubs == NULL || value == NULL || derivatives == NULL) { + return -1; + } + + if (dim >= table->ndims) { + *value = NDTable_get_value_subs(table, nsubs); + return 0; + } + + // find the right function: + if(table->dims[dim] < 2) { + func = interp_hold; + } else if (t[dim] < 0.0 || t[dim] > 1.0) { + // extrapolate + switch (extrap_method) { + case NDTABLE_EXTRAP_HOLD: + func = extrap_hold; + break; + case NDTABLE_EXTRAP_LINEAR: + switch (interp_method) { + case NDTABLE_INTERP_AKIMA: func = interp_akima; break; + case NDTABLE_INTERP_FRITSCH_BUTLAND: func = interp_fritsch_butland; break; + default: func = extrap_linear; break; + } + break; + default: + NDTable_set_error_message("Requested value is outside data range"); + return -1; + } + } else { + // interpolate + switch (interp_method) { + case NDTABLE_INTERP_HOLD: func = interp_hold; break; + case NDTABLE_INTERP_NEAREST: func = interp_nearest; break; + case NDTABLE_INTERP_LINEAR: func = interp_linear; break; + case NDTABLE_INTERP_AKIMA: func = interp_akima; break; + case NDTABLE_INTERP_FRITSCH_BUTLAND: func = interp_fritsch_butland; break; + case NDTABLE_INTERP_STEFFEN: func = interp_steffen; break; + default: return -1; // TODO: set error message + } + } + + return (*func)(table, t, subs, nsubs, dim, interp_method, extrap_method, value, derivatives); +} + +static int interp_hold(const NDTable_h table, const double *t, const int *subs, int *nsubs, int dim, NDTable_InterpMethod_t interp_method, NDTable_ExtrapMethod_t extrap_method, double *value, double der_values[]) { + nsubs[dim] = subs[dim]; // always take the left sample value + der_values[dim] = 0; + return NDTable_evaluate_internal(table, t, subs, nsubs, dim + 1, interp_method, extrap_method, value, der_values); +} + +static int interp_nearest(const NDTable_h table, const double *t, const int *subs, int *nsubs, int dim, NDTable_InterpMethod_t interp_method, NDTable_ExtrapMethod_t extrap_method, double *value, double der_values[]) { + int err; + nsubs[dim] = t[dim] < 0.5 ? subs[dim] : subs[dim] + 1; + der_values[dim] = 0; + + if ((err = NDTable_evaluate_internal(table, t, subs, nsubs, dim + 1, interp_method, extrap_method, value, der_values)) != 0) { + return err; + } + + // if the value is not finite return NAN + if (!ISFINITE(*value)) { + *value = NAN; + der_values[dim] = NAN; + } + + return 0; +} + +static int interp_linear(const NDTable_h table, const double *t, const int *subs, int *nsubs, int dim, NDTable_InterpMethod_t interp_method, NDTable_ExtrapMethod_t extrap_method, double *value, double der_values[]) { + int err; + double a, b; + + // get the left value + nsubs[dim] = subs[dim]; + if ((err = NDTable_evaluate_internal(table, t, subs, nsubs, dim + 1, interp_method, extrap_method, &a, der_values)) != 0) { + return err; + } + + // get the right value + nsubs[dim] = subs[dim] + 1; + if ((err = NDTable_evaluate_internal(table, t, subs, nsubs, dim + 1, interp_method, extrap_method, &b, der_values)) != 0) { + return err; + } + + // if any of the values is not finite return NAN + if (!ISFINITE(a) || !ISFINITE(b)) { + *value = NAN; + der_values[dim] = NAN; + return 0; + } + + // calculate the interpolated value + *value = (1 - t[dim]) * a + t[dim] * b; + + // calculate the derivative + der_values[dim] = (b - a) / (table->scales[dim][subs[dim] + 1] - table->scales[dim][subs[dim]]); + + return 0; +} + +static void cubic_hermite_spline(const double x0, const double x1, const double y0, const double y1, const double t, const double c[4], double *value, double *derivative) { + + double v; + + if (t < 0) { // extrapolate left + + *value = y0 + c[2] * ((x1 - x0) * t); + *derivative = c[2]; + + } else if (t <= 1) { // interpolate + + v = (x1 - x0) * t; + *value = ((c[0] * v + c[1]) * v + c[2]) * v + c[3]; + *derivative = (3 * c[0] * v + (2 * c[1])) * v + c[2]; + + } else { // extrapolate right + + v = x1 - x0; + *value = y1 + ((3 * c[0] * v + 2 * c[1]) * v + c[2]) * (v * (t - 1)); + *derivative = (3 * c[0] * v + 2 * c[1]) * v + c[2]; + } +} + +static int interp_akima(const NDTable_h table, const double *t, const int *subs, int *nsubs, int dim, NDTable_InterpMethod_t interp_method, NDTable_ExtrapMethod_t extrap_method, double *value, double der_values[]) { + + double x[6] = { 0, 0, 0, 0, 0, 0}; + double y[6] = { 0, 0, 0, 0, 0, 0}; + double c[4] = { 0, 0, 0, 0 }; // spline coefficients + double d[5] = { 0, 0, 0, 0, 0 }; // divided differences + double c2 = 0; + double dx = 0; + double a = 0; + + int n = table->dims[dim]; // extent of the current dimension + int sub = subs[dim]; // subscript of current dimension + int err, i, idx; + + for (i = 0; i < 6; i++) { + idx = sub - 2 + i; + + if (idx >= 0 && idx < n) { + x[i] = table->scales[dim][idx]; + + nsubs[dim] = idx; + if ((err = NDTable_evaluate_internal(table, t, subs, nsubs, dim + 1, interp_method, extrap_method, &y[i], der_values)) != 0) { + return err; + } + } + } + + // if any of the values is not finite return NAN + for (i = 0; i < 6; i++) { + if (!ISFINITE(y[i])) { + *value = NAN; + der_values[dim] = NAN; + return 0; + } + } + + // calculate the divided differences + for (i = MAX(0, 2 - sub); i < MIN(5, 1 + n - sub); i++) { + d[i] = (y[i + 1] - y[i]) / (x[i + 1] - x[i]); + } + + // pad left + if (sub < 2) { + if (sub < 1) { + d[1] = 2.0 * d[2] - d[3]; + } + d[0] = 2.0 * d[1] - d[2]; + } + + // pad right + if (sub > n - 4) { + if (sub > n - 3) { + d[3] = 2.0 * d[2] - d[1]; + } + d[4] = 2.0 * d[3] - d[2]; + } + + // initialize the left boundary slope + c2 = fabs(d[3] - d[2]) + fabs(d[1] - d[0]); + + if (c2 > 0) { + a = fabs(d[1] - d[0]) / c2; + c2 = (1 - a) * d[1] + a * d[2]; + } else { + c2 = 0.5 * d[1] + 0.5 * d[2]; + } + + // calculate the coefficients + dx = x[3] - x[2]; + + c[2] = c2; + c2 = fabs(d[4] - d[3]) + fabs(d[2] - d[1]); + + if (c2 > 0) { + a = fabs(d[2] - d[1]) / c2; + c2 = (1 - a) * d[2] + a * d[3]; + } else { + c2 = 0.5 * d[2] + 0.5 * d[3]; + } + + c[1] = (3 * d[2] - 2 * c[2] - c2) / dx; + c[0] = (c[2] + c2 - 2 * d[2]) / (dx * dx); + + c[3] = y[2]; + + cubic_hermite_spline(x[2], x[3], y[2], y[3], t[dim], c, value, &der_values[dim]); + + return 0; +} + +static int interp_fritsch_butland(const NDTable_h table, const double *t, const int *subs, int *nsubs, int dim, NDTable_InterpMethod_t interp_method, NDTable_ExtrapMethod_t extrap_method, double *value, double der_values[]) { + + double x [4] = { 0, 0, 0, 0 }; + double y [4] = { 0, 0, 0, 0 }; + double dx[3] = { 0, 0, 0 }; + double d [3] = { 0, 0, 0 }; // divided differences + double c [4] = { 0, 0, 0, 0 }; // spline coefficients + double c2 = 0; + + int n = table->dims[dim]; // extent of the current dimension + int sub = subs[dim]; // subscript of current dimension + int err, i, idx; + + for (i = 0; i < 4; i++) { + idx = sub - 1 + i; + + if (idx >= 0 && idx < n) { + x[i] = table->scales[dim][idx]; + + nsubs[dim] = idx; + if ((err = NDTable_evaluate_internal(table, t, subs, nsubs, dim + 1, interp_method, extrap_method, &y[i], der_values)) != 0) { + return err; + } + } + } + + // if any of the values is not finite return NAN + for (i = 0; i < 4; i++) { + if (!ISFINITE(y[i])) { + *value = NAN; + der_values[dim] = NAN; + return 0; + } + } + + // calculate the divided differences + //for (i = MAX(0, 1 - sub); i < MIN(3, n - 1 - sub); i++) { + for (i = 0; i < 3; i++) { + dx[i] = x[i + 1] - x[i]; + d[i] = (y[i + 1] - y[i]) / dx[i]; + } + + // initialize the left boundary slope + + // calculate the coefficients + + if (sub == 0) { + c2 = d[1]; + } else if (d[0] == 0 || d[1] == 0 || (d[0] < 0 && d[1] > 0) || (d[0] > 0 && d[1] < 0)) { + c2 = 0; + } else { + c2 = 3 * (dx[0] + dx[1]) / ((dx[0] + 2 * dx[1]) / d[0] + (dx[1] + 2 * dx[0]) / d[1]); + } + + c[2] = c2; + + if (sub == n - 2) { + c2 = d[1]; + } else if (d[1] == 0 || d[2] == 0 || (d[1] < 0 && d[2] > 0) || (d[1] > 0 && d[2] < 0)) { + c2 = 0; + } else { + c2 = 3 * (dx[1] + dx[2]) / ((dx[1] + 2 * dx[2]) / d[1] + (dx[2] + 2 * dx[1]) / d[2]); + } + + c[1] = (3 * d[1] - 2 * c[2] - c2) / dx[1]; + c[0] = (c[2] + c2 - 2 * d[1]) / (dx[1] * dx[1]); + + c[3] = y[1]; + + cubic_hermite_spline(x[1], x[2], y[1], y[2], t[dim], c, value, &der_values[dim]); + + return 0; +} + +static int interp_steffen(const NDTable_h table, const double *t, const int *subs, int *nsubs, int dim, NDTable_InterpMethod_t interp_method, NDTable_ExtrapMethod_t extrap_method, double *value, double der_values[]) { + + double x [4] = { 0, 0, 0, 0 }; + double y [4] = { 0, 0, 0, 0 }; + double dx[3] = { 0, 0, 0 }; + double d [3] = { 0, 0, 0 }; // divided differences + double c [4] = { 0, 0, 0, 0 }; // spline coefficients + double c2 = 0; + + const int n = table->dims[dim]; // extent of the current dimension + const int sub = subs[dim]; // subscript of current dimension + int err, i, idx; + + for (i = 0; i < 4; i++) { + idx = sub - 1 + i; + + if (idx >= 0 && idx < n) { + x[i] = table->scales[dim][idx]; + + nsubs[dim] = idx; + if ((err = NDTable_evaluate_internal(table, t, subs, nsubs, dim + 1, interp_method, extrap_method, &y[i], der_values)) != 0) { + return err; + } + } + } + + // if any of the values is not finite return NAN + for (i = 0; i < 4; i++) { + if (!ISFINITE(y[i])) { + *value = NAN; + der_values[dim] = NAN; + return 0; + } + } + + // calculate the divided differences + for (i = 0; i < 3; i++) { + dx[i] = x[i + 1] - x[i]; + d[i] = (y[i + 1] - y[i]) / dx[i]; + } + + // calculate the coefficients + if (sub == 0) { + c2 = d[1]; + } else if (d[0] == 0 || d[1] == 0 || (d[0] < 0 && d[1] > 0) || (d[0] > 0 && d[1] < 0)) { + c2 = 0; + } else { + double half_abs_c2, abs_di, abs_di1; + c2 = (d[0] * dx[1] + d[1] * dx[0]) / (dx[0] + dx[1]); + half_abs_c2 = 0.5 * fabs(c2); + abs_di = fabs(d[0]); + abs_di1 = fabs(d[1]); + if (half_abs_c2 > abs_di || half_abs_c2 > abs_di1) { + const double two_a = d[0] > 0 ? 2 : -2; + c2 = two_a*(abs_di < abs_di1 ? abs_di : abs_di1); + } + } + + c[2] = c2; + + if (sub == n - 2) { + c2 = d[1]; + } else if (d[1] == 0 || d[2] == 0 || (d[1] < 0 && d[2] > 0) || (d[1] > 0 && d[2] < 0)) { + c2 = 0; + } else { + double half_abs_c2, abs_di, abs_di1; + c2 = (d[1] * dx[2] + d[2] * dx[1]) / (dx[1] + dx[2]); + half_abs_c2 = 0.5 * fabs(c2); + abs_di = fabs(d[1]); + abs_di1 = fabs(d[2]); + if (half_abs_c2 > abs_di || half_abs_c2 > abs_di1) { + const double two_a = d[1] > 0 ? 2 : -2; + c2 = two_a*(abs_di < abs_di1 ? abs_di : abs_di1); + } + } + + c[1] = (3 * d[1] - 2 * c[2] - c2) / dx[1]; + c[0] = (c[2] + c2 - 2 * d[1]) / (dx[1] * dx[1]); + c[3] = y[1]; + + cubic_hermite_spline(x[1], x[2], y[1], y[2], t[dim], c, value, &der_values[dim]); + + return 0; +} + +static int extrap_hold(const NDTable_h table, const double *t, const int *subs, int *nsubs, int dim, NDTable_InterpMethod_t interp_method, NDTable_ExtrapMethod_t extrap_method, double *value, double der_values[]) { + int err; + nsubs[dim] = t[dim] < 0.0 ? subs[dim] : subs[dim] + 1; + der_values[dim] = 0; + + if ((err = NDTable_evaluate_internal(table, t, subs, nsubs, dim + 1, interp_method, extrap_method, value, der_values)) != 0) { + return err; + } + + // if the value is not finite return NAN + if (!ISFINITE(*value)) { + *value = NAN; + der_values[dim] = NAN; + } + + return 0; +} + +static int extrap_linear(const NDTable_h table, const double *t, const int *subs, int *nsubs, int dim, NDTable_InterpMethod_t interp_method, NDTable_ExtrapMethod_t extrap_method, double *value, double der_values[]) { + int err; + double a, b; + + nsubs[dim] = subs[dim]; + if ((err = NDTable_evaluate_internal(table, t, subs, nsubs, dim + 1, interp_method, extrap_method, &a, der_values)) != 0) { + return err; + } + + nsubs[dim] = subs[dim] + 1; + if ((err = NDTable_evaluate_internal(table, t, subs, nsubs, dim + 1, interp_method, extrap_method, &b, der_values)) != 0) { + return err; + } + + // if any of the values is not finite return NAN + if (!ISFINITE(a) || !ISFINITE(b)) { + *value = NAN; + der_values[dim] = NAN; + return 0; + } + + // calculate the extrapolated value + *value = (1 - t[dim]) * a + t[dim] * b; + + // calculate the derivative + der_values[dim] = (b - a) / (table->scales[dim][subs[dim] + 1] - table->scales[dim][subs[dim]]); + + return 0; +} diff --git a/test/CoreTests.cpp b/test/CoreTests.cpp new file mode 100644 index 0000000..49af4f0 --- /dev/null +++ b/test/CoreTests.cpp @@ -0,0 +1,72 @@ +#include +#include + +#include "NDTable.h" + + +TEST_CASE("Core") { + + SECTION("Set error message") { + NDTable_set_error_message("%d plus %.1f equals %s", 1, 1.5, "two point five"); + auto message = NDTable_get_error_message(); + REQUIRE_THAT(message, Catch::Matchers::Equals("1 plus 1.5 equals two point five")); + } + + SECTION("Convert index to subscripts with legal index") { + int subs[3] = { 0 }; + auto ds = NDTable_alloc_table(); + + ds->ndims = 3; + ds->offs[0] = 12; + ds->offs[1] = 4; + ds->offs[2] = 1; + + NDTable_ind2sub(17, ds, subs); + + REQUIRE(1 == subs[0]); + REQUIRE(1 == subs[1]); + REQUIRE(1 == subs[2]); + } + + SECTION("Convert subscripts to index with legal subscripts") { + + auto ds = NDTable_alloc_table(); + + ds->ndims = 3; + ds->dims[0] = 2; + ds->dims[1] = 3; + ds->dims[2] = 4; + + int subs[3] = { 1, 1, 1 }; + int index = 0; + + NDTable_sub2ind(subs, ds, &index); + + REQUIRE(17 == index); + } + + SECTION("Allocate and free") { + + // de-allocate with null pointer + NDTable_free_table(nullptr); + + auto ds = NDTable_alloc_table(); + + // de-allocate with no data + NDTable_free_table(ds); + + // get a new dataset + ds = NDTable_alloc_table(); + + // allocate some dummy space + ds->data = static_cast(malloc(1)); + + for(int i = 0; i < MAX_NDIMS; i++) { + ds->scales[i] = (double *)malloc(1); + } + + // de-allocate with data + NDTable_free_table(ds); + } + +} diff --git a/test/InterpolationTests.cpp b/test/InterpolationTests.cpp new file mode 100644 index 0000000..16e5122 --- /dev/null +++ b/test/InterpolationTests.cpp @@ -0,0 +1,388 @@ +#include + +#define _USE_MATH_DEFINES // for M_PI +#include +#include + +#include "NDTable.h" + + +//#ifdef _MSC_VER +//static const unsigned long __nan[2] = { 0xffffffff, 0x7fffffff }; +//#define NAN (*(const float *) __nan) +//#define INFINITY (DBL_MAX + DBL_MAX) +//#endif + +#define N_NON_FINITE 3 +static const double non_finite[N_NON_FINITE] = { NAN, INFINITY, -INFINITY }; + +TEST_CASE("Interplation") { + + SECTION("Find index") { + double value = 1.2; + double values[4] = { 0, 1, 2, 3 }; + int num_values = 4; + int index = -1; + double t = 0; + + NDTable_find_index(value, num_values, values, &index, &t, NDTABLE_EXTRAP_HOLD); + + REQUIRE(index == 1); + } + + SECTION("Interpolate nearest") { + double x[2] = { 0, 1 }; + double y[2] = { 2, 3 }; + + NDTable_t ds; + ds.data = y; + ds.scales[0] = x; + ds.dims[0] = 2; + ds.numel = 2; + ds.ndims = 1; + ds.offs[0] = 1; + + int subs[1] = { 0 }; + int nsubs[1]; + double t[1]; + + double value; + double derivatives[1]; + + // interpolate left boundary + t[0] = 0.0; + REQUIRE(0 == NDTable_evaluate_internal(&ds, t, subs, nsubs, 0, NDTABLE_INTERP_NEAREST, NDTABLE_EXTRAP_NONE, &value, derivatives)); + REQUIRE(2.0 == value); + REQUIRE(0.0 == derivatives[0]); + + // interpolate before step + t[0] = 0.49; + REQUIRE(0 == NDTable_evaluate_internal(&ds, t, subs, nsubs, 0, NDTABLE_INTERP_NEAREST, NDTABLE_EXTRAP_NONE, &value, derivatives)); + REQUIRE(2.0 == value); + REQUIRE(0.0 == derivatives[0]); + + // interpolate after step + t[0] = 0.51; + REQUIRE(0 == NDTable_evaluate_internal(&ds, t, subs, nsubs, 0, NDTABLE_INTERP_NEAREST, NDTABLE_EXTRAP_NONE, &value, derivatives)); + REQUIRE(3.0 == value); + REQUIRE(0.0 == derivatives[0]); + + // interpolate right boundary + t[0] = 1.0; + REQUIRE(0 == NDTable_evaluate_internal(&ds, t, subs, nsubs, 0, NDTABLE_INTERP_NEAREST, NDTABLE_EXTRAP_NONE, &value, derivatives)); + REQUIRE(3.0 == value); + REQUIRE(0.0 == derivatives[0]); + + // non-finite values + for(int i = 0; i < N_NON_FINITE; i++) { + // left + y[0] = non_finite[i]; + y[1] = 3; + t[0] = 0.0; + REQUIRE(0 == NDTable_evaluate_internal(&ds, t, subs, nsubs, 0, NDTABLE_INTERP_NEAREST, NDTABLE_EXTRAP_LINEAR, &value, derivatives)); + REQUIRE(1 == std::isnan(value)); + REQUIRE(1 == std::isnan(derivatives[0])); + + // right boundary + y[0] = 3; + y[1] = non_finite[i]; + t[0] = 1.0; + REQUIRE(0 == NDTable_evaluate_internal(&ds, t, subs, nsubs, 0, NDTABLE_INTERP_NEAREST, NDTABLE_EXTRAP_LINEAR, &value, derivatives)); + REQUIRE(1 == std::isnan(value)); + REQUIRE(1 == std::isnan(derivatives[0])); + } + } + + SECTION("Interplate linear") { + double x[2] = { 0, 1 }; + double y[2] = { 2, 3 }; + + NDTable_t ds; + ds.data = y; + ds.scales[0] = x; + ds.dims[0] = 2; + ds.numel = 2; + ds.ndims = 1; + ds.offs[0] = 1; + + int subs[1] = { 0 }; + int nsubs[1]; + double t[1] = { 0.3 }; + + double value; + double derivatives[1]; + + REQUIRE(0 == NDTable_evaluate_internal(&ds, t, subs, nsubs, 0, NDTABLE_INTERP_LINEAR, NDTABLE_EXTRAP_NONE, &value, derivatives)); + REQUIRE(2.3 == value); + REQUIRE(1.0 == derivatives[0]); + + // non-finite values + for(int i = 0; i < N_NON_FINITE; i++) { + // left value non-finite + y[0] = non_finite[i]; + y[1] = 2; + REQUIRE(0 == NDTable_evaluate_internal(&ds, t, subs, nsubs, 0, NDTABLE_INTERP_LINEAR, NDTABLE_EXTRAP_NONE, &value, derivatives)); + REQUIRE(1 == std::isnan(value)); + REQUIRE(1 == std::isnan(derivatives[0])); + + // right value non-finite + y[0] = non_finite[i]; + y[1] = 3; + REQUIRE(0 == NDTable_evaluate_internal(&ds, t, subs, nsubs, 0, NDTABLE_INTERP_LINEAR, NDTABLE_EXTRAP_NONE, &value, derivatives)); + REQUIRE(1 == std::isnan(value)); + REQUIRE(1 == std::isnan(derivatives[0])); + } + } + + SECTION("Interpolate akima") { + #define N 8 + + int i, j, k; + + double x[N]; + double y[N]; + + for (i = 0; i < N; i++) { + x[i] = i * ((double)N) / ((double) N-1); + y[i] = sin(x[i] * M_PI); + } + + NDTable_t ds; + ds.data = y; + ds.scales[0] = x; + ds.dims[0] = N; + ds.numel = N; + ds.ndims = 1; + ds.offs[0] = 1; + + int subs[1] = { 1 }; + int nsubs[1]; + double t[1] = { 0.99 }; + + double value; + double derivatives[1]; + + REQUIRE(0 == NDTable_evaluate_internal(&ds, t, subs, nsubs, 0, NDTABLE_INTERP_AKIMA, NDTABLE_EXTRAP_NONE, &value, derivatives)); + + // non-finite values + subs[0] = 2; + t[0] = 0.5; + for(i = 0; i < N_NON_FINITE; i++) { + for(j = 0; j < N; j++) { + for (k = 0; k < N; k++) { + y[i] = sin(x[i] * M_PI); + } + y[j] = non_finite[i]; + REQUIRE(0 == NDTable_evaluate_internal(&ds, t, subs, nsubs, 0, NDTABLE_INTERP_AKIMA, NDTABLE_EXTRAP_NONE, &value, derivatives)); + if (!std::isnan(value)) { + i = i; + } + REQUIRE(1 == std::isnan(value)); + REQUIRE(1 == std::isnan(derivatives[0])); + } + } + + } + + SECTION("Extrapolate hold") { + double x[2] = { 0, 1 }; + double y[2] = { 2, 3 }; + + NDTable_t ds; + ds.data = y; + ds.scales[0] = x; + ds.dims[0] = 2; + ds.numel = 2; + ds.ndims = 1; + ds.offs[0] = 1; + + int subs[1] = { 0 }; + int nsubs[1]; + double t[1]; + + double value; + double derivatives[1]; + + int i; + + // extrapolate left + t[0] = -0.1; + REQUIRE(0 == NDTable_evaluate_internal(&ds, t, subs, nsubs, 0, NDTABLE_INTERP_NEAREST, NDTABLE_EXTRAP_HOLD, &value, derivatives)); + REQUIRE(2.0 == value); + REQUIRE(0.0 == derivatives[0]); + + // extrapolate right + t[0] = 1.1; + REQUIRE(0 == NDTable_evaluate_internal(&ds, t, subs, nsubs, 0, NDTABLE_INTERP_NEAREST, NDTABLE_EXTRAP_HOLD, &value, derivatives)); + REQUIRE(3.0 == value); + REQUIRE(0.0 == derivatives[0]); + + // non-finite values + for(i = 0; i < N_NON_FINITE; i++) { + y[0] = y[1] = non_finite[i]; + + // left + t[0] = -0.1; + REQUIRE(0 == NDTable_evaluate_internal(&ds, t, subs, nsubs, 0, NDTABLE_INTERP_NEAREST, NDTABLE_EXTRAP_HOLD, &value, derivatives)); + REQUIRE(1 == std::isnan(value)); + REQUIRE(1 == std::isnan(derivatives[0])); + + // right + t[0] = 1.1; + REQUIRE(0 == NDTable_evaluate_internal(&ds, t, subs, nsubs, 0, NDTABLE_INTERP_NEAREST, NDTABLE_EXTRAP_HOLD, &value, derivatives)); + REQUIRE(1 == std::isnan(value)); + REQUIRE(1 == std::isnan(derivatives[0])); + } + } + + SECTION("Extrapolate linear") { + double x[2] = { 0, 1 }; + double y[2] = { 2, 3 }; + + NDTable_t ds; + ds.data = y; + ds.scales[0] = x; + ds.dims[0] = 2; + ds.numel = 2; + ds.ndims = 1; + ds.offs[0] = 1; + + int subs[1] = { 0 }; + int nsubs[1]; + double t[1]; + + double value; + double derivatives[1]; + + int i, j; + + // extrapolate left + t[0] = -1.0; + REQUIRE(0 == NDTable_evaluate_internal(&ds, t, subs, nsubs, 0, NDTABLE_INTERP_NEAREST, NDTABLE_EXTRAP_LINEAR, &value, derivatives)); + REQUIRE(1.0 == value); + REQUIRE(1.0 == derivatives[0]); + + // extrapolate right + t[0] = 2.0; + REQUIRE(0 == NDTable_evaluate_internal(&ds, t, subs, nsubs, 0, NDTABLE_INTERP_NEAREST, NDTABLE_EXTRAP_LINEAR, &value, derivatives)); + REQUIRE(4.0 == value); + REQUIRE(1.0 == derivatives[0]); + + // non-finite values + for(i = 0; i < N_NON_FINITE; i++) { + for(j = 0; j < 2; j++) { + y[0] = 2; y[1] = 3; + y[j] = non_finite[i]; + REQUIRE(0 == NDTable_evaluate_internal(&ds, t, subs, nsubs, 0, NDTABLE_INTERP_NEAREST, NDTABLE_EXTRAP_LINEAR, &value, derivatives)); + REQUIRE(1 == std::isnan(value)); + REQUIRE(1 == std::isnan(derivatives[0])); + } + } + } + + SECTION("Extrapolate none") { + double x[2] = { 0, 1 }; + double y[2] = { 2, 3 }; + + NDTable_t ds; + ds.data = y; + ds.scales[0] = x; + ds.dims[0] = 2; + ds.numel = 2; + ds.ndims = 1; + ds.offs[0] = 1; + + int subs[1] = { 0 }; + int nsubs[1]; + double t[1]; + + double value; + double derivatives[MAX_NDIMS]; + + // extrapolate left + t[0] = -1; + REQUIRE(-1 == NDTable_evaluate_internal(&ds, t, subs, nsubs, 0, NDTABLE_INTERP_NEAREST, NDTABLE_EXTRAP_NONE, &value, derivatives)); + + // extrapolate right + t[0] = 2.0; + REQUIRE(-1 == NDTable_evaluate_internal(&ds, t, subs, nsubs, 0, NDTABLE_INTERP_NEAREST, NDTABLE_EXTRAP_NONE, &value, derivatives)); + } + + SECTION("Dummy dimension 1-d") { + // this tests checks the direct return of the sample value + // without interpolation for dimensions with extent < 2 + + // create a 1-d dataset with length 1 + double x[1] = { 0 }; + double y[1] = { 1.1 }; + + NDTable_t ds; + ds.data = y; + ds.scales[0] = x; + ds.dims[0] = 1; + ds.numel = 1; + ds.ndims = 1; + ds.offs[0] = 1; + + int subs[1] = { 0 }; + int nsubs[1] = { 0 }; // not used + double t[1] = { 0.0 }; + + double value; + double derivative; + + REQUIRE(0 == NDTable_evaluate_internal(&ds, t, subs, nsubs, 0, NDTABLE_INTERP_NEAREST, NDTABLE_EXTRAP_NONE, &value, &derivative)); + REQUIRE(y[0] == value); + REQUIRE(0.0 == derivative); + } + + SECTION("Dummy dimension 2-d") { + // this tests checks the direct return of the sample value + // without interpolation for dimensions with extent < 2 + + // create a 2-d dataset with size (2,1) + double x[2] = { 0, 1 }; + double y[1] = { 1 }; + double z[2] = { 1, 2 }; + + NDTable_t ds; + ds.data = z; + ds.scales[0] = x; + ds.scales[1] = y; + ds.dims[0] = 2; + ds.dims[1] = 1; + ds.numel = 2; + ds.ndims = 2; + ds.offs[0] = 1; + ds.offs[1] = 2; + + int subs[2] = { 0, 0 }; + int nsubs[2]; + double t[2]; + + double value; + double der_values[2]; + + t[0] = -0.1; + t[1] = -0.1; + REQUIRE(0 == NDTable_evaluate_internal(&ds, t, subs, nsubs, 0, NDTABLE_INTERP_LINEAR, NDTABLE_EXTRAP_HOLD, &value, der_values)); + REQUIRE(1.0 == value); + REQUIRE(0.0 == der_values[0]); + REQUIRE(0.0 == der_values[1]); + + t[0] = 0.5; + t[1] = 0.5; + REQUIRE(0 == NDTable_evaluate_internal(&ds, t, subs, nsubs, 0, NDTABLE_INTERP_LINEAR, NDTABLE_EXTRAP_HOLD, &value, der_values)); + REQUIRE(1.5 == value); + REQUIRE(1.0 == der_values[0]); + REQUIRE(0.0 == der_values[1]); + + t[0] = 1.1; + t[1] = 1.1; + REQUIRE(0 == NDTable_evaluate_internal(&ds, t, subs, nsubs, 0, NDTABLE_INTERP_LINEAR, NDTABLE_EXTRAP_HOLD, &value, der_values)); + REQUIRE(2.0 == value); + REQUIRE(0.0 == der_values[0]); + REQUIRE(0.0 == der_values[1]); + } + +}