diff --git a/.github/workflows/build.yml b/.github/workflows/build.yml index 3552565..464f5e2 100644 --- a/.github/workflows/build.yml +++ b/.github/workflows/build.yml @@ -21,12 +21,9 @@ jobs: - name: linux64 image: ubuntu-22.04 cmake-args: -G "Unix Makefiles" - - name: darwin64 - image: macos-13 - cmake-args: -G "Unix Makefiles" -D CMAKE_OSX_ARCHITECTURES=x86_64 runs-on: ${{ matrix.image }} steps: - - uses: actions/checkout@v4.1.7 + - uses: actions/checkout@v7 with: submodules: recursive - run: | @@ -39,7 +36,7 @@ jobs: cmake ${{ matrix.cmake-args }} -B ${{ matrix.name }} . cmake --build ${{ matrix.name }} --config Release ctest --test-dir ${{ matrix.name }} --build-config Release - - uses: actions/upload-artifact@v4.4.0 + - uses: actions/upload-artifact@v7 with: name: ${{ matrix.name }} path: SDF @@ -48,25 +45,21 @@ jobs: runs-on: ubuntu-22.04 needs: [build] steps: - - uses: actions/checkout@v4.1.7 - - uses: actions/download-artifact@v4.1.8 + - uses: actions/checkout@v7 + - uses: actions/download-artifact@v8 with: name: win64 path: SDF - - uses: actions/download-artifact@v4.1.8 + - uses: actions/download-artifact@v8 with: name: win32 path: SDF - - uses: actions/download-artifact@v4.1.8 + - uses: actions/download-artifact@v8 with: name: linux64 path: SDF - - uses: actions/download-artifact@v4.1.8 - with: - name: darwin64 - path: SDF - run: find -type f -name .gitignore -delete - - uses: actions/upload-artifact@v4.4.0 + - uses: actions/upload-artifact@v7 with: name: SDF-Modelica path: | diff --git a/CMakeLists.txt b/CMakeLists.txt index 32a9e18..9a2b64d 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -1,4 +1,4 @@ -cmake_minimum_required (VERSION 3.2) +cmake_minimum_required (VERSION 3.20) set (CMAKE_CXX_STANDARD 14) @@ -99,3 +99,15 @@ set_property(TARGET ModelicaSDF_Test PROPERTY VS_DEBUGGER_WORKING_DIRECTORY "${C enable_testing() add_test(NAME ModelicaSDF_Test COMMAND $ WORKING_DIRECTORY "${CMAKE_CURRENT_SOURCE_DIR}/C/tests") + +add_library(NDTable SHARED + SDF/Resources/C-Sources/NDTable.h + SDF/Resources/C-Sources/ModelicaNDTable.c + SDF/Resources/C-Sources/NDTable.h + ThirdParty/ModelicaUtilities.h + ThirdParty/ModelicaUtilities.c +) + +target_include_directories(NDTable PUBLIC + ThirdParty +) diff --git a/SDF/Examples/InterpolationMethods.mo b/SDF/Examples/InterpolationMethods.mo index 41f0d04..a7aef64 100644 --- a/SDF/Examples/InterpolationMethods.mo +++ b/SDF/Examples/InterpolationMethods.mo @@ -1,39 +1,49 @@ within SDF.Examples; -model InterpolationMethods "Interpolate a sampled signal using different interpolation methods" +model InterpolationMethods "Interpolate a sampled signal using different interpolation and extrapolation methods" extends Modelica.Icons.Example; - parameter Real x[:] = linspace(0,2,6) * Modelica.Constants.pi; + parameter Real x[:] = linspace(0, 2, 6) * Modelica.Constants.pi; parameter Real y[size(x, 1)] = Modelica.Math.sin(x); parameter Real data[:] = cat(1, {1}, {size(x, 1)}, x, y); parameter String filename = Modelica.Utilities.Files.loadResource("modelica://SDF/Resources/Data/Examples/sine.sdf"); parameter String dataset = "/y"; - NDTable linear( + Modelica.Blocks.Sources.Ramp ramp( + duration=1, + offset=-1, + height=8.3) + annotation (Placement(transformation(extent={{-80,-30},{-60,-10}}))); + + NDTable hold( nin=1, readFromFile=false, data=data, - interpMethod=SDF.Types.InterpolationMethod.Linear, + interpMethod=SDF.Types.InterpolationMethod.Hold, extrapMethod=SDF.Types.ExtrapolationMethod.Linear, filename=filename, dataset=dataset) - annotation (Placement(transformation(extent={{-10,10},{10,30}}))); + annotation (Placement(transformation(extent={{-40,40},{-20,60}}))); - Modelica.Blocks.Sources.Ramp ramp( - duration=1, - offset=-1, - height=8.3) - annotation (Placement(transformation(extent={{-70,-10},{-50,10}}))); + NDTable nearest( + nin=1, + readFromFile=false, + data=data, + interpMethod=SDF.Types.InterpolationMethod.Nearest, + extrapMethod=SDF.Types.ExtrapolationMethod.Linear, + filename=filename, + dataset=dataset) + annotation (Placement(transformation(extent={{0,40},{20,60}}))); - NDTable hold( + NDTable linear( nin=1, readFromFile=false, data=data, - interpMethod=SDF.Types.InterpolationMethod.Hold, - extrapMethod=SDF.Types.ExtrapolationMethod.Hold, + interpMethod=SDF.Types.InterpolationMethod.Linear, + extrapMethod=SDF.Types.ExtrapolationMethod.Linear, filename=filename, dataset=dataset) - annotation (Placement(transformation(extent={{-10,90},{10,110}}))); + annotation (Placement(transformation(extent={{40,40},{60,60}}))); NDTable akima( nin=1, @@ -43,56 +53,75 @@ model InterpolationMethods "Interpolate a sampled signal using different interpo extrapMethod=SDF.Types.ExtrapolationMethod.Linear, filename=filename, dataset=dataset) - annotation (Placement(transformation(extent={{-10,-30},{10,-10}}))); + annotation (Placement(transformation(extent={{-40,0},{-20,20}}))); - NDTable nearest( + NDTable fritsch_butland( nin=1, readFromFile=false, data=data, - interpMethod=SDF.Types.InterpolationMethod.Nearest, - extrapMethod=SDF.Types.ExtrapolationMethod.Hold, + interpMethod=SDF.Types.InterpolationMethod.FritschButland, + extrapMethod=SDF.Types.ExtrapolationMethod.Linear, filename=filename, dataset=dataset) - annotation (Placement(transformation(extent={{-10,50},{10,70}}))); + annotation (Placement(transformation(extent={{0,0},{20,20}}))); - NDTable fritsch_butland( + NDTable steffen( nin=1, readFromFile=false, data=data, + interpMethod=SDF.Types.InterpolationMethod.Steffen, extrapMethod=SDF.Types.ExtrapolationMethod.Linear, - interpMethod=SDF.Types.InterpolationMethod.FritschButland, filename=filename, dataset=dataset) - annotation (Placement(transformation(extent={{-10,-70},{10,-50}}))); + annotation (Placement(transformation(extent={{40,0},{60,20}}))); - NDTable steffen( + NDTable extrap_hold( + nin=1, + readFromFile=false, + data=data, + interpMethod=SDF.Types.InterpolationMethod.Linear, + extrapMethod=SDF.Types.ExtrapolationMethod.Hold, + filename=filename, + dataset=dataset) + annotation (Placement(transformation(extent={{-20,-60},{0,-40}}))); + + NDTable extrap_linear( nin=1, readFromFile=false, data=data, + interpMethod=SDF.Types.InterpolationMethod.Linear, extrapMethod=SDF.Types.ExtrapolationMethod.Linear, - interpMethod=SDF.Types.InterpolationMethod.Steffen, filename=filename, dataset=dataset) - annotation (Placement(transformation(extent={{-10,-110},{10,-90}}))); + annotation (Placement(transformation(extent={{20,-60},{40,-40}}))); equation - connect(ramp.y, hold.u[1]) annotation (Line(points={{-49,0},{-30,0},{-30,100}, - {-12,100}}, - color={0,0,127})); - connect(ramp.y, akima.u[1]) annotation (Line(points={{-49,0},{-30,0},{-30,-20}, - {-12,-20}}, color={0,0,127})); - connect(ramp.y, linear.u[1]) annotation (Line(points={{-49,0},{-30,0},{-30,20}, - {-12,20}}, color={0,0,127})); - connect(ramp.y, nearest.u[1]) annotation (Line(points={{-49,0},{-30,0},{-30, - 60},{-12,60}}, + connect(ramp.y, hold.u[1]) annotation (Line(points={{-59,-20},{-50,-20},{-50, + 50},{-42,50}}, color={0,0,127})); - connect(ramp.y, fritsch_butland.u[1]) annotation (Line(points={{-49,0},{-30,0}, - {-30,-60},{-12,-60}}, color={0,0,127})); - connect(ramp.y, steffen.u[1]) annotation (Line(points={{-49,0},{-30,0},{-30, - -100},{-12,-100}}, - color={0,0,127})); + connect(akima.u[1], ramp.y) annotation (Line(points={{-42,10},{-50,10},{-50, + -20},{-59,-20}}, + color={0,0,127})); + connect(nearest.u[1], ramp.y) annotation (Line(points={{-2,50},{-10,50},{-10, + -20},{-59,-20}}, + color={0,0,127})); + connect(fritsch_butland.u[1], ramp.y) annotation (Line(points={{-2,10},{-10, + 10},{-10,-20},{-59,-20}}, + color={0,0,127})); + connect(linear.u[1], ramp.y) annotation (Line(points={{38,50},{30,50},{30,-20}, + {-59,-20}}, + color={0,0,127})); + connect(steffen.u[1], ramp.y) annotation (Line(points={{38,10},{30,10},{30, + -20},{-59,-20}}, + color={0,0,127})); + connect(extrap_hold.u[1], ramp.y) annotation (Line(points={{-22,-50},{-30,-50}, + {-30,-20},{-59,-20}}, + color={0,0,127})); + connect(extrap_linear.u[1], ramp.y) annotation (Line(points={{18,-50},{10,-50}, + {10,-20},{-59,-20}}, + color={0,0,127})); annotation ( - Diagram(coordinateSystem(preserveAspectRatio=false, extent={{-100,-120},{ - 100,120}})), - Icon(coordinateSystem(extent={{-100,-80},{100,100}}))); + Diagram(coordinateSystem(preserveAspectRatio=false, extent={{-100,-100},{ + 100,100}})), + Icon(coordinateSystem(extent={{-100,-100},{100,120}}))); end InterpolationMethods; diff --git a/SDF/NDTable.mo b/SDF/NDTable.mo index e815d7b..e637a11 100644 --- a/SDF/NDTable.mo +++ b/SDF/NDTable.mo @@ -12,10 +12,22 @@ parameter SDF.Types.InterpolationMethod interpMethod=SDF.Types.InterpolationMeth "Interpolation method"; parameter SDF.Types.ExtrapolationMethod extrapMethod=SDF.Types.ExtrapolationMethod.None "Extrapolation method"; +parameter Boolean enableDy_du=false "Enable partial derivative output dy_du" + annotation(Evaluate=true); - parameter Real data[:] = { 0} "Table data (as returned by readTableData())" annotation(Dialog(enable=not readFromFile), Evaluate=true); +parameter Real data[:] = { 0} "Table data (as returned by readTableData())" annotation(Dialog(enable=not readFromFile), Evaluate=true); + +Modelica.Blocks.Interfaces.RealOutput dy_du[nin] if enableDy_du + "Partial derivatives of y with respect to each input u[i]" + annotation (Placement( + transformation(extent={{-10,-10},{10,10}}, + rotation=0, + origin={110,-60}), + visible=enableDy_du)); protected + Modelica.Blocks.Interfaces.RealInput local_dy_du[nin] annotation(HideResult=true); + function evaluate input SDF.Types.ExternalNDTable table; input Real[:] params; @@ -25,8 +37,48 @@ protected external "C" value = ModelicaNDTable_evaluate(table, size(params, 1), params, interpMethod, extrapMethod) annotation ( Include="#include ", IncludeDirectory="modelica://SDF/Resources/C-Sources"); + annotation( + derivative(zeroDerivative = table, + zeroDerivative = interpMethod, + zeroDerivative = extrapMethod) = evaluate_der); end evaluate; + function evaluate_der + input SDF.Types.ExternalNDTable table; + input Real[:] params; + input SDF.Types.InterpolationMethod interpMethod; + input SDF.Types.ExtrapolationMethod extrapMethod; + input Real[size(params, 1)] der_params; + output Real der_value; + external "C" der_value = ModelicaNDTable_evaluate_derivative( + table, + size(params, 1), + params, + interpMethod, + extrapMethod, + der_params) + annotation ( + Include="#include ", + IncludeDirectory="modelica://SDF/Resources/C-Sources"); + end evaluate_der; + + function evaluatePartialDerivatives + input SDF.Types.ExternalNDTable table; + input Real[:] params; + input SDF.Types.InterpolationMethod interpMethod; + input SDF.Types.ExtrapolationMethod extrapMethod; + output Real derivatives[size(params, 1)]; + external "C" ModelicaNDTable_evaluate_partial_derivatives( + table, + size(params, 1), + params, + interpMethod, + extrapMethod, + derivatives) annotation ( + Include="#include ", + IncludeDirectory="modelica://SDF/Resources/C-Sources"); + end evaluatePartialDerivatives; + SDF.Types.ExternalNDTable externalTable=SDF.Types.ExternalNDTable(nin, if readFromFile then SDF.Functions.readTableData( Modelica.Utilities.Files.loadResource(filename), dataset, @@ -34,12 +86,24 @@ protected scaleUnits) else data); equation - y = evaluate( + y = evaluate( externalTable, u, interpMethod, extrapMethod); + if enableDy_du then + connect(dy_du, local_dy_du); + local_dy_du = evaluatePartialDerivatives( + externalTable, + u, + interpMethod, + extrapMethod); + else + local_dy_du = zeros(nin); + end if; + + annotation (Documentation(info="

The NDTable block is a multi-dimensional lookup-table (up to 32 dimensions) that supports various inter- and extrapolation methods.

diff --git a/SDF/Resources/C-Sources/ModelicaNDTable.c b/SDF/Resources/C-Sources/ModelicaNDTable.c index ed204ac..b0e42c6 100644 --- a/SDF/Resources/C-Sources/ModelicaNDTable.c +++ b/SDF/Resources/C-Sources/ModelicaNDTable.c @@ -5,10 +5,20 @@ #include "ModelicaUtilities.h" #include "NDTable.h" -#include "NDTable.c" +#include "NDTableCore.c" #define NDTABLE_INTERPSTATUS_OK 0 -#include "Interpolation.c" +#include "NDTableInterpolation.c" + +#define MODELICA_ASSERT_NOT_NULL(argument) \ + if ((argument) == NULL) { \ + ModelicaFormatError("Argument '%s' must not be NULL", #argument); \ + } + +#define MODELICA_ASSERT_INT_ARG(actual, expected) \ + if ((actual) != (expected)) { \ + ModelicaFormatError("Argument %s=%d must be equal to %s=%d", #actual, actual, #expected, expected); \ + } NDTable_h ModelicaNDTable_open(const int ndims, const double *data, const int size) { @@ -92,14 +102,19 @@ void ModelicaNDTable_close(NDTable_h externalTable) { double ModelicaNDTable_evaluate( NDTable_h table, - int nparams, - const double params[], + int nparams, + const double params[], NDTable_InterpMethod_t interp_method, - NDTable_ExtrapMethod_t extrap_method) { + NDTable_ExtrapMethod_t extrap_method +) { + MODELICA_ASSERT_NOT_NULL(table); + MODELICA_ASSERT_INT_ARG(nparams, table->ndims); + MODELICA_ASSERT_NOT_NULL(params); - double value; + double value = 0.0; + double derivative = 0.0; - if (NDTable_evaluate(table, nparams, params, interp_method, extrap_method, &value)) { + if (NDTable_evaluate(table, params, 0, interp_method, extrap_method, &value, &derivative)) { ModelicaError(NDTable_get_error_message()); } @@ -108,19 +123,64 @@ double ModelicaNDTable_evaluate( double ModelicaNDTable_evaluate_derivative( NDTable_h table, - int nparams, + int nparams, const double params[], NDTable_InterpMethod_t interp_method, - NDTable_ExtrapMethod_t extrap_method, + NDTable_ExtrapMethod_t extrap_method, const double delta_params[]) { - double value; + MODELICA_ASSERT_NOT_NULL(table); + MODELICA_ASSERT_INT_ARG(nparams, table->ndims); + MODELICA_ASSERT_NOT_NULL(params); + MODELICA_ASSERT_NOT_NULL(delta_params); - if(NDTable_evaluate_derivative(table, nparams, params, delta_params, interp_method, extrap_method, &value)) { - ModelicaError(NDTable_get_error_message()); + double total_derivative = 0.0; + + for (int i = 0; i < nparams; i++) { + + if (delta_params[i] == 0.0) { + continue; + } + + double value = 0.0; + double derivative = 0.0; + + if (NDTable_evaluate(table, params, i, interp_method, extrap_method, &value, &derivative)) { + ModelicaError(NDTable_get_error_message()); + } + + total_derivative += derivative * delta_params[i]; } - return value; + return total_derivative; } +void ModelicaNDTable_evaluate_partial_derivatives( + NDTable_h table, + int nparams, + const double params[], + NDTable_InterpMethod_t interp_method, + NDTable_ExtrapMethod_t extrap_method, + double derivatives[] +) { + MODELICA_ASSERT_NOT_NULL(table); + MODELICA_ASSERT_INT_ARG(nparams, table->ndims); + MODELICA_ASSERT_NOT_NULL(params); + MODELICA_ASSERT_NOT_NULL(derivatives); + + for (int i = 0; i < nparams; i++) { + double value = 0.0; + double derivative = 0.0; + + if (NDTable_evaluate(table, params, i, interp_method, extrap_method, &value, &derivative)) { + ModelicaError(NDTable_get_error_message()); + } + + derivatives[i] = derivative; + } +} + +#undef MODELICA_ASSERT_NOT_NULL +#undef MODELICA_ASSERT_INT_ARG + #endif // MODELICA_NDTABLE_C diff --git a/SDF/Resources/C-Sources/NDTable.h b/SDF/Resources/C-Sources/NDTable.h index f0574b2..3497931 100644 --- a/SDF/Resources/C-Sources/NDTable.h +++ b/SDF/Resources/C-Sources/NDTable.h @@ -55,30 +55,25 @@ 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 - * + * @param [in] table table handle + * @param [in] params sample point + * @param [in] start_dim start dimension point + * @param [in] interp_method interpolation method (all dimensions) + * @param [in] extrap_method extrapolation method (all dimensions) + * @param [out] value value at the sample point + * @param [out] derivative derivative w.r.t start dimension + * * @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); +int NDTable_evaluate( + NDTable_h table, + const double params[], + int start_dim, + NDTable_InterpMethod_t interp_method, + NDTable_ExtrapMethod_t extrap_method, + double *value, + double *derivative +); /*! The maximum length of an error message */ #define MAX_MESSAGE_LENGTH 256 @@ -126,9 +121,20 @@ double NDTable_get_value_subs(const NDTable_h table, const int subs[]); * * @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); +void NDTable_find_index(double value, int num_values, const double values[], int *index, double *t); + +int NDTable_evaluate_internal( + const NDTable_h table, + const double *t, + const int *subs, + int *nsubs, + int dim, + int dim0, + 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); diff --git a/SDF/Resources/C-Sources/NDTable.c b/SDF/Resources/C-Sources/NDTableCore.c similarity index 71% rename from SDF/Resources/C-Sources/NDTable.c rename to SDF/Resources/C-Sources/NDTableCore.c index 7b81b2d..d74a836 100644 --- a/SDF/Resources/C-Sources/NDTable.c +++ b/SDF/Resources/C-Sources/NDTableCore.c @@ -15,18 +15,6 @@ static char error_message[MAX_MESSAGE_LENGTH] = ""; - -/* A string comparison that returns 1 if a is NULL, empty ("") or equal b and 0 otherwise */ -static int streq(const char *a, const char *b) { - if(a == NULL || strlen(a) == 0) - return 1; - - if(b == NULL || strlen(b) == 0) - return 0; - - return strcmp(a, b) == 0; -} - void NDTable_set_error_message(const char *msg, ...) { va_list vargs; va_start(vargs, msg); @@ -56,17 +44,17 @@ const char * NDTable_get_error_message() { return error_message; } -void NDTable_calculate_offsets(int ndims, const int dims[], int *offs) { +void NDTable_calculate_offsets(int ndims, const int dims[], int offs[]) { int i; if(ndims < 1) { return; } - offs[ndims-1] = 1; + offs[ndims - 1] = 1; for(i = ndims-2; i >= 0; i--) { - offs[i] = offs[i+1] * dims[i+1]; + offs[i] = offs[i + 1] * dims[i + 1]; } } @@ -76,67 +64,67 @@ int NDTable_validate_table(NDTable_h table) { // 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); + NDTable_set_error_message("The rank of must be in the range [0;32] but was %d", 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]); + NDTable_set_error_message("Extent of dimension %d of must be >=0 but was %d", i, 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"); + NDTable_set_error_message("The size of does not match its extent"); 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]); + NDTable_set_error_message("The offset[%d] of must be %d but was %d", i, 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"); + NDTable_set_error_message("Scale for dimension %d of is not set", i); 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); + NDTable_set_error_message("Scale for dimension %d of is not strictly monotonic increasing at index %d", i, 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]); + NDTable_set_error_message("The offset[%d] of must be %d but was %d", i, 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"); + NDTable_set_error_message("The data value at index %d of is not finite", i); return -1; } } @@ -154,7 +142,7 @@ int NDTable_calculate_numel(int ndims, const int dims[]) { return numel; } -void NDTable_ind2sub(const int index, const NDTable_h table, int *subs) { +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++) { diff --git a/SDF/Resources/C-Sources/Interpolation.c b/SDF/Resources/C-Sources/NDTableInterpolation.c similarity index 53% rename from SDF/Resources/C-Sources/Interpolation.c rename to SDF/Resources/C-Sources/NDTableInterpolation.c index b547e17..b98492e 100644 --- a/SDF/Resources/C-Sources/Interpolation.c +++ b/SDF/Resources/C-Sources/NDTableInterpolation.c @@ -6,6 +6,14 @@ #include "NDTable.h" +#define NDTABLE_ASSERT_NOT_NULL(argument) \ + do { \ + if ((argument) == NULL) { \ + NDTable_set_error_message("Argument '%s' must not be NULL", #argument); \ + return -1; \ + } \ + } while (0) + #ifndef MAX #define MAX(a,b) (((a) > (b)) ? (a) : (b)) #endif @@ -25,35 +33,56 @@ static const unsigned long __nan[2] = { 0xffffffff, 0x7fffffff }; #define ISFINITE(x) isfinite(x) #endif +/** + * Calculate the dimension for the current iteration + */ +static int calc_dim(const NDTable_h table, int depth, int start_dim) { + return (depth + start_dim) % table->ndims; +} + /** 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 subs [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 +@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 neighboring sample point +@param depth [in] index of the current dimension +@param start_dim [in] start dimension +@param interp_method [in] interpolation method (all dimensions) +@param extrap_method [in] extrapolation method (all dimensions) +@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[]); +typedef int interp_fun_t( + const NDTable_h table, + const double t[], + const int subs[], + int nsubs[], + int depth, + int start_dim, + NDTable_InterpMethod_t interp_method, + NDTable_ExtrapMethod_t extrap_method, + double *value, + double *derivative +); + +// Interpolation function handle +typedef interp_fun_t *interp_fun_h; // 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) { +static interp_fun_t interp_hold; +static interp_fun_t interp_nearest; +static interp_fun_t interp_linear; +static interp_fun_t interp_akima; +static interp_fun_t interp_fritsch_butland; +static interp_fun_t interp_steffen; +static interp_fun_t extrap_hold; +static interp_fun_t extrap_linear; + +void NDTable_find_index(double value, int nvalues, const double values[], int *index, double *t) { int i; double a, b; double min = values[0]; @@ -83,80 +112,105 @@ void NDTable_find_index(double value, int nvalues, const double *values, int *in *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) { +int NDTable_evaluate( + NDTable_h table, + const double params[], + int start_dim, + NDTable_InterpMethod_t interp_method, + NDTable_ExtrapMethod_t extrap_method, + double *value, + double *derivative +) { + int i; + double t[MAX_NDIMS] = { 0.0 }; // the weights for the interpolation + int subs[MAX_NDIMS] = { 0 }; // the subscripts + int nsubs[MAX_NDIMS] = { 0 }; // the neighboring subscripts + + NDTABLE_ASSERT_NOT_NULL(table); + NDTABLE_ASSERT_NOT_NULL(params); + NDTABLE_ASSERT_NOT_NULL(value); + NDTABLE_ASSERT_NOT_NULL(derivative); + + if (table->ndims == 0) { // dataset is scalar *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); + NDTable_find_index(params[i], table->dims[i], table->scales[i], &subs[i], &t[i]); } - return NDTable_evaluate_internal(table, t, subs, nsubs, 0, interp_method, extrap_method, value, derivatives); + return NDTable_evaluate_internal(table, t, subs, nsubs, 0, start_dim, interp_method, extrap_method, value, derivative); } -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]; +/*! + * @param [in] table NDTable handle + * @param [in] t normalized interpolation factors + * @param [in] subs subscripts + * @param [in] nsubs next subscripts + * @param [in] depth recursion depth + * @param [in] dim_start start dimension + * @param [in] interp_method interpolation method (all dimensions) + * @param [in] extrap_method extrapolation method (all dimensions) + * @param [out] value interpolated value + * @param [out] derivative the derivative value w.r.t the last dimension + */ +int NDTable_evaluate_internal( + const NDTable_h table, + const double t[], + const int subs[], + int nsubs[], + int depth, + int start_dim, + NDTable_InterpMethod_t interp_method, + NDTable_ExtrapMethod_t extrap_method, + double *value, + double *derivative +) { + interp_fun_h func; + + NDTABLE_ASSERT_NOT_NULL(table); + NDTABLE_ASSERT_NOT_NULL(t); + NDTABLE_ASSERT_NOT_NULL(subs); + NDTABLE_ASSERT_NOT_NULL(nsubs); + NDTABLE_ASSERT_NOT_NULL(value); + NDTABLE_ASSERT_NOT_NULL(derivative); + + if (depth < 0 || depth > table->ndims) { + NDTable_set_error_message("Argument 'depth' must be in the range [0, %d] but was %d", table->ndims, depth); + return -1; } - 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) { + if (start_dim < 0 || start_dim >= table->ndims) { + NDTable_set_error_message("Argument 'start_dim' must be in the range [0, %d] but was %d", table->ndims - 1, depth); return -1; } - if (dim >= table->ndims) { + if (depth == table->ndims) { // final dimension reached *value = NDTable_get_value_subs(table, nsubs); return 0; } + int dim = calc_dim(table, depth, start_dim); + // find the right function: - if(table->dims[dim] < 2) { + if(table->dims[dim] < 2) { // dimension has less than two values func = interp_hold; - } else if (t[dim] < 0.0 || t[dim] > 1.0) { + } else if (t[dim] >= 0.0 && t[dim] <= 1.0) { // t inside [0; 1] + // 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: + NDTable_set_error_message("Unknown interpolation method"); + return NDTABLE_INTERPSTATUS_UNKNOWN_METHOD; + } + } else { // t outside [0; 1] // extrapolate switch (extrap_method) { case NDTABLE_EXTRAP_HOLD: @@ -169,70 +223,69 @@ int NDTable_evaluate_internal(const NDTable_h table, const double *t, const int default: func = extrap_linear; break; } break; - default: + case NDTABLE_EXTRAP_NONE: 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 NDTABLE_INTERPSTATUS_OUTOFBOUNS; + default: + NDTable_set_error_message("Unknown extrapolation method"); + return NDTABLE_INTERPSTATUS_UNKNOWN_METHOD; } } - return (*func)(table, t, subs, nsubs, dim, interp_method, extrap_method, value, derivatives); + return (*func)(table, t, subs, nsubs, depth, start_dim, interp_method, extrap_method, value, derivative); } -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_hold(const NDTable_h table, const double t[], const int subs[], int nsubs[], int depth, int start_dim, NDTable_InterpMethod_t interp_method, NDTable_ExtrapMethod_t extrap_method, double *value, double *derivative) { + int dim = calc_dim(table, depth, start_dim); + + nsubs[dim] = subs[dim]; // take the left sample value + *derivative = 0.0; + + return NDTable_evaluate_internal(table, t, subs, nsubs, dim + 1, 0, interp_method, extrap_method, value, derivative); } -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[]) { +static int interp_nearest(const NDTable_h table, const double t[], const int subs[], int nsubs[], int depth, int start_dim, NDTable_InterpMethod_t interp_method, NDTable_ExtrapMethod_t extrap_method, double *value, double *derivative) { int err; + int dim = calc_dim(table, depth, start_dim); + nsubs[dim] = t[dim] < 0.5 ? subs[dim] : subs[dim] + 1; - der_values[dim] = 0; + *derivative = 0.0; - if ((err = NDTable_evaluate_internal(table, t, subs, nsubs, dim + 1, interp_method, extrap_method, value, der_values)) != 0) { + if ((err = NDTable_evaluate_internal(table, t, subs, nsubs, dim + 1, 0, interp_method, extrap_method, value, derivative)) != 0) { return err; } // if the value is not finite return NAN if (!ISFINITE(*value)) { *value = NAN; - der_values[dim] = NAN; + *derivative = 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[]) { +static int interp_linear(const NDTable_h table, const double t[], const int subs[], int nsubs[], int depth, int start_dim, NDTable_InterpMethod_t interp_method, NDTable_ExtrapMethod_t extrap_method, double *value, double *derivative) { int err; double a, b; + int dim = calc_dim(table, depth, start_dim); + // 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) { + if ((err = NDTable_evaluate_internal(table, t, subs, nsubs, depth + 1, start_dim, interp_method, extrap_method, &a, derivative)) != 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) { + if ((err = NDTable_evaluate_internal(table, t, subs, nsubs, depth + 1, start_dim, interp_method, extrap_method, &b, derivative)) != 0) { return err; } // if any of the values is not finite return NAN if (!ISFINITE(a) || !ISFINITE(b)) { *value = NAN; - der_values[dim] = NAN; + *derivative = NAN; return 0; } @@ -240,7 +293,7 @@ static int interp_linear(const NDTable_h table, const double *t, const int *subs *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]]); + *derivative = (b - a) / (table->scales[dim][subs[dim] + 1] - table->scales[dim][subs[dim]]); return 0; } @@ -268,17 +321,17 @@ static void cubic_hermite_spline(const double x0, const double x1, const double } } -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[]) { +static int interp_akima(const NDTable_h table, const double t[], const int subs[], int nsubs[], int depth, int start_dim, NDTable_InterpMethod_t interp_method, NDTable_ExtrapMethod_t extrap_method, double *value, double *derivative) { - 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; - double v = 0; - + double x[6] = { 0.0 }; + double y[6] = { 0.0 }; + double c[4] = { 0.0 }; // spline coefficients + double d[5] = { 0.0 }; // divided differences + double c2 = 0.0; + double dx = 0.0; + double a = 0.0; + + int dim = calc_dim(table, depth, start_dim); int n = table->dims[dim]; // extent of the current dimension int sub = subs[dim]; // subscript of current dimension int err, i, idx; @@ -290,7 +343,7 @@ static int interp_akima(const NDTable_h table, const double *t, const int *subs, 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) { + if ((err = NDTable_evaluate_internal(table, t, subs, nsubs, dim + 1, 0, interp_method, extrap_method, &y[i], derivative)) != 0) { return err; } } @@ -300,7 +353,7 @@ static int interp_akima(const NDTable_h table, const double *t, const int *subs, for (i = 0; i < 6; i++) { if (!ISFINITE(y[i])) { *value = NAN; - der_values[dim] = NAN; + *derivative = NAN; return 0; } } @@ -354,21 +407,21 @@ static int interp_akima(const NDTable_h table, const double *t, const int *subs, c[3] = y[2]; - cubic_hermite_spline(x[2], x[3], y[2], y[3], t[dim], c, value, &der_values[dim]); + cubic_hermite_spline(x[2], x[3], y[2], y[3], t[dim], c, value, derivative); 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[]) { +static int interp_fritsch_butland(const NDTable_h table, const double t[], const int subs[], int nsubs[], int depth, int start_dim, NDTable_InterpMethod_t interp_method, NDTable_ExtrapMethod_t extrap_method, double *value, double *derivative) { - 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; - double v = 0; - + double x [4] = { 0.0 }; + double y [4] = { 0.0 }; + double dx[3] = { 0.0 }; + double d [3] = { 0.0 }; // divided differences + double c [4] = { 0.0 }; // spline coefficients + double c2 = 0.0; + + int dim = calc_dim(table, depth, start_dim); int n = table->dims[dim]; // extent of the current dimension int sub = subs[dim]; // subscript of current dimension int err, i, idx; @@ -380,7 +433,7 @@ static int interp_fritsch_butland(const NDTable_h table, const double *t, const 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) { + if ((err = NDTable_evaluate_internal(table, t, subs, nsubs, dim + 1, 0, interp_method, extrap_method, &y[i], derivative)) != 0) { return err; } } @@ -390,7 +443,7 @@ static int interp_fritsch_butland(const NDTable_h table, const double *t, const for (i = 0; i < 4; i++) { if (!ISFINITE(y[i])) { *value = NAN; - der_values[dim] = NAN; + derivative[dim] = NAN; return 0; } } @@ -429,21 +482,21 @@ static int interp_fritsch_butland(const NDTable_h table, const double *t, const c[3] = y[1]; - cubic_hermite_spline(x[1], x[2], y[1], y[2], t[dim], c, value, &der_values[dim]); + cubic_hermite_spline(x[1], x[2], y[1], y[2], t[dim], c, value, derivative); 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[]) { +static int interp_steffen(const NDTable_h table, const double t[], const int subs[], int nsubs[], int depth, int start_dim, NDTable_InterpMethod_t interp_method, NDTable_ExtrapMethod_t extrap_method, double *value, double *derivative) { - 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 x [4] = { 0.0 }; + double y [4] = { 0.0 }; + double dx[3] = { 0 }; + double d [3] = { 0 }; // divided differences + double c [4] = { 0 }; // spline coefficients double c2 = 0; - double v = 0; + int dim = calc_dim(table, depth, start_dim); const int n = table->dims[dim]; // extent of the current dimension const int sub = subs[dim]; // subscript of current dimension int err, i, idx; @@ -455,7 +508,7 @@ static int interp_steffen(const NDTable_h table, const double *t, const int *sub 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) { + if ((err = NDTable_evaluate_internal(table, t, subs, nsubs, dim + 1, 0, interp_method, extrap_method, &y[i], derivative)) != 0) { return err; } } @@ -465,7 +518,7 @@ static int interp_steffen(const NDTable_h table, const double *t, const int *sub for (i = 0; i < 4; i++) { if (!ISFINITE(y[i])) { *value = NAN; - der_values[dim] = NAN; + *derivative = NAN; return 0; } } @@ -515,47 +568,49 @@ static int interp_steffen(const NDTable_h table, const double *t, const int *sub 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]); + cubic_hermite_spline(x[1], x[2], y[1], y[2], t[dim], c, value, derivative); 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[]) { +static int extrap_hold(const NDTable_h table, const double t[], const int subs[], int nsubs[], int depth, int start_dim, NDTable_InterpMethod_t interp_method, NDTable_ExtrapMethod_t extrap_method, double *value, double *derivative) { int err; + int dim = calc_dim(table, depth, start_dim); nsubs[dim] = t[dim] < 0.0 ? subs[dim] : subs[dim] + 1; - der_values[dim] = 0; + *derivative = 0; - if ((err = NDTable_evaluate_internal(table, t, subs, nsubs, dim + 1, interp_method, extrap_method, value, der_values)) != 0) { + if ((err = NDTable_evaluate_internal(table, t, subs, nsubs, dim + 1, 0, interp_method, extrap_method, value, derivative)) != 0) { return err; } // if the value is not finite return NAN if (!ISFINITE(*value)) { *value = NAN; - der_values[dim] = NAN; + *derivative = 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[]) { +static int extrap_linear(const NDTable_h table, const double t[], const int subs[], int nsubs[], int depth, int start_dim, NDTable_InterpMethod_t interp_method, NDTable_ExtrapMethod_t extrap_method, double *value, double *derivative) { int err; double a, b; - + int dim = calc_dim(table, depth, start_dim); nsubs[dim] = subs[dim]; - if ((err = NDTable_evaluate_internal(table, t, subs, nsubs, dim + 1, interp_method, extrap_method, &a, der_values)) != 0) { + + if ((err = NDTable_evaluate_internal(table, t, subs, nsubs, dim + 1, 0, interp_method, extrap_method, &a, derivative)) != 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) { + if ((err = NDTable_evaluate_internal(table, t, subs, nsubs, dim + 1, 0, interp_method, extrap_method, &b, derivative)) != 0) { return err; } // if any of the values is not finite return NAN if (!ISFINITE(a) || !ISFINITE(b)) { *value = NAN; - der_values[dim] = NAN; + *derivative = NAN; return 0; } @@ -563,7 +618,7 @@ static int extrap_linear(const NDTable_h table, const double *t, const int *subs *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]]); + *derivative = (b - a) / (table->scales[dim][subs[dim] + 1] - table->scales[dim][subs[dim]]); return 0; } diff --git a/ThirdParty/ModelicaUtilities.c b/ThirdParty/ModelicaUtilities.c new file mode 100644 index 0000000..916b770 --- /dev/null +++ b/ThirdParty/ModelicaUtilities.c @@ -0,0 +1,36 @@ +#include "ModelicaUtilities.h" +#include +#include + +//void ModelicaMessage(const char *string) {} +// +//void ModelicaFormatMessage(const char *format, ...) {} +// +//void ModelicaVFormatMessage(const char *format, va_list args) {} + +void ModelicaError(const char *string) { + printf(string); +} + +//void ModelicaWarning(const char *string) {} +// +//void ModelicaFormatWarning(const char *format, ...) MODELICA_FORMATATTR_PRINTF; +// +//void ModelicaVFormatWarning(const char *format, va_list args) MODELICA_FORMATATTR_VPRINTF; + +void ModelicaFormatError(const char* format, ...) { + va_list args; + va_start(args, format); + vprintf(format, args); + va_end(args); +} + +//MODELICA_NORETURN void ModelicaVFormatError(const char *format, va_list args) MODELICA_FORMATATTR_VPRINTF; +// +//char* ModelicaAllocateString(size_t len); +// +//char* ModelicaAllocateStringWithErrorReturn(size_t len); +// +//char* ModelicaDuplicateString(const char* str); +// +//char* ModelicaDuplicateStringWithErrorReturn(const char* str); diff --git a/ThirdParty/ModelicaUtilities.h b/ThirdParty/ModelicaUtilities.h new file mode 100644 index 0000000..e966529 --- /dev/null +++ b/ThirdParty/ModelicaUtilities.h @@ -0,0 +1,189 @@ +/* ModelicaUtilities.h - External utility functions header + + Copyright (C) 2010-2025, Modelica Association and contributors + All rights reserved. + + Redistribution and use in source and binary forms, with or without + modification, are permitted provided that the following conditions are met: + + 1. Redistributions of source code must retain the above copyright notice, + this list of conditions and the following disclaimer. + + 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. + + 3. 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. + + THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND + ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED + WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE + DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE + FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL + DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR + SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER + CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, + OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE + OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. +*/ + +/* Utility functions which can be called by external Modelica functions. + + These functions are defined in section 12.9.6 of the Modelica Specification 3.6. + + A generic C-implementation of these functions cannot be given, + because it is tool dependent how strings are output in a + console or window of the respective simulation tool. +*/ + +#ifndef MODELICA_UTILITIES_H +#define MODELICA_UTILITIES_H + +#include +#include + +#if defined(__cplusplus) +extern "C" { +#endif + +/* + Some of the functions never return to the caller. In order to compile + external Modelica C-code in most compilers, noreturn attributes need to + be present to avoid warnings or errors. + + The following macro handles the noreturn attribute according to the + C23/C11/C++11 standard with fallback to the MSVC/Borland extension. +*/ +#undef MODELICA_NORETURN +#if (defined(__STDC_VERSION__) && __STDC_VERSION__ >= 202311L) || \ + (defined(__cplusplus) && __cplusplus >= 201103L && \ + (!defined(__GNUC__) || __GNUC__ >= 5 || (__GNUC__ == 4 && defined(__GNUC_MINOR__) && __GNUC_MINOR__ >= 8))) +#define MODELICA_NORETURN [[noreturn]] +#elif defined(__STDC_VERSION__) && __STDC_VERSION__ >= 201112L +#define MODELICA_NORETURN _Noreturn +#elif (defined(_MSC_VER) && _MSC_VER >= 1200) || defined(__BORLANDC__) +#define MODELICA_NORETURN __declspec(noreturn) +#else +#define MODELICA_NORETURN +#endif + +/* + The following macros handle format attributes for type-checks against a + format string. +*/ + +#if defined(__clang__) +/* Encapsulated for Clang since GCC fails to process __has_attribute */ +#if __has_attribute(format) +#define MODELICA_FORMATATTR_PRINTF __attribute__((format(printf, 1, 2))) +#define MODELICA_FORMATATTR_VPRINTF __attribute__((format(printf, 1, 0))) +#else +#define MODELICA_FORMATATTR_PRINTF +#define MODELICA_FORMATATTR_VPRINTF +#endif +#elif defined(__GNUC__) && __GNUC__ >= 3 +#define MODELICA_FORMATATTR_PRINTF __attribute__((format(printf, 1, 2))) +#define MODELICA_FORMATATTR_VPRINTF __attribute__((format(printf, 1, 0))) +#else +#define MODELICA_FORMATATTR_PRINTF +#define MODELICA_FORMATATTR_VPRINTF +#endif + +void ModelicaMessage(const char *string); +/* +Output the message string (no format control). +*/ + +void ModelicaFormatMessage(const char *format, ...) MODELICA_FORMATATTR_PRINTF; +/* +Output the message under the same format control as the C-function printf. +*/ + +void ModelicaVFormatMessage(const char *format, va_list args) MODELICA_FORMATATTR_VPRINTF; +/* +Output the message under the same format control as the C-function vprintf. +*/ + +MODELICA_NORETURN void ModelicaError(const char *string); +/* +Output the error message string (no format control). This function +never returns to the calling function, but handles the error +similarly to an assert in the Modelica code. +*/ + +void ModelicaWarning(const char *string); +/* +Output the warning message string (no format control). +*/ + +void ModelicaFormatWarning(const char *format, ...) MODELICA_FORMATATTR_PRINTF; +/* +Output the warning message under the same format control as the C-function printf. +*/ + +void ModelicaVFormatWarning(const char *format, va_list args) MODELICA_FORMATATTR_VPRINTF; +/* +Output the warning message under the same format control as the C-function vprintf. +*/ + +MODELICA_NORETURN void ModelicaFormatError(const char *format, ...) MODELICA_FORMATATTR_PRINTF; +/* +Output the error message under the same format control as the C-function +printf. This function never returns to the calling function, +but handles the error similarly to an assert in the Modelica code. +*/ + +MODELICA_NORETURN void ModelicaVFormatError(const char *format, va_list args) MODELICA_FORMATATTR_VPRINTF; +/* +Output the error message under the same format control as the C-function +vprintf. This function never returns to the calling function, +but handles the error similarly to an assert in the Modelica code. +*/ + +char* ModelicaAllocateString(size_t len); +/* +Allocate memory for a Modelica string which is used as return +argument of an external Modelica function. Note, that the storage +for string arrays (= pointer to string array) is still provided by the +calling program, as for any other array. If an error occurs, this +function does not return, but calls "ModelicaError". +*/ + +char* ModelicaAllocateStringWithErrorReturn(size_t len); +/* +Same as ModelicaAllocateString, except that in case of error, the +function returns 0. This allows the external function to close files +and free other open resources in case of error. After cleaning up +resources use ModelicaError or ModelicaFormatError to signal +the error. +*/ + +char* ModelicaDuplicateString(const char* str); +/* +Duplicate (= allocate memory and deep copy) a Modelica string which +is used as return argument of an external Modelica function. Note, +that the storage for string arrays (= pointer to string array) is still +provided by the calling program, as for any other array. If an error +occurs, this function does not return, but calls "ModelicaError". +*/ + +char* ModelicaDuplicateStringWithErrorReturn(const char* str); +/* +Same as ModelicaDuplicateString, except that in case of error, the +function returns 0. This allows the external function to close files +and free other open resources in case of error. After cleaning up +resources use, ModelicaError or ModelicaFormatError to signal +the error. +*/ + +#undef MODELICA_NORETURN +#undef MODELICA_FORMATATTR_PRINTF +#undef MODELICA_FORMATATTR_VPRINTF + +#if defined(__cplusplus) +} +#endif + +#endif