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 / 155
Functions: 0.0% 0 / 0 / 18
Branches: 0.0% 0 / 0 / 170
OMCompiler/SimulationRuntime/c/moo/evaluations.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 <base/block_sparsity.h>
29
30 #include "hessian_finite_diff.h"
31
32 #include "evaluations.h"
33
34 namespace OpenModelica {
35
36 ✗ void init_eval(InfoGDOP& info, GDOP::FullSweepLayout& layout_lfg, GDOP::BoundarySweepLayout& layout_mr) {
37 ✗ init_eval_lfg(info, layout_lfg);
38 ✗ init_eval_mr(info, layout_mr);
39 ✗ }
40
41 /* just enumerate them in order: L -> f -> g. The correct placement will be handled in eval */
42 ✗ void init_eval_lfg(InfoGDOP& info, GDOP::FullSweepLayout& layout_lfg) {
43 int buf_index = 0;
44
45 ✗ if (layout_lfg.L) {
46 ✗ layout_lfg.L->buf_index = buf_index++;
47 }
48
49 ✗ for (auto& f : layout_lfg.f) {
50 ✗ f.buf_index = buf_index++;
51 }
52
53 ✗ for (auto& g : layout_lfg.g) {
54 ✗ g.buf_index = buf_index++;
55 }
56 ✗ }
57
58 /* just enumerate them in order: M -> r. The correct placement will be handled in eval */
59 ✗ void init_eval_mr(InfoGDOP& info, GDOP::BoundarySweepLayout& layout_mr) {
60 int buf_index = 0;
61
62 ✗ if (layout_mr.M) {
63 ✗ layout_mr.M->buf_index = buf_index++;
64 }
65
66 ✗ for (auto& r : layout_mr.r) {
67 ✗ r.buf_index = buf_index++;
68 }
69 ✗ }
70
71 ✗ void init_jac(InfoGDOP& info, GDOP::FullSweepLayout& layout_lfg, GDOP::BoundarySweepLayout& layout_mr) {
72 ✗ init_jac_lfg(info, layout_lfg);
73 ✗ init_jac_mr(info, layout_mr);
74 ✗ }
75
76 // rows are sorted as fLg (as in OpenModelica)
77 ✗ void init_jac_lfg(InfoGDOP& info, GDOP::FullSweepLayout& layout_lfg) {
78 /* full B Jacobian */
79 ✗ for (int nz = 0; nz < info.exc_jac->B.sparsity.nnz; nz++) {
80 ✗ int row = info.exc_jac->B.sparsity.row[nz]; // OpenModelica B matrix row
81 ✗ int col = info.exc_jac->B.sparsity.col[nz]; // OpenModelica B matrix col
82 int csc_buffer_entry_B = info.exc_jac->B.sparsity.coo_to_csc(nz); // OpenModelica B matrix CSC buffer index
83 ✗ FunctionLFG& fn = access_fLg_from_row(layout_lfg, row); // get function corresponding to the OM row
84 ✗ if (col < info.x_size) {
85 ✗ fn.jac.dx.push_back(JacobianSparsity{col, csc_buffer_entry_B});
86 }
87 ✗ else if (col < info.xu_size) {
88 ✗ fn.jac.du.push_back(JacobianSparsity{col - info.x_size, csc_buffer_entry_B});
89 }
90 else {
91 ✗ fn.jac.dp.push_back(JacobianSparsity{col - info.xu_size, csc_buffer_entry_B});
92 }
93 }
94 ✗ }
95
96 ✗ void init_jac_mr(InfoGDOP& info, GDOP::BoundarySweepLayout& layout_mr) {
97 /* M (first row) in C(COO) Jacobian */
98 int nz_C = 0;
99 ✗ if (info.mayer_exists) {
100 ✗ assert(layout_mr.M);
101
102 ✗ while (info.exc_jac->C.sparsity.row[nz_C] == 0) {
103 ✗ int col = info.exc_jac->C.sparsity.col[nz_C];
104
105 /* for now only final states: xf! no parameters, no dx0 */
106 ✗ if (col < info.x_size) {
107 /* just point to nz_C, since for 1 row, CSC == COO */
108 ✗ layout_mr.M->jac.dxf.push_back(JacobianSparsity{col, nz_C});
109 }
110 ✗ else if (col < info.xu_size) {
111 /* just point to nz_C, since for 1 row, CSC == COO */
112 ✗ layout_mr.M->jac.duf.push_back(JacobianSparsity{col - info.x_size, nz_C});
113 }
114
115 ✗ nz_C++;
116 }
117 }
118
119 /* r in D Jacobian */
120 ✗ for (int nz_D = 0; nz_D < info.exc_jac->D.sparsity.nnz; nz_D++) {
121 ✗ int row = info.exc_jac->D.sparsity.row[nz_D];
122 ✗ int col = info.exc_jac->D.sparsity.col[nz_D];
123 int csc_buffer_entry_D = info.exc_jac->D.sparsity.coo_to_csc(nz_D); // jac_buffer == OpenModelica D CSC buffer!
124 ✗ auto& fn = layout_mr.r[row];
125 ✗ if (col < info.x_size) {
126 /* add the Mayer offset, since the values f64* is [M, r] */
127 // Attention: this offset only works if D contains just r!!
128 ✗ fn.jac.dxf.push_back(JacobianSparsity{col, info.exc_jac->D.sparsity.nnz_offset + csc_buffer_entry_D});
129 }
130 ✗ else if (col < info.xu_size) {
131 ✗ fn.jac.duf.push_back(JacobianSparsity{col - info.x_size, info.exc_jac->D.sparsity.nnz_offset + csc_buffer_entry_D});
132 }
133 }
134 ✗ }
135
136 ✗ void init_hes(InfoGDOP& info, GDOP::FullSweepLayout& layout_lfg, GDOP::BoundarySweepLayout& layout_mr) {
137 ✗ init_hes_lfg(info, layout_lfg);
138 ✗ init_hes_mr(info, layout_mr);
139 ✗ }
140
141 ✗ void init_hes_lfg(InfoGDOP& info, GDOP::FullSweepLayout& layout_lfg) {
142 /* TODO: PARAMETERS include hes_lfg_pp and make it threaded (~numberThreads buffers with fancy pooling) */
143 auto& hes = layout_lfg.hes;
144
145 ✗ HESSIAN_PATTERN* hes_b = info.exc_hes->B.hessian;
146
147 ✗ for (int lnz = 0; lnz < hes_b->lnnz; lnz++) {
148 ✗ int row = hes_b->row[lnz];
149 ✗ int col = hes_b->col[lnz];
150
151 ✗ if (row < info.x_size && col < info.x_size) {
152 ✗ hes.dx_dx.push_back({row, col, lnz});
153 }
154 ✗ else if (row >= info.x_size && col < info.x_size) {
155 ✗ hes.du_dx.push_back({row - info.x_size, col, lnz});
156 }
157 else {
158 ✗ hes.du_du.push_back({row - info.x_size, col - info.x_size, lnz});
159 }
160 }
161 ✗ }
162
163 ✗ void init_hes_mr(InfoGDOP& info, GDOP::BoundarySweepLayout& layout_mr) {
164 auto& hes_mr = layout_mr.hes;
165
166 ✗ HESSIAN_PATTERN* hes_c = info.exc_hes->C.hessian;
167 ✗ HESSIAN_PATTERN* hes_d = info.exc_hes->D.hessian;
168
169 OrderedIndexSet hessian_mr;
170 std::vector<HessianSparsity> M_sparsities;
171 std::vector<int> M_cols;
172 ✗ if (info.mayer_exists) {
173 ✗ assert(layout_mr.M);
174
175 auto& mayer_term = *layout_mr.M;
176
177 // ignore x0 and p for now
178 ✗ for (auto& M_dxf : mayer_term.jac.dxf) M_cols.push_back(M_dxf.col);
179 ✗ for (auto& M_duf : mayer_term.jac.duf) M_cols.push_back(info.x_size + M_duf.col);
180
181 // estimate lower triangle sparsity pattern
182 ✗ for (size_t i = 0; i < M_cols.size(); i++) {
183 ✗ for (size_t j = 0; j <= i; j++) {
184 ✗ M_sparsities.push_back({M_cols[i], M_cols[j], 0});
185 }
186 }
187 ✗ hessian_mr.insert_sparsity(M_sparsities, 0, 0);
188 }
189
190 ✗ if (hes_d != nullptr) {
191 ✗ for (int nz = 0; nz < hes_d->lnnz; nz++) {
192 ✗ hessian_mr.set.insert({hes_d->row[nz], hes_d->col[nz]});
193 }
194 }
195
196 // create sparsity pattern struct(H(M) + H(r))
197 int lnz = 0;
198 std::map<std::pair<int, int>, int> sparsity_to_lnz;
199 ✗ for (auto pair : hessian_mr.set) {
200 ✗ auto [row, col] = pair;
201 ✗ sparsity_to_lnz[pair] = lnz;
202
203 ✗ if (row < info.x_size && col < info.x_size) {
204 ✗ hes_mr.dxf_dxf.push_back({row, col, lnz});
205 }
206 ✗ else if (row >= info.x_size && col < info.x_size) {
207 ✗ hes_mr.duf_dxf.push_back({row - info.x_size, col, lnz});
208 }
209 else {
210 ✗ hes_mr.duf_duf.push_back({row - info.x_size, col - info.x_size, lnz});
211 }
212 ✗ lnz++;
213 }
214
215 int c_index = 0;
216 ✗ info.exc_hes->C_to_Mr_buffer = FixedVector<std::pair<int, int>>(info.mayer_exists ? M_sparsities.size() : 0);
217
218 ✗ if (hes_c) {
219 int hes_c_index = 0;
220 ✗ for (const auto& mayer_hess : M_sparsities) {
221 ✗ while (hes_c_index < hes_c->lnnz &&
222 ✗ (hes_c->row[hes_c_index] < mayer_hess.row ||
223 ✗ (hes_c->row[hes_c_index] == mayer_hess.row && hes_c->col[hes_c_index] < mayer_hess.col))) {
224 ✗ hes_c_index++;
225 }
226 ✗ auto it = sparsity_to_lnz.find({mayer_hess.row, mayer_hess.col});
227 ✗ if (it != sparsity_to_lnz.end()) {
228 ✗ info.exc_hes->C_to_Mr_buffer[c_index++] = {hes_c_index, it->second};
229 ✗ hes_c_index++;
230 }
231 else {
232 ✗ Log::error("Hessian entry row = {}, col = {} from hes_d not found in pattern!", mayer_hess.row, mayer_hess.col);
233 ✗ abort();
234 }
235 }
236 }
237
238 ✗ info.exc_hes->D_to_Mr_buffer = FixedVector<std::pair<int, int>>(hes_d ? hes_d->lnnz : 0);
239 ✗ if (hes_d != nullptr) {
240 ✗ for (int i = 0; i < hes_d->lnnz; i++) {
241 ✗ int row = hes_d->row[i];
242 ✗ int col = hes_d->col[i];
243 ✗ auto it = sparsity_to_lnz.find({row, col});
244 ✗ if (it != sparsity_to_lnz.end()) {
245 info.exc_hes->D_to_Mr_buffer[i] = {i, it->second};
246 } else {
247 ✗ Log::error("Hessian entry row = {}, col = {} from hes_d not found in pattern!", row, col);
248 ✗ abort();
249 }
250 }
251 }
252 ✗ }
253
254 /* TODO: PARAMETERS add me */
255 ✗ void set_parameters(InfoGDOP& info, const f64* p) {
256 ✗ return;
257 }
258
259 ✗ void set_states(InfoGDOP& info, const f64* x_ij) {
260 ✗ std::memcpy(
261 ✗ info.data->localData[0]->realVars + info.index_x_real_vars,
262 x_ij,
263 ✗ info.x_size * sizeof(f64)
264 );
265 ✗ }
266
267 ✗ void set_inputs(InfoGDOP& info, const f64* u_ij) {
268 ✗ for (int u = 0; u < info.u_size; u++) {
269 ✗ info.data->localData[0]->realVars[info.u_indices_real_vars[u]] = u_ij[u];
270 }
271 ✗ }
272
273 ✗ void set_states_inputs(InfoGDOP& info, const f64* xu_ij) {
274 ✗ set_states(info, xu_ij);
275 ✗ set_inputs(info, xu_ij + info.x_size);
276 ✗ }
277
278 ✗ void set_time(InfoGDOP& info, const f64 t_ij) {
279 /* move time horizon to Modelica model time */
280 ✗ info.data->localData[0]->timeValue = t_ij;
281 ✗ }
282
283 ✗ void eval_ode_write(InfoGDOP& info, f64* eval_ode_buffer) {
284 /* f */
285 ✗ for (int der_x = 0; der_x < info.f_size; der_x++) {
286 ✗ eval_ode_buffer[der_x] = info.data->localData[0]->realVars[info.index_der_x_real_vars + der_x];
287 }
288 ✗ }
289
290 ✗ void eval_lfg_write(InfoGDOP& info, f64* eval_lfg_buffer) {
291 int nz = 0;
292 /* L */
293 ✗ if (info.lagrange_exists) {
294 ✗ eval_lfg_buffer[nz++] = info.data->localData[0]->realVars[info.index_lagrange_real_vars];
295 }
296 /* f */
297 ✗ for (int der_x = 0; der_x < info.f_size; der_x++) {
298 ✗ eval_lfg_buffer[nz++] = info.data->localData[0]->realVars[info.index_der_x_real_vars + der_x];
299 }
300 /* g */
301 ✗ for (int g = 0; g < info.g_size; g++) {
302 ✗ eval_lfg_buffer[nz++] = info.data->localData[0]->realVars[info.index_g_real_vars + g];
303 }
304 ✗ }
305
306 ✗ void eval_mr_write(InfoGDOP& info, f64* eval_mr_buffer) {
307 int nz = 0;
308 /* M */
309 ✗ if (info.mayer_exists) {
310 ✗ eval_mr_buffer[nz++] = info.data->localData[0]->realVars[info.index_mayer_real_vars];
311 }
312 /* r */
313 ✗ for (int r = 0; r < info.r_size; r++) {
314 ✗ eval_mr_buffer[nz++] = info.data->localData[0]->realVars[info.index_r_real_vars + r];
315 }
316 ✗ }
317
318 ✗ void jac_eval_write_first_row_as_csc(InfoGDOP& info, JACOBIAN* jacobian, f64* full_buffer,
319 f64* eval_jac_buffer, CscToCoo& exc) {
320 ✗ assert(jacobian && jacobian->sparsePattern);
321 ✗ evalJacobian(info.data, info.threadData, jacobian, NULL, full_buffer, FALSE);
322
323 ✗ for (int nz = 0; nz < exc.nnz_moved_row; nz++) {
324 ✗ eval_jac_buffer[nz] = full_buffer[exc.coo_to_csc(nz)];
325 }
326 ✗ }
327
328 } // namespace OpenModelica
329