Linux GNU 11.4.0 Code Coverage Report

Directory: ./
Coverage: low: ≥ 0% medium: ≥ 75.0% high: ≥ 90.0%
Coverage Exec / Excl / Total
Lines: 0.0% 0 / 0 / 316
Functions: 0.0% 0 / 0 / 9
Branches: 0.0% 0 / 0 / 248
OMCompiler/SimulationRuntime/c/moo/hessian_finite_diff.cpp
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 #include <stdlib.h>
29 #include <sstream>
30 #include <iomanip>
31 #include <vector>
32
33 #include "simulation/arrayIndex.h"
34
35 #include "hessian_finite_diff.h"
36
37 /**
38 * @brief Constructs a compressed Hessian sparsity pattern C struct using Jacobian coloring.
39 *
40 * Given a sparse Jacobian J(x) of F: R^n → R^m, this builds the structure of the Hessian
41 * of a scalar adjoint G(x) = Σ λ[i]·F[i](x), based on co-occurrence of variables in J(x).
42 *
43 * Variable pairs (i,j) are collected if they appear together in any function row f.
44 * These define the nonzero Hessian structure (i.e. where ∂²G/∂xi∂xj != 0).
45 *
46 * A second coloring is induced: if variable x_i ∈ color c₁ and x_j ∈ color c₂, then (i,j) is assigned to the color pair (c₁, c₂).
47 * This allows evaluating Hessian entries H[i,j] via directional finite differences: perturb x along seed vector s₁ (color c₁),
48 * and apply a Jacobian-vector product with seed vector s₂ (color c₂).
49 *
50 * For each color pair (c₁, c₂), only a subset of Hessian entries is affected. Each entry H[i,j] receives contributions
51 * only from those function rows r where both variables x_i and x_j appear (i.e. where ∂f/∂x_i and ∂f/∂x_j are nonzero).
52 * The directional second derivative is given using:
53 *
54 * H[i,j] = ∑_r λ[r] · ((J(x + h · s_1) - J(x)) · s_2 / h)[r]
55 *
56 * where λ ∈ ℝᵐ is the adjoint vector. This corresponds to evaluating the contraction λᵗ · ∇²F(x) · v without forming the full Hessian.
57 *
58 * The function makes heavy use of the STL, but the result is a `HESSIAN_PATTERN` pure C struct that includes:
59 * - COO row/col index arrays for Hessian nonzeros (lower triangle).
60 * - A lookup from (color₁, color₂) to `ColorPair`, listing variable pairs and contributing rows.
61 *
62 * @param jac [in] Pointer to a `JACOBIAN` struct (sparsity pattern and coloring).
63 * @return [out] Pointer to newly allocated `HESSIAN_PATTERN` struct.
64 */
65 ✗ HESSIAN_PATTERN* generate_hessian_pattern(JACOBIAN* jac) {
66 ✗ if (jac == nullptr || jac->sparsePattern == nullptr) { return nullptr; }
67
68 ✗ int numVars = jac->sizeCols;
69 ✗ int numFuncs = jac->sizeRows;
70 SPARSE_PATTERN* sp = jac->sparsePattern;
71 ✗ int numColors = sp->maxColors;
72
73 // 1. build adjacency list: which variables affect which functions
74 ✗ std::vector<std::vector<int>> adj(numFuncs);
75 ✗ for (int col = 0; col < numVars; col++) {
76 ✗ for (unsigned int nz = sp->leadindex[col]; nz < sp->leadindex[col + 1]; nz++) {
77 ✗ int row = sp->index[nz];
78 ✗ adj[row].push_back(col);
79 }
80 }
81
82 // 2. build M[v1, v2] = list of function rows where both variables appear
83 std::map<std::pair<int, int>, std::vector<int>> M;
84 ✗ for (int f = 0; f < numFuncs; f++) {
85 ✗ const auto& vars = adj[f];
86 ✗ for (size_t i = 0; i < vars.size(); i++) {
87 ✗ for (size_t j = 0; j <= i; j++) {
88 ✗ int v1 = vars[i];
89 ✗ int v2 = vars[j];
90 ✗ if (v1 < v2) std::swap(v1, v2);
91 ✗ M[{v1, v2}].push_back(f);
92 }
93 }
94 }
95
96 // 3. assign flat indices (lower nnz) directly from sorted M keys
97 std::map<std::pair<int, int>, int> cooMap;
98 int lnnz = 0;
99 ✗ for (const auto& [pair, _] : M) {
100 ✗ cooMap[pair] = lnnz++;
101 }
102
103 // 4. build color groups :: TODO: implement this in OpenModelica for the JACOBIAN
104 ✗ std::vector<std::vector<int>> colorCols(numColors);
105 ✗ for (int col = 0; col < numVars; col++) {
106 ✗ int c = sp->colorCols[col];
107 ✗ if (c > 0) {
108 ✗ colorCols[c - 1].push_back(col);
109 }
110 }
111
112 // 5. allocate pattern
113 ✗ HESSIAN_PATTERN* hes_pattern = (HESSIAN_PATTERN*)malloc(sizeof(HESSIAN_PATTERN));
114 ✗ hes_pattern->colorPairs = (ColorPair**)calloc(numColors * (numColors + 1) / 2, sizeof(ColorPair*));
115 ✗ hes_pattern->row = (int*)malloc(lnnz * sizeof(int));
116 ✗ hes_pattern->col = (int*)malloc(lnnz * sizeof(int));
117 ✗ hes_pattern->colsForColor = (int**)malloc(numColors * sizeof(int*));
118 ✗ hes_pattern->colorSizes = (int*)malloc(numColors * sizeof(int));
119 ✗ hes_pattern->numColors = numColors;
120 ✗ hes_pattern->numFuncs = numFuncs;
121 ✗ hes_pattern->size = numVars;
122 ✗ hes_pattern->lnnz = lnnz;
123 ✗ hes_pattern->jac = jac;
124
125 // workspace memory
126 ✗ hes_pattern->ws_oldX = (modelica_real*)malloc(numVars * sizeof(modelica_real));
127 ✗ hes_pattern->ws_h = (modelica_real*)malloc(numVars * sizeof(modelica_real));
128 ✗ hes_pattern->ws_baseJac = (modelica_real**)malloc(numFuncs * sizeof(modelica_real*));
129 ✗ for (int row = 0; row < numFuncs; row++) {
130 ✗ hes_pattern->ws_baseJac[row] = (modelica_real*)calloc(numColors, sizeof(modelica_real));
131 }
132
133 // 6. remember columns in each color
134 ✗ for (int i = 0; i < numColors; i++) {
135 ✗ int size = colorCols[i].size();
136 ✗ hes_pattern->colorSizes[i] = size;
137 ✗ hes_pattern->colsForColor[i] = (int*)malloc(size * sizeof(int));
138 memcpy(hes_pattern->colsForColor[i], colorCols[i].data(), size * sizeof(int));
139 }
140
141 // 7. set mapping from Jacobian[row][color] -> Jacobian CSC index
142 ✗ hes_pattern->cscJacIndexFromRowColor = (int**)malloc(numFuncs * sizeof(int*));
143 ✗ for (int row = 0; row < numFuncs; row++) {
144 ✗ hes_pattern->cscJacIndexFromRowColor[row] = (int*)malloc(numColors * sizeof(int));
145 ✗ for (int color = 0; color < numColors; color++) {
146 ✗ hes_pattern->cscJacIndexFromRowColor[row][color] = -1;
147 }
148 }
149
150 ✗ for (int color = 0; color < numColors; color++) {
151 ✗ const int* cols = hes_pattern->colsForColor[color];
152 ✗ for (int colIdx = 0; colIdx < hes_pattern->colorSizes[color]; colIdx++) {
153 ✗ int col = cols[colIdx];
154 ✗ for (unsigned int nz = sp->leadindex[col]; nz < sp->leadindex[col + 1]; nz++) {
155 ✗ int row = sp->index[nz];
156 ✗ hes_pattern->cscJacIndexFromRowColor[row][color] = nz;
157 }
158 }
159 }
160
161 // 8. fill the coordinate format sparsity
162 ✗ for (const auto& coo : cooMap) {
163 ✗ int var_row = coo.first.first;
164 ✗ int var_col = coo.first.second;
165 ✗ int nz = coo.second;
166
167 ✗ hes_pattern->row[nz] = var_row;
168 ✗ hes_pattern->col[nz] = var_col;
169 }
170
171 ColorPair* colorPair;
172
173 // 9. fill HESSIAN_PATTERN.colorPairs[c1][c2] -> ColorPair
174 ✗ for (int c1 = 0; c1 < numColors; c1++) {
175 ✗ for (int c2 = 0; c2 <= c1; c2++) {
176 std::vector<std::vector<int>> rowsVec;
177 std::vector<int> nnzIndicesVec;
178 std::vector<VarPair> pairVec;
179
180 ✗ for (int i1 : colorCols[c1]) {
181 ✗ for (int i2 : colorCols[c2]) {
182 // copy and swap if needed
183 int v1 = i1;
184 int v2 = i2;
185 ✗ if (v1 < v2){
186 std::swap(v1, v2);
187 }
188
189 ✗ auto v_pair = std::make_pair(v1, v2);
190 auto it = M.find(v_pair);
191 ✗ if (it == M.end()) continue;
192
193 auto cooIt = cooMap.find(v_pair);
194 ✗ if (cooIt == cooMap.end()) continue;
195
196 ✗ rowsVec.push_back(it->second); // function rows
197 ✗ nnzIndicesVec.push_back(cooIt->second); // flat Hessian index, nz index
198 ✗ pairVec.push_back({v1, v2}); // variable pair
199 }
200 }
201
202 // create and allocate ColorPair
203 ✗ int variablePairCount = rowsVec.size();
204 ✗ if (variablePairCount == 0) {
205 colorPair = nullptr;
206 }
207 else {
208 ✗ colorPair = (ColorPair*)malloc(sizeof(ColorPair));
209 ✗ colorPair->size = variablePairCount;
210 ✗ colorPair->contributingRows = (int**)malloc(variablePairCount * sizeof(int*));
211 ✗ colorPair->numContributingRows = (int*)malloc(variablePairCount * sizeof(int));
212 ✗ colorPair->lnnzIndices = (int*)malloc(variablePairCount * sizeof(int));
213 ✗ colorPair->varPairs = (VarPair*)malloc(variablePairCount * sizeof(VarPair));
214
215 ✗ for (int i = 0; i < variablePairCount; i++) {
216 ✗ int sz = rowsVec[i].size();
217 ✗ colorPair->contributingRows[i] = (int*)malloc(sz * sizeof(int));
218 memcpy(colorPair->contributingRows[i], rowsVec[i].data(), sz * sizeof(int));
219 ✗ colorPair->numContributingRows[i] = sz;
220 ✗ colorPair->lnnzIndices[i] = nnzIndicesVec[i];
221 ✗ colorPair->varPairs[i] = pairVec[i];
222 }
223 }
224
225 ✗ hes_pattern->colorPairs[get_color_pair_index(c1, c2)] = colorPair;
226 ✗ }
227 }
228
229 return hes_pattern;
230 ✗ }
231
232 /**
233 * @brief Compute Hessian-vector product λᵗH(x) using forward finite differences of the Jacobian.
234 *
235 * @note I found out that this is very shady and only works well for well-scaled / posed problems, as
236 * h can not be taken differently within a color. Therefore, we choose the geometric mean of the
237 * nominal h as the h for a given color.
238 * => use eval_hessian_fwd_differences_safe for safe version
239 *
240 * Approximates the entries of the Hessian matrix H(x) using first-order directional derivatives.
241 * The method uses seed vector coloring for efficient evaluation and exploits sparse Hessian structure.
242 * Assumes the current point x has all controls and states set in `data->localData[0]->realVars`.
243 * For a more detailed explanation of the algorithm (see generate_hessian_pattern).
244 *
245 * Runtime: O(#colors * (#colors + 1) / 2 * T_{JVP} + #colors * T_{JVP} + #funcs_{avg} * nnz(Hessian)),
246 * where T_{JVP} is the time of one Jacobian column evaluation and funcs_{avg} is the average
247 * number of functions for each variable pair
248 *
249 * Driving term: O(#colors * (#colors + 1) / 2 * T_{JVP}, since #colors * T_{JVP} will be precomputed
250 * for the Jacobian anyway and #funcs_{avg} <= #funcs, thus comparably insignificant
251 * => just (#colors + 1) / 2 times the time for the Jacobian evaluation
252 *
253 * @param[in] data Runtime simulation data structure.
254 * @param[in] threadData Thread-local data.
255 * @param[in] hes_pattern Precomputed sparsity and coloring pattern for Hessian and Jacobian.
256 * @param[in] h Perturbation step size (pre-scaling).
257 * @param[in] lambda Adjoint vector (size = number of functions).
258 * @param[in] u_indices Indices of the input variables. For Optimization, these can be obtained by calling data->callback->getInputVarIndicesInOptimization(). (I hate it that this is an arg; it should be somewhere in DATA or so.)
259 * @param[in] jac_csc (Optional) Jacobian values in CSC format, used to speed up Hessian calculation. NULL -> compute from scratch.
260 * @param[out] hes Output sparse Hessian values (COO format of hes_pattern, length = hes_pattern->nnz).
261 */
262 ✗ void eval_hessian_fwd_differences_fast(
263 DATA* data,
264 threadData_t* threadData,
265 HESSIAN_PATTERN* hes_pattern,
266 modelica_real h,
267 int* u_indices,
268 const modelica_real* lambda,
269 modelica_real* jac_csc,
270 modelica_real* hes)
271 {
272 /* 0. retrieve pointers */
273 ✗ JACOBIAN* jacobian = hes_pattern->jac;
274 ✗ modelica_real** ws_baseJac = hes_pattern->ws_baseJac;
275 ✗ modelica_real* ws_oldX = hes_pattern->ws_oldX;
276 modelica_real* ws_h = hes_pattern->ws_h;
277 ✗ modelica_real* seeds = jacobian->seedVars;
278 ✗ modelica_real* jvp = jacobian->resultVars;
279 ✗ unsigned int* jacLeadIndex = jacobian->sparsePattern->leadindex;
280 ✗ unsigned int* jacIndex = jacobian->sparsePattern->index;
281
282 ✗ int nStates = data->modelData->nStates;
283
284 /* 1. compute standard Jacobian, if jac_csc is NULL, else use the jac_csc as precomputed Jacobian */
285 ✗ if (!jac_csc) {
286 /* 1.a. evaluate base system (needed for Jacobian columns) */
287 ✗ data->callback->functionDAE(data, threadData);
288
289 /* 1.b. evaluate all JVPs J(x) * s_{c} of the current point x */
290 ✗ for (int color = 0; color < hes_pattern->numColors; color++) {
291 ✗ set_seed_vector(hes_pattern->colorSizes[color], hes_pattern->colsForColor[color], 1, seeds);
292 ✗ jacobian->evalColumn(data, threadData, jacobian, NULL);
293
294 ✗ for (int colIndex = 0; colIndex < hes_pattern->colorSizes[color]; colIndex++) {
295 ✗ int col = hes_pattern->colsForColor[color][colIndex];
296 ✗ for (unsigned int nz = jacLeadIndex[col]; nz < jacLeadIndex[col + 1]; nz++) {
297 ✗ int row = jacIndex[nz];
298 ✗ ws_baseJac[row][color] = jvp[row];
299 }
300 }
301
302 ✗ set_seed_vector(hes_pattern->colorSizes[color], hes_pattern->colsForColor[color], 0, seeds);
303 }
304 }
305
306 /* 2. loop over all colors c1 */
307 ✗ for (int c1 = 0; c1 < hes_pattern->numColors; c1++) {
308 /* 3. define seed vector s_{c_1} with all cols in c_1 active (implicitly) */
309 /* 4. peturbate current x_{c_1} := x + h * s_{c_1} */
310 modelica_real c1_h = 0;
311 ✗ for (int columnIndex = 0; columnIndex < hes_pattern->colorSizes[c1]; columnIndex++) {
312 ✗ int col = hes_pattern->colsForColor[c1][columnIndex];
313 ✗ int realVarsIndex = (col < nStates ? col : u_indices[col - nStates]);
314
315 /* create perturbation size based on nominals and current entry */
316 ✗ const modelica_real nom = getNominalFromScalarIdx(data->simulationInfo, data->modelData, VAR_KIND_VARIABLE, realVarsIndex);
317 ✗ c1_h += log(1.0 + fmax(ws_oldX[col], nom));
318 }
319
320 // geometric mean of h
321 ✗ c1_h = h / (1 + hes_pattern->colorSizes[c1]) * exp(c1_h);
322
323 ✗ for (int columnIndex = 0; columnIndex < hes_pattern->colorSizes[c1]; columnIndex++) {
324 ✗ int col = hes_pattern->colsForColor[c1][columnIndex];
325 ✗ int realVarsIndex = (col < nStates ? col : u_indices[col - nStates]);
326 /* remember the current realVars (to be perturbed) and perturbate */
327 ✗ ws_oldX[col] = data->localData[0]->realVars[realVarsIndex];
328 ✗ data->localData[0]->realVars[realVarsIndex] += c1_h;
329 }
330
331 /* evaluate perturbed system (needed for Jacobian columns) */
332 ✗ data->callback->functionDAE(data, threadData);
333
334 /* 5. loop over all colors c2 with index less or equal to c_1 */
335 ✗ for (int c2 = 0; c2 <= c1; c2++) {
336 /* 6. define seed vector s_{c_2} with all cols in c_2 active */
337 ✗ set_seed_vector(hes_pattern->colorSizes[c2], hes_pattern->colsForColor[c2], 1, seeds);
338
339 /* 7. evaluate JVP J(x_{c_1}) * s_{c_2}: writes column to jvp = jacobian->resultVars */
340 ✗ jacobian->evalColumn(data, threadData, jacobian, NULL);
341
342 /* 8. retrieve Hessian approximation */
343 ✗ ColorPair* colorPair = hes_pattern->colorPairs[get_color_pair_index(c1, c2)];
344 ✗ if (colorPair) {
345 ✗ for (int varPairIdx = 0; varPairIdx < colorPair->size; varPairIdx++) {
346 /* nz index in flattened Hessian array (COO format) */
347 ✗ int nz = colorPair->lnnzIndices[varPairIdx];
348
349 /* rows (functions) where both ∂f/∂xi and ∂f/∂xj are nonzero */
350 ✗ int* contributingRows = colorPair->contributingRows[varPairIdx];
351 ✗ int numContributingRows = colorPair->numContributingRows[varPairIdx];
352
353 /* second derivative eval at nz index */
354 modelica_real der = 0.0;
355
356 /* 10. Approximate directional second derivative:
357 * (1/h) ∑_{f ∈ rows} λ[f] · (J(x + h·s_{c₁})[s_{c₂}][f] - J(x)[s_{c₂}][f])
358 * where:
359 * - f / fnRow indexes function rows where both ∂f/∂xᵢ and ∂f/∂xⱼ are nonzero */
360 ✗ for (int fIdx = 0; fIdx < numContributingRows; fIdx++) {
361 ✗ int fnRow = contributingRows[fIdx];
362 ✗ modelica_real J_fnRow_c2 = (jac_csc ? jac_csc[hes_pattern->cscJacIndexFromRowColor[fnRow][c2]] : ws_baseJac[fnRow][c2]);
363 ✗ der += lambda[fnRow] * (jvp[fnRow] - J_fnRow_c2);
364 }
365
366 /* store and divide by step size, retrieve step size via nz col index / same as for the perturbation (step 4) */
367 ✗ int col = hes_pattern->col[nz];
368 ✗ hes[nz] = der / c1_h;
369 }
370 }
371
372 /* 11. reset s_{c_2} */
373 ✗ set_seed_vector(hes_pattern->colorSizes[c2], hes_pattern->colsForColor[c2], 0.0, seeds);
374 }
375
376 /* 12. reset perturbation in x */
377 ✗ for (int columnIndex = 0; columnIndex < hes_pattern->colorSizes[c1]; columnIndex++) {
378 ✗ int col = hes_pattern->colsForColor[c1][columnIndex];
379 ✗ int realVarsIndex = (col < nStates ? col : u_indices[col - nStates]);
380 ✗ data->localData[0]->realVars[realVarsIndex] = ws_oldX[col];
381 }
382 }
383 ✗ }
384
385 /**
386 * @brief Compute Hessian-vector product λᵗH(x) using forward finite differences of the Jacobian.
387 *
388 * This version iterates through colors c1, but inside that loop, it iterates through every
389 * variable belonging to c1 individually. This prevents step-size scaling issues by allowing
390 * a fixed h for each variable.
391 *
392 * Runtime: O(1/2 * #vars * #colors * T_{JVP}) - Significantly slower than the colored version,
393 * but numerically more robust for bad scaling.
394 *
395 * @param[in] data Runtime simulation data structure.
396 * @param[in] threadData Thread-local data.
397 * @param[in] hes_pattern Precomputed sparsity and coloring pattern.
398 * @param[in] h (Ignored in this version, uses fixed 1e-6).
399 * @param[in] lambda Adjoint vector.
400 * @param[in] u_indices Indices of input variables.
401 * @param[in] jac_csc (Optional) Jacobian values in CSC format.
402 * @param[out] hes Output sparse Hessian values.
403 */
404 ✗ void eval_hessian_fwd_differences(
405 DATA* data,
406 threadData_t* threadData,
407 HESSIAN_PATTERN* hes_pattern,
408 modelica_real h,
409 int* u_indices,
410 const modelica_real* lambda,
411 modelica_real* jac_csc,
412 modelica_real* hes)
413 {
414 /* 0. retrieve pointers */
415 ✗ JACOBIAN* jacobian = hes_pattern->jac;
416 ✗ modelica_real** ws_baseJac = hes_pattern->ws_baseJac;
417 ✗ modelica_real* seeds = jacobian->seedVars;
418 ✗ modelica_real* jvp = jacobian->resultVars;
419 ✗ unsigned int* jacLeadIndex = jacobian->sparsePattern->leadindex;
420 ✗ unsigned int* jacIndex = jacobian->sparsePattern->index;
421
422 ✗ int nStates = data->modelData->nStates;
423
424 /* 1. compute standard Jacobian, if jac_csc is NULL */
425 ✗ if (!jac_csc) {
426 /* 1.a. evaluate base system */
427 ✗ data->callback->functionDAE(data, threadData);
428
429 /* 1.b. evaluate all JVPs J(x) * s_{c} */
430 ✗ for (int color = 0; color < hes_pattern->numColors; color++) {
431 ✗ set_seed_vector(hes_pattern->colorSizes[color], hes_pattern->colsForColor[color], 1, seeds);
432 ✗ jacobian->evalColumn(data, threadData, jacobian, NULL);
433
434 ✗ for (int colIndex = 0; colIndex < hes_pattern->colorSizes[color]; colIndex++) {
435 ✗ int col = hes_pattern->colsForColor[color][colIndex];
436 ✗ for (unsigned int nz = jacLeadIndex[col]; nz < jacLeadIndex[col + 1]; nz++) {
437 ✗ int row = jacIndex[nz];
438 ✗ ws_baseJac[row][color] = jvp[row];
439 }
440 }
441 ✗ set_seed_vector(hes_pattern->colorSizes[color], hes_pattern->colsForColor[color], 0, seeds);
442 }
443 }
444
445 /* 2. Loop over all colors c1 */
446 ✗ for (int c1 = 0; c1 < hes_pattern->numColors; c1++) {
447
448 /* 3. Loop over each variable in color c1 individually */
449 ✗ for (int columnIndex = 0; columnIndex < hes_pattern->colorSizes[c1]; columnIndex++) {
450
451 /* Identify the specific variable to perturb */
452 ✗ int col = hes_pattern->colsForColor[c1][columnIndex];
453 ✗ int realVarsIndex = (col < nStates ? col : u_indices[col - nStates]);
454
455 /* 4. Perturb current x_{col} := x + h */
456 ✗ modelica_real oldVal = data->localData[0]->realVars[realVarsIndex];
457 ✗ modelica_real h_col = h * (1 + 1e-5 * std::abs(oldVal));
458 ✗ data->localData[0]->realVars[realVarsIndex] += h_col;
459
460 /* Evaluate perturbed system */
461 ✗ data->callback->functionDAE(data, threadData);
462
463 /* 5. Loop over all colors c2 with index less or equal to c_1 */
464 ✗ for (int c2 = 0; c2 <= c1; c2++) {
465 /* 6. Define seed vector s_{c_2} with all cols in c_2 active */
466 ✗ set_seed_vector(hes_pattern->colorSizes[c2], hes_pattern->colsForColor[c2], 1, seeds);
467
468 /* 7. Evaluate JVP J(x + h*e_{col}) * s_{c_2} */
469 ✗ jacobian->evalColumn(data, threadData, jacobian, NULL);
470
471 /* 8. Retrieve Hessian approximation */
472 ✗ ColorPair* colorPair = hes_pattern->colorPairs[get_color_pair_index(c1, c2)];
473 ✗ if (colorPair) {
474 ✗ for (int varPairIdx = 0; varPairIdx < colorPair->size; varPairIdx++) {
475 /* skip unrelated variables (not perturbated pairs) */
476 ✗ if (!(colorPair->varPairs[varPairIdx].i == col || colorPair->varPairs[varPairIdx].j == col)) continue;
477
478 /* nz index in flattened Hessian array */
479 ✗ int nz = colorPair->lnnzIndices[varPairIdx];
480
481 ✗ int* contributingRows = colorPair->contributingRows[varPairIdx];
482 ✗ int numContributingRows = colorPair->numContributingRows[varPairIdx];
483 modelica_real der = 0.0;
484
485 /* 10. Approximate directional second derivative */
486 ✗ for (int fIdx = 0; fIdx < numContributingRows; fIdx++) {
487 ✗ int fnRow = contributingRows[fIdx];
488 ✗ modelica_real J_fnRow_c2 = (jac_csc ? jac_csc[hes_pattern->cscJacIndexFromRowColor[fnRow][c2]] : ws_baseJac[fnRow][c2]);
489
490 ✗ der += lambda[fnRow] * (jvp[fnRow] - J_fnRow_c2);
491 }
492
493 ✗ hes[nz] = der / h_col;
494 }
495 }
496
497 /* 11. Reset s_{c_2} */
498 ✗ set_seed_vector(hes_pattern->colorSizes[c2], hes_pattern->colsForColor[c2], 0.0, seeds);
499 }
500
501 /* 12. Reset perturbation in x */
502 ✗ data->localData[0]->realVars[realVarsIndex] = oldVal;
503 }
504 }
505 ✗ }
506
507 ✗ void print_hessian_pattern(const HESSIAN_PATTERN* hes_pattern) {
508 ✗ if (!hes_pattern) {
509 ✗ errorStreamPrint(OMC_LOG_MOO, 0, "Hessian pattern is NULL.");
510 ✗ return;
511 }
512
513 ✗ infoStreamPrint(OMC_LOG_MOO, 0, "\n=== HESSIAN SPARSITY INFO ===");
514 ✗ infoStreamPrint(OMC_LOG_MOO, 0, "Matrix size: %d x %d", hes_pattern->size, hes_pattern->size);
515 ✗ infoStreamPrint(OMC_LOG_MOO, 0, "Lower triangle NNZ: %d", hes_pattern->lnnz);
516 ✗ infoStreamPrint(OMC_LOG_MOO, 0, "Number of colors: %d", hes_pattern->numColors);
517
518 ✗ infoStreamPrint(OMC_LOG_MOO, 0, "\nBase Jacobian Colors:");
519 ✗ for (int c = 0; c < hes_pattern->numColors; c++) {
520 ✗ std::ostringstream oss;
521 ✗ oss << " Color " << c << " (size " << hes_pattern->colorSizes[c] << "): ";
522 ✗ for (int j = 0; j < hes_pattern->colorSizes[c]; j++) {
523 ✗ oss << hes_pattern->colsForColor[c][j] << " ";
524 }
525 ✗ infoStreamPrint(OMC_LOG_MOO, 0, "%s", oss.str().c_str());
526 ✗ }
527
528 ✗ infoStreamPrint(OMC_LOG_MOO, 0, "\nCoordinate Format (COO, lower triangle):");
529 ✗ infoStreamPrint(OMC_LOG_MOO, 0, " lnnz | Row | Col");
530 ✗ infoStreamPrint(OMC_LOG_MOO, 0, "------------------");
531 ✗ for (int i = 0; i < hes_pattern->lnnz; i++) {
532 ✗ infoStreamPrint(OMC_LOG_MOO, 0, " %3d | %3d | %3d", i, hes_pattern->row[i], hes_pattern->col[i]);
533 }
534
535 ✗ infoStreamPrint(OMC_LOG_MOO, 0, "\nColor Pair Entries:");
536 ✗ for (int c1 = 0; c1 < hes_pattern->numColors; c1++) {
537 ✗ for (int c2 = 0; c2 <= c1; c2++) { // symmetric lower triangle
538 int idx = get_color_pair_index(c1, c2);
539 ✗ ColorPair* colorPair = hes_pattern->colorPairs[idx];
540 ✗ if (!colorPair) continue;
541
542 ✗ infoStreamPrint(OMC_LOG_MOO, 0, " Color pair (%d, %d): %d variable pairs",
543 c1, c2, colorPair->size);
544
545 ✗ for (int i = 0; i < colorPair->size; i++) {
546 ✗ int nnzIdx = colorPair->lnnzIndices[i];
547
548 ✗ std::ostringstream oss;
549 ✗ oss << " VarPair: (" << hes_pattern->row[nnzIdx]
550 ✗ << ", " << hes_pattern->col[nnzIdx]
551 ✗ << "), nnz_index = " << nnzIdx << ", Functions = [";
552
553 ✗ for (int j = 0; j < colorPair->numContributingRows[i]; j++) {
554 ✗ oss << colorPair->contributingRows[i][j];
555 ✗ if (j + 1 < colorPair->numContributingRows[i]) oss << ", ";
556 }
557 ✗ oss << "]";
558 ✗ infoStreamPrint(OMC_LOG_MOO, 0, "%s", oss.str().c_str());
559 ✗ }
560
561 ✗ infoStreamPrint(OMC_LOG_MOO, 0, "----------------------------------------------------------------------------");
562 }
563 }
564
565 ✗ int n = hes_pattern->size;
566 {
567 ✗ infoStreamPrint(OMC_LOG_MOO, 0, "\n=== HESSIAN SPARSITY PLOT (λᵗ·∇²F) ===");
568
569 ✗ std::ostringstream oss;
570 ✗ oss << " ";
571 ✗ for (int j = 0; j < n; j++) oss << j;
572 ✗ infoStreamPrint(OMC_LOG_MOO, 0, "%s", oss.str().c_str());
573 ✗ }
574
575 ✗ char* sparsity = (char*)calloc(n * n, sizeof(char));
576 ✗ for (int lnz = 0; lnz < hes_pattern->lnnz; lnz++) {
577 ✗ int i = hes_pattern->row[lnz];
578 ✗ int j = hes_pattern->col[lnz];
579 ✗ sparsity[i + n * j] = 1;
580 ✗ sparsity[j + n * i] = 1; // symmetric for display
581 }
582
583 ✗ for (int i = 0; i < n; i++) {
584 ✗ std::ostringstream oss;
585 ✗ oss << std::setw(2) << i << ": ";
586 ✗ for (int j = 0; j < n; j++) {
587 ✗ oss << (sparsity[i + n * j] ? '*' : ' ');
588 }
589 ✗ infoStreamPrint(OMC_LOG_MOO, 0, "%s", oss.str().c_str());
590 ✗ }
591
592 ✗ free(sparsity);
593 ✗ infoStreamPrint(OMC_LOG_MOO, 0, "=====================================");
594 }
595
596 ✗ void free_hessian_pattern(HESSIAN_PATTERN* hes_pattern) {
597 ✗ if (!hes_pattern) return;
598
599 ✗ int numColorPairs = hes_pattern->numColors * (hes_pattern->numColors + 1) / 2;
600 ✗ for (int i = 0; i < numColorPairs; i++) {
601 ✗ ColorPair* colorPair = hes_pattern->colorPairs[i];
602 ✗ if (!colorPair) continue;
603
604 ✗ for (int j = 0; j < colorPair->size; j++) {
605 ✗ free(colorPair->contributingRows[j]);
606 }
607
608 ✗ free(colorPair->contributingRows);
609 ✗ free(colorPair->numContributingRows);
610 ✗ free(colorPair->lnnzIndices);
611 ✗ free(colorPair->varPairs);
612 ✗ free(colorPair);
613 }
614
615 ✗ free(hes_pattern->colorPairs);
616 ✗ free(hes_pattern->row);
617 ✗ free(hes_pattern->col);
618
619 ✗ if (hes_pattern->colsForColor) {
620 ✗ for (int i = 0; i < hes_pattern->numColors; i++) {
621 ✗ free(hes_pattern->colsForColor[i]);
622 }
623 ✗ free(hes_pattern->colsForColor);
624 }
625 ✗ free(hes_pattern->colorSizes);
626
627 ✗ if (hes_pattern->cscJacIndexFromRowColor) {
628 ✗ for (int row = 0; row < hes_pattern->numFuncs; row++) {
629 ✗ free(hes_pattern->cscJacIndexFromRowColor[row]);
630 }
631 ✗ free(hes_pattern->cscJacIndexFromRowColor);
632 }
633
634 ✗ for (int row = 0; row < hes_pattern->numFuncs; row++) {
635 ✗ free(hes_pattern->ws_baseJac[row]);
636 }
637 ✗ free(hes_pattern->ws_baseJac);
638 ✗ free(hes_pattern->ws_oldX);
639 ✗ free(hes_pattern->ws_h);
640
641 ✗ free(hes_pattern);
642 }
643
644 // ====== EXTRAPOLATION ======
645
646 /**
647 * @brief Allocate and initialize internal workspace for Richardson extrapolation.
648 *
649 * @param[in] resultSize Number of result values computed by `fn` (length of result array).
650 * @param[in] maxSteps Maximum number of extrapolation steps that may be used.
651 * @return Pointer to an initialized ExtrapolationData struct.
652 */
653 ✗ ExtrapolationData* init_extrapolation_data(int resultSize, int maxSteps) {
654 ✗ ExtrapolationData* extrData = (ExtrapolationData*)malloc(sizeof(ExtrapolationData));
655 ✗ extrData->resultSize = resultSize;
656 ✗ extrData->maxSteps = maxSteps;
657 ✗ extrData->ws_results = (modelica_real**)malloc(maxSteps * sizeof(modelica_real*));
658 ✗ for (int i = 0; i < maxSteps; i++) {
659 ✗ extrData->ws_results[i] = (modelica_real*)malloc(resultSize * sizeof(modelica_real));
660 }
661 ✗ return extrData;
662 }
663
664 ✗ void free_extrapolation_data(ExtrapolationData* extrData) {
665 ✗ for (int i = 0; i < extrData->maxSteps; i++) {
666 ✗ free(extrData->ws_results[i]);
667 }
668 ✗ free(extrData->ws_results);
669 ✗ free(extrData);
670 ✗ }
671
672 /**
673 * @brief Apply in-place Richardson extrapolation using a generic computation function.
674 *
675 * Accepts a function of the form `f(args, h, result)`, evaluated at decreasing step sizes.
676 * Performs in-place extrapolation to increase accuracy. `steps <= 5` recommended to limit roundoff error.
677 *
678 * @param[in] extrData Workspace from init_extrapolation_data.
679 * @param[in] fn Function pointer: computes result := f(args, h).
680 * @param[in] args User data passed to fn.
681 * @param[in] h0 Initial step size.
682 * @param[in] steps Number of extrapolation steps (1 means no extrapolation!).
683 * @param[in] stepDivisor Step reduction factor (e.g. 2, then h_{i+1} = h_i / 2).
684 * @param[in] methodOrder Order of the underlying method (e.g. 1 for Forward Differences).
685 * @param[out] result Final extrapolated result.
686 */
687 ✗ void richardson_extrapolation(ExtrapolationData* extrData, computation_fn_t fn, void* args, modelica_real h0,
688 int steps, modelica_real stepDivisor, int methodOrder, modelica_real* result) {
689 /* call fn_ptr if no extrapolation is executed */
690 ✗ if (steps <= 1) {
691 ✗ fn(args, h0, result);
692 ✗ return;
693 }
694 ✗ else if (steps > extrData->maxSteps) {
695 ✗ warningStreamPrint(OMC_LOG_MOO, 0, "Requested extrapolation steps '%d' exceed maximum '%d', set in init_extrapolation_data. Using '%d' instead.\n",
696 steps, extrData->maxSteps, extrData->maxSteps);
697 ✗ steps = extrData->maxSteps;
698 }
699
700 /* compute all stages for extrapolation */
701 ✗ for (int i = 0; i < steps; i++) {
702 ✗ modelica_real h = h0 / pow(stepDivisor, i);
703 ✗ fn(args, h, extrData->ws_results[i]);
704 }
705
706 /* perform extrapolation: cancel taylor terms, in-place */
707 ✗ for (int j = 0; j < extrData->resultSize; j++) {
708 ✗ for (int k = 1; k < steps; k++) {
709 ✗ for (int i = steps - 1; i >= k; i--) {
710 ✗ modelica_real factor = pow(stepDivisor, methodOrder * k);
711 ✗ extrData->ws_results[i][j] = (factor * extrData->ws_results[i][j] - extrData->ws_results[i - 1][j]) / (factor - 1);
712 }
713 }
714 ✗ result[j] = extrData->ws_results[steps - 1][j];
715 }
716 }
717
718 /* wrapper for eval_hessian_fwd_differences */
719 ✗ void hessian_fwd_differences_wrapper(void* args, modelica_real h, modelica_real* result) {
720 HessianFiniteDiffArgs* hessianArgs = (HessianFiniteDiffArgs*)args;
721 ✗ eval_hessian_fwd_differences(hessianArgs->data, hessianArgs->threadData, hessianArgs->hes_pattern, h,
722 hessianArgs->u_indices, hessianArgs->lambda, hessianArgs->jac_csc, result);
723 ✗ }
724