OMCompiler/SimulationRuntime/c/moo/hessian_finite_diff.h
| Line | Branch | Exec | Source |
|---|---|---|---|
| 1 | /* | ||
| 2 | * This file belongs to the OpenModelica Run-Time System | ||
| 3 | * | ||
| 4 | * Copyright (c) 1998-2026, Open Source Modelica Consortium (OSMC), c/o Linköpings | ||
| 5 | * universitet, Department of Computer and Information Science, SE-58183 Linköping, Sweden. All rights | ||
| 6 | * reserved. | ||
| 7 | * | ||
| 8 | * THIS PROGRAM IS PROVIDED UNDER THE TERMS OF THE BSD NEW LICENSE OR THE | ||
| 9 | * AGPL VERSION 3 LICENSE OR THE OSMC PUBLIC LICENSE (OSMC-PL) VERSION 1.8. ANY | ||
| 10 | * USE, REPRODUCTION OR DISTRIBUTION OF THIS PROGRAM CONSTITUTES RECIPIENT'S | ||
| 11 | * ACCEPTANCE OF THE BSD NEW LICENSE OR THE OSMC PUBLIC LICENSE OR THE AGPL | ||
| 12 | * VERSION 3, ACCORDING TO RECIPIENTS CHOICE. | ||
| 13 | * | ||
| 14 | * The OpenModelica software and the OSMC (Open Source Modelica Consortium) Public License | ||
| 15 | * (OSMC-PL) are obtained from OSMC, either from the above address, from the URLs: | ||
| 16 | * http://www.openmodelica.org or https://github.com/OpenModelica/ or | ||
| 17 | * http://www.ida.liu.se/projects/OpenModelica, and in the OpenModelica distribution. GNU | ||
| 18 | * AGPL version 3 is obtained from: https://www.gnu.org/licenses/licenses.html#GPL. The BSD NEW | ||
| 19 | * License is obtained from: http://www.opensource.org/licenses/BSD-3-Clause. | ||
| 20 | * | ||
| 21 | * This program is distributed WITHOUT ANY WARRANTY; without even the implied warranty of | ||
| 22 | * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE, EXCEPT AS EXPRESSLY | ||
| 23 | * SET FORTH IN THE BY RECIPIENT SELECTED SUBSIDIARY LICENSE CONDITIONS OF | ||
| 24 | * OSMC-PL. | ||
| 25 | * | ||
| 26 | */ | ||
| 27 | |||
| 28 | #ifndef MOO_OM_HESSIAN_FINITE_DIFFERENCES_H | ||
| 29 | #define MOO_OM_HESSIAN_FINITE_DIFFERENCES_H | ||
| 30 | |||
| 31 | // TODO: maybe split this into the C++ part with pattern generation and the C part of evaluation and structs | ||
| 32 | |||
| 33 | #include <map> | ||
| 34 | |||
| 35 | #include "simulation_data.h" | ||
| 36 | |||
| 37 | typedef struct { | ||
| 38 | int i; // first variable index | ||
| 39 | int j; // second variable index | ||
| 40 | } VarPair; | ||
| 41 | |||
| 42 | /* Maps a (color1, color2) pair to all (i, j) variable pairs sharing these colors. | ||
| 43 | * For each (i, j), stores the list of function rows f where both ∂f/∂xi and ∂f/∂xj are nonzero (overestimate). | ||
| 44 | * Also stores the flat COO nz index for (i, j). */ | ||
| 45 | typedef struct { | ||
| 46 | // actual variable pairs, i.e. varPair[k] == (v1, v2); can also be accessed via HESSIAN->(row, col)[lnnzIndices[k]] | ||
| 47 | VarPair* varPairs; // is the variable pair contributing to the functions contributingRows[k] | ||
| 48 | int** contributingRows; // contributingRows[k] = functions affecting the varPair[k] | ||
| 49 | int* numContributingRows; // number of rows for each pair | ||
| 50 | int* lnnzIndices; // mapping from variable pair to Hessian COO index | ||
| 51 | int size; // number of variable pairs in this color group | ||
| 52 | } ColorPair; | ||
| 53 | |||
| 54 | /* Holds the compressed Hessian structure derived from a Jacobian. | ||
| 55 | * COO format row/col lists lower-triangular nonzeros (∂²G/∂xi∂xj). | ||
| 56 | * Variable pairs are grouped by (color1, color2) inside ColorPair blocks. */ | ||
| 57 | typedef struct { | ||
| 58 | /* this is an array of ptrs to ColorPair, is NULL if (c1, c2) is not contained */ | ||
| 59 | ColorPair** colorPairs; // get_color_pair_index(c1, c2) with c1 >= c2 -> variable pairs for color pair | ||
| 60 | int* row; // flat COO row indices (i) | ||
| 61 | int* col; // flat COO column indices (j) | ||
| 62 | int size; // number of variables (Hessian is size × size) | ||
| 63 | int numFuncs; // number of functions in the Hessian | ||
| 64 | int lnnz; // number of lower triangular nonzeros | ||
| 65 | int** colsForColor; // colsForColor[c] is an array of column indices in color c | ||
| 66 | int* colorSizes; // colorSizes[c] is the number of columns in colorCols[c] | ||
| 67 | int numColors; // number of seed vector colors | ||
| 68 | JACOBIAN* jac; // input Jacobian with sparsity + coloring | ||
| 69 | int** cscJacIndexFromRowColor; // mapping of J[function / row][color] -> index in flat Jacobian CSC buffer | ||
| 70 | modelica_real* ws_oldX; // workspace array to remember old x values and seed vector for JVPs | size = #vars | ||
| 71 | modelica_real* ws_h; // workspace array to remember perturbation for variables during numerical Hessian eval | size = #vars | ||
| 72 | modelica_real** ws_baseJac; // workspace stores all rows x colors of the base Jacobian J(x) | ||
| 73 | } HESSIAN_PATTERN; | ||
| 74 | |||
| 75 | /* always use this if accessing HESSIAN_PATTERN.colorPairs | ||
| 76 | * returns the index of a colorPair (c1, c2) in the HESSIAN_PATTERN.colorPairs */ | ||
| 77 | static inline int get_color_pair_index(int c1, int c2) { | ||
| 78 | ✗ | if (c1 >= c2) return c1 * (c1 + 1) / 2 + c2; | |
| 79 | else return c2 * (c2 + 1) / 2 + c1; | ||
| 80 | } | ||
| 81 | |||
| 82 | static inline void set_seed_vector(int size, const int* cols, modelica_real value, modelica_real* seeds) { | ||
| 83 | ✗ | for (int i = 0; i < size; i++) { seeds[cols[i]] = value; } | |
| 84 | } | ||
| 85 | |||
| 86 | HESSIAN_PATTERN* generate_hessian_pattern(JACOBIAN* jac); | ||
| 87 | |||
| 88 | void print_hessian_pattern(const HESSIAN_PATTERN* hes_pattern); | ||
| 89 | |||
| 90 | void free_hessian_pattern(HESSIAN_PATTERN* hes_pattern); | ||
| 91 | |||
| 92 | void eval_hessian_fwd_differences(DATA* data, threadData_t* threadData, HESSIAN_PATTERN* hes_pattern, modelica_real h, | ||
| 93 | int* u_indices, const modelica_real* lambda, modelica_real* jac_csc, modelica_real* hes); | ||
| 94 | |||
| 95 | void hessian_fwd_differences_wrapper(void* args, modelica_real h, modelica_real* result); | ||
| 96 | |||
| 97 | // ===== EXTRAPOLATION ===== | ||
| 98 | |||
| 99 | /* generic computation function of the form "result := f(args, h0)" */ | ||
| 100 | typedef void (*computation_fn_t)(void* args, modelica_real h0, modelica_real* result); | ||
| 101 | |||
| 102 | /** | ||
| 103 | * @brief Workspace for Richardson extrapolation. | ||
| 104 | * Stores intermediate results and metadata for extrapolation steps. | ||
| 105 | */ | ||
| 106 | typedef struct { | ||
| 107 | modelica_real** ws_results; | ||
| 108 | int resultSize; | ||
| 109 | int maxSteps; | ||
| 110 | } ExtrapolationData; | ||
| 111 | |||
| 112 | /* Hessian structure for richardson extrapolation scheme */ | ||
| 113 | typedef struct { | ||
| 114 | DATA* data; | ||
| 115 | threadData_t* threadData; | ||
| 116 | HESSIAN_PATTERN* hes_pattern; // Hessian structure | ||
| 117 | int* u_indices; // indices of the inputs in realVars (REMOVE ME!) | ||
| 118 | const modelica_real* lambda; // dual variables for each function | ||
| 119 | modelica_real* jac_csc; // precomputed Jacobian entries in CSC format (can be NULL) | ||
| 120 | } HessianFiniteDiffArgs; | ||
| 121 | |||
| 122 | ExtrapolationData* init_extrapolation_data(int resultSize, int maxSteps); | ||
| 123 | |||
| 124 | void free_extrapolation_data(ExtrapolationData* extrData); | ||
| 125 | |||
| 126 | void richardson_extrapolation(ExtrapolationData* extrData, computation_fn_t fn, void* args, modelica_real h0, | ||
| 127 | int steps, modelica_real stepDivisor, int methodOrder, modelica_real* result); | ||
| 128 | |||
| 129 | #endif // MOO_OM_HESSIAN_FINITE_DIFFERENCES_H | ||
| 130 |