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 / 1378
Functions: 0.0% 0 / 1 / 55
Branches: 0.0% 0 / 0 / 2440
OMCompiler/SimulationRuntime/c/dataReconciliation/dataReconciliation.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 "util/omc_error.h"
29 #include "simulation_data.h"
30 #include "openmodelica_func.h"
31 #include "simulation/solver/external_input.h"
32 #include "simulation/options.h"
33 #include "simulation/solver/model_help.h"
34 #include "simulation/jacobian_util.h"
35 #include <iostream>
36 #include <sstream>
37 #include <string>
38 #include <fstream>
39 #include <vector>
40 #include <algorithm>
41 #include <iomanip>
42 #include <stdlib.h>
43 #include <math.h>
44 #include <ctime>
45 #include <regex>
46 #include "omc_config.h"
47 #include "../util/omc_file.h"
48 #include <cmath>
49 #include "dataReconciliation.h"
50 using namespace std;
51
52 extern "C"
53 {
54 int dgesv_(int *n, int *nrhs, double *a, int *lda, int *ipiv, double *b, int *ldb, int *info);
55 int dgemm_(char *transa, char *transb, int *m, int *n, int *k, double *alpha, double *a, int *lda,
56 double *b, int *ldb, double *beta, double *c, int *ldc);
57 int dgetrf_(int *m, int *n, double *a, int *lda, int *ipiv, int *info);
58 int dgetri_(int *n, double *a, int *lda, int *ipiv, double *work, int *lwork, int *info);
59 int dscal_(int *n, double *da, double *dx, int *incx);
60 int dcopy_(int *n, double *dx, int *incx, double *dy, int *incy);
61 }
62
63 // only 200 values of chisquared x^2 values are added with degree of freedom
64 static double chisquaredvalue[200] = {3.84146, 5.99146, 7.81473, 9.48773, 11.0705, 12.5916, 14.0671, 15.5073, 16.919, 18.307, 19.6751, 21.0261, 22.362, 23.6848, 24.9958, 26.2962, 27.5871, 28.8693, 30.1435, 31.4104, 32.6706, 33.9244, 35.1725, 36.415, 37.6525, 38.8851, 40.1133, 41.3371, 42.557, 43.773, 44.9853, 46.1943, 47.3999, 48.6024, 49.8018, 50.9985, 52.1923, 53.3835, 54.5722, 55.7585, 56.9424, 58.124, 59.3035, 60.4809, 61.6562, 62.8296, 64.0011, 65.1708, 66.3386, 67.5048, 68.6693, 69.8322, 70.9935, 72.1532, 73.3115, 74.4683, 75.6237, 76.7778, 77.9305, 79.0819, 80.2321, 81.381, 82.5287, 83.6753, 84.8206, 85.9649, 87.1081, 88.2502, 89.3912, 90.5312, 91.6702, 92.8083, 93.9453, 95.0815, 96.2167, 97.351, 98.4844, 99.6169, 100.749, 101.879, 103.01, 104.139, 105.267, 106.395, 107.522, 108.648, 109.773, 110.898, 112.022, 113.145, 114.268, 115.39, 116.511, 117.632, 118.752, 119.871, 120.99, 122.108, 123.225, 124.342, 125.458, 126.574, 127.689, 128.804, 129.918, 131.031, 132.144, 133.257, 134.369, 135.48, 136.591, 137.701, 138.811, 139.921, 141.03, 142.138, 143.246, 144.354, 145.461, 146.567, 147.674, 148.779, 149.885, 150.989, 152.094, 153.198, 154.302, 155.405, 156.508, 157.61, 158.712, 159.814, 160.915, 162.016, 163.116, 164.216, 165.316, 166.415, 167.514, 168.613, 169.711, 170.809, 171.907, 173.004, 174.101, 175.198, 176.294, 177.39, 178.485, 179.581, 180.676, 181.77, 182.865, 183.959, 185.052, 186.146, 187.239, 188.332, 189.424, 190.516, 191.608, 192.7, 193.791, 194.883, 195.973, 197.064, 198.154, 199.244, 200.334, 201.423, 202.513, 203.602, 204.69, 205.779, 206.867, 207.955, 209.042, 210.13, 211.217, 212.304, 213.391, 214.477, 215.563, 216.649, 217.735, 218.82, 219.906, 220.991, 222.076, 223.16, 224.245, 225.329, 226.413, 227.496, 228.58, 229.663, 230.746, 231.829, 232.912};
65
66 struct csvData
67 {
68 int linecount;
69 int rowcount;
70 int columncount;
71 vector<double> xdata;
72 vector<double> sxdata;
73 vector<string> headers;
74 vector< vector<string> > rx;
75 };
76
77 struct errorData
78 {
79 string name;
80 string x;
81 string sx;
82 };
83
84 struct correlationData
85 {
86 vector<double> data;
87 vector<string> rowHeaders;
88 vector<string> columnHeaders;
89 };
90
91 struct correlationDataWarning
92 {
93 vector<string> diagonalEntry;
94 vector<string> aboveDiagonalEntry;
95 vector<errorData> warningInfo; // to warn user about entries closer to 1
96 };
97
98 struct matrixData
99 {
100 int rows;
101 int column;
102 double * data;
103 };
104
105 ✗ struct inputData
106 {
107 int rows;
108 int column;
109 double * data;
110 vector<int> index;
111 };
112
113 ✗ struct dataReconciliationData
114 {
115 csvData csvinputs;
116 matrixData xdiag;
117 matrixData reconciled_X;
118 matrixData reconciled_SX; // full covariance matrix
119 matrixData copyreconSx_diag; // only diagonal elements
120 double *newX;
121 double eps;
122 int iterationcount;
123 double value;
124 double J;
125 correlationDataWarning warning;
126 };
127
128 ✗ struct boundaryConditionData
129 {
130 std::vector<std::string> boundaryConditionVars;
131 double *boundaryConditionVarsResults;
132 double *reconSt_diag;
133 };
134
135 ✗ void copyReferenceFile(DATA * data, const std::string & filename)
136 {
137 ✗ std::string outputPath = std::string(omc_flagValue[FLAG_OUTPUT_PATH]) + "/" + std::string(data->modelData->modelFilePrefix) + filename;
138
139 // Read the reference file from -inputPath when it's given -- that's what
140 // OMEdit passes, since the model isn't necessarily run from its build
141 // directory -- falling back to the current directory otherwise.
142 ✗ std::string inputDir = omc_flag[FLAG_INPUT_PATH] ? std::string(omc_flagValue[FLAG_INPUT_PATH]) : std::string(".");
143 ✗ std::string referenceFile = inputDir + "/" + std::string(data->modelData->modelFilePrefix) + filename;
144
145 // If the output directory is the same as the input directory, the file is
146 // already there -- skip the copy. ofstream truncates on open, so opening
147 // the same file for both reading and writing here would wipe it out
148 // instead of copying it.
149 char resolvedOutputDir[PATH_MAX];
150 char resolvedInputDir[PATH_MAX];
151 ✗ if (realpath(omc_flagValue[FLAG_OUTPUT_PATH], resolvedOutputDir) && realpath(inputDir.c_str(), resolvedInputDir)
152 ✗ && std::string(resolvedOutputDir) == std::string(resolvedInputDir))
153 {
154 return;
155 }
156
157 ✗ ifstream ifstreamfile;
158 ✗ ifstreamfile.open(referenceFile);
159
160 ✗ if (ifstreamfile.good())
161 {
162 ✗ ofstream ofstreamfile;
163 ✗ ofstreamfile.open(outputPath);
164 ✗ ofstreamfile << ifstreamfile.rdbuf();
165 ✗ ofstreamfile.close();
166 ✗ ifstreamfile.close();
167 ✗ }
168 ✗ }
169
170 ✗ int getRelatedBoundaryConditions(DATA * data)
171 {
172 // check for _relatedBoundaryConditionsEquations.txt file exists to map the nonReconciled Vars failing with condition-2 of extraction algorithm
173 ✗ std::string relatedBoundaryConditionsFilename = string(data->modelData->modelFilePrefix) + "_relatedBoundaryConditionsEquations.html";
174
175 ✗ if (omc_flag[FLAG_OUTPUT_PATH])
176 {
177 ✗ relatedBoundaryConditionsFilename = string(omc_flagValue[FLAG_OUTPUT_PATH]) + "/" + relatedBoundaryConditionsFilename;
178 ✗ copyReferenceFile(data, "_relatedBoundaryConditionsEquations.html");
179 }
180
181 ✗ ifstream relatedBoundaryConditionsFilenameip(relatedBoundaryConditionsFilename);
182 string line;
183 int count = 0;
184 ✗ if (relatedBoundaryConditionsFilenameip.good())
185 {
186 ✗ while (relatedBoundaryConditionsFilenameip.good())
187 {
188 ✗ getline(relatedBoundaryConditionsFilenameip, line);
189 ✗ if (!line.empty())
190 {
191 ✗ count = count + 1;
192 }
193 }
194 ✗ relatedBoundaryConditionsFilenameip.close();
195 //omc_unlink(relatedBoundaryConditionsFilename.c_str());
196 }
197 ✗ return count;
198 ✗ }
199
200 /*
201 * create html report with error logs for D.1
202 */
203 ✗ void createErrorHtmlReport(DATA * data, int status = 0)
204 {
205 // create HTML Report with Error Logs
206 ✗ ofstream myfile;
207 ✗ time_t now = time(0);
208 ✗ std::stringstream htmlfile;
209 ✗ if (omc_flag[FLAG_OUTPUT_PATH])
210 {
211 ✗ htmlfile << string(omc_flagValue[FLAG_OUTPUT_PATH]) << "/" << data->modelData->modelFilePrefix << ".html";
212 }
213 else
214 {
215 ✗ htmlfile << data->modelData->modelFilePrefix << ".html";
216 }
217 string html = htmlfile.str();
218 ✗ myfile.open(html.c_str());
219
220 /* Add Overview Data */
221 ✗ myfile << "<!DOCTYPE html><html>\n <head> <h1> Data Reconciliation Report</h1></head> \n <body> \n ";
222 ✗ myfile << "<h2> Overview: </h2>\n";
223 ✗ myfile << "<table> \n";
224 ✗ myfile << "<tr> \n" << "<th align=right> Model file: </th> \n" << "<td>" << data->modelData->modelFileName << "</td> </tr>\n";
225 ✗ myfile << "<tr> \n" << "<th align=right> Model name: </th> \n" << "<td>" << data->modelData->modelName << "</td> </tr>\n";
226 ✗ myfile << "<tr> \n" << "<th align=right> Model directory: </th> \n" << "<td>" << data->modelData->modelDir << "</td> </tr>\n";
227 ✗ if (omc_flagValue[FLAG_DATA_RECONCILE_Sx])
228 {
229 ✗ myfile << "<tr> \n" << "<th align=right> Measurement input file: </th> \n" << "<td>" << omc_flagValue[FLAG_DATA_RECONCILE_Sx] << "</td> </tr>\n";
230 }
231 else
232 {
233 ✗ myfile << "<tr> \n" << "<th align=right> Measurement input file: </th> \n" << "<td style=color:red>" << "no file provided" << "</td> </tr>\n";
234 }
235 ✗ myfile << "<tr> \n" << "<th align=right> Correlation matrix input file: </th> \n" << "<td>" << "no file provided" << "</td> </tr>\n";
236 ✗ myfile << "<tr> \n" << "<th align=right> Generated: </th> \n" << "<td>" << ctime(&now) << " by "<< "<b>" << CONFIG_VERSION << "</b>" << "</td> </tr>\n";
237 ✗ myfile << "</table>\n";
238
239 /* add analysis section */
240 ✗ myfile << "<h2> Analysis: </h2>\n";
241 ✗ myfile << "<table> \n";
242 ✗ myfile << "<tr> \n" << "<th align=right> Number of auxiliary conditions: </th> \n" << "<td>" << data->modelData->nSetcVars << "</td> </tr>\n";
243 ✗ myfile << "<tr> \n" << "<th align=right> Number of measured variables: </th> \n" << "<td>" << data->modelData->ndataReconVars << "</td> </tr>\n";
244 ✗ myfile << "<tr> \n" << "<th align=right> Number of unmeasured variables: </th> \n" << "<td>" << data->modelData->nSetbVars << "</td> </tr>\n";
245 ✗ myfile << "<tr> \n" << "<th align=right> Number of related boundary conditions: </th> \n" << "<td>" << data->modelData->nRelatedBoundaryConditions << "</td> </tr>\n";
246 ✗ myfile << "</table> \n";
247
248 // Auxiliary Conditions
249 ✗ myfile << "<h3> <a href=" << data->modelData->modelFilePrefix << "_AuxiliaryConditions.html" << " target=_blank> Auxiliary conditions </a> </h3>\n";
250 // Intermediate Conditions
251 ✗ myfile << "<h3> <a href=" << data->modelData->modelFilePrefix << "_IntermediateEquations.html" << " target=_blank> Intermediate equations </a> </h3>\n";
252
253 ✗ if (data->modelData->nSetbVars > 0)
254 {
255 // Boundary Conditions
256 //myfile << "<h3> <a href=" << data->modelData->modelFilePrefix << "_BoundaryConditionsEquations.html" << " target=_blank> Boundary conditions </a> </h3>\n";
257 // Intermediate Conditions
258 ✗ myfile << "<h3> <a href=" << data->modelData->modelFilePrefix << "_BoundaryConditionIntermediateEquations.html" << " target=_blank> Intermediate equations for unmeasured variables </a> </h3>\n";
259 }
260
261 ✗ if (data->modelData->nRelatedBoundaryConditions > 0)
262 ✗ myfile << "<h3> <a href=" << data->modelData->modelFilePrefix << "_relatedBoundaryConditionsEquations.html" << " target=_blank> Related boundary conditions </a> </h3>\n";
263
264 // Error log
265 ✗ myfile << "<h2> <a href=" << data->modelData->modelFilePrefix << ".log" << " target=_blank> Errors </a> </h2>\n";
266 // copy the error log to output path
267 ✗ if (omc_flag[FLAG_OUTPUT_PATH])
268 {
269 ✗ copyReferenceFile(data, ".log");
270 }
271
272 // iteration vars file
273 ✗ if (status == 0)
274 {
275 ✗ myfile << "<h2> <a href=" << data->modelData->modelFilePrefix << "_iterationVars.txt" << " target=_blank> Iteration vars </a> </h2>\n";
276 }
277 // debug log
278 ✗ if (status == 0)
279 {
280 ✗ myfile << "<h2> <a href=" << data->modelData->modelFilePrefix << "_debug.txt" << " target=_blank> Debug log </a> </h2>\n";
281 }
282
283 ✗ myfile << "</table>\n";
284 ✗ myfile << "</body>\n</html>";
285 ✗ myfile.flush();
286 ✗ myfile.close();
287 ✗ }
288
289 /*
290 * create warning report for correlation coefficients data
291 */
292 ✗ void createCorrelationWarningReport(DATA * data, ofstream &myfile, correlationDataWarning &warningCorrelationData)
293 {
294 // create a warning log for correlation input file
295 ✗ if (!warningCorrelationData.aboveDiagonalEntry.empty() || !warningCorrelationData.diagonalEntry.empty() || !warningCorrelationData.warningInfo.empty())
296 {
297 ✗ myfile << "<h3> <a href=" << data->modelData->modelFilePrefix << "_warning.txt" << " target=_blank> Warnings </a> </h3>\n";
298 /* create a warning log file */
299 ✗ ofstream warningfile;
300 ✗ std::stringstream warning_file;
301 ✗ if (omc_flag[FLAG_OUTPUT_PATH])
302 {
303 ✗ warning_file << string(omc_flagValue[FLAG_OUTPUT_PATH]) << "/" << data->modelData->modelFilePrefix << "_warning.txt";
304 }
305 else
306 {
307 ✗ warning_file << data->modelData->modelFilePrefix << "_warning.txt";
308 }
309
310 ✗ warningfile.open(warning_file.str().c_str());
311 // user warning #1 : Diagonal entry for variable of interest <variable name> in correlation input file <input file name> is ignored
312 ✗ for (const auto &index : warningCorrelationData.diagonalEntry)
313 {
314 ✗ warningfile << "| warning | " << "Diagonal entry for variable of interest " << index << " in correlation input file " << omc_flagValue[FLAG_DATA_RECONCILE_Cx] << " is ignored" << "\n";
315 }
316 // user warning #2 : Above diagonal entry for variable of interest <variable name> in correlation input file <input file name> is ignored
317 ✗ for (const auto &index : warningCorrelationData.aboveDiagonalEntry)
318 {
319 ✗ warningfile << "| warning | " << "Above diagonal entry for variable of interest " << index << " in correlation input file " << omc_flagValue[FLAG_DATA_RECONCILE_Cx] << " is ignored" << "\n";
320 }
321 // user warning #3 : Entries in correlation matrix closer to 1 for variable of interest <variable name> in correlation input file <input file name>
322 ✗ for (const auto &info : warningCorrelationData.warningInfo)
323 {
324 ✗ warningfile << "| warning | " << "Entry for variable of interest " << info.name << " and variable of interest " << info.x << " in correlation input file " << omc_flagValue[FLAG_DATA_RECONCILE_Cx] << " is closer to 1: " << "[" << info.sx << "]" <<"\n";
325 }
326 ✗ warningfile.close();
327 ✗ }
328 ✗ }
329
330
331 /*
332 * create html report for data Reconciliation D.1
333 */
334 ✗ void createHtmlReportFordataReconciliation(DATA *data, csvData &csvinputs, matrixData &xdiag, matrixData &reconciled_X, matrixData &copyreconSx_diag, double *newX, double &eps, int &iterationcount, double &value, double &J, correlationDataWarning &warningCorrelationData, boundaryConditionData& boundaryconditiondata)
335 {
336 ✗ ofstream myfile;
337 ✗ time_t now = time(0);
338 ✗ std::stringstream htmlfile;
339 ✗ if (omc_flag[FLAG_OUTPUT_PATH])
340 {
341 ✗ htmlfile << string(omc_flagValue[FLAG_OUTPUT_PATH]) << "/" << data->modelData->modelFilePrefix << ".html";
342 }
343 else
344 {
345 ✗ htmlfile << data->modelData->modelFilePrefix << ".html";
346 }
347 string html = htmlfile.str();
348 ✗ myfile.open(html.c_str());
349
350 /* create a csv file */
351 ✗ ofstream csvfile;
352 ✗ std::stringstream csv_file;
353 ✗ if (omc_flag[FLAG_OUTPUT_PATH])
354 {
355 ✗ csv_file << string(omc_flagValue[FLAG_OUTPUT_PATH]) << "/" << data->modelData->modelFilePrefix << "_Outputs.csv";
356 }
357 else
358 {
359 ✗ csv_file << data->modelData->modelFilePrefix << "_Outputs.csv";
360 }
361
362 string tmpcsv = csv_file.str();
363 ✗ csvfile.open(tmpcsv.c_str());
364
365 // check for nonReconciledVars.txt file exists to map the nonReconciled Vars failing with condition-2 of extraction algorithm
366 std::string nonReconciledVarsFilename;
367 ✗ if (omc_flag[FLAG_OUTPUT_PATH])
368 {
369 ✗ nonReconciledVarsFilename = std::string(omc_flagValue[FLAG_OUTPUT_PATH]) + "/" + std::string(data->modelData->modelFilePrefix) + "_NonReconcilcedVars.txt";
370 ✗ copyReferenceFile(data, "_NonReconcilcedVars.txt");
371 }
372 else
373 {
374 ✗ nonReconciledVarsFilename = string(data->modelData->modelFilePrefix) + "_NonReconcilcedVars.txt";
375 }
376 vector<std::string> nonReconciledVars;
377
378 ✗ ifstream nonreconcilevarsip(nonReconciledVarsFilename);
379 string line;
380 ✗ if (nonreconcilevarsip.good())
381 {
382 ✗ while (nonreconcilevarsip.good())
383 {
384 ✗ getline(nonreconcilevarsip, line);
385 ✗ if (!line.empty())
386 {
387 //std::cout << "\n reading nonVariables of interest : " << line;
388 ✗ nonReconciledVars.push_back(line);
389 }
390 }
391 ✗ nonreconcilevarsip.close();
392 //omc_unlink(nonReconciledVarsFilename.c_str());
393 }
394
395 /* Add Overview Data */
396 ✗ myfile << "<!DOCTYPE html><html>\n <head> <h1> Data Reconciliation Report</h1></head> \n <body> \n ";
397 ✗ myfile << "<h2> Overview: </h2>\n";
398 ✗ myfile << "<table> \n";
399 ✗ myfile << "<tr> \n" << "<th align=right> Model file: </th> \n" << "<td>" << data->modelData->modelFileName << "</td> </tr>\n";
400 ✗ myfile << "<tr> \n" << "<th align=right> Model name: </th> \n" << "<td>" << data->modelData->modelName << "</td> </tr>\n";
401 ✗ myfile << "<tr> \n" << "<th align=right> Model directory: </th> \n" << "<td>" << data->modelData->modelDir << "</td> </tr>\n";
402 ✗ myfile << "<tr> \n" << "<th align=right> Measurement input file: </th> \n" << "<td>" << omc_flagValue[FLAG_DATA_RECONCILE_Sx] << "</td> </tr>\n";
403 ✗ if (omc_flagValue[FLAG_DATA_RECONCILE_Cx])
404 {
405 ✗ myfile << "<tr> \n" << "<th align=right> Correlation matrix input file: </th> \n" << "<td>" << omc_flagValue[FLAG_DATA_RECONCILE_Cx] << "</td> </tr>\n";
406 }
407 else
408 {
409 ✗ myfile << "<tr> \n" << "<th align=right> Correlation matrix input file: </th> \n" << "<td>" << "no file provided" << "</td> </tr>\n";
410 }
411 ✗ myfile << "<tr> \n" << "<th align=right> Generated: </th> \n" << "<td>" << ctime(&now) << " by "<< "<b>" << CONFIG_VERSION << "</b>" << "</td> </tr>\n";
412 ✗ myfile << "</table>\n";
413
414 /* Add Analysis data */
415 ✗ myfile << "<h2> Analysis: </h2>\n";
416 ✗ myfile << "<table> \n";
417 ✗ myfile << "<tr> \n" << "<th align=right> Number of auxiliary conditions: </th> \n" << "<td>" << data->modelData->nSetcVars << "</td> </tr>\n";
418 ✗ myfile << "<tr> \n" << "<th align=right> Number of measured variables: </th> \n" << "<td>" << data->modelData->ndataReconVars << "</td> </tr>\n";
419 ✗ myfile << "<tr> \n" << "<th align=right> Number of unmeasured variables: </th> \n" << "<td>" << data->modelData->nSetbVars << "</td> </tr>\n";
420 ✗ myfile << "<tr> \n" << "<th align=right> Number of related boundary conditions: </th> \n" << "<td>" << data->modelData->nRelatedBoundaryConditions << "</td> </tr>\n";
421 ✗ myfile << "<tr> \n" << "<th align=right> Number of iterations to convergence: </th> \n" << "<td>" << iterationcount << "</td> </tr>\n";
422 ✗ myfile << "<tr> \n" << "<th align=right> Final value of (J*/r) : </th> \n" << "<td>" << value << "</td> </tr>\n";
423 ✗ myfile << "<tr> \n" << "<th align=right> Epsilon : </th> \n" << "<td>" << eps << "</td> </tr>\n";
424 ✗ myfile << "<tr> \n" << "<th align=right> Final value of the objective function (J) : </th> \n" << "<td>" << J << "</td> </tr>\n";
425 //myfile << "<tr> \n" << "<th align=right> Chi-square value : </th> \n" << "<td>" << quantile(complement(chi_squared(data->modelData->nSetcVars), 0.05)) << "</td> </tr>\n";
426
427 ✗ if (data->modelData->nSetcVars > 200)
428 {
429 ✗ myfile << "<tr> \n" << "<th align=right> Chi-square value : </th> \n" << "<td>" << "NOT Available for equations > 200 in setC" << "</td> </tr>\n";
430 }
431 else
432 {
433 ✗ myfile << "<tr> \n" << "<th align=right> Chi-square value : </th> \n" << "<td>" << chisquaredvalue[data->modelData->nSetcVars - 1] << "</td> </tr>\n";
434 }
435 ✗ if (J <= chisquaredvalue[data->modelData->nSetcVars - 1])
436 {
437 ✗ myfile << "<tr> \n" << "<th align=right> Result of global test : </th> \n" << "<td>" << "TRUE" << "</td> </tr>\n";
438 }
439 else
440 {
441 ✗ myfile << "<tr> \n" << "<th align=right> Result of global test : </th> \n" << "<td>" << "FALSE" << "</td> </tr>\n";
442 }
443 ✗ myfile << "<tr> \n" << "<th align=right> Quality (J/Chi-square) : </th> \n" << "<td>" << J/chisquaredvalue[data->modelData->nSetcVars - 1] << "</td> </tr>\n";
444 ✗ myfile << "</table>\n";
445
446 // Auxiliary Conditions
447 ✗ myfile << "<h3> <a href=" << data->modelData->modelFilePrefix << "_AuxiliaryConditions.html" << " target=_blank> Auxiliary conditions </a> </h3>\n";
448
449 // Intermediate Conditions
450 ✗ myfile << "<h3> <a href=" << data->modelData->modelFilePrefix << "_IntermediateEquations.html" << " target=_blank> Intermediate equations for measured variables </a> </h3>\n";
451
452 ✗ if (data->modelData->nSetbVars > 0)
453 {
454 // Boundary Conditions
455 //myfile << "<h3> <a href=" << data->modelData->modelFilePrefix << "_BoundaryConditionsEquations.html" << " target=_blank> Boundary conditions </a> </h3>\n";
456 // Intermediate Conditions
457 ✗ myfile << "<h3> <a href=" << data->modelData->modelFilePrefix << "_BoundaryConditionIntermediateEquations.html" << " target=_blank> Intermediate equations for unmeasured variables </a> </h3>\n";
458 }
459
460 ✗ if (data->modelData->nRelatedBoundaryConditions > 0)
461 ✗ myfile << "<h3> <a href=" << data->modelData->modelFilePrefix << "_relatedBoundaryConditionsEquations.html" << " target=_blank> Related boundary conditions </a> </h3>\n";
462
463 // iteration vars file
464 ✗ myfile << "<h3> <a href=" << data->modelData->modelFilePrefix << "_iterationVars.txt" << " target=_blank> Iteration vars </a> </h3>\n";
465
466 // Debug log
467 ✗ myfile << "<h3> <a href=" << data->modelData->modelFilePrefix << "_debug.txt" << " target=_blank> Debug log </a> </h3>\n";
468
469 // create a warning log for correlation input file
470 ✗ createCorrelationWarningReport(data, myfile, warningCorrelationData);
471
472
473 /* Add Results data */
474 ✗ myfile << "<h2> Results: </h2>\n";
475 ✗ myfile << "<table border=2>\n";
476 ✗ myfile << "<tr>\n" << "<th> Variable to be Estimated </th>\n" << "<th> Unit </th>\n" << "<th> Description </th>\n" << "<th> Initial Measured Value </th>\n" << "<th> Estimated Value </th>\n" << "<th> Initial Uncertainty </th>\n" <<"<th> Estimated Uncertainty </th>\n";
477 ✗ csvfile << "Variable to be Estimated ," << "Initial Measured Value ," << "Estimated Value ," << "Initial Uncertainty ," << "Estimated Uncertainty,";
478 ✗ myfile << "<th> Result of Local Test </th>\n" << "<th> Local Quality </th>\n" << "<th> Comment </th>\n" << "</tr>\n";
479 ✗ csvfile << "Result of Local Test ," << "Local Quality ," << "\n";
480
481 // collect units and description
482
483 std::vector<std::string> unitString, unitStringUnMeasuredVariables;
484 std::vector<std::string> description, descriptionUnMeasuredVariables;
485 ✗ for (unsigned int r = 0; r < csvinputs.headers.size(); r++)
486 {
487 ✗ for (int i = 0; i < data->modelData->nVariablesReal; i++)
488 {
489 ✗ if (strcmp(data->modelData->realVarsData[i].info.name, csvinputs.headers[r].c_str())==0)
490 {
491 ✗ char *unitStr = omc_string_data(data->modelData->realVarsData[i].attribute.displayUnit);
492 ✗ unitString.push_back(unitStr);
493 ✗ description.push_back(data->modelData->realVarsData[i].info.comment);
494 }
495 else
496 {
497 // check units and description for unmeasured variables of interest
498 ✗ auto it = find(boundaryconditiondata.boundaryConditionVars.begin(), boundaryconditiondata.boundaryConditionVars.end(), data->modelData->realVarsData[i].info.name);
499 ✗ if (it != boundaryconditiondata.boundaryConditionVars.end())
500 {
501 ✗ char *unitStr = omc_string_data(data->modelData->realVarsData[i].attribute.displayUnit);
502 ✗ unitStringUnMeasuredVariables.push_back(unitStr);
503 ✗ descriptionUnMeasuredVariables.push_back(data->modelData->realVarsData[i].info.comment);
504 }
505 }
506 }
507 }
508
509 ✗ for (unsigned int r = 0; r < csvinputs.headers.size(); r++)
510 {
511 bool reconciled = true;
512 ✗ if (!nonReconciledVars.empty())
513 {
514 ✗ auto nonReconciledVar = std::find(nonReconciledVars.begin(), nonReconciledVars.end(), csvinputs.headers[r]);
515 ✗ if (nonReconciledVar != nonReconciledVars.end())
516 {
517 reconciled = false;
518 }
519 }
520
521 ✗ myfile << "<tr>\n";
522 // variables of interest
523 ✗ myfile << "<td>" << csvinputs.headers[r] << "</td>\n";
524 ✗ csvfile << csvinputs.headers[r] << ",";
525
526 // Unit
527 ✗ myfile << "<td>" << unitString[r] << "</td>\n";
528 //csvfile << unitString[r] << ",";
529
530 // description
531 ✗ myfile << "<td>" << description[r] << "</td>\n";
532 //csvfile << description[r] << ",";
533
534 // Initial Measured Values
535 ✗ myfile << "<td>" << xdiag.data[r] << "</td>\n";
536 ✗ csvfile << xdiag.data[r] << ",";
537
538 // Reconciled Values
539 ✗ myfile << "<td>" << reconciled_X.data[r] << "</td>\n";
540 ✗ csvfile << reconciled_X.data[r] << ",";
541
542 // Initial Uncertainty Values
543 ✗ myfile << "<td>" << csvinputs.sxdata[r] << "</td>\n";
544 ✗ csvfile << csvinputs.sxdata[r] << ",";
545
546 // Reconciled Uncertainty Values
547 ✗ myfile << "<td>" << copyreconSx_diag.data[r] << "</td>\n";
548 ✗ csvfile << copyreconSx_diag.data[r] << ",";
549
550 // Results of Local Tests
551 ✗ if (newX[r] < 1.96)
552 {
553 ✗ myfile << "<td>" << "TRUE" << "</td>\n";
554 ✗ csvfile << "TRUE" << ",";
555 }
556 else
557 {
558 ✗ myfile << "<td>" << "FALSE" << "</td>\n";
559 ✗ csvfile << "FALSE" << ",";
560 }
561
562 // Values of Local Tests
563 ✗ myfile << "<td>" << newX[r]/1.96 << "</td>\n";
564 ✗ csvfile << newX[r]/1.96 << ",\n";
565
566 // // Margin to Correctness(distance from 1.96)
567 // myfile << "<td>" << (1.96 - newX[r]) << "</td>\n";
568 // csvfile << (1.96 - newX[r]) << ",";
569
570 // comments
571 ✗ if (reconciled)
572 {
573 ✗ myfile << "<td>" << "" << "</td>\n";
574 //csvfile << "" << ",\n";
575 }
576 else
577 {
578 ✗ myfile << "<td style=color:red>" << "Not reconciled" << "</td>\n";
579 //csvfile << "Not reconciled" << ",\n";
580 }
581 ✗ myfile << "</tr>\n";
582 }
583
584 // check for boundary condition vars and unmeasured vars
585 ✗ if (!boundaryconditiondata.boundaryConditionVars.empty())
586 {
587 ✗ for (int i=0; i<boundaryconditiondata.boundaryConditionVars.size(); i++)
588 {
589 ✗ myfile << "<tr>\n";
590 ✗ myfile << "<td>" << boundaryconditiondata.boundaryConditionVars[i] << "</td>\n";
591 ✗ csvfile << boundaryconditiondata.boundaryConditionVars[i] << ",";
592
593 // unit
594 ✗ myfile << "<td>" << unitStringUnMeasuredVariables[i] << "</td>\n";
595 //csvfile << unitStringUnMeasuredVariables[i] << ",";
596
597 // description
598 ✗ myfile << "<td>"<< descriptionUnMeasuredVariables[i] << "</td>\n";
599 //csvfile << descriptionUnMeasuredVariables[i] << ",";
600
601 ✗ myfile << "<td> </td>\n";
602 ✗ csvfile << "" << ",";
603
604 ✗ myfile << "<td>" << boundaryconditiondata.boundaryConditionVarsResults[i] << "</td>\n";
605 ✗ csvfile << boundaryconditiondata.boundaryConditionVarsResults[i] << ",";
606
607 // Initial Uncertainty Values
608 ✗ myfile << "<td>" << "" << "</td>\n";
609 ✗ csvfile << "" << ",";
610
611 // Reconciled Uncertainty Values
612 ✗ myfile << "<td>" << boundaryconditiondata.reconSt_diag[i] << "</td>\n";
613 ✗ csvfile << boundaryconditiondata.reconSt_diag[i] << ",";
614
615 ✗ myfile << "<td>" << "" << "</td>\n";
616 ✗ csvfile << "" << ",";
617
618 ✗ myfile << "<td>" << "" << "</td>\n";
619 ✗ csvfile << "" << ",\n";
620
621 ✗ myfile << "<td>" << "" << "</td>\n";
622 //csvfile << "" << ",\n";
623 ✗ myfile << "</tr>\n";
624 }
625 }
626 ✗ csvfile.flush();
627 ✗ csvfile.close();
628 ✗ myfile << "</table>\n";
629 ✗ myfile << "</body>\n</html>";
630 ✗ myfile.flush();
631 ✗ myfile.close();
632 ✗ }
633
634 /*
635 * create html report with error logs for Boundary conditions D.2
636 */
637 ✗ void createErrorHtmlReportForBoundaryConditions(DATA * data, int status = 0)
638 {
639 // create HTML Report with Error Logs
640 ✗ ofstream myfile;
641 ✗ time_t now = time(0);
642 ✗ std::stringstream htmlfile;
643 ✗ if (omc_flag[FLAG_OUTPUT_PATH])
644 {
645 ✗ htmlfile << string(omc_flagValue[FLAG_OUTPUT_PATH]) << "/" << data->modelData->modelFilePrefix << "_BoundaryConditions.html";
646 }
647 else
648 {
649 ✗ htmlfile << data->modelData->modelFilePrefix << "_BoundaryConditions.html";
650 }
651 string html = htmlfile.str();
652 ✗ myfile.open(html.c_str());
653
654 /* Add Overview Data */
655 ✗ myfile << "<!DOCTYPE html><html>\n <head> <h1> Boundary Conditions Report </h1></head> \n <body> \n ";
656 ✗ myfile << "<h2> Overview: </h2>\n";
657 ✗ myfile << "<table> \n";
658 ✗ myfile << "<tr> \n" << "<th align=right> Model file: </th> \n" << "<td>" << data->modelData->modelFileName << "</td> </tr>\n";
659 ✗ myfile << "<tr> \n" << "<th align=right> Model name: </th> \n" << "<td>" << data->modelData->modelName << "</td> </tr>\n";
660 ✗ myfile << "<tr> \n" << "<th align=right> Model directory: </th> \n" << "<td>" << data->modelData->modelDir << "</td> </tr>\n";
661 // Sx input file
662 ✗ if (omc_flagValue[FLAG_DATA_RECONCILE_Sx])
663 {
664 ✗ myfile << "<tr> \n" << "<th align=right> Reconciled values input file: </th> \n" << "<td>" << omc_flagValue[FLAG_DATA_RECONCILE_Sx] << "</td> </tr>\n";
665 }
666 else
667 {
668 ✗ myfile << "<tr> \n" << "<th align=right> Reconciled values input file: </th> \n" << "<td style=color:red>" << "no file provided" << "</td> </tr>\n";
669 }
670 // Cx input file
671 ✗ if (omc_flagValue[FLAG_DATA_RECONCILE_Cx])
672 {
673 ✗ myfile << "<tr> \n" << "<th align=right> Reconciled covariance matrix input file: </th> \n" << "<td>" << omc_flagValue[FLAG_DATA_RECONCILE_Cx] << "</td> </tr>\n";
674 }
675 else
676 {
677 ✗ myfile << "<tr> \n" << "<th align=right> Reconciled covariance matrix input file: </th> \n" << "<td style=color:red>" << "no file provided" << "</td> </tr>\n";
678 }
679 ✗ myfile << "<tr> \n" << "<th align=right> Generated: </th> \n" << "<td>" << ctime(&now) << " by "<< "<b>" << CONFIG_VERSION << "</b>" << "</td> </tr>\n";
680 ✗ myfile << "</table>\n";
681
682 /* add analysis section */
683 ✗ myfile << "<h2> Analysis: </h2>\n";
684 ✗ myfile << "<table> \n";
685 ✗ myfile << "<tr> \n" << "<th align=right> Number of boundary conditions: </th> \n" << "<td>" << data->modelData->nSetcVars << "</td> </tr>\n";
686 ✗ myfile << "<tr> \n" << "<th align=right> Number of variables to be reconciled: </th> \n" << "<td>" << data->modelData->ndataReconVars << "</td> </tr>\n";
687 ✗ myfile << "</table> \n";
688
689 // Boundary Conditions
690 ✗ myfile << "<h3> <a href=" << data->modelData->modelFilePrefix << "_BoundaryConditionsEquations.html" << " target=_blank> Boundary conditions </a> </h3>\n";
691 // Intermediate Conditions
692 ✗ myfile << "<h3> <a href=" << data->modelData->modelFilePrefix << "_BoundaryConditionIntermediateEquations.html" << " target=_blank> Intermediate equations </a> </h3>\n";
693
694 // Error log
695 ✗ myfile << "<h2> <a href=" << data->modelData->modelFilePrefix << ".log" << " target=_blank> Errors </a> </h2>\n";
696 // copy the error log to output path
697 ✗ if (omc_flag[FLAG_OUTPUT_PATH])
698 {
699 ✗ copyReferenceFile(data, ".log");
700 }
701
702 // iteration vars file
703 ✗ if (status == 0)
704 {
705 ✗ myfile << "<h2> <a href=" << data->modelData->modelFilePrefix << "_iterationVars.txt" << " target=_blank> Iteration vars </a> </h2>\n";
706 }
707 // debug log
708 ✗ if (status == 0)
709 {
710 ✗ myfile << "<h2> <a href=" << data->modelData->modelFilePrefix << "_BoundaryConditions_debug.txt" << " target=_blank> Debug log </a> </h2>\n";
711 }
712
713 ✗ myfile << "</table>\n";
714 ✗ myfile << "</body>\n</html>";
715 ✗ myfile.flush();
716 ✗ myfile.close();
717 ✗ }
718
719 /*
720 * create HTML Report for Boundary Conditions D.2
721 */
722 ✗ void createHtmlReportForBoundaryConditions(DATA * data, std::vector<std::string> & boundaryConditionVars, double* values, double* uncertaintyValues, correlationDataWarning &warningCorrelationData)
723 {
724 ✗ ofstream myfile;
725 ✗ time_t now = time(0);
726 ✗ std::stringstream htmlfile;
727 ✗ if (omc_flag[FLAG_OUTPUT_PATH])
728 {
729 ✗ htmlfile << string(omc_flagValue[FLAG_OUTPUT_PATH]) << "/" << data->modelData->modelFilePrefix << "_BoundaryConditions.html";
730 }
731 else
732 {
733 ✗ htmlfile << data->modelData->modelFilePrefix << "_BoundaryConditions.html";
734 }
735 string html = htmlfile.str();
736 ✗ myfile.open(html.c_str());
737
738 /* create a csv file */
739 ✗ ofstream csvfile;
740 ✗ std::stringstream csv_file;
741 ✗ if (omc_flag[FLAG_OUTPUT_PATH])
742 {
743 ✗ csv_file << string(omc_flagValue[FLAG_OUTPUT_PATH]) << "/" << data->modelData->modelFilePrefix << "_BoundaryConditions_Outputs.csv";
744 }
745 else
746 {
747 ✗ csv_file << data->modelData->modelFilePrefix << "_BoundaryConditions_Outputs.csv";
748 }
749
750 string tmpcsv = csv_file.str();
751 ✗ csvfile.open(tmpcsv.c_str());
752
753 /* Add Overview Data */
754 ✗ myfile << "<!DOCTYPE html><html>\n <head> <h1> Boundary Conditions Report </h1></head> \n <body> \n ";
755 ✗ myfile << "<h2> Overview: </h2>\n";
756 ✗ myfile << "<table> \n";
757 ✗ myfile << "<tr> \n" << "<th align=right> Model file: </th> \n" << "<td>" << data->modelData->modelFileName << "</td> </tr>\n";
758 ✗ myfile << "<tr> \n" << "<th align=right> Model name: </th> \n" << "<td>" << data->modelData->modelName << "</td> </tr>\n";
759 ✗ myfile << "<tr> \n" << "<th align=right> Model directory: </th> \n" << "<td>" << data->modelData->modelDir << "</td> </tr>\n";
760 ✗ myfile << "<tr> \n" << "<th align=right> Reconciled values input file: </th> \n" << "<td>" << omc_flagValue[FLAG_DATA_RECONCILE_Sx] << "</td> </tr>\n";
761 ✗ if (omc_flagValue[FLAG_DATA_RECONCILE_Cx])
762 {
763 ✗ myfile << "<tr> \n" << "<th align=right> Reconciled covariance matrix input file: </th> \n" << "<td>" << omc_flagValue[FLAG_DATA_RECONCILE_Cx] << "</td> </tr>\n";
764 }
765 else
766 {
767 ✗ myfile << "<tr> \n" << "<th align=right> Correlation matrix input file: </th> \n" << "<td>" << "no file provided" << "</td> </tr>\n";
768 }
769 ✗ myfile << "<tr> \n" << "<th align=right> Generated: </th> \n" << "<td>" << ctime(&now) << " by "<< "<b>" << CONFIG_VERSION << "</b>" << "</td> </tr>\n";
770 ✗ myfile << "</table>\n";
771
772 /* Add Analysis data */
773 ✗ myfile << "<h2> Analysis: </h2>\n";
774 ✗ myfile << "<table> \n";
775 ✗ myfile << "<tr> \n" << "<th align=right> Number of boundary conditions: </th> \n" << "<td>" << data->modelData->nSetcVars << "</td> </tr>\n";
776 ✗ myfile << "<tr> \n" << "<th align=right> Number of variables to be reconciled: </th> \n" << "<td>" << data->modelData->ndataReconVars << "</td> </tr>\n";
777 ✗ myfile << "</table>\n";
778
779 // Boundary Conditions
780 ✗ myfile << "<h3> <a href=" << data->modelData->modelFilePrefix << "_BoundaryConditionsEquations.html" << " target=_blank> Boundary conditions </a> </h3>\n";
781
782 // Intermediate Conditions
783 ✗ myfile << "<h3> <a href=" << data->modelData->modelFilePrefix << "_BoundaryConditionIntermediateEquations.html" << " target=_blank> Intermediate equations </a> </h3>\n";
784
785 // iteration vars file
786 ✗ myfile << "<h3> <a href=" << data->modelData->modelFilePrefix << "_iterationVars.txt" << " target=_blank> Iteration vars </a> </h3>\n";
787
788 // Debug log
789 ✗ myfile << "<h3> <a href=" << data->modelData->modelFilePrefix << "_BoundaryConditions_debug.txt" << " target=_blank> Debug log </a> </h3>\n";
790
791 // create a warning log for correlation input file
792 ✗ createCorrelationWarningReport(data, myfile, warningCorrelationData);
793
794 /* Add Results data */
795 ✗ myfile << "<h2> Results: </h2>\n";
796 ✗ myfile << "<table border=2>\n";
797 ✗ myfile << "<tr>\n" << "<th> Boundary conditions </th>\n" << "<th> Values </th>\n" << "<th> Reconciled Half-width Confidence Intervals </th> </tr>\n";
798 ✗ csvfile << "Boundary conditions ," << "Values ," << "Reconciled Half-width Confidence Intervals," << "\n";
799
800 ✗ for (unsigned int r = 0; r < boundaryConditionVars.size(); r++)
801 {
802 ✗ myfile << "<tr>\n";
803 // Boundary Conditions
804 ✗ myfile << "<td>" << boundaryConditionVars[r] << "</td>\n";
805 ✗ csvfile << boundaryConditionVars[r] << ",";
806
807 // simulation Values
808 ✗ myfile << "<td>" << values[r] << "</td>\n";
809 ✗ csvfile << values[r] << ",";
810
811 // uncertainty Values
812 ✗ myfile << "<td>" << uncertaintyValues[r] << "</td>\n";
813 ✗ myfile << "</tr>\n";
814 ✗ csvfile << uncertaintyValues[r] << "," << "\n";
815 }
816
817 ✗ myfile << "</table>\n</html>";
818 ✗ myfile.close();
819 ✗ csvfile.close();
820 ✗ }
821
822 /*
823 * function which returns the index pos
824 * of input variables
825 */
826 ✗ int getVariableIndex(vector<string> headers, string name, ofstream & logfile, DATA * data)
827 {
828 int pos = -1;
829 ✗ for (unsigned int i = 0; i < headers.size(); i++)
830 {
831 //logfile << "founded headers " << headers[i] << i << "\n";
832 ✗ if (strcmp(headers[i].c_str(), name.c_str()) == 0)
833 {
834 ✗ pos = i;
835 ✗ break;
836 }
837 }
838 //logfile << "founded pos " << name << ": " << pos << "\n";
839 ✗ if (pos == -1)
840 {
841 //logfile << "Variable Name not Matched :" << name;
842 ✗ logfile << "| error | " << "CoRelation-Coefficient Variable Name not Matched: " << name << " ,getVariableIndex() failed!" << "\n";
843 ✗ logfile.close();
844 ✗ exit(1);
845 createErrorHtmlReport(data);
846 }
847 ✗ return pos;
848 }
849
850 /*
851 * check string is a valid double
852 */
853 ✗ bool isStringValidDouble(std::string &cref)
854 {
855 ✗ return std::regex_match(cref, std::regex("[-+]?[0-9]*\\.?[0-9]+([eE][-+]?[0-9]+)?"));
856 }
857
858 /*
859 * check string is a empty, (i.e) contains only "," in csv input
860 * also ignore lines starting with c comments //
861 */
862 ✗ bool isLineEmptyData(std::string &cref)
863 {
864 ✗ return std::regex_match(cref, std::regex("^[,|/]+.*"));
865 }
866
867 /*
868 * function which checks whether a variable
869 * is unmeasured (i.e) uncertain=Uncertainty.propagate
870 */
871 ✗ bool isUnmeasuredVariables(DATA* data, const char* name)
872 {
873 ✗ char **unmeasuredvariable = (char**) malloc(data->modelData->nSetbVars * sizeof(char*));
874 ✗ data->callback->dataReconciliationUnmeasuredVariables(data, unmeasuredvariable);
875 bool found = false;
876 // check for unmeasured variables
877 ✗ for (int i = 0; i < data->modelData->nSetbVars; i++)
878 {
879 ✗ if (strcmp(unmeasuredvariable[i], name) == 0)
880 {
881 found = true;
882 break;
883 }
884 }
885 ✗ free(unmeasuredvariable);
886 ✗ return found;
887 }
888
889 //----------------------------------------------
890 // Helper: Update Reconciled.mo with reconciled values
891 //----------------------------------------------
892 ✗ void updateReconciledMo(DATA * data, threadData_t * threadData, vector<string> headers, double * reconciled_X, ofstream & logfile)
893 {
894 ✗ std::string modelPrefix(data->modelData->modelFilePrefix);
895 std::replace(modelPrefix.begin(), modelPrefix.end(), '.', '_');
896
897 // check for reconciled.mo file to update with reconciled values
898 std::string reconciledMoFile, reconciledValuesMoFile;
899 ✗ if (omc_flag[FLAG_OUTPUT_PATH])
900 {
901 ✗ reconciledMoFile = std::string(omc_flagValue[FLAG_OUTPUT_PATH]) + "/" + data->modelData->modelFilePrefix +"_Reconciled_tmp.mo";
902 ✗ reconciledValuesMoFile = std::string(omc_flagValue[FLAG_OUTPUT_PATH]) + "/" + "Reconciled_" + modelPrefix + ".mo";
903 ✗ copyReferenceFile(data, "_Reconciled_tmp.mo");
904 }
905 else
906 {
907 ✗ reconciledMoFile = std::string(data->modelData->modelFilePrefix) + "_Reconciled_tmp.mo";
908 ✗ reconciledValuesMoFile = "Reconciled_" + modelPrefix + ".mo";
909 }
910 ✗ std::ifstream infile(reconciledMoFile);
911 ✗ if (!infile.is_open())
912 {
913 // just give a warning, if file not found as this is optional and user may not want to update the file with reconciled values
914 ✗ warningStreamPrint(OMC_LOG_STDOUT, 0, "Reconciled modelica file path not found %s.", reconciledMoFile.c_str());
915 ✗ logfile << "| warning | " << "Measurement input file path not found " << reconciledMoFile << "\n";
916 }
917
918 ✗ std::ofstream outfile(reconciledValuesMoFile);
919 ✗ if (!outfile.is_open())
920 {
921 // just give a warning, if file not found as this is optional and user may not want to update the file with reconciled values
922 ✗ warningStreamPrint(OMC_LOG_STDOUT, 0, "Cannot open reconciled values output file %s.", reconciledValuesMoFile.c_str());
923 ✗ logfile << "| warning | " << "Cannot open reconciled values output file " << reconciledValuesMoFile << "\n";
924 }
925
926 std::string line;
927 int count = 1;
928 int varCount = 1; // counter for variables of interest
929 ✗ while (std::getline(infile, line))
930 {
931 ✗ if (count > 3 && varCount <= data->modelData->ndataReconVars)
932 {
933 ✗ std::string variableName(headers[varCount-1]);
934 std::replace(variableName.begin(), variableName.end(), '.', '_');
935 ✗ outfile << " parameter Real " << variableName << " = " << reconciled_X[varCount-1] << ";\n";
936 ✗ varCount++;
937 ✗ }
938 else
939 {
940 // copy other lines as it is
941 ✗ count++;
942 ✗ outfile << line << "\n";
943 }
944 }
945 ✗ infile.close();
946 ✗ outfile.close();
947 //omc_unlink(reconciledMoFile.c_str());
948 ✗ logfile << "| info | " << "Reconciled modelica file updated successfully " << reconciledValuesMoFile << "\n";
949 ✗ }
950
951 /*
952 * Function which reads the csv file
953 * and stores the initial measured value X and HalfWidth confidence
954 * interval Wx and also the input variable names
955 */
956 ✗ csvData readMeasurementInputFile(ofstream & logfile, DATA * data, threadData_t * threadData, bool boundaryConditions = false)
957 {
958 char * filename = NULL;
959 ✗ filename = (char*) omc_flagValue[FLAG_DATA_RECONCILE_Sx];
960
961 ✗ if (filename == NULL && !boundaryConditions)
962 {
963 ✗ errorStreamPrint(OMC_LOG_STDOUT, 0, "Measurement input file not provided (eg:-sx=filename.csv), DataReconciliation cannot be computed!.");
964 ✗ logfile << "| error | " << "Measurement input file not provided (eg:-sx=filename.csv), DataReconciliation cannot be computed!.\n";
965 ✗ logfile.close();
966 ✗ createErrorHtmlReport(data);
967 ✗ exit(1);
968 }
969
970 ✗ if (filename == NULL && boundaryConditions)
971 {
972 ✗ errorStreamPrint(OMC_LOG_STDOUT, 0, "Reconciled values input file not provided (eg:-sx=filename.csv), Boundary conditions cannot be computed!.");
973 ✗ logfile << "| error | " << "Reconciled values input file not provided (eg:-sx=filename.csv), Boundary conditions cannot be computed!.\n";
974 ✗ logfile.close();
975 ✗ createErrorHtmlReportForBoundaryConditions(data);
976 ✗ exit(1);
977 }
978 /*
979 * fix issue https://github.com/OpenModelica/OpenModelica/issues/13797
980 * check if filepath uses uri format and convert it to absolute path
981 * eg: modelica://Modelica/Resources/Files/filename.csv => /absolute/path/to/filename.csv
982 */
983 ✗ modelica_string uri = OpenModelica_uriToFilename(omc_string_new(filename));
984 ✗ std::string filenameStr = omc_string_data(uri);
985
986 ✗ ifstream ip(filenameStr.c_str());
987 string line;
988 vector<double> xdata;
989 vector<double> sxdata;
990 vector<string> names;
991 //vector<double> rx_ik;
992 vector< vector<string> > rx;
993 int Sxrowcount = 0;
994 ✗ int linecount = 1;
995 int Sxcolscount = 0;
996 bool flag = false;
997 vector<errorData> errorInfo;
998 vector<int> errorInfoHeaders;
999
1000 ✗ if (!ip.good() && !boundaryConditions)
1001 {
1002 ✗ errorStreamPrint(OMC_LOG_STDOUT, 0, "Measurement input file path not found %s.",filename);
1003 ✗ logfile << "| error | " << "Measurement input file path not found " << filename << "\n";
1004 ✗ logfile.close();
1005 ✗ createErrorHtmlReport(data);
1006 ✗ exit(1);
1007 }
1008
1009 ✗ if (!ip.good() && boundaryConditions)
1010 {
1011 ✗ errorStreamPrint(OMC_LOG_STDOUT, 0, "Reconciled values input file path not found %s.", filename);
1012 ✗ logfile << "| error | " << "Reconciled values input file path not found " << filename << "\n";
1013 ✗ logfile.close();
1014 ✗ createErrorHtmlReportForBoundaryConditions(data);
1015 ✗ exit(1);
1016 }
1017
1018 ✗ while (ip.good())
1019 {
1020 ✗ getline(ip, line);
1021 vector<string> t1;
1022
1023 // allow comments on the line#1 until finding the headers
1024 ✗ if (linecount == 1 && isLineEmptyData(line))
1025 {
1026 continue;
1027 }
1028
1029 ✗ if (linecount > 1 && !line.empty() && !isLineEmptyData(line))
1030 {
1031 //std::cout << "\nline info:" << line;
1032 std::replace(line.begin(), line.end(), ';', ',');
1033 // remove whitespace in a string
1034 ✗ line.erase(std::remove_if(line.begin(), line.end(), ::isspace), line.end());
1035
1036 ✗ stringstream ss(line);
1037 string temp;
1038 int columnCount = 0;
1039 bool col0 = false, col1 = false, col2 = false;
1040 ✗ while (getline(ss, temp, ','))
1041 {
1042 // ignore unmeasured variables of interest
1043 ✗ if (columnCount == 0 && isUnmeasuredVariables(data, temp.c_str()))
1044 {
1045 col0 = true;
1046 col1 = true;
1047 col2 = true;
1048 break;
1049 }
1050
1051 ✗ if (columnCount == 0)
1052 {
1053 // // error : no variable of interest is provided by user at column #1
1054 ✗ if (temp.empty())
1055 {
1056 ✗ errorInfoHeaders.push_back(linecount);
1057 }
1058 col0 = true;
1059 ✗ names.push_back(temp.c_str());
1060 ✗ Sxrowcount++;
1061 ✗ if (!flag)
1062 {
1063 ✗ Sxcolscount++;
1064 }
1065 }
1066 ✗ if (columnCount == 1 && !temp.empty() && isStringValidDouble(temp))
1067 {
1068 //std::cout << "\n x: " << temp << " type : "<< isStringValidDouble(temp);
1069 //logfile << "xdata" << temp << " double" << atof(temp.c_str()) <<"\n";
1070 col1 = true;
1071 ✗ xdata.push_back(atof(temp.c_str()));
1072 ✗ if (!flag)
1073 {
1074 ✗ Sxcolscount++;
1075 }
1076 }
1077 ✗ if (columnCount == 2 && !temp.empty() && isStringValidDouble(temp))
1078 {
1079 //std::cout << "\n sxdata: " << temp << " valid type " << isStringValidDouble(temp);
1080 //logfile << "sxdata" << temp << " double" << atof(temp.c_str()) <<"\n";
1081 col2 = true;
1082 ✗ sxdata.push_back(atof(temp.c_str()));
1083 ✗ if (!flag)
1084 {
1085 ✗ Sxcolscount++;
1086 }
1087 }
1088 // ignore columns greater than 3
1089 ✗ if (columnCount > 2)
1090 {
1091 break;
1092 }
1093 ✗ columnCount++;
1094 }
1095 flag = true;
1096
1097 ✗ if (!col0 || !col1 || !col2)
1098 {
1099 ✗ std::string column1 = "(no-Value/wrong-Type)", column2 = "(no-Value/wrong-Type)", column3 = "(no-Value/wrong-Type)";
1100 ✗ if (col0)
1101 {
1102 column1 = names.back();
1103 }
1104 ✗ if (col1)
1105 {
1106 ✗ column2 = std::to_string(xdata.back());
1107 }
1108 ✗ if (col2)
1109 {
1110 ✗ column3 = std::to_string(sxdata.back());
1111 }
1112 errorData info = {column1, column2, column3};
1113 ✗ errorInfo.push_back(info);
1114 ✗ }
1115 ✗ }
1116 ✗ linecount++;
1117 ✗ }
1118
1119 // user error : variable of interest is missing in column #1
1120 ✗ if (!errorInfoHeaders.empty())
1121 {
1122 ✗ for (const auto &line : errorInfoHeaders)
1123 {
1124 ✗ errorStreamPrint(OMC_LOG_STDOUT, 0, "the name of the variable of interest in measurement input file %s is missing in line #%d ", filename, line);
1125 ✗ logfile << "| error | " << "the name of the variable of interest in measurement input file " << filename << " is missing in line #" << line << "\n";
1126 }
1127 ✗ logfile.close();
1128 ✗ createErrorHtmlReport(data);
1129 ✗ exit(1);
1130 }
1131
1132 // user error #6: Entry for variable of interest <variable name> in measurement input file <input file name> is incorrect: <reason (no value or incorrect value type)>
1133 ✗ if (!errorInfo.empty())
1134 {
1135 ✗ for (const auto & info : errorInfo)
1136 {
1137 ✗ errorStreamPrint(OMC_LOG_STDOUT, 0, "Entry for variable of interest %s in measurement input file %s is incorrect because of (no-Value/wrong-Type), with following data: [%s, %s, %s] ", info.name.c_str(), filename, info.name.c_str(), info.x.c_str(), info.sx.c_str());
1138 ✗ logfile << "| error | " << "Entry for variable of interest " << info.name << " in measurement input file " << filename << " is incorrect because of (no-Value/wrong-Type), with following data: " << "[" << info.name << ", " << info.x << ", " << info.sx << "]" <<"\n";
1139 }
1140 ✗ logfile.close();
1141 ✗ createErrorHtmlReport(data);
1142 ✗ exit(1);
1143 }
1144
1145 //logfile << "csvdata header:" << "header length: " << names.size() << " " << names[0] << names[1] << names[2] << "" << "\n";
1146 //logfile << "linecount:" << linecount << " " << "rowcount :" << Sxrowcount << " " << "colscount:" << Sxcolscount << "\n";
1147
1148 ✗ csvData csvData = {linecount, Sxrowcount, Sxcolscount, xdata, sxdata, names, rx};
1149 ✗ return csvData;
1150 ✗ }
1151
1152 /*
1153 * Function which arranges the elements in column major
1154 */
1155 ✗ void initColumnMatrix(vector<double> data, int rows, int cols, double * tempSx)
1156 {
1157 ✗ for (int i = 0; i < rows; i++)
1158 {
1159 ✗ for (int j = 0; j < cols; j++)
1160 {
1161 // store the matrix in column order
1162 ✗ tempSx[j + i * rows] = data[i + j * rows];
1163 }
1164 }
1165 ✗ }
1166
1167 /*
1168 * Function to print and debug whether the matrices are stored in column major
1169 */
1170 ✗ void printColumnAlginment(double * matrix, int rows, int cols, string name)
1171 {
1172 ✗ cout << "\n" << "************ " << name << " **********" << "\n";
1173 ✗ for (int i = 0; i < rows * cols; i++)
1174 {
1175 ✗ cout << matrix[i] << " ";
1176 }
1177 ✗ cout << "\n";
1178 ✗ }
1179
1180 /*
1181 * Function to Print the matrix in row based format
1182 */
1183 ✗ void printMatrix(double * matrix, int rows, int cols, string name, ofstream & logfile)
1184 {
1185 ✗ logfile << "\n" << "************ " << name << " **********" << "\n";
1186 ✗ for (int i = 0; i < rows; i++)
1187 {
1188 ✗ for (int j = 0; j < cols; j++)
1189 {
1190 //cout << setprecision(5);
1191 ✗ logfile << std::right << setw(15) << matrix[i + j * rows];
1192 ✗ logfile.flush();
1193 }
1194 ✗ logfile << "\n";
1195 }
1196 ✗ logfile << "\n";
1197 ✗ }
1198
1199 /*
1200 * Function to Print the matrix in row based format
1201 */
1202 ✗ void printMatrixModelicaFormat(double * matrix, int rows, int cols, string name, ofstream & logfile)
1203 {
1204 ✗ logfile << "\n" << "************ " << name << " **********" << "\n";
1205 ✗ logfile << "\n[";
1206 ✗ for (int i = 0; i < rows; i++)
1207 {
1208 ✗ for (int j = 0; j < cols; j++)
1209 {
1210 //cout << setprecision(5);
1211 ✗ if (j == cols - 1)
1212 {
1213 ✗ logfile << std::right << setw(15) << matrix[i + j * rows] << ";\n";
1214 }
1215 else
1216 {
1217 ✗ logfile << std::right << setw(15) << matrix[i + j * rows] << ",";
1218 }
1219
1220 ✗ logfile.flush();
1221 }
1222 //logfile << ";\n";
1223 }
1224 ✗ logfile << "\n";
1225 ✗ }
1226
1227 /*
1228 *
1229 Function to Print the matrix in row based format with headers
1230 */
1231 ✗ void printMatrixWithHeaders(double * matrix, int rows, int cols, vector<string> headers, string name, ofstream & logfile)
1232 {
1233 ✗ logfile << "\n" << "************ " << name << " **********" << "\n";
1234 ✗ for (int i = 0; i < rows; i++)
1235 {
1236 ✗ logfile << std::right << setw(10) << headers[i];
1237 ✗ for (int j = 0; j < cols; j++)
1238 {
1239 //cout << setprecision(5);
1240 ✗ logfile << std::right << setw(15) << matrix[i + j * rows];
1241 ✗ logfile.flush();
1242 //printf("% .5e ", matrix[i+j*rows]);
1243 }
1244 ✗ logfile << "\n";
1245 }
1246 ✗ logfile << "\n";
1247 ✗ }
1248
1249 /*
1250 *
1251 Function to Print the matrix in row based format with headers
1252 */
1253 ✗ void printBoundaryConditionsResults(double * matrixA, double * matrixB, int rows, int cols, vector<string> headers, string name, ofstream & logfile)
1254 {
1255 ✗ logfile << "\n" << "************ " << name << " **********" << "\n";
1256 ✗ logfile << "\n Boundary conditions" << setw(20) << "Values" << setw(45) << "Half-width Confidence Interval" << "\n";
1257 ✗ for (int i = 0; i < rows; i++)
1258 {
1259 ✗ logfile << std::right << setw(20) << headers[i];
1260 ✗ for (int j = 0; j < cols; j++)
1261 {
1262 //cout << setprecision(5);
1263 ✗ logfile << std::right << setw(20) << matrixA[i + j * rows] << setw(25) << matrixB[i + j * rows];
1264 ✗ logfile.flush();
1265 //printf("% .5e ", matrix[i+j*rows]);
1266 }
1267 ✗ logfile << "\n";
1268 }
1269 ✗ logfile << "\n";
1270 ✗ }
1271
1272 ✗ void dumpReconciledSxToCSV(double * matrix, int rows, int cols, vector<string> headers, DATA * data)
1273 {
1274 /* create a csv file */
1275 ✗ ofstream csvfile;
1276 ✗ std::stringstream csv_file;
1277 ✗ if (omc_flag[FLAG_OUTPUT_PATH])
1278 {
1279 ✗ csv_file << string(omc_flagValue[FLAG_OUTPUT_PATH]) << "/" << data->modelData->modelFilePrefix << "_Reconciled_Sx.csv";
1280 }
1281 else
1282 {
1283 ✗ csv_file << data->modelData->modelFilePrefix << "_Reconciled_Sx.csv";
1284 }
1285
1286 string tmpcsv = csv_file.str();
1287 ✗ csvfile.open(tmpcsv.c_str());
1288
1289 ✗ csvfile << "Sxij" << ",";
1290 ✗ for (auto it : headers)
1291 {
1292 //std::cout << "headers : " << it << "\n";
1293 ✗ csvfile << it << ",";
1294 }
1295 ✗ csvfile << "\n";
1296
1297 ✗ for (int i = 0; i < rows; i++)
1298 {
1299 ✗ csvfile << headers[i] << ",";
1300 ✗ for (int j = 0; j < cols; j++)
1301 {
1302 //cout << setprecision(5);
1303 ✗ csvfile << matrix[i + j * rows] << ",";
1304 //csvfile.flush();
1305 //printf("% .5e ", matrix[i+j*rows]);
1306 }
1307 ✗ csvfile << "\n";
1308 }
1309 //csvfile << "\n";
1310 ✗ csvfile.flush();
1311 ✗ csvfile.close();
1312 ✗ }
1313
1314 /*
1315 *Function to Print the vecomatrix in row based format with headers
1316 *based on vector arrays
1317 */
1318 ✗ void printVectorMatrixWithHeaders(vector<double> matrix, int rows, int cols, vector<string> headers, string name, ofstream & logfile)
1319 {
1320 ✗ logfile << "\n" << "************ " << name << " **********" << "\n";
1321 ✗ for (int i = 0; i < rows; i++)
1322 {
1323 ✗ logfile << std::right << setw(10) << headers[i];
1324 ✗ for (int j = 0; j < cols; j++)
1325 {
1326 //cout << setprecision(5);
1327 ✗ logfile << std::right << setw(15) << matrix[i + j * rows];
1328 ✗ logfile.flush();
1329 //printf("% .5e ", matrix[i+j*rows]);
1330 }
1331 ✗ logfile << "\n";
1332 }
1333 ✗ logfile << "\n";
1334 ✗ }
1335
1336 /*
1337 Function to Print the corelation matrix in row based format with headers
1338 */
1339 ✗ void printCorelationMatrix(vector<double> cx_data, vector<string> rowHeaders, vector<string> columnHeaders, string name, ofstream & logfile, correlationDataWarning & warningCorrelationData)
1340 {
1341 ✗ if (cx_data.empty())
1342 {
1343 return;
1344 }
1345
1346 ✗ logfile << "\n" << "************ " << name << " **********" << "\n";
1347 ✗ for (int i = 0; i < rowHeaders.size(); i++)
1348 {
1349 logfile << std::right << setw(10) << rowHeaders[i];
1350 ✗ for (int j = 0; j < columnHeaders.size(); j++)
1351 {
1352 ✗ if (i == j && cx_data[columnHeaders.size() * i + j] != 0)
1353 {
1354 ✗ warningCorrelationData.diagonalEntry.push_back(rowHeaders[i]);
1355 }
1356 ✗ else if (j > i && cx_data[columnHeaders.size() * i + j] != 0)
1357 {
1358 ✗ warningCorrelationData.aboveDiagonalEntry.push_back(rowHeaders[i]);
1359 }
1360 ✗ logfile << std::right << setw(15) << cx_data[columnHeaders.size() * i + j];
1361 }
1362 ✗ logfile << "\n";
1363 }
1364 ✗ logfile << "\n";
1365 }
1366
1367 /*
1368 *
1369 Function Which gets the diagonal elements of the matrix
1370 */
1371 ✗ void getDiagonalElements(double * matrix, int rows, int cols, double * result)
1372 {
1373 int k = 0;
1374 ✗ for (int i = 0; i < rows; i++)
1375 {
1376 ✗ for (int j = 0; j < cols; j++)
1377 {
1378 ✗ if (i == j)
1379 {
1380 ✗ result[k++] = matrix[i + j * rows];
1381 }
1382 }
1383 }
1384 ✗ }
1385
1386 /*
1387 * Function to transpose the Matrix
1388 */
1389 ✗ void transposeMatrix(double * jacF, double * jacFT, int rows, int cols)
1390 {
1391 ✗ for (int i = 0; i < rows; i++)
1392 {
1393 ✗ for (int j = 0; j < cols; j++)
1394 {
1395 // Perform matrix transpose store the elements in column major
1396 ✗ jacFT[i * cols + j] = jacF[i + j * rows];
1397 }
1398 }
1399 ✗ }
1400
1401
1402 /*
1403 * Matrix Multiplication using dgemm LaPack routine
1404 */
1405 ✗ void solveMatrixMultiplication(double *matrixA, double *matrixB, int rowsa, int colsa, int rowsb, int colsb, double *matrixC, ofstream &logfile, DATA * data)
1406 {
1407 ✗ char trans = 'N';
1408 ✗ double one = 1.0, zero = 0.0;
1409 ✗ int rowsA = rowsa;
1410 int colsA = colsa;
1411 int rowsB = rowsb;
1412 ✗ int colsB = colsb;
1413 ✗ int common = colsa;
1414
1415 ✗ if (colsA != rowsB)
1416 {
1417 //cout << "\n Error: Column of First Matrix not equal to Rows of Second Matrix \n ";
1418 ✗ errorStreamPrint(OMC_LOG_STDOUT, 0, "solveMatrixMultiplication() Failed!, Column of First Matrix not equal to Rows of Second Matrix %i != %i.",colsA,rowsB);
1419 ✗ logfile << "| error | " << "solveMatrixMultiplication() Failed!, Column of First Matrix not equal to Rows of Second Matrix " << colsA << " != " << rowsB << "\n";
1420 ✗ logfile.close();
1421 ✗ createErrorHtmlReport(data);
1422 ✗ exit(1);
1423 }
1424 // solve matrix multiplication using dgemm_ LAPACK routine
1425 ✗ dgemm_(&trans, &trans, &rowsA, &colsB, &common, &one, matrixA, &rowsA, matrixB, &common, &zero, matrixC, &rowsA);
1426 ✗ }
1427
1428 /*
1429 * Solve the Linear System A*x=b using LAPACK Solver routine dgesv_
1430 */
1431 ✗ void solveSystemFstar(int n, int nhrs, double *tmpMatrixD, double *tmpMatrixC, ofstream &logfile, DATA * data)
1432 {
1433 ✗ int N = n; // number of rows of Matrix A
1434 ✗ int NRHS = nhrs; // number of columns of Matrix B
1435 ✗ int LDA = N;
1436 ✗ int LDB = N;
1437 ✗ int* ipiv = new int[N];
1438 int info;
1439 // call the external function
1440 ✗ dgesv_(&N, &NRHS, tmpMatrixD, &LDA, ipiv, tmpMatrixC, &LDB, &info);
1441 ✗ delete[] ipiv;
1442
1443 ✗ if (info > 0)
1444 {
1445 //cout << "The solution could not be computed, The info satus is : " << info;
1446 ✗ errorStreamPrint(OMC_LOG_STDOUT, 0, "solveSystemFstar() Failed !, The solution could not be computed, The info satus is %i", info);
1447 ✗ logfile << "| error | " << "solveSystemFstar() Failed !, The solution could not be computed, The info satus is " << info << "\n";
1448 ✗ logfile.close();
1449 ✗ createErrorHtmlReport(data);
1450 ✗ exit(1);
1451 }
1452 ✗ }
1453
1454 /*
1455 * Solve the matrix Subtraction of two matrices
1456 */
1457 ✗ void solveMatrixSubtraction(matrixData A, matrixData B, double *result, ofstream &logfile, DATA * data)
1458 {
1459 ✗ if (A.rows != B.rows && A.column != B.column)
1460 {
1461 //cout << "The Matrix Dimensions are not equal to Compute ! \n";
1462 ✗ errorStreamPrint(OMC_LOG_STDOUT, 0, "solveMatrixSubtraction() Failed !, The Matrix Dimensions are not equal to Compute ! %i != %i.", A.rows,B.rows);
1463 ✗ logfile << "| error | " << "solveMatrixSubtraction() Failed !, The Matrix Dimensions are not equal to Compute" << A.rows << " != " << B.rows << "\n";
1464 ✗ logfile.close();
1465 ✗ createErrorHtmlReport(data);
1466 ✗ exit(1);
1467 }
1468
1469 //printColumnAlginment(A.data,A.rows,A.column,"A-Matrix");
1470 //printColumnAlginment(B.data,B.rows,B.column,"B-Matrix");
1471
1472 // subtract elements in cloumn major
1473 ✗ for (int i = 0; i < A.rows * A.column; i++)
1474 {
1475 ✗ result[i] = A.data[i] - B.data[i];
1476 }
1477 ✗ }
1478
1479 /*
1480 * Solve the matrix addition of two matrices
1481 */
1482 ✗ matrixData solveMatrixAddition(matrixData A, matrixData B, ofstream &logfile, DATA * data)
1483 {
1484 ✗ double *result = (double*) calloc(A.rows * A.column, sizeof(double));
1485 ✗ if (A.rows != B.rows && A.column != B.column)
1486 {
1487 //cout << "The Matrix Dimensions are not equal to Compute ! \n";
1488 ✗ errorStreamPrint(OMC_LOG_STDOUT, 0, "solveMatrixAddition() Failed !, The Matrix Dimensions are not equal to Compute ! %i != %i.", A.rows,B.rows);
1489 ✗ logfile << "| error | " << "solveMatrixAddition() Failed !, The Matrix Dimensions are not equal to Compute" << A.rows << " != " << B.rows << "\n";
1490 ✗ logfile.close();
1491 ✗ createErrorHtmlReport(data);
1492 ✗ exit(1);
1493 }
1494
1495 //printColumnAlginment(A.data,A.rows,A.column,"A-Matrix");
1496 //printColumnAlginment(B.data,B.rows,B.column,"B-Matrix");
1497
1498 // Add the elements in cloumn major
1499 ✗ for (int i = 0; i < A.rows * A.column; i++)
1500 {
1501 ✗ result[i] = A.data[i] + B.data[i];
1502 }
1503 matrixData tmpadd_a_b = {A.rows, A.column, result};
1504 ✗ return tmpadd_a_b;
1505 }
1506
1507 /*
1508 * Function which Calculates the Matrix Multiplication
1509 * of (Sx*Ft)*Fstar
1510 */
1511 ✗ matrixData Calculate_Sx_Ft_Fstar(matrixData Sx, matrixData Ft, matrixData Fstar, ofstream &logfile, DATA * data)
1512 {
1513 // Sx*Ft
1514 ✗ double *tmpMatrixA = (double*) calloc(Sx.rows * Ft.column, sizeof(double));
1515 ✗ solveMatrixMultiplication(Sx.data, Ft.data, Sx.rows, Sx.column, Ft.rows, Ft.column, tmpMatrixA, logfile, data);
1516 //printMatrix1(tmpMatrixA,Sx.rows,Ft.column,"Reconciled-(Sx*Ft)");
1517 //printMatrix1(Fstar.data,Fstar.rows,Fstar.column,"REconciled-FStar");
1518
1519 //(Sx*Ft)*Fstar
1520 ✗ double *tmpMatrixB = (double*) calloc(Sx.rows * Fstar.column, sizeof(double));
1521 ✗ solveMatrixMultiplication(tmpMatrixA, Fstar.data, Sx.rows, Ft.column, Fstar.rows, Fstar.column, tmpMatrixB, logfile, data);
1522 matrixData rhsdata = {Sx.rows, Fstar.column, tmpMatrixB};
1523
1524 ✗ free(tmpMatrixA);
1525 ✗ free(tmpMatrixB);
1526 ✗ return rhsdata;
1527 }
1528
1529 /*
1530 * Solves the system
1531 * recon_x = x - (Sx*Ft*fstar)
1532 */
1533 ✗ matrixData solveReconciledX(matrixData x, matrixData Sx, matrixData Ft, matrixData Fstar, ofstream &logfile, DATA * data)
1534 {
1535 // Sx*Ft
1536 ✗ double *tmpMatrixAf = (double*) calloc(Sx.rows * Ft.column, sizeof(double));
1537 ✗ solveMatrixMultiplication(Sx.data, Ft.data, Sx.rows, Sx.column, Ft.rows, Ft.column, tmpMatrixAf, logfile, data);
1538 //printMatrix(tmpMatrixAf,Sx.rows,Ft.column,"Sx*Ft");
1539
1540 //(Sx*Ft)*fstar
1541 ✗ double *tmpMatrixBf = (double*) calloc(Sx.rows * Fstar.column, sizeof(double));
1542 ✗ solveMatrixMultiplication(tmpMatrixAf, Fstar.data, Sx.rows, Ft.column, Fstar.rows, Fstar.column, tmpMatrixBf, logfile, data);
1543 //printMatrix(tmpMatrixBf,Sx.rows,Fstar.column,"(Sx*Ft*fstar)");
1544
1545 ✗ matrixData rhs = {Sx.rows, Fstar.column, tmpMatrixBf};
1546 //matrixData rhs = Calculate_Sx_Ft_Fstar(Sx, Ft, Fstar, data);
1547
1548 ✗ double *reconciledX = (double*) calloc(x.rows * x.column, sizeof(double));
1549 ✗ solveMatrixSubtraction(x, rhs, reconciledX, logfile, data);
1550 //printMatrix(reconciledX,x.rows,x.column,"reconciled X^cap ===> (x - (Sx*Ft*fstar))");
1551
1552 ✗ if (OMC_ACTIVE_STREAM(OMC_LOG_JAC))
1553 {
1554 ✗ logfile << "Calculations of Reconciled_x ==> (x - (Sx*Ft*f*))" << "\n";
1555 ✗ logfile << "====================================================";
1556 ✗ printMatrix(tmpMatrixAf, Sx.rows, Ft.column, "Sx*Ft", logfile);
1557 ✗ printMatrix(tmpMatrixBf, Sx.rows, Fstar.column, "(Sx*Ft*f*)", logfile);
1558 ✗ printMatrix(reconciledX, x.rows, x.column, "x - (Sx*Ft*f*))", logfile);
1559 ✗ logfile << "***** Completed ****** \n\n";
1560 }
1561 matrixData recon_x = {x.rows, x.column, reconciledX};
1562 //free(reconciledX);
1563 ✗ free(tmpMatrixAf);
1564 ✗ free(tmpMatrixBf);
1565 ✗ return recon_x;
1566 }
1567
1568 /*
1569 * Solves the system
1570 * recon_Sx = Sx - (Sx*Ft*Fstar)
1571 */
1572 ✗ matrixData solveReconciledSx(matrixData Sx, matrixData Ft, matrixData Fstar, ofstream &logfile, DATA * data)
1573 {
1574 // Sx*Ft
1575 ✗ double *tmpMatrixA = (double*) calloc(Sx.rows * Ft.column, sizeof(double));
1576 ✗ solveMatrixMultiplication(Sx.data, Ft.data, Sx.rows, Sx.column, Ft.rows, Ft.column, tmpMatrixA, logfile, data);
1577 //printMatrix(tmpMatrixA,Sx.rows,Ft.column,"Reconciled-(Sx*Ft)");
1578 //printMatrix(Fstar.data,Fstar.rows,Fstar.column,"REconciled-FStar");
1579
1580 //(Sx*Ft)*Fstar
1581 ✗ double *tmpMatrixB = (double*) calloc(Sx.rows * Fstar.column, sizeof(double));
1582 ✗ solveMatrixMultiplication(tmpMatrixA, Fstar.data, Sx.rows, Ft.column, Fstar.rows, Fstar.column, tmpMatrixB, logfile, data);
1583 //printMatrix(tmpMatrixB,Sx.rows,Fstar.column,"Reconciled-(Sx*Ft*Fstar)");
1584
1585 ✗ matrixData rhs = {Sx.rows, Fstar.column, tmpMatrixB};
1586 //matrixData rhs = Calculate_Sx_Ft_Fstar(Sx, Ft, Fstar, data);
1587
1588 ✗ double *reconciledSx = (double*) calloc(Sx.rows * Sx.column, sizeof(double));
1589 ✗ solveMatrixSubtraction(Sx, rhs, reconciledSx, logfile, data);
1590 //printMatrix(reconciledSx,Sx.rows,Sx.column,"reconciled Sx ===> (Sx - (Sx*Ft*Fstar))");
1591
1592 ✗ if (OMC_ACTIVE_STREAM(OMC_LOG_JAC))
1593 {
1594 ✗ logfile << "Calculations of Reconciled_Sx ===> (Sx - (Sx*Ft*F*))" << "\n";
1595 ✗ logfile << "============================================";
1596 ✗ printMatrix(tmpMatrixA, Sx.rows, Ft.column, "(Sx*Ft)", logfile);
1597 ✗ printMatrix(tmpMatrixB, Sx.rows, Fstar.column, "(Sx*Ft*F*)", logfile);
1598 ✗ printMatrix(reconciledSx, Sx.rows, Sx.column, "Sx - (Sx*Ft*F*))", logfile);
1599 ✗ logfile << "***** Completed ****** \n\n";
1600 }
1601 matrixData recon_sx = {Sx.rows, Sx.column, reconciledSx};
1602 //free(reconciledSx);
1603 ✗ free(tmpMatrixA);
1604 ✗ free(tmpMatrixB);
1605 ✗ return recon_sx;
1606 }
1607
1608 /*
1609 * Function Which Computes the
1610 * Jacobian Matrix F
1611 */
1612 ✗ matrixData getJacobianMatrixF(DATA * data, threadData_t * threadData, ofstream & logfile, bool boundaryConditions = false)
1613 {
1614 // initialize the jacobian call
1615 ✗ const int index = data->callback->INDEX_JAC_F;
1616 ✗ JACOBIAN *jacobian = &(data->simulationInfo->analyticJacobians[index]);
1617 ✗ data->callback->initialAnalyticJacobianF(data, threadData, jacobian);
1618 ✗ int cols = jacobian->sizeCols;
1619 ✗ int rows = jacobian->sizeRows;
1620 ✗ if (cols == 0)
1621 {
1622 ✗ errorStreamPrint(OMC_LOG_STDOUT, 0, "Cannot Compute Jacobian Matrix F");
1623 ✗ logfile << "| error | " << "Cannot Compute Jacobian Matrix F" << "\n";
1624 ✗ logfile.close();
1625 ✗ if (!boundaryConditions)
1626 {
1627 ✗ createErrorHtmlReport(data);
1628 }
1629 else
1630 {
1631 ✗ createErrorHtmlReportForBoundaryConditions(data);
1632 }
1633 ✗ exit(1);
1634 }
1635 ✗ double *jacF = (double*) calloc(rows * cols, sizeof(double)); // allocate for Matrix F
1636 int k = 0;
1637 ✗ for (int x = 0; x < cols; x++)
1638 {
1639 ✗ jacobian->seedVars[x] = 1.0;
1640 ✗ data->callback->functionJacF_column(data, threadData, jacobian, NULL);
1641 //cout << "Calculate one column\n:";
1642 ✗ for (int y = 0; y < rows; y++)
1643 {
1644 ✗ jacF[k++] = jacobian->resultVars[y];
1645 }
1646 ✗ jacobian->seedVars[x] = 0.0;
1647 }
1648 matrixData Fdata = {rows, cols, jacF};
1649 // free the allocated seed and result variables in jacobian struct after computing the jacobian matrix F to avoid memory leak
1650 ✗ freeJacobian(jacobian);
1651 ✗ return Fdata;
1652 }
1653
1654 /*
1655 * Function Which Computes the
1656 * Jacobian Matrix H
1657 */
1658 ✗ matrixData getJacobianMatrixH(DATA * data, threadData_t * threadData, ofstream & logfile, bool boundaryConditions = false)
1659 {
1660 // initialize the jacobian call
1661 ✗ const int index = data->callback->INDEX_JAC_H;
1662 ✗ JACOBIAN *jacobian = &(data->simulationInfo->analyticJacobians[index]);
1663 ✗ data->callback->initialAnalyticJacobianH(data, threadData, jacobian);
1664 ✗ int cols = jacobian->sizeCols;
1665 ✗ int rows = jacobian->sizeRows;
1666 //std::cout << "\n check jacobian H :" << rows << "==>" << cols;
1667 ✗ if (cols == 0)
1668 {
1669 ✗ errorStreamPrint(OMC_LOG_STDOUT, 0, "Cannot Compute Jacobian Matrix H");
1670 ✗ logfile << "| error | " << "Cannot Compute Jacobian Matrix H" << "\n";
1671 ✗ logfile.close();
1672 ✗ if (!boundaryConditions)
1673 {
1674 ✗ createErrorHtmlReport(data);
1675 }
1676 else
1677 {
1678 ✗ createErrorHtmlReportForBoundaryConditions(data);
1679 }
1680 ✗ exit(1);
1681 }
1682 ✗ double *jacF = (double*) calloc(rows * cols, sizeof(double)); // allocate for Matrix F
1683 int k = 0;
1684 ✗ for (int x = 0; x < cols; x++)
1685 {
1686 ✗ jacobian->seedVars[x] = 1.0;
1687 ✗ data->callback->functionJacH_column(data, threadData, jacobian, NULL);
1688 //cout << "Calculate one column\n:";
1689 ✗ for (int y = 0; y < rows; y++)
1690 {
1691 ✗ jacF[k++] = jacobian->resultVars[y];
1692 }
1693 ✗ jacobian->seedVars[x] = 0.0;
1694 }
1695 matrixData Fdata = {rows, cols, jacF};
1696 // free the allocated seed and result variables in jacobian struct after computing the jacobian matrix F to avoid memory leak
1697 ✗ freeJacobian(jacobian);
1698 ✗ return Fdata;
1699 }
1700
1701 /*
1702 * Function Which Computes the
1703 * Transpose of Jacobian Matrix FT
1704 */
1705 ✗ matrixData getTransposeMatrix(matrixData jacF)
1706 {
1707 int rows = jacF.column;
1708 int cols = jacF.rows;
1709 ✗ double *jacFT = (double*) calloc(rows * cols, sizeof(double)); // allocate for Matrix F-transpose
1710 int k = 0;
1711 ✗ for (int i = 0; i < jacF.rows; i++)
1712 {
1713 ✗ for (int j = 0; j < jacF.column; j++)
1714 {
1715 // Perform matrix transpose store the elements in column major
1716 //cout << (i1*jacF.rows+j1) << " index :" << (i1+j1*jacF.rows) << " value is: " << jacF.data[i1+j1*jacF.rows] << "\n";
1717 ✗ jacFT[k++] = jacF.data[i + j * jacF.rows];
1718
1719 }
1720 }
1721 matrixData Ft_data = {rows, cols, jacFT};
1722 ✗ return Ft_data;
1723 }
1724
1725 /*
1726 * Function which reads the vector
1727 * and assign to c pointer arrays
1728 */
1729 ✗ matrixData getCovarianceMatrixSx(csvData Sx_result, DATA *data, threadData_t *threadData)
1730 {
1731 ✗ double *tempSx = (double*) calloc(Sx_result.rowcount * Sx_result.columncount, sizeof(double));
1732 ✗ initColumnMatrix(Sx_result.sxdata, Sx_result.rowcount, Sx_result.columncount, tempSx);
1733 ✗ matrixData Sx_data = {Sx_result.rowcount, Sx_result.columncount, tempSx};
1734 ✗ return Sx_data;
1735 }
1736
1737 /*
1738 * function which validates corelation inputs
1739 * and displays error messages for
1740 * user error #7: variable of interest has multiple entry in measurement input file
1741 * user error #8: variable of interest in correlation input file does not correspond to variable of interest
1742 */
1743 ✗ void validateCorelationInputs(csvData Sx_result, DATA * data, ofstream &logfile, vector<string> headers, string comments, bool boundaryConditions = false)
1744 {
1745 vector<string> noEntry, multipleEntry, entry;
1746 ✗ for (int i = 0; i < headers.size(); i++)
1747 {
1748 bool flag = false;
1749 ✗ for (int j = 0; j < Sx_result.headers.size(); j++)
1750 {
1751 ✗ if (strcmp(headers[i].c_str(), Sx_result.headers[j].c_str()) == 0)
1752 {
1753 //std::cout << "\n matched variables of interest: "<< headers[i] << " => " << Sx_result.headers[j] << " pos: " << j << "\n";
1754 flag = true;
1755 ✗ auto it = find(entry.begin(), entry.end(), headers[i]);
1756 ✗ if (it != entry.end())
1757 {
1758 //std::cout << "Element found in myvector: " << headers[i] << "\n";
1759 ✗ multipleEntry.push_back(headers[i]);
1760 }
1761 else
1762 {
1763 ✗ entry.push_back(headers[i]);
1764 }
1765 }
1766 }
1767
1768 // user error variable of interest in correlation input file does not correspond to variable of interest
1769 ✗ if (!flag)
1770 {
1771 ✗ noEntry.push_back(headers[i]);
1772 }
1773 }
1774
1775 // dump user error #7: variable of interest , has multiple entry in measurement input file
1776 ✗ for (int i = 0; i < multipleEntry.size(); i++)
1777 {
1778 ✗ if (!boundaryConditions)
1779 {
1780 ✗ errorStreamPrint(OMC_LOG_STDOUT, 0, "variable of interest %s, at %s has multiple entries in correlation input file %s ", multipleEntry[i].c_str(), comments.c_str(), omc_flagValue[FLAG_DATA_RECONCILE_Cx]);
1781 ✗ logfile << "| error | " << "variable of interest " << multipleEntry[i] << " at " << comments << " has multiple entries in correlation input file " << omc_flagValue[FLAG_DATA_RECONCILE_Cx] << "\n";
1782 }
1783 else
1784 {
1785 ✗ errorStreamPrint(OMC_LOG_STDOUT, 0, "variable of interest %s, at %s has multiple entries in reconciled covariance matrix input file %s ", multipleEntry[i].c_str(), comments.c_str(), omc_flagValue[FLAG_DATA_RECONCILE_Cx]);
1786 ✗ logfile << "| error | " << "variable of interest " << multipleEntry[i] << " at " << comments << " has multiple entries in reconciled covariance matrix input file " << omc_flagValue[FLAG_DATA_RECONCILE_Cx] << "\n";
1787 }
1788 }
1789
1790 // dump user error #8: variable of interest in correlation input file does not correspond to variable of interest
1791 ✗ for (int i = 0; i < noEntry.size(); i++)
1792 {
1793 ✗ if (!boundaryConditions)
1794 {
1795 ✗ errorStreamPrint(OMC_LOG_STDOUT, 0, "variable of interest %s, at %s entry in correlation input file %s does not correspond to a variable of interest ", noEntry[i].c_str(), comments.c_str(), omc_flagValue[FLAG_DATA_RECONCILE_Cx]);
1796 ✗ logfile << "| error | " << "variable of interest " << noEntry[i] << ", at " << comments << " entry in correlation input file " << omc_flagValue[FLAG_DATA_RECONCILE_Cx] << " does not correspond to a variable of interest" << "\n";
1797 }
1798 else
1799 {
1800 ✗ errorStreamPrint(OMC_LOG_STDOUT, 0, "variable of interest %s, at %s entry in reconciled covariance matrix input file %s does not correspond to a variable of interest ", noEntry[i].c_str(), comments.c_str(), omc_flagValue[FLAG_DATA_RECONCILE_Cx]);
1801 ✗ logfile << "| error | " << "variable of interest " << noEntry[i] << ", at " << comments << " entry in reconciled covariance matrix input file " << omc_flagValue[FLAG_DATA_RECONCILE_Cx] << " does not correspond to a variable of interest" << "\n";
1802 }
1803 }
1804
1805 ✗ if (!noEntry.empty() || !multipleEntry.empty())
1806 {
1807 ✗ logfile.close();
1808 ✗ if (!boundaryConditions)
1809 {
1810 ✗ createErrorHtmlReport(data);
1811 }
1812 else
1813 {
1814 ✗ createErrorHtmlReportForBoundaryConditions(data);
1815 }
1816 ✗ exit(1);
1817 }
1818 ✗ }
1819
1820 /*
1821 * function which validates corelation inputs is
1822 * square matrix or not and displays error messages for
1823 * user error #10: Lines and columns are in different orders
1824 */
1825 ✗ void validateCorelationInputsSquareMatrix(DATA * data, ofstream &logfile, vector<string> rowHeaders, vector<string> columnHeaders, bool boundaryConditions = false)
1826 {
1827 ✗ if (rowHeaders != columnHeaders)
1828 {
1829 ✗ if (!boundaryConditions)
1830 {
1831 ✗ errorStreamPrint(OMC_LOG_STDOUT, 0, "Lines and columns of correlation matrix in correlation input file %s, do not have identical names in the same order.", omc_flagValue[FLAG_DATA_RECONCILE_Cx]);
1832 ✗ logfile << "| error | " << "Lines and columns of correlation matrix in correlation input file " << omc_flagValue[FLAG_DATA_RECONCILE_Cx] << " do not have identical names in the same order." << "\n";
1833 }
1834 else
1835 {
1836 ✗ errorStreamPrint(OMC_LOG_STDOUT, 0, "Lines and columns of covariance matrix in reconciled covariance matrix input file %s, do not have identical names in the same order.", omc_flagValue[FLAG_DATA_RECONCILE_Cx]);
1837 ✗ logfile << "| error | " << "Lines and columns of covariance matrix in reconciled covariance matrix input file " << omc_flagValue[FLAG_DATA_RECONCILE_Cx] << " do not have identical names in the same order." << "\n";
1838 }
1839
1840 // user error #10: missing line headers
1841 ✗ for (const auto & index : columnHeaders)
1842 {
1843 ✗ auto it = std::find(rowHeaders.begin(), rowHeaders.end(), index);
1844 ✗ if (it == rowHeaders.end())
1845 {
1846 ✗ errorStreamPrint(OMC_LOG_STDOUT, 0, "Line %s is missing", index.c_str());
1847 ✗ logfile << "| error | " << "Line " << index << " is missing " << "\n";
1848 }
1849 }
1850
1851 // user error #10 : missing column headers
1852 ✗ for (const auto & index : rowHeaders)
1853 {
1854 ✗ auto it = find(columnHeaders.begin(), columnHeaders.end(), index);
1855 ✗ if (it == columnHeaders.end())
1856 {
1857 ✗ errorStreamPrint(OMC_LOG_STDOUT, 0, "Column %s is missing", index.c_str());
1858 ✗ logfile << "| error | " << "Column " << index << " is missing " << "\n";
1859 }
1860 }
1861
1862 // user error #10 : missing line Vs column
1863 ✗ for (int i = 0; i < rowHeaders.size(); i++)
1864 {
1865 //std::cout << "\n " << i << rowHeaders[i] << " Vs " << columnHeaders[i];
1866 ✗ if (rowHeaders[i] != columnHeaders[i])
1867 {
1868 ✗ errorStreamPrint(OMC_LOG_STDOUT, 0, "Lines and columns are in different orders %s Vs %s", rowHeaders[i].c_str(), columnHeaders[i].c_str());
1869 ✗ logfile << "| error | " << "Lines and columns are in different orders " << rowHeaders[i] << " Vs " << columnHeaders[i] << "\n";
1870 }
1871 }
1872
1873 ✗ logfile.close();
1874 ✗ if (!boundaryConditions)
1875 {
1876 ✗ createErrorHtmlReport(data);
1877 }
1878 else
1879 {
1880 ✗ createErrorHtmlReportForBoundaryConditions(data);
1881 }
1882 ✗ exit(1);
1883 }
1884 ✗ }
1885
1886 /*
1887 * Function which reads the correlation coefficient input file
1888 * and stores the correlation coefficient matrix Cx for DataReconciliation
1889 */
1890 ✗ correlationData readCorrelationCoefficientFile(csvData Sx_result, ofstream & logfile, DATA * data, threadData_t * threadData, correlationDataWarning & warningCorrelationData, bool boundaryConditions = false)
1891 {
1892 char * filename = NULL;
1893 ✗ filename = (char*) omc_flagValue[FLAG_DATA_RECONCILE_Cx];
1894
1895 vector<string> columnHeaders, rowHeaders;
1896 vector<double> cx_data;
1897 vector<errorData> errorInfo;
1898 correlationData Cx;
1899 vector<int> errorInfoHeaders;
1900
1901 ✗ if (filename == NULL && !boundaryConditions)
1902 {
1903 ✗ Cx = {cx_data, rowHeaders, columnHeaders};
1904 ✗ return Cx;
1905 }
1906
1907 ✗ if (filename == NULL && boundaryConditions)
1908 {
1909 ✗ errorStreamPrint(OMC_LOG_STDOUT, 0, "Reconciled covariance matrix input file not provided (eg:-cx=filename.csv), Boundary conditions cannot be computed!.");
1910 ✗ logfile << "| error | " << "Reconciled covariance matrix input file not provided (eg:-cx=filename.csv), Boundary conditions cannot be computed!.\n";
1911 ✗ logfile.close();
1912 ✗ createErrorHtmlReportForBoundaryConditions(data);
1913 ✗ exit(1);
1914 }
1915
1916 /*
1917 * fix issue https://github.com/OpenModelica/OpenModelica/issues/13797
1918 * check if filepath uses uri format and convert it to absolute path
1919 * eg: modelica://Modelica/Resources/Files/filename.csv => /absolute/path/to/filename.csv
1920 */
1921 ✗ modelica_string uri = OpenModelica_uriToFilename(omc_string_new(filename));
1922 ✗ std::string filenameStr = omc_string_data(uri);
1923
1924 // read the file
1925 ✗ ifstream ip(filenameStr.c_str());
1926 string line;
1927 ✗ if (!ip.good() && !boundaryConditions)
1928 {
1929 ✗ errorStreamPrint(OMC_LOG_STDOUT, 0, "correlation coefficient input file path not found %s.", filename);
1930 ✗ logfile << "| error | " << "correlation coefficient input file path not found " << filename << "\n";
1931 ✗ logfile.close();
1932 ✗ createErrorHtmlReport(data);
1933 ✗ exit(1);
1934 }
1935
1936 ✗ if (!ip.good() && boundaryConditions)
1937 {
1938 ✗ errorStreamPrint(OMC_LOG_STDOUT, 0, "Reconciled covariance matrix input file path not found %s.", filename);
1939 ✗ logfile << "| error | " << "Reconciled covariance matrix input file path not found " << filename << "\n";
1940 ✗ logfile.close();
1941 ✗ createErrorHtmlReportForBoundaryConditions(data);
1942 ✗ exit(1);
1943 }
1944
1945 ✗ int linecount = 1;
1946 ✗ while (ip.good())
1947 {
1948 ✗ getline(ip, line);
1949 vector<string> t1;
1950 std::replace(line.begin(), line.end(), ';', ',');
1951 // remove whitespace in a string
1952 ✗ line.erase(std::remove_if(line.begin(), line.end(), ::isspace), line.end());
1953 ✗ stringstream ss(line);
1954 string temp;
1955
1956 // allow comments on the line#1 until finding the headers
1957 ✗ if (linecount == 1 && isLineEmptyData(line))
1958 {
1959 continue;
1960 }
1961
1962 ✗ if (linecount == 1 && !line.empty())
1963 {
1964 int columnCount = 1;
1965 ✗ while (getline(ss, temp, ','))
1966 {
1967 ✗ if (columnCount > 1)
1968 {
1969 //std::cout << "\nreading Column headers : " << temp;
1970 ✗ columnHeaders.push_back(temp);
1971 }
1972 ✗ columnCount++;
1973 }
1974 }
1975 ✗ else if (linecount > 1 && !line.empty() && !isLineEmptyData(line))
1976 {
1977 //std::cout << "\nreading covariance matrix : " << line << " size : " << columnHeaders.size();
1978 int columnCount = 1;
1979 ✗ while (getline(ss, temp, ','))
1980 {
1981 ✗ if (columnCount == 1)
1982 {
1983 // error : no variable of interest is provided by user at column #1
1984 ✗ if (temp.empty())
1985 {
1986 ✗ errorInfoHeaders.push_back(linecount);
1987 }
1988 ✗ rowHeaders.push_back(temp);
1989 }
1990 else
1991 {
1992 //std::cout << "\nreading Column values : " << temp << " : " << columnCount << " value= " << atof(temp.c_str());
1993 ✗ if (temp.empty())
1994 {
1995 ✗ cx_data.push_back(0);
1996 }
1997 ✗ else if (!isStringValidDouble(temp))
1998 {
1999 //std::cout << "\n check wrong type : " << columnCount << " = " << columnHeaders[columnCount-2];
2000 ✗ errorData info = {rowHeaders.back(), columnHeaders[columnCount-2], temp};
2001 ✗ errorInfo.push_back(info);
2002 ✗ }
2003 // warn users if value is closer to 1
2004 ✗ else if (atof(temp.c_str()) > 0.99 && atof(temp.c_str()) < 1.01)
2005 {
2006 ✗ errorData info = {rowHeaders.back(), columnHeaders[columnCount-2], temp};
2007 ✗ warningCorrelationData.warningInfo.push_back(info);
2008 ✗ cx_data.push_back(atof(temp.c_str()));
2009 ✗ }
2010 else
2011 {
2012 ✗ cx_data.push_back(atof(temp.c_str()));
2013 }
2014 }
2015 ✗ columnCount++;
2016 }
2017 ✗ if (columnCount - 2 == columnHeaders.size())
2018 {
2019 // rows and columns have values correctly entered
2020 }
2021 else
2022 {
2023 // just in case, fill the empty rows with zeros
2024 ✗ int count = columnHeaders.size() - (columnCount -2);
2025 ✗ for (int i = 0; i < count; ++i)
2026 {
2027 ✗ cx_data.push_back(0);
2028 }
2029 }
2030 }
2031 ✗ linecount++;
2032 ✗ }
2033
2034 // user error : variable of interest is missing in column #1
2035 ✗ if (!errorInfoHeaders.empty())
2036 {
2037 ✗ for (const auto &line : errorInfoHeaders)
2038 {
2039 ✗ if (!boundaryConditions)
2040 {
2041 ✗ errorStreamPrint(OMC_LOG_STDOUT, 0, "the name of the variable of interest in correlation input file %s is missing in line #%d ", filename, line);
2042 ✗ logfile << "| error | " << "the name of the variable of interest in correlation input file " << filename << " is missing in line #" << line << "\n";
2043 }
2044 else
2045 {
2046 ✗ errorStreamPrint(OMC_LOG_STDOUT, 0, "the name of the variable of interest in reconciled covariance matrix input file %s is missing in line #%d ", filename, line);
2047 ✗ logfile << "| error | " << "the name of the variable of interest in reconciled covariance matrix input file " << filename << " is missing in line #" << line << "\n";
2048 }
2049 }
2050 ✗ logfile.close();
2051 ✗ if (!boundaryConditions)
2052 {
2053 ✗ createErrorHtmlReport(data);
2054 }
2055 else
2056 {
2057 ✗ createErrorHtmlReportForBoundaryConditions(data);
2058 }
2059 ✗ exit(1);
2060 }
2061
2062 // user error #9: Entry for variable of interest <variable name> and variable of interest <variable name> in correlation input file <input file name> is incorrect: incorrect value type.
2063 ✗ if (!errorInfo.empty())
2064 {
2065 ✗ for (const auto & info : errorInfo)
2066 {
2067 ✗ if (!boundaryConditions)
2068 {
2069 ✗ errorStreamPrint(OMC_LOG_STDOUT, 0, "Entry for variable of interest %s and variable of interest %s in correlation input file %s is incorrect because of wrong-Type: [%s] ", info.name.c_str(), info.x.c_str(), filename, info.sx.c_str());
2070 ✗ logfile << "| error | " << "Entry for variable of interest " << info.name << " and variable of interest " << info.x << " in correlation input file " << filename << " is incorrect because of wrong-Type: " << "[" << info.sx << "]" <<"\n";
2071 }
2072 else
2073 {
2074 ✗ errorStreamPrint(OMC_LOG_STDOUT, 0, "Entry for variable of interest %s and variable of interest %s in reconciled covariance matrix input file %s is incorrect because of wrong-Type: [%s] ", info.name.c_str(), info.x.c_str(), filename, info.sx.c_str());
2075 ✗ logfile << "| error | " << "Entry for variable of interest " << info.name << " and variable of interest " << info.x << " in reconciled covariance matrix input file " << filename << " is incorrect because of wrong-Type: " << "[" << info.sx << "]" <<"\n";
2076 }
2077 }
2078 ✗ logfile.close();
2079 ✗ if (!boundaryConditions)
2080 {
2081 ✗ createErrorHtmlReport(data);
2082 }
2083 else
2084 {
2085 ✗ createErrorHtmlReportForBoundaryConditions(data);
2086 }
2087 ✗ exit(1);
2088 }
2089
2090 // user warning : Entry for variable of interest <variable name> and variable of interest <variable name> in correlation input file <input file name> is closer to 1.
2091 ✗ if (!warningCorrelationData.warningInfo.empty())
2092 {
2093 ✗ for (const auto & info : warningCorrelationData.warningInfo)
2094 {
2095 ✗ warningStreamPrint(OMC_LOG_STDOUT, 0, "Entry for variable of interest %s and variable of interest %s in correlation input file %s is closer to 1: [%s] ", info.name.c_str(), info.x.c_str(), filename, info.sx.c_str());
2096 }
2097 }
2098
2099 // validate correlation input column headers
2100 ✗ validateCorelationInputs(Sx_result, data, logfile, columnHeaders, "column headers", boundaryConditions);
2101
2102 // validate correlation input column headers
2103 ✗ validateCorelationInputs(Sx_result, data, logfile, rowHeaders, "row headers", boundaryConditions);
2104
2105 // check for square matrix
2106 ✗ validateCorelationInputsSquareMatrix(data, logfile, rowHeaders, columnHeaders, boundaryConditions);
2107
2108 ✗ Cx = {cx_data, rowHeaders, columnHeaders};
2109
2110 return Cx;
2111 ✗ }
2112
2113
2114 /*
2115 * Function which Computes
2116 * covariance matrix Sx based on
2117 * Half width confidence interval provided by user (i.e) csvData.sxdata
2118 * Sx=(Wxi/1.96)^2
2119 */
2120 ✗ matrixData computeCovarianceMatrixSx(csvData Sx_result, correlationData Cx_data, ofstream &logfile, DATA * data)
2121 {
2122 ✗ double *tempSx = (double*) calloc(Sx_result.sxdata.size() * Sx_result.sxdata.size(), sizeof(double));
2123 vector<double> tmpdata;
2124 int k = 0;
2125
2126 ✗ for (unsigned int i = 0; i < Sx_result.sxdata.size(); i++)
2127 {
2128 ✗ double data = pow(Sx_result.sxdata[k] / 1.96, 2);
2129 ✗ for (unsigned int j = 0; j < Sx_result.sxdata.size(); j++)
2130 {
2131 ✗ if (i == j)
2132 {
2133 //tmpdata.push_back(pow(Sx_result.sxdata[k]/1.96,2));
2134 //k++;
2135 ✗ tmpdata.push_back(data);
2136 }
2137 else
2138 {
2139 ✗ tmpdata.push_back(0);
2140 }
2141 // logfile << " data " << count << "=="<< tmpdata[count++] << "\n";
2142 }
2143 ✗ k++;
2144 }
2145
2146 // check for correlation coefficient Cx_data is not empty and recompute the covariance matrix
2147 ✗ if (! Cx_data.data.empty())
2148 {
2149 ✗ for (int i = 0; i < Cx_data.rowHeaders.size(); i++)
2150 {
2151 ✗ for (int j = 0; j < Cx_data.columnHeaders.size(); j++)
2152 {
2153 // consider the values which are strictly below the diagonal entry
2154 ✗ if (j < i && Cx_data.data[Cx_data.columnHeaders.size() * i + j] != 0)
2155 {
2156 //std::cout << "\n value : " << Cx_data.rowHeaders[i] << "=" << Cx_data.data[columnHeaders.size() * i + j];
2157 ✗ int rowpos = getVariableIndex(Sx_result.headers, Cx_data.rowHeaders[i], logfile, data);
2158 ✗ int colpos = getVariableIndex(Sx_result.headers, Cx_data.columnHeaders[j], logfile, data);
2159
2160 ✗ double xi = tmpdata[(Sx_result.rowcount * rowpos) + rowpos];
2161 ✗ double xk = tmpdata[(Sx_result.rowcount * colpos) + colpos];
2162
2163 //std::cout << "\n row pos : " << rowpos << " col pos: " << colpos << " xi " << xi;
2164 ✗ double tmprx = Cx_data.data[Cx_data.columnHeaders.size() * i + j] * sqrt(xi) * sqrt(xk);
2165
2166 // find the symmetric position and insert the elements
2167 ✗ tmpdata[(Sx_result.rowcount * rowpos) + colpos] = tmprx;
2168 ✗ tmpdata[(Sx_result.rowcount * colpos) + rowpos] = tmprx;
2169 }
2170 }
2171 }
2172 }
2173
2174 ✗ initColumnMatrix(tmpdata, Sx_result.rowcount, Sx_result.rowcount, tempSx);
2175 ✗ matrixData Sx_data = {Sx_result.rowcount, Sx_result.rowcount, tempSx};
2176 ✗ return Sx_data;
2177 }
2178
2179 /*
2180 * function which validates the inputs read from measurement input file
2181 * and validates (i.e) all the variables of interest in input file = all variables of interest in model
2182 * and displays error message if the above condition fails.
2183 * Also it perform internal mapping to sort the variable of interest in input file to match with variables of interest in model
2184 * which is very important for jacobians and also to get correct numerical results
2185 */
2186
2187 ✗ csvData validateMeasurementInputs(csvData Sx_result, DATA * data, ofstream &logfile, bool boundaryConditions = false)
2188 {
2189 ✗ if (data->modelData->ndataReconVars != Sx_result.headers.size())
2190 {
2191 ✗ errorStreamPrint(OMC_LOG_STDOUT, 0, "invalid input file %s, number of variable of interest(%li) != (%zu)number of variables in measurement input file", omc_flagValue[FLAG_DATA_RECONCILE_Sx], data->modelData->ndataReconVars, Sx_result.headers.size());
2192 ✗ logfile << "| error | " << "invalid input file "<< omc_flagValue[FLAG_DATA_RECONCILE_Sx] << ", number of variable of interest(" << data->modelData->ndataReconVars << ")" << " != " << "(" << Sx_result.headers.size() << ")" << "number of variables in measurement input file";
2193 ✗ logfile.close();
2194 ✗ if (!boundaryConditions)
2195 {
2196 ✗ createErrorHtmlReport(data);
2197 }
2198 else
2199 {
2200 ✗ createErrorHtmlReportForBoundaryConditions(data);
2201 }
2202 ✗ exit(1);
2203 }
2204
2205 ✗ char **knowns = (char**) malloc(data->modelData->ndataReconVars * sizeof(char*));
2206 ✗ data->callback->dataReconciliationInputNames(data, knowns);
2207
2208 vector<string> noEntry, multipleEntry;
2209 vector<int> mapindex;
2210 ✗ for (int i = 0; i < data->modelData->ndataReconVars; i++)
2211 {
2212 bool flag = false;
2213 int count = 0;
2214 ✗ for (int j = 0; j < Sx_result.headers.size(); j++)
2215 {
2216 ✗ if (strcmp(knowns[i], Sx_result.headers[j].c_str()) == 0)
2217 {
2218 //std::cout << "\n matched variables of interest: "<< knowns[i] << " => " << Sx_result.headers[j] << " pos: " << j << " size : " << Sx_result.headers.size() << "\n";
2219 ✗ mapindex.push_back(j);
2220 flag = true;
2221 ✗ count ++;
2222 }
2223 }
2224
2225 // user error variable of interest , has no entry in measurement input file
2226 ✗ if (!flag)
2227 {
2228 ✗ noEntry.push_back(knowns[i]);
2229 }
2230 // user error variable of interest , has multiple entry in measurement input file
2231 ✗ if (count > 1)
2232 {
2233 ✗ multipleEntry.push_back(knowns[i]);
2234 }
2235 }
2236
2237 // dump user error #3: variable of interest, has no entry in measurement input file
2238 ✗ for (int i = 0; i < noEntry.size(); i++)
2239 {
2240 ✗ errorStreamPrint(OMC_LOG_STDOUT, 0, "variable of interest %s, has no entry in measurement input file %s ", noEntry[i].c_str(), omc_flagValue[FLAG_DATA_RECONCILE_Sx]);
2241 ✗ logfile << "| error | " << "variable of interest " << noEntry[i] << ", has no entry in measurement input file" << omc_flagValue[FLAG_DATA_RECONCILE_Sx] << "\n";
2242 }
2243
2244 // dump user error #4: variable of interest , has multiple entry in measurement input file
2245 ✗ for (int i = 0; i < multipleEntry.size(); i++)
2246 {
2247 ✗ errorStreamPrint(OMC_LOG_STDOUT, 0, "variable of interest %s, has multiple entries in measurement input file %s ", multipleEntry[i].c_str(), omc_flagValue[FLAG_DATA_RECONCILE_Sx]);
2248 ✗ logfile << "| error | " << "variable of interest " << multipleEntry[i] << ", has multiple entries in measurement input file " << omc_flagValue[FLAG_DATA_RECONCILE_Sx] << "\n";
2249 }
2250
2251 // dump user error #5 entry in measurement input file does not correspond to a variable of interest
2252 bool userError5 = false;
2253 ✗ for (int i = 0; i < data->modelData->ndataReconVars; i++)
2254 {
2255 ✗ auto it = find(mapindex.begin(), mapindex.end(), i);
2256 ✗ if (it == mapindex.end())
2257 {
2258 //std::cout << "\n not found" << i;
2259 userError5 = true;
2260 ✗ errorStreamPrint(OMC_LOG_STDOUT, 0, "variable of interest %s, entry in measurement input file %s does not correspond to a variable of interest ", Sx_result.headers[i].c_str(), omc_flagValue[FLAG_DATA_RECONCILE_Sx]);
2261 ✗ logfile << "| error | " << "variable of interest " << Sx_result.headers[i] << ", entry in measurement input file " << omc_flagValue[FLAG_DATA_RECONCILE_Sx] << " does not correspond to a variable of interest" << "\n";
2262 }
2263 }
2264
2265 ✗ if (!noEntry.empty() || !multipleEntry.empty() || userError5)
2266 {
2267 ✗ logfile.close();
2268 ✗ free(knowns);
2269 ✗ if (!boundaryConditions)
2270 {
2271 ✗ createErrorHtmlReport(data);
2272 }
2273 else
2274 {
2275 ✗ createErrorHtmlReportForBoundaryConditions(data);
2276 }
2277 ✗ exit(1);
2278 }
2279
2280 // map csv inputs with order of variable of interest in the model
2281 vector<double> mapped_xdata, mapped_Sxdata;
2282 vector<string> mappedHeader;
2283 ✗ for (const auto & index : mapindex)
2284 {
2285 ✗ mapped_xdata.push_back(Sx_result.xdata[index]);
2286 ✗ mapped_Sxdata.push_back(Sx_result.sxdata[index]);
2287 ✗ mappedHeader.push_back(Sx_result.headers[index]);
2288 }
2289
2290 // assign the mapped order
2291 ✗ Sx_result.xdata = mapped_xdata;
2292 ✗ Sx_result.sxdata = mapped_Sxdata;
2293 ✗ Sx_result.headers = mappedHeader;
2294
2295 ✗ free(knowns);
2296 ✗ return Sx_result;
2297 ✗ }
2298
2299
2300 /*
2301 * Function which reads the input data from csvData.xdata
2302 * and stores the input values in double array
2303 */
2304 ✗ inputData getInputData(csvData Sx_result, ofstream & logfile)
2305 {
2306 ✗ double * tempx = (double*) calloc(Sx_result.rowcount, sizeof(double));
2307 vector<int> index;
2308
2309 ✗ for (int i = 0; i < Sx_result.headers.size(); i++)
2310 {
2311 ✗ tempx[i] = Sx_result.xdata[i];
2312 }
2313
2314 ✗ inputData x_data = {Sx_result.rowcount, 1, tempx, index};
2315 ✗ return x_data;
2316 }
2317
2318 /*
2319 * Function which reads the input data from csvData.xdata
2320 * and stores the input values in double array
2321 * Note: this function takes the output.csv file from data reconciliation
2322 * as input and get the reconciled_X values from third column
2323 */
2324 ✗ inputData getReconciledX(csvData Sx_result, ofstream & logfile)
2325 {
2326 ✗ double * tempx = (double*) calloc(Sx_result.rowcount, sizeof(double));
2327 vector<int> index;
2328
2329 ✗ for (int i = 0; i < Sx_result.headers.size(); i++)
2330 {
2331 ✗ tempx[i] = Sx_result.sxdata[i];
2332 }
2333
2334 ✗ inputData x_data = {Sx_result.rowcount, 1, tempx, index};
2335 ✗ return x_data;
2336 }
2337
2338 /*
2339 * Function which Copy Matrix
2340 * using dcopy_ LAPACK routine
2341 * this is mostly used when LAPACK routines override arrays
2342 */
2343 ✗ matrixData copyMatrix(matrixData matdata)
2344 {
2345 ✗ double *tmpcopymatrix = (double*) calloc(matdata.rows * matdata.column, sizeof(double));
2346 ✗ int n = matdata.rows * matdata.column;
2347 ✗ int inc = 1;
2348 ✗ dcopy_(&n, matdata.data, &inc, tmpcopymatrix, &inc);
2349
2350 // for (int i=0; i < matdata.rows*matdata.column; i++)
2351 // {
2352 // tmpcopymatrix[i]=matdata.data[i];
2353 // }
2354 matrixData tmpcopymatrixdata = {matdata.rows, matdata.column, tmpcopymatrix};
2355 ✗ return tmpcopymatrixdata;
2356 }
2357
2358 /*
2359 * Function which scales the MAtrix with constant
2360 * dscal_ LAPACK_routine and result is updated in data
2361 * eg: alpha = 2, data=[2,4,5]
2362 * result = [4,8,10]
2363 */
2364 ✗ void scaleVector(int rows, int cols, double alpha, double *data)
2365 {
2366 ✗ int n = rows * cols;
2367 ✗ int inc = 1;
2368 ✗ dscal_(&n, &alpha, data, &inc);
2369 ✗ }
2370
2371 /*
2372 * Function which calculates the square root of elements
2373 * eg : a=[1,2,3,4]
2374 * result a = [srt(1),sqrt(2).......]
2375 */
2376 ✗ void calculateSquareRoot(double *data, int length)
2377 {
2378 ✗ for (int i = 0; i < length; i++)
2379 {
2380 ✗ data[i] = sqrt(data[i]);
2381 }
2382 ✗ }
2383
2384 /*
2385 * Function which calculates
2386 * J*=(recon_x-x)T*(Sx^-1)*(recon_x-x)+2.[f+F*(recon_x-x)]T*fstar
2387 * where T= transpose of matrix
2388 * and returns the converged value
2389 */
2390 ✗ double solveConvergence(DATA *data, matrixData conv_recon_x, matrixData conv_recon_sx, inputData conv_x, matrixData conv_sx, matrixData conv_jacF, matrixData conv_vector_c, matrixData conv_fstar, ofstream &logfile)
2391 {
2392
2393 //printMatrix(conv_vector_c.data,conv_vector_c.rows,conv_vector_c.column,"Convergence_C(x,y)");
2394 //printMatrix(conv_fstar.data,conv_fstar.rows,conv_fstar.column,"Convergence_f*");
2395 //printMatrix(conv_recon_x.data,conv_recon_x.rows,conv_recon_x.column,"check_recon_x*");
2396
2397 // calculate(recon_x-x)
2398 ✗ double *conv_data1 = (double*) calloc(conv_x.rows * conv_x.column, sizeof(double));
2399 ✗ matrixData conv_inputs = {conv_x.rows, conv_x.column, conv_x.data};
2400 ✗ solveMatrixSubtraction(conv_recon_x, conv_inputs, conv_data1, logfile, data);
2401 ✗ matrixData conv_data1result = {conv_x.rows, conv_x.column, conv_data1};
2402 ✗ matrixData copy_reconx_x = copyMatrix(conv_data1result);
2403
2404 //printMatrix(conv_inputs.data,conv_inputs.rows,conv_inputs.column,"check_inputs");
2405 //printMatrix(conv_data1result.data,conv_data1result.rows,conv_data1result.column,"(recon_X - X)");
2406
2407 // calculate Transpose_(recon_x-x)
2408 ✗ matrixData conv_data1Transpose = getTransposeMatrix(conv_data1result);
2409 //printMatrix(conv_data1Transpose.data,conv_data1Transpose.rows,conv_data1Transpose.column,"Transpose(recon_X - X)");
2410
2411 /* solves (Sx^-1)*(recon_x-x)
2412 * Solve the inverse of matrix Sx using linear form
2413 * Ax=b
2414 * where A=Sx and b= (recon_x-x) to avoid inversion of Sx which is
2415 * expensive
2416 */
2417 ✗ solveSystemFstar(conv_sx.rows, 1, conv_sx.data, conv_data1result.data, logfile, data);
2418 //printMatrix(conv_data1result.data,conv_sx.rows,conv_data1result.column,"inverse multiplication_without inverse");
2419
2420 ✗ double *conv_tmpmatrixlhs = (double*) calloc(conv_data1Transpose.rows * conv_data1result.column, sizeof(double));
2421 /*
2422 * Solve (recon_x-x)T*(Sx^-1)*(recon_x-x)
2423 */
2424 ✗ solveMatrixMultiplication(conv_data1Transpose.data, conv_data1result.data, conv_data1Transpose.rows, conv_data1Transpose.column, conv_data1result.rows, conv_data1result.column, conv_tmpmatrixlhs, logfile, data);
2425 //printMatrix(conv_tmpmatrixlhs,conv_data1Transpose.rows,conv_data1result.column,"(recon_x-x)T*(Sx^-1)*(recon_x-x)");
2426 ✗ matrixData struct_conv_tmpmatrixlhs = {conv_data1Transpose.rows, conv_data1result.column, conv_tmpmatrixlhs};
2427
2428 /*
2429 * Solve rhs = 2.[f+F*(recon_x-x)]T*fstar
2430 *
2431 */
2432 // Calculate F*(recon_x-x)
2433 ✗ double *tmp_F_recon_x_x = (double*) calloc(conv_jacF.rows * copy_reconx_x.column, sizeof(double));
2434 ✗ solveMatrixMultiplication(conv_jacF.data, copy_reconx_x.data, conv_jacF.rows, conv_jacF.column, copy_reconx_x.rows, copy_reconx_x.column, tmp_F_recon_x_x, logfile, data);
2435 //printMatrix(tmp_F_recon_x_x,conv_jacF.rows,copy_reconx_x.column,"F*(recon_x-x)");
2436 ✗ matrixData mult_F_recon_x_x = {conv_jacF.rows, copy_reconx_x.column, tmp_F_recon_x_x};
2437
2438 // Calculate f + F*(recon_x-x)
2439 ✗ matrixData add_f_F_recon_x_x = solveMatrixAddition(conv_vector_c, mult_F_recon_x_x, logfile, data);
2440 //printMatrix(add_f_F_recon_x_x.data,add_f_F_recon_x_x.rows,add_f_F_recon_x_x.column,"f + F*(recon_x-x)");
2441
2442 ✗ matrixData transpose_add_f_F_recon_x_x = getTransposeMatrix(add_f_F_recon_x_x);
2443 //printMatrix(transpose_add_f_F_recon_x_x.data,transpose_add_f_F_recon_x_x.rows,transpose_add_f_F_recon_x_x.column,"transpose-[f + F*(recon_x-x)]");
2444
2445 // calculate [f + F*(recon_x-x)]T*fstar
2446 ✗ double *conv_tmpmatrixrhs = (double*) calloc(transpose_add_f_F_recon_x_x.rows * conv_fstar.column, sizeof(double));
2447 ✗ solveMatrixMultiplication(transpose_add_f_F_recon_x_x.data, conv_fstar.data, transpose_add_f_F_recon_x_x.rows, transpose_add_f_F_recon_x_x.column, conv_fstar.rows, conv_fstar.column, conv_tmpmatrixrhs, logfile, data);
2448 //printMatrix(conv_tmpmatrixrhs, transpose_add_f_F_recon_x_x.rows, conv_fstar.column,"[f + F*(recon_x-x)]*fstar");
2449
2450 // scale the matrix with 2*[f + F*(recon_x-x)]T*fstar
2451 ✗ scaleVector(transpose_add_f_F_recon_x_x.rows, conv_fstar.column, 2.0, conv_tmpmatrixrhs);
2452 //printMatrix(conv_tmpmatrixrhs, transpose_add_f_F_recon_x_x.rows, conv_fstar.column,"2*[f + F*(recon_x-x)]*fstar");
2453 ✗ matrixData struct_conv_tmpmatrixrhs = {transpose_add_f_F_recon_x_x.rows, conv_fstar.column, conv_tmpmatrixrhs};
2454
2455 /*
2456 * solve the final J*=J*=(recon_x-x)T*(Sx^-1)*(recon_x-x)+2.[f+F*(recon_x-x)]T*fstar
2457 * J*=_struct_conv_tmpmatrixlhs + struct_conv_tmpmatrixrhs
2458 */
2459 ✗ matrixData struct_Jstar = solveMatrixAddition(struct_conv_tmpmatrixlhs, struct_conv_tmpmatrixrhs, logfile, data);
2460 //printMatrix(struct_Jstar.data,struct_Jstar.rows,struct_Jstar.column,"J*",logfile);
2461
2462 ✗ int r = data->modelData->nSetcVars; // number of setc equations
2463 ✗ double val = 1.0 / r;
2464
2465 /*
2466 * calculate J/r < epselon
2467 */
2468 ✗ scaleVector(struct_Jstar.rows, struct_Jstar.column, val, struct_Jstar.data);
2469 //printMatrix(struct_Jstar.data,struct_Jstar.rows,struct_Jstar.column,"J*/r ");
2470 ✗ double convergedvalue = struct_Jstar.data[0];
2471
2472 // free the tmp matrices
2473 ✗ free(conv_data1);
2474 ✗ free(conv_tmpmatrixlhs);
2475 ✗ free(tmp_F_recon_x_x);
2476 ✗ free(conv_tmpmatrixrhs);
2477 ✗ free(add_f_F_recon_x_x.data);
2478 ✗ free(struct_Jstar.data);
2479 ✗ free(conv_data1Transpose.data);
2480 ✗ free(transpose_add_f_F_recon_x_x.data);
2481 //free(conv_recon_x.data);
2482 ✗ free(copy_reconx_x.data);
2483
2484 ✗ return convergedvalue;
2485 }
2486
2487 /*
2488 * function which calculates qualityValue J = transpose (x_reconciled – x_measured)*Sx^-1*(x_reconciled – x_measured)
2489 */
2490 ✗ double calculateQualityValue(matrixData reconciledX, matrixData Sx, csvData measuredX, ofstream & logfile, DATA* data)
2491 {
2492 ✗ logfile << "Calculations of Quality Value (J) " << "\n";
2493 ✗ logfile << "=================================\n";
2494
2495 ✗ printMatrix(reconciledX.data, reconciledX.rows, reconciledX.column, "reconciled_x", logfile);
2496 ✗ inputData measured_x = getInputData(measuredX, logfile);
2497 ✗ printMatrix(measured_x.data, measured_x.rows, measured_x.column, "measured_X", logfile);
2498 ✗ printMatrix(Sx.data, Sx.rows, Sx.column, "Sx", logfile);
2499 ✗ matrixData measured_xCopy = {measured_x.rows, measured_x.column, measured_x.data};
2500 ✗ double *newX = (double*) calloc (measured_x.rows * 1, sizeof(double));
2501 ✗ solveMatrixSubtraction(reconciledX, measured_xCopy, newX, logfile, data);
2502 ✗ printMatrix(newX, measured_x.rows, measured_x.column, "x_reconciled - measured_X", logfile);
2503
2504 /* solves (Sx^-1)*(x_reconciled – x_measured)
2505 * Solve the inverse of matrix Sx using linear form
2506 * Ax=b
2507 * where A=Sx and b= (recon_x-x_measured) to avoid inversion of Sx which is
2508 * expensive
2509 */
2510 //solveSystemFstar(conv_sx.rows, 1, conv_sx.data, conv_data1result.data, logfile, data);
2511 ✗ matrixData sub = {measured_x.rows, measured_x.column, newX};
2512 ✗ matrixData subCopy = copyMatrix(sub);
2513 ✗ solveSystemFstar(Sx.rows, 1, Sx.data, newX, logfile, data);
2514 ✗ printMatrix(newX, measured_x.rows, measured_x.column, "Sx-inverse", logfile);
2515
2516 ✗ matrixData subCopyTranspose = getTransposeMatrix(subCopy);
2517
2518 ✗ double *J = (double*) calloc(subCopyTranspose.rows * measured_x.column, sizeof(double));
2519 /*
2520 * Solve transpose (x_reconciled – x_measured)*Sx^-1*(x_reconciled – x_measured)
2521 */
2522 ✗ solveMatrixMultiplication(subCopyTranspose.data, newX, subCopyTranspose.rows, subCopyTranspose.column, measured_x.rows, measured_x.column, J, logfile, data);
2523 ✗ printMatrix(J, subCopyTranspose.rows, measured_x.column, "J", logfile);
2524
2525 ✗ double J_value = J[0];
2526
2527 // --- Free all temporary allocations ---
2528 ✗ free(newX);
2529 ✗ free(subCopy.data);
2530 ✗ free(subCopyTranspose.data);
2531 ✗ free(measured_x.data);
2532 ✗ free(J);
2533
2534 ✗ return J_value;
2535 }
2536
2537 /*
2538 * Example Function which performs matrix inverse
2539 * using dgetri_ and dgetrf_ LAPACK routine
2540 * which is expensive one and not recommended
2541 * use it when no other way to compute it
2542 */
2543 ✗ void checkExpensiveMatrixInverse()
2544 {
2545 ✗ double newval[3*3]={3,2,0,
2546 0,0,1,
2547 2,-2,1};
2548
2549 ✗ int N = 3;
2550 int LDA = N;
2551 int LDB = N;
2552 int ipiv[3];
2553 ✗ int info = 1;
2554 ✗ int LWORK = N;
2555 ✗ double *WORK = (double*) calloc(LWORK, sizeof(double));
2556
2557 ✗ dgetrf_(&N, &N, newval, &N, ipiv, &info);
2558 ✗ dgetri_(&N, newval, &N, ipiv, WORK, &LWORK, &info);
2559 //printMatrix(newval,3,3,"Expensive_Matrix_Inverse");
2560 ✗ }
2561
2562 /*
2563 * Function which performs matrix inverse without performing
2564 * actual matrix inverse, Instead use the dgesv to get result
2565 * Ax=b where matrix mutiplication of x=bA gives the inversed
2566 * mutiplication result b with A inverse
2567 */
2568 ✗ void checkInExpensiveMatrixInverse(ofstream & logfile, DATA * data)
2569 {
2570 ✗ double newchecksx[3*3]={1,1,1,
2571 0,0.95,0,
2572 0,0,0.95};
2573 ✗ double checksx[3 * 1] = {-0.028, 0.026, -0.004};
2574 ✗ solveSystemFstar(3, 1, newchecksx, checksx, logfile, data);
2575 //printMatrix(checksx,3,1,"InExpensive_Matrix_Inverse");
2576 ✗ }
2577
2578 ✗ int RunReconciliation(DATA *data, threadData_t *threadData, inputData x, matrixData Sx, matrixData tmpjacF, matrixData tmpjacFt, double eps, int iterationcount, csvData csvinputs, matrixData xdiag, matrixData sxdiag, ofstream &logfile, correlationDataWarning & warningCorrelationData, dataReconciliationData& datareconciliationdata)
2579 {
2580 // set the inputs from csv file to simulationInfo datainputVars
2581 ✗ for (int i = 0; i < x.rows * x.column; i++)
2582 {
2583 ✗ data->simulationInfo->datainputVars[i] = x.data[i];
2584 //logfile << "input data:" << x.data[i]<<"\n";
2585 }
2586
2587 /* set the inputs via this special function generated for dataReconciliation
2588 * which also sets inputs for models not involving top level inputs
2589 */
2590 ✗ data->callback->data_function(data, threadData);
2591 //data->callback->input_function(data, threadData);
2592 ✗ data->callback->functionDAE(data, threadData);
2593 //data->callback->functionODE(data,threadData);
2594 ✗ data->callback->setc_function(data, threadData);
2595
2596 //data->callback->setb_function(data, threadData);
2597
2598 ✗ matrixData jacF = getJacobianMatrixF(data, threadData, logfile);
2599 ✗ matrixData jacFt = getTransposeMatrix(jacF);
2600
2601 ✗ printMatrix(jacF.data, jacF.rows, jacF.column, "F", logfile);
2602 ✗ printMatrix(jacFt.data, jacFt.rows, jacFt.column, "Ft", logfile);
2603
2604 // allocate data for setc array
2605 ✗ double *setc = (double*) calloc (data->modelData->nSetcVars, sizeof(double));
2606
2607 // store the setc data to compute for convergence as setc will be overriddeen with new values
2608 ✗ double *tmpsetc = (double*) calloc (data->modelData->nSetcVars, sizeof(double));
2609
2610 /* loop to store the data C(x,y) rhs side, get the elements in reverse order */
2611 int t = 0;
2612 ✗ for (int i = data->modelData->nSetcVars; i > 0; i--)
2613 {
2614 ✗ setc[t] = data->simulationInfo->setcVars[i-1];
2615 ✗ tmpsetc[t] = data->simulationInfo->setcVars[i-1];
2616 ✗ t++;
2617 //cout << "array_setc_vars:=>" << t << ":" << data->simulationInfo->setcVars[i-1] << "\n";
2618 }
2619
2620 int nsetcvars = data->modelData->nSetcVars;
2621 ✗ matrixData vector_c = {nsetcvars, 1, tmpsetc};
2622
2623 //allocate data for matrix multiplication F*Sx
2624 ✗ double *tmpmatrixC = (double*) calloc (jacF.rows * Sx.column, sizeof(double));
2625 ✗ solveMatrixMultiplication (jacF.data, Sx.data, jacF.rows, jacF.column,Sx.rows, Sx.column, tmpmatrixC, logfile, data);
2626 //printMatrix(tmpmatrixC,jacF.rows,Sx.column,"F*Sx");
2627
2628 //allocate data for matrix multiplication (F*Sx)*Ftranspose
2629 ✗ double *tmpmatrixD = (double*) calloc (jacF.rows * jacFt.column, sizeof(double));
2630 ✗ solveMatrixMultiplication (tmpmatrixC, jacFt.data, jacF.rows, Sx.column, jacFt.rows, jacFt.column, tmpmatrixD, logfile, data);
2631
2632 //printMatrix(tmpmatrixD,jacF.rows,jacFt.column,"F*Sx*Ft");
2633 //printMatrix(setc,nsetcvars,1,"c(x,y)");
2634
2635 /*
2636 * Copy tmpmatrixC and tmpmatrixD to avoid loss of data
2637 * when calculating F*
2638 */
2639 ✗ matrixData cpytmpmatrixC = {jacF.rows, Sx.column, tmpmatrixC};
2640 ✗ matrixData cpytmpmatrixD = {jacF.rows, jacFt.column, tmpmatrixD};
2641 ✗ matrixData tmpmatrixC1 = copyMatrix(cpytmpmatrixC);
2642 ✗ matrixData tmpmatrixD1 = copyMatrix(cpytmpmatrixD);
2643
2644 ✗ if(OMC_ACTIVE_STREAM(OMC_LOG_JAC))
2645 {
2646 ✗ logfile << "Calculations of Matrix (F*Sx*Ft) f* = c(x,y) " << "\n";
2647 ✗ logfile << "============================================\n";
2648 ✗ printMatrix(tmpmatrixC, jacF.rows, Sx.column, "F*Sx", logfile);
2649 ✗ printMatrix(tmpmatrixD, jacF.rows, jacFt.column, "F*Sx*Ft", logfile);
2650 ✗ printMatrix(setc, nsetcvars, 1, "c(x,y)", logfile);
2651 }
2652
2653 /*
2654 * calculate f* for covariance matrix (F*Sx*Ftranspose).F*= c(x,y)
2655 * matrix setc will be overridden with new values which is the output
2656 * for the calculation A *x =B
2657 * A = tmpmatrixD
2658 * B = setc
2659 */
2660 ✗ solveSystemFstar(jacF.rows, 1, tmpmatrixD, setc, logfile, data);
2661
2662 ✗ if(OMC_ACTIVE_STREAM(OMC_LOG_JAC))
2663 {
2664 ✗ printMatrix(setc, jacF.rows, 1, "f*", logfile);
2665 ✗ logfile << "***** Completed ****** \n\n";
2666 }
2667
2668 ✗ matrixData tmpxcap = {x.rows, 1, x.data};
2669 ✗ matrixData tmpfstar = {jacF.rows, 1, setc};
2670 ✗ matrixData reconciled_X = solveReconciledX(tmpxcap, Sx, jacFt, tmpfstar, logfile, data);
2671 //printMatrix(reconciled_X.data,reconciled_X.rows,reconciled_X.column,"reconciled_X ===> (x - (Sx*Ft*fstar))");
2672
2673 ✗ if (OMC_ACTIVE_STREAM(OMC_LOG_JAC))
2674 {
2675 ✗ logfile << "Calculations of Matrix (F*Sx*Ft) F* = F*Sx " << "\n";
2676 ✗ logfile << "===============================================\n";
2677 //printMatrix(tmpmatrixC,jacF.rows,Sx.column,"F*Sx");
2678 //printMatrix(tmpmatrixD,jacF.rows,jacFt.column,"F*Sx*Ft");
2679 ✗ printMatrix(tmpmatrixC1.data, tmpmatrixC1.rows, tmpmatrixC1.column, "F*Sx", logfile);
2680 ✗ printMatrix(tmpmatrixD1.data, tmpmatrixD1.rows, tmpmatrixD1.column, "F*Sx*Ft", logfile);
2681 }
2682
2683 /*
2684 * calculate F* for covariance matrix (F*Sx*Ftranspose).F*= (F*Sx)
2685 * tmpmatrixC1 will be overridden with new values which is the output
2686 * for the calculation A *x =B
2687 * A = tmpmatrixD
2688 * B = tmpmatrixC
2689 */
2690 ✗ solveSystemFstar(jacF.rows, Sx.column, tmpmatrixD1.data, tmpmatrixC1.data, logfile, data);
2691
2692 ✗ if (OMC_ACTIVE_STREAM(OMC_LOG_JAC))
2693 {
2694 ✗ printMatrix(tmpmatrixC1.data, jacF.rows, Sx.column, "F*", logfile);
2695 ✗ logfile << "***** Completed ****** \n\n";
2696 }
2697
2698 ✗ matrixData tmpFstar = {jacF.rows, Sx.column, tmpmatrixC1.data};
2699 ✗ matrixData reconciled_Sx = solveReconciledSx(Sx, jacFt, tmpFstar, logfile, data);
2700 //printMatrix(reconciled_Sx.data,reconciled_Sx.rows,reconciled_Sx.column,"reconciled Sx ===> (Sx - (Sx*Ft*Fstar))");
2701
2702 ✗ matrixData copySx = copyMatrix(Sx);
2703 ✗ double value = solveConvergence(data, reconciled_X, reconciled_Sx, x, copySx, jacF, vector_c, tmpfstar, logfile);
2704 ✗ if (value > eps)
2705 {
2706 ✗ logfile << "J*/r" << "(" << value << ")" << " > " << eps << ", Value not Converged \n";
2707 ✗ logfile << "==========================================\n\n";
2708 }
2709 ✗ if (value > eps)
2710 {
2711 ✗ logfile << "Running Convergence iteration: " << iterationcount << " with the following reconciled values:" << "\n";
2712 ✗ logfile << "========================================================================" << "\n";
2713 //cout << "J*/r :=" << value << "\n";
2714 //printMatrix(jacF.data,jacF.rows,jacF.column,"F");
2715 //printMatrix(jacFt.data,jacFt.rows,jacFt.column,"Ft");
2716 //printMatrix(setc,jacF.rows,1,"f*");
2717 //printMatrix(tmpmatrixC,jacF.rows,Sx.column,"F*");
2718 ✗ printMatrixWithHeaders(reconciled_X.data, reconciled_X.rows, reconciled_X.column, csvinputs.headers, "reconciled_X ===> (x - (Sx*Ft*fstar))", logfile);
2719 ✗ printMatrixWithHeaders(reconciled_Sx.data, reconciled_Sx.rows, reconciled_Sx.column, csvinputs.headers, "reconciled_Sx ===> (Sx - (Sx*Ft*Fstar))", logfile);
2720 //free(x.data);
2721 //free(Sx.data);
2722 //x.data = reconciled_X.data;
2723 ✗ for (int i = 0; i < reconciled_X.rows * reconciled_X.column; i++)
2724 {
2725 ✗ x.data[i] = reconciled_X.data[i];
2726 }
2727
2728 ✗ free(reconciled_X.data);
2729 ✗ free(jacF.data);
2730 ✗ free(jacFt.data);
2731 ✗ free(tmpmatrixC);
2732 ✗ free(tmpmatrixD);
2733 ✗ free(tmpmatrixC1.data);
2734 ✗ free(tmpmatrixD1.data);
2735 ✗ free(setc);
2736 ✗ free(tmpsetc);
2737 ✗ free(copySx.data);
2738 ✗ free(reconciled_Sx.data);
2739 ✗ iterationcount++;
2740 ✗ return RunReconciliation(data, threadData, x, Sx, jacF, jacFt, eps, iterationcount, csvinputs, xdiag, sxdiag, logfile, warningCorrelationData, datareconciliationdata);
2741 }
2742
2743 ✗ if (value < eps && iterationcount == 1)
2744 {
2745 ✗ logfile << "J*/r" << "(" << value << ")" << " > " << eps << ", Convergence iteration not required \n\n";
2746 }
2747 else
2748 {
2749 ✗ logfile << "***** Value Converged, Convergence Completed******* \n\n";
2750 }
2751
2752 ✗ double J = calculateQualityValue(reconciled_X, Sx, csvinputs, logfile, data);
2753
2754 ✗ logfile << "Final Results:\n";
2755 ✗ logfile << "=============\n";
2756 ✗ logfile << "Total Iteration to Converge : " << iterationcount << "\n";
2757 ✗ logfile << "Final Converged Value(J*/r) : " << value << "\n";
2758 ✗ logfile << "Final value of the objective function (J) : " << J << "\n";
2759 ✗ logfile << "Epsilon : " << eps << "\n";
2760 ✗ printMatrixWithHeaders(reconciled_X.data, reconciled_X.rows, reconciled_X.column, csvinputs.headers, "reconciled_X ===> (x - (Sx*Ft*fstar))", logfile);
2761 ✗ printMatrixWithHeaders(reconciled_Sx.data, reconciled_Sx.rows, reconciled_Sx.column, csvinputs.headers, "reconciled_Sx ===> (Sx - (Sx*Ft*Fstar))", logfile);
2762
2763 ✗ dumpReconciledSxToCSV(reconciled_Sx.data, reconciled_Sx.rows, reconciled_Sx.column, csvinputs.headers, data);
2764
2765 // copy the reconciledSx matrix for state Estimation
2766 ✗ matrixData copyReconciledSx = copyMatrix(reconciled_Sx);
2767
2768 /*
2769 * Calculate half width Confidence interval
2770 * W=lambda*sqrt(Sx)
2771 * where lamba = 1.96 and
2772 * Sx - diagonal elements of reconciled_Sx
2773 */
2774 ✗ double *reconSx_diag = (double*) calloc (reconciled_Sx.rows * 1, sizeof(double));
2775 ✗ getDiagonalElements(reconciled_Sx.data, reconciled_Sx.rows, reconciled_Sx.column, reconSx_diag);
2776
2777 ✗ matrixData copyreconSx_diag = {reconciled_Sx.rows, 1, reconSx_diag};
2778 ✗ matrixData tmpcopyreconSx_diag = copyMatrix(copyreconSx_diag);
2779
2780 ✗ if (OMC_ACTIVE_STREAM(OMC_LOG_JAC))
2781 {
2782 ✗ logfile << "Calculations of HalfWidth Confidence Interval " << "\n";
2783 ✗ logfile << "===============================================\n";
2784 ✗ printMatrix(copyreconSx_diag.data, reconciled_Sx.rows, 1, "reconciled-Sx_Diagonal", logfile);
2785 }
2786
2787 ✗ calculateSquareRoot(copyreconSx_diag.data, reconciled_Sx.rows);
2788
2789 ✗ if (OMC_ACTIVE_STREAM(OMC_LOG_JAC))
2790 {
2791 ✗ printMatrix(copyreconSx_diag.data, reconciled_Sx.rows, 1, "reconciled-Sx_SquareRoot", logfile);
2792 ✗ logfile << "*****Completed***********\n";
2793 }
2794
2795 ✗ scaleVector(reconciled_Sx.rows, 1, 1.96, copyreconSx_diag.data);
2796 ✗ printMatrixWithHeaders(copyreconSx_diag.data, reconciled_Sx.rows, 1, csvinputs.headers, "Wx-HalfWidth-Interval-(1.96)*sqrt(Sx_diagonal)", logfile);
2797
2798 /*
2799 * Calculate individual tests
2800 * (recon_x - x)/sqrt(Sx-recon_Sx)
2801 */
2802 ✗ double *newSx_diag = (double*) calloc (reconciled_Sx.rows * 1, sizeof(double));
2803 ✗ solveMatrixSubtraction(sxdiag, tmpcopyreconSx_diag, newSx_diag, logfile, data);
2804
2805 ✗ if (OMC_ACTIVE_STREAM(OMC_LOG_JAC))
2806 {
2807 ✗ logfile << "Calculations of Individual Tests " << "\n";
2808 ✗ logfile << "===============================================\n";
2809 ✗ printMatrix(newSx_diag, sxdiag.rows, sxdiag.column, "Sx-recon_Sx", logfile);
2810 }
2811
2812 ✗ calculateSquareRoot(newSx_diag, reconciled_Sx.rows);
2813
2814 ✗ if (OMC_ACTIVE_STREAM(OMC_LOG_JAC))
2815 {
2816 ✗ printMatrix (newSx_diag, sxdiag.rows, sxdiag.column, "squareroot-newSx", logfile);
2817 }
2818
2819 ✗ double *newX = (double*) calloc (xdiag.rows * 1, sizeof(double));
2820 ✗ solveMatrixSubtraction(reconciled_X, xdiag, newX, logfile, data);
2821
2822 // calculate absolute value for this numeric analysis
2823 ✗ for (int a = 0; a < xdiag.rows; a++)
2824 {
2825 ✗ newX[a] = fabs (newX[a]);
2826 }
2827
2828 ✗ if (OMC_ACTIVE_STREAM(OMC_LOG_JAC))
2829 {
2830 ✗ printMatrix(newX, xdiag.rows, xdiag.column, "recon_X - X", logfile);
2831 ✗ logfile << "*********Completed***********\n";
2832 }
2833
2834 ✗ for (int val = 0; val < xdiag.rows; val++)
2835 {
2836 ✗ newX[val] = newX[val] / max(newSx_diag[val], sqrt(sxdiag.data[val] / 10));
2837 }
2838
2839 ✗ printMatrixWithHeaders(newX, xdiag.rows, xdiag.column, csvinputs.headers, "IndividualTests_Value- (recon_x-x)/sqrt(Sx_diag)", logfile);
2840
2841 // copy the outputs for state Estimation
2842 ✗ if (omc_flag[FLAG_DATA_RECONCILE_STATE])
2843 ✗ datareconciliationdata = {csvinputs, xdiag, reconciled_X, copyReconciledSx, copyreconSx_diag, newX, eps, iterationcount, value, J, warningCorrelationData};
2844
2845 // read and update reconciled mo file with new values
2846
2847 ✗ updateReconciledMo(data, threadData, csvinputs.headers, reconciled_X.data, logfile);
2848
2849 boundaryConditionData boundaryconditiondata;
2850 // create HTML Report for D.1
2851 ✗ if (omc_flag[FLAG_DATA_RECONCILE])
2852 {
2853 ✗ createHtmlReportFordataReconciliation(data, csvinputs, xdiag, reconciled_X, copyreconSx_diag, newX, eps, iterationcount, value, J, warningCorrelationData, boundaryconditiondata);
2854 // free the memory for data Reconciliation
2855 ✗ free(jacF.data);
2856 ✗ free(jacFt.data);
2857 ✗ free(tmpmatrixC);
2858 ✗ free(tmpmatrixD);
2859 ✗ free(tmpmatrixC1.data);
2860 ✗ free(tmpmatrixD1.data);
2861 ✗ free(setc);
2862 ✗ free(tmpsetc);
2863 ✗ free(copySx.data);
2864 ✗ free(copyReconciledSx.data);
2865 ✗ free(reconciled_X.data);
2866 ✗ free(reconciled_Sx.data);
2867 ✗ free(newX);
2868 ✗ free(newSx_diag);
2869 ✗ free(copyreconSx_diag.data);
2870 ✗ free(tmpcopyreconSx_diag.data);
2871 }
2872 else
2873 {
2874 // state estimation, do not free the commented data as it is used in state estimation report
2875 ✗ free(jacF.data);
2876 ✗ free(jacFt.data);
2877 ✗ free(tmpmatrixC);
2878 ✗ free(tmpmatrixD);
2879 ✗ free(tmpmatrixC1.data);
2880 ✗ free(tmpmatrixD1.data);
2881 ✗ free(setc);
2882 ✗ free(tmpsetc);
2883 ✗ free(copySx.data);
2884 //free(copyReconciledSx.data);
2885 //free(reconciled_X.data);
2886 ✗ free(reconciled_Sx.data);
2887 //free(newX);
2888 ✗ free(newSx_diag);
2889 //free(copyreconSx_diag.data);
2890 ✗ free(tmpcopyreconSx_diag.data);
2891 }
2892
2893 return 0;
2894 }
2895
2896 /*
2897 * Runs the numerical procedure to compute Boundary conditions (D.2)
2898 */
2899 ✗ int reconcileBoundaryConditions(DATA * data, threadData_t * threadData, inputData reconciled_x, matrixData reconciled_Sx, boundaryConditionData& boundaryconditiondata, correlationDataWarning& warningCorrelationData, ofstream& logfile)
2900 {
2901 // set the inputs from csv file to simulationInfo datainputVars
2902 ✗ for (int i = 0; i < reconciled_x.rows * reconciled_x.column; i++)
2903 {
2904 ✗ data->simulationInfo->datainputVars[i] = reconciled_x.data[i];
2905 }
2906
2907 /* set the inputs via this special function generated for dataReconciliation
2908 * which also sets inputs for models not involving top level inputs
2909 */
2910 ✗ data->callback->data_function(data, threadData);
2911 // solve the system with reconciled input values got from D.1
2912 ✗ data->callback->functionDAE(data, threadData);
2913 // call the setc function which stores the results of boundary conditions variable
2914 ✗ if (omc_flag[FLAG_DATA_RECONCILE_BOUNDARY])
2915 ✗ data->callback->setc_function(data, threadData);
2916 // call the setb function for state estimation
2917 ✗ if (omc_flag[FLAG_DATA_RECONCILE_STATE])
2918 ✗ data->callback->setb_function(data, threadData);
2919
2920 // Compute the Jacobian Matrix F or H
2921 matrixData jacF;
2922 ✗ if (omc_flag[FLAG_DATA_RECONCILE_BOUNDARY])
2923 {
2924 ✗ jacF = getJacobianMatrixF(data, threadData, logfile, true);
2925 }
2926 else
2927 {
2928 ✗ jacF = getJacobianMatrixH(data, threadData, logfile, true);
2929 }
2930
2931 ✗ printMatrix(jacF.data, jacF.rows, jacF.column, "F", logfile);
2932
2933 // Compute the Transpose of jacobian Matrix F
2934 ✗ matrixData jacFt = getTransposeMatrix(jacF);
2935 ✗ printMatrix(jacFt.data, jacFt.rows, jacFt.column, "Ft", logfile);
2936
2937 /*
2938 * Compute St = jacF*reconciles_Sx*jacFt
2939 */
2940 // F*reconciledSx
2941 ✗ double *tmpMatrixAf = (double *)calloc(jacF.rows * reconciled_Sx.column, sizeof(double));
2942 ✗ solveMatrixMultiplication(jacF.data, reconciled_Sx.data, jacF.rows, jacF.column, reconciled_Sx.rows, reconciled_Sx.column, tmpMatrixAf, logfile, data);
2943 ✗ printMatrix(tmpMatrixAf, jacF.rows, reconciled_Sx.column, "F*reconciled_Sx", logfile);
2944
2945 //(F*reconciledSx)*ftranspose
2946 ✗ double *S_t = (double*) calloc(jacF.rows * jacFt.column, sizeof(double));
2947 ✗ solveMatrixMultiplication(tmpMatrixAf, jacFt.data, jacF.rows, jacF.column, jacFt.rows, jacFt.column, S_t, logfile, data);
2948 ✗ printMatrix(S_t, jacF.rows, jacFt.column, "(s_t = F*reconciled_Sx*Ft)", logfile);
2949
2950 /*
2951 * Calculate half width Confidence interval
2952 * W=lambda*sqrt(S_t)
2953 * where lamba = 1.96 and
2954 * S_t - diagonal elements of S_t
2955 */
2956 ✗ double *reconSt_diag = (double*) calloc (jacF.rows * 1, sizeof(double));
2957 ✗ getDiagonalElements(S_t, jacF.rows, jacFt.column, reconSt_diag);
2958
2959 ✗ if (OMC_ACTIVE_STREAM(OMC_LOG_JAC))
2960 {
2961 ✗ logfile << "Calculations of half-width confidence interval" << "\n";
2962 ✗ logfile << "===============================================\n";
2963 ✗ printMatrix(reconSt_diag, jacF.rows, 1, "S_t_Diagonal", logfile);
2964 }
2965
2966 ✗ calculateSquareRoot(reconSt_diag, jacF.rows);
2967
2968 ✗ if (OMC_ACTIVE_STREAM(OMC_LOG_JAC))
2969 {
2970 ✗ printMatrix(reconSt_diag, jacF.rows, 1, "S_t_SquareRoot", logfile);
2971 }
2972
2973 ✗ scaleVector(jacF.rows, 1, 1.96, reconSt_diag);
2974
2975 // check for BoundaryConditionVars.txt file exists to generate the html report
2976 std::string boundaryConditionsVarsFilename;
2977 ✗ if (omc_flag[FLAG_OUTPUT_PATH])
2978 {
2979 ✗ boundaryConditionsVarsFilename = std::string(omc_flagValue[FLAG_OUTPUT_PATH]) + "/" + std::string(data->modelData->modelFilePrefix) + "_BoundaryConditionVars.txt";
2980 ✗ copyReferenceFile(data, "_BoundaryConditionVars.txt");
2981 }
2982 else
2983 {
2984 ✗ boundaryConditionsVarsFilename = std::string(data->modelData->modelFilePrefix) + "_BoundaryConditionVars.txt";
2985 }
2986
2987 vector<std::string> boundaryConditionVars;
2988
2989 ✗ ifstream boundaryConditionVarsip(boundaryConditionsVarsFilename);
2990 std::string line;
2991 ✗ if (boundaryConditionVarsip.good())
2992 {
2993 ✗ while (boundaryConditionVarsip.good())
2994 {
2995 ✗ getline(boundaryConditionVarsip, line);
2996 ✗ if (!line.empty())
2997 {
2998 //std::cout << "\n reading nonVariables of interest : " << line;
2999 ✗ boundaryConditionVars.push_back(line);
3000 }
3001 }
3002 ✗ boundaryConditionVarsip.close();
3003 //omc_unlink(boundaryConditionsVarsFilename.c_str());
3004 }
3005 else
3006 {
3007 ✗ errorStreamPrint(OMC_LOG_STDOUT, 0, "Boundary conditions vars filename not found: %s.", boundaryConditionsVarsFilename.c_str());
3008 ✗ logfile << "| error | " << "Boundary conditions vars filename not found: " << boundaryConditionsVarsFilename << "\n";
3009 ✗ logfile.close();
3010 ✗ createErrorHtmlReportForBoundaryConditions(data);
3011 ✗ exit(1);
3012 }
3013
3014 ✗ printMatrixWithHeaders(reconSt_diag, jacF.rows, 1, boundaryConditionVars, "Half-width Confidence Interval(1.96*S_t_SquareRoot)", logfile);
3015
3016 // allocate data for boundaryconditions vars simulation results
3017 double *boundaryConditionVarsResults;
3018 ✗ if (omc_flag[FLAG_DATA_RECONCILE_BOUNDARY])
3019 {
3020 ✗ boundaryConditionVarsResults = (double *)calloc(data->modelData->nSetcVars, sizeof(double));
3021 int t = 0;
3022 ✗ for (int i = data->modelData->nSetcVars; i > 0; i--)
3023 {
3024 ✗ boundaryConditionVarsResults[t] = data->simulationInfo->setcVars[i - 1];
3025 ✗ t++;
3026 }
3027 }
3028 else
3029 {
3030 ✗ boundaryConditionVarsResults = (double *)calloc(data->modelData->nSetbVars, sizeof(double));
3031 int t = 0;
3032 ✗ for (int i = data->modelData->nSetbVars; i > 0; i--)
3033 {
3034 ✗ boundaryConditionVarsResults[t] = data->simulationInfo->setbVars[i - 1];
3035 ✗ t++;
3036 }
3037 }
3038
3039 ✗ printBoundaryConditionsResults(boundaryConditionVarsResults, reconSt_diag, jacF.rows, 1, boundaryConditionVars, "Final Results", logfile);
3040
3041 // create html report for boundary conditions
3042 ✗ if (omc_flag[FLAG_DATA_RECONCILE_BOUNDARY])
3043 ✗ createHtmlReportForBoundaryConditions(data, boundaryConditionVars, boundaryConditionVarsResults, reconSt_diag, warningCorrelationData);
3044
3045 // copy the results for state estimation
3046 ✗ if (omc_flag[FLAG_DATA_RECONCILE_STATE])
3047 ✗ boundaryconditiondata = {boundaryConditionVars, boundaryConditionVarsResults, reconSt_diag};
3048
3049 // free the memory
3050 ✗ if (omc_flag[FLAG_DATA_RECONCILE_BOUNDARY])
3051 {
3052 ✗ free(reconciled_Sx.data);
3053 ✗ free(reconciled_x.data);
3054 ✗ free(tmpMatrixAf);
3055 ✗ free(S_t);
3056 ✗ free(jacF.data);
3057 ✗ free(jacFt.data);
3058 ✗ free(reconSt_diag);
3059 ✗ free(boundaryConditionVarsResults);
3060 }
3061 else
3062 {
3063 // state Estimation
3064 ✗ free(tmpMatrixAf);
3065 ✗ free(S_t);
3066 ✗ free(jacF.data);
3067 ✗ free(jacFt.data);
3068 }
3069
3070 ✗ return 0;
3071 ✗ }
3072
3073 /*
3074 * State Estimation numerical computation is a combination of
3075 * Data Reconciliation and boundary condition computation
3076 */
3077
3078 ✗ int stateEstimation(DATA *data, threadData_t *threadData, inputData x, matrixData Sx, matrixData tmpjacF, matrixData tmpjacFt, double eps, int iterationcount, csvData csvinputs, matrixData xdiag, matrixData sxdiag, ofstream &logfile, correlationDataWarning & warningCorrelationData)
3079 {
3080 // run the data Reconciliation
3081 dataReconciliationData datareconciliationdata;
3082 ✗ RunReconciliation(data, threadData, x, Sx, tmpjacF, tmpjacFt, eps, 1, csvinputs, xdiag, sxdiag, logfile, warningCorrelationData, datareconciliationdata);
3083
3084 //printMatrixWithHeaders(datareconciliationdata.reconciled_X.data, datareconciliationdata.reconciled_X.rows, datareconciliationdata.reconciled_X.column, csvinputs.headers, "ARRRRRreconciled_X ===> (x - (Sx*Ft*fstar))", logfile);
3085 //printMatrixWithHeaders(datareconciliationdata.copyreconSx_diag.data, datareconciliationdata.copyreconSx_diag.rows, datareconciliationdata.copyreconSx_diag.column, csvinputs.headers, "ARRRRRRreconciled_Sx ===> (Sx - (Sx*Ft*Fstar))", logfile);
3086 //printMatrixWithHeaders(datareconciliationdata.reconciled_SX.data, datareconciliationdata.reconciled_SX.rows, datareconciliationdata.reconciled_SX.column, csvinputs.headers, "NovakreconciledS_X ===> (x - (Sx*Ft*fstar))", logfile);
3087
3088 // Compute Boundary conditions only if unmeasured variables exist
3089 boundaryConditionData boundaryconditiondata;
3090 ✗ if (data->modelData->nSetbVars > 0)
3091 {
3092 // copy the reference files "BoundaryConditionsEquations.html and _BoundaryConditionIntermediateEquations.html to output path"
3093 ✗ if (omc_flag[FLAG_OUTPUT_PATH])
3094 {
3095 ✗ copyReferenceFile(data, "_BoundaryConditionIntermediateEquations.html");
3096 }
3097 // pepare data to compute boundary condition
3098 ✗ inputData reconciled_x = {datareconciliationdata.reconciled_X.rows, datareconciliationdata.reconciled_X.column, datareconciliationdata.reconciled_X.data, {}};
3099 ✗ matrixData reconciled_Sx = {datareconciliationdata.reconciled_SX.rows, datareconciliationdata.reconciled_SX.column, datareconciliationdata.reconciled_SX.data};
3100
3101 ✗ logfile << "\n\nCalculation of Boundary condition \n" << "====================================\n";
3102 ✗ reconcileBoundaryConditions(data, threadData, reconciled_x, reconciled_Sx, boundaryconditiondata, warningCorrelationData, logfile);
3103 // printBoundaryConditionsResults(boundaryconditiondata.boundaryConditionVarsResults, boundaryconditiondata.reconSt_diag, boundaryconditiondata.boundaryConditionVars.size(), 1, boundaryconditiondata.boundaryConditionVars, "Final Results Copied", logfile);
3104 }
3105
3106 ✗ createHtmlReportFordataReconciliation(data, datareconciliationdata.csvinputs, datareconciliationdata.xdiag, datareconciliationdata.reconciled_X, datareconciliationdata.copyreconSx_diag, datareconciliationdata.newX, eps, datareconciliationdata.iterationcount, datareconciliationdata.value, datareconciliationdata.J, warningCorrelationData, boundaryconditiondata);
3107
3108 // free data Reconciliation data
3109 ✗ free(datareconciliationdata.reconciled_SX.data);
3110 ✗ free(datareconciliationdata.reconciled_X.data);
3111 ✗ free(datareconciliationdata.copyreconSx_diag.data);
3112 ✗ free(datareconciliationdata.newX);
3113
3114 // free boundaryCondition data
3115 ✗ if (data->modelData->nSetbVars > 0)
3116 {
3117 ✗ free(boundaryconditiondata.boundaryConditionVarsResults);
3118 ✗ free(boundaryconditiondata.reconSt_diag);
3119 }
3120 ✗ return 0;
3121 }
3122
3123 /*
3124 * Runs the numerical procedure to compute constraint equation (D.1)
3125 */
3126 ✗ int dataReconciliation(DATA * data, threadData_t * threadData, int status)
3127 {
3128 // copy the reference files "AuxiliaryConditions and IntermediateEquations.html to output path"
3129 ✗ if (omc_flag[FLAG_OUTPUT_PATH])
3130 {
3131 ✗ copyReferenceFile(data, "_AuxiliaryConditions.html");
3132 ✗ copyReferenceFile(data, "_IntermediateEquations.html");
3133 ✗ copyReferenceFile(data, "_relatedBoundaryConditionsEquations.html");
3134 ✗ copyReferenceFile(data, "_iterationVars.txt");
3135 }
3136
3137 // report run time initialization and non linear convergence error to html
3138 ✗ if (status != 0)
3139 {
3140 ✗ createErrorHtmlReport(data, status);
3141 ✗ exit(1);
3142 }
3143
3144 const char * epselon = NULL;
3145 ✗ epselon = (char*) omc_flagValue[FLAG_DATA_RECONCILE_Eps];
3146
3147 // create a debug log file
3148 ✗ ofstream logfile;
3149 ✗ std::stringstream logfilename;
3150 ✗ if (omc_flag[FLAG_OUTPUT_PATH])
3151 {
3152 ✗ logfilename << omc_flagValue[FLAG_OUTPUT_PATH] << "/" << data->modelData->modelFilePrefix << "_debug.txt";
3153 }
3154 else
3155 {
3156 ✗ logfilename << data->modelData->modelFilePrefix << "_debug.txt";
3157 }
3158
3159 string tmplogfilename = logfilename.str();
3160 ✗ logfile.open(tmplogfilename.c_str());
3161
3162 ✗ if (omc_flag[FLAG_DATA_RECONCILE])
3163 {
3164 ✗ logfile << "| info | " << "DataReconciliation Starting!\n";
3165 ✗ logfile << "| info | " << data->modelData->modelName << "\n";
3166 }
3167 ✗ if (omc_flag[FLAG_DATA_RECONCILE_STATE])
3168 {
3169 ✗ logfile << "| info | " << "State Estimation Starting!\n";
3170 ✗ logfile << "| info | " << data->modelData->modelName << "\n";
3171 }
3172
3173 // set default value (epselon = 1.e-10), if no value provided by user
3174 ✗ if (epselon == NULL)
3175 {
3176 epselon = "0.0000000001";
3177 }
3178
3179 // read the measurement input data provide by user
3180 ✗ csvData csvdata = readMeasurementInputFile(logfile, data, threadData);
3181
3182 // validate the input data read from measurement input file
3183 ✗ csvData Sx_data = validateMeasurementInputs(csvdata, data, logfile);
3184
3185 // extracts the input data (x) from csvData
3186 ✗ inputData x = getInputData(Sx_data, logfile);
3187
3188 // read the correlation coefficient input data provide by user
3189 correlationDataWarning warningCorrelationData;
3190
3191 ✗ correlationData Cx_data = readCorrelationCoefficientFile(Sx_data, logfile, data, threadData, warningCorrelationData);
3192
3193 // Compute the covariance matrix (Sx) from csvData
3194 ✗ matrixData Sx = computeCovarianceMatrixSx(Sx_data, Cx_data, logfile, data);
3195
3196 // Compute the Jacobian Matrix F
3197 ✗ matrixData jacF = getJacobianMatrixF(data, threadData, logfile);
3198
3199 // Compute the Transpose of jacobian Matrix F
3200 ✗ matrixData jacFt = getTransposeMatrix(jacF);
3201
3202 ✗ double * Sx_diag = (double*) calloc(Sx.rows * 1, sizeof(double));
3203 ✗ getDiagonalElements(Sx.data, Sx.rows, Sx.column, Sx_diag);
3204 ✗ matrixData tmpSx_diag = {Sx.rows, 1, Sx_diag};
3205
3206 ✗ matrixData tmp_x = {x.rows, x.column, x.data};
3207 ✗ matrixData x_diag = copyMatrix(tmp_x);
3208
3209 // Print the initial information
3210 ✗ logfile << "\n\nInitial Data \n" << "=============\n";
3211 ✗ printMatrixWithHeaders(x.data, x.rows, x.column, Sx_data.headers, "X", logfile);
3212 ✗ printVectorMatrixWithHeaders(Sx_data.sxdata, Sx_data.rowcount, 1, Sx_data.headers, "Half-WidthConfidenceInterval", logfile);
3213 ✗ printCorelationMatrix(Cx_data.data, Cx_data.rowHeaders, Cx_data.columnHeaders, "Co-Relation_Coefficient", logfile, warningCorrelationData);
3214 ✗ printMatrixWithHeaders(Sx.data, Sx.rows, Sx.column, Sx_data.headers, "Sx", logfile);
3215
3216 // Start the Algorithm
3217 ✗ if (omc_flag[FLAG_DATA_RECONCILE])
3218 {
3219 dataReconciliationData datareconciliationdata;
3220 ✗ RunReconciliation(data, threadData, x, Sx, jacF, jacFt, atof(epselon), 1, Sx_data, x_diag, tmpSx_diag, logfile, warningCorrelationData, datareconciliationdata);
3221 ✗ logfile << "| info | " << "DataReconciliation Completed! \n";
3222 }
3223 ✗ if (omc_flag[FLAG_DATA_RECONCILE_STATE])
3224 {
3225 ✗ stateEstimation(data, threadData, x, Sx, jacF, jacFt, atof(epselon), 1, Sx_data, x_diag, tmpSx_diag, logfile, warningCorrelationData);
3226 ✗ logfile << "| info | " << "state estimation Completed! \n";
3227 }
3228 ✗ logfile.flush();
3229 ✗ logfile.close();
3230 ✗ free(Sx.data);
3231 ✗ free(x.data);
3232 ✗ free(jacF.data);
3233 ✗ free(jacFt.data);
3234 ✗ free(tmpSx_diag.data);
3235 ✗ free(x_diag.data);
3236 ✗ return 0;
3237 ✗ }
3238
3239 ✗ int boundaryConditions(DATA * data, threadData_t * threadData, int status)
3240 {
3241 // copy the reference files "BoundaryConditionsEquations.html and _BoundaryConditionIntermediateEquations.html to output path"
3242 ✗ if (omc_flag[FLAG_OUTPUT_PATH])
3243 {
3244 ✗ copyReferenceFile(data, "_BoundaryConditionsEquations.html");
3245 ✗ copyReferenceFile(data, "_BoundaryConditionIntermediateEquations.html");
3246 ✗ copyReferenceFile(data, "_iterationVars.txt");
3247 }
3248
3249 // report run time initialization and non linear convergence error to html
3250 ✗ if (status != 0)
3251 {
3252 ✗ createErrorHtmlReportForBoundaryConditions(data, status);
3253 ✗ exit(1);
3254 }
3255
3256 // create a debug log file
3257 ✗ ofstream logfile;
3258 ✗ std::stringstream logfilename;
3259 ✗ if (omc_flag[FLAG_OUTPUT_PATH])
3260 {
3261 ✗ logfilename << omc_flagValue[FLAG_OUTPUT_PATH] << "/" << data->modelData->modelFilePrefix << "_BoundaryConditions_debug.txt";
3262 }
3263 else
3264 {
3265 ✗ logfilename << data->modelData->modelFilePrefix << "_BoundaryConditions_debug.txt";
3266 }
3267
3268 string tmplogfilename = logfilename.str();
3269 ✗ logfile.open(tmplogfilename.c_str());
3270 ✗ logfile << "| info | " << "Reconcile Boundary Conditions Starting!\n";
3271 ✗ logfile << "| info | " << data->modelData->modelName << "\n";
3272
3273 // read the measurement input data provide by user
3274 ✗ csvData csvdata = readMeasurementInputFile(logfile, data, threadData, true);
3275
3276 // validate the input data read from measurement input file
3277 ✗ csvData Sx_data = validateMeasurementInputs(csvdata, data, logfile, true);
3278
3279 // extracts the input data (x) from csvData
3280 ✗ inputData reconciled_x = getReconciledX(Sx_data, logfile);
3281
3282 // read the reconciled covariance matrix input file provided by user
3283 correlationDataWarning warningCorrelationData;
3284 ✗ correlationData cx_data = readCorrelationCoefficientFile(Sx_data, logfile, data, threadData, warningCorrelationData, true);
3285
3286 // create the column matrix from the covariance matrix
3287 ✗ int rowsize = cx_data.rowHeaders.size();
3288 ✗ int colsize = cx_data.columnHeaders.size();
3289 ✗ double *tempSx = (double*) calloc(rowsize * colsize, sizeof(double));
3290 ✗ initColumnMatrix(cx_data.data, rowsize, colsize, tempSx);
3291 matrixData reconciled_Sx = {rowsize, colsize, tempSx};
3292
3293
3294 ✗ logfile << "\n\nInitial Data \n" << "=============\n";
3295 ✗ printMatrixWithHeaders(reconciled_x.data, reconciled_x.rows, reconciled_x.column, Sx_data.headers, "Reconciled_X", logfile);
3296 //printCorelationMatrix(reconciled_Sx.data, reconciled_Sx.rowHeaders, reconciled_Sx.columnHeaders, "Reconciled_Sx", logfile, warningCorrelationData);
3297 ✗ printMatrixWithHeaders(reconciled_Sx.data, reconciled_Sx.rows, reconciled_Sx.column, Sx_data.headers, "Reconciled_Sx", logfile);
3298 boundaryConditionData boundaryconditiondata;
3299 ✗ reconcileBoundaryConditions(data, threadData, reconciled_x, reconciled_Sx, boundaryconditiondata, warningCorrelationData, logfile);
3300 ✗ updateReconciledMo(data, threadData, Sx_data.headers, reconciled_x.data, logfile);
3301
3302 ✗ logfile << "*****Completed***********\n";
3303 ✗ logfile << "| info | " << "Reconcile Boundary Conditions Completed! \n";
3304 ✗ logfile.flush();
3305 ✗ logfile.close();
3306
3307 //free(reconciled_Sx.data);
3308 //free(reconciled_x.data);
3309 ✗ return 0;
3310 ✗ }
3311