Bug Summary

File:OMCompiler/SimulationRuntime/c/./simulation/../dataReconciliation/dataReconciliation.cpp
Warning:line 435, column 2
Declared variable-length array (VLA) has zero size

Annotated Source Code

[?] Use j/k keys for keyboard navigation

1/*
2 * This file is part of OpenModelica.
3 *
4 * Copyright (c) 1998-2010, Linköpings University,
5 * Department of Computer and Information Science,
6 * SE-58183 Linköping, Sweden.
7 *
8 * All rights reserved.
9 *
10 * THIS PROGRAM IS PROVIDED UNDER THE TERMS OF THIS OSMC PUBLIC
11 * LICENSE (OSMC-PL). ANY USE, REPRODUCTION OR DISTRIBUTION OF
12 * THIS PROGRAM CONSTITUTES RECIPIENT'S ACCEPTANCE OF THE OSMC
13 * PUBLIC LICENSE.
14 *
15 * The OpenModelica software and the Open Source Modelica
16 * Consortium (OSMC) Public License (OSMC-PL) are obtained
17 * from Linköpings University, either from the above address,
18 * from the URL: http://www.ida.liu.se/projects/OpenModelica
19 * and in the OpenModelica distribution.
20 *
21 * This program is distributed WITHOUT ANY WARRANTY; without
22 * even the implied warranty of MERCHANTABILITY or FITNESS
23 * FOR A PARTICULAR PURPOSE, EXCEPT AS EXPRESSLY SET FORTH
24 * IN THE BY RECIPIENT SELECTED SUBSIDIARY LICENSE CONDITIONS
25 * OF OSMC-PL.
26 *
27 * See the full OSMC Public License conditions for more details.
28 *
29 */
30
31#include "util/omc_error.h"
32#include "simulation_data.h"
33#include "openmodelica_func.h"
34#include "simulation/solver/external_input.h"
35#include "simulation/options.h"
36#include "simulation/solver/model_help.h"
37#include <iostream>
38#include <sstream>
39#include <string>
40#include <fstream>
41#include <vector>
42#include <algorithm>
43#include <iomanip>
44#include <stdlib.h>
45#include <math.h>
46#include <ctime>
47#include "omc_config.h"
48#include <cmath>
49#include "dataReconciliation.h"
50
51extern "C"
52{
53int dgesv_(int *n, int *nrhs, double *a, int *lda, int *ipiv, double *b, int *ldb, int *info);
54int dgemm_(char *transa, char *transb, int *m, int *n, int *k, double *alpha, double *a, int *lda,
55 double *b, int *ldb, double *beta, double *c, int *ldc);
56int dgetrf_(int *m, int *n, double *a, int *lda, int *ipiv, int *info);
57int dgetri_(int *n, double *a, int *lda, int *ipiv, double *work, int *lwork, int *info);
58int dscal_(int *n, double *da, double *dx, int *incx);
59int dcopy_(int *n, double *dx, int *incx, double *dy, int *incy);
60}
61
62
63using namespace std;
64
65extern "C" {
66
67// only 200 values of chisquared x^2 values are added with degree of freedom
68static 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};
69
70struct csvData {
71 int linecount;
72 int rowcount;
73 int columncount;
74 vector<double> xdata;
75 vector<double> sxdata;
76 vector<string> headers;
77 vector< vector<string> > rx;
78};
79
80struct matrixData {
81 int rows;
82 int column;
83 double * data;
84};
85
86struct inputData {
87 int rows;
88 int column;
89 double * data;
90 vector<int> index;
91};
92
93/*
94 * Function which reads the csv file
95 * and stores the covariance matrix Sx for DataReconciliation
96 */
97csvData readcsvfiles(const char * filename, ofstream & logfile)
98{
99 ifstream ip(filename);
100 string line;
101 vector<double> xdata;
102 vector<double> vals;
103 vector<string> names;
104 vector< vector<string> > rx;
105 int Sxrowcount=0;
106 int linecount=1;
107 int Sxcolscount=0;
108 bool flag=false;
109 int myarraycount=0;
110 if(!ip.good())
111 {
112 //errorStreamPrint(LOG_STDOUT, 0, "file name not found %s.",filename);
113 logfile << "| error | " << "file name not found " << filename << "\n";
114 logfile.close();
115 exit(1);
116 }
117 while(ip.good())
118 {
119 getline(ip,line);
120 if(linecount>1 && !line.empty())
121 {
122 //cout << "array info:" << line << "\n";
123 std::replace(line.begin(), line.end(), ';', ' ');
124 std::replace(line.begin(), line.end(), ',', ' ');
125 stringstream ss(line);
126 string temp;
127 int skip=0;
128 while(ss >> temp){
129 if(skip==0)
130 {
131 names.push_back(temp.c_str());
132 Sxrowcount++;
133 }
134 if(skip>0){
135 //cout << "check temp:" << temp << " double" << atof(temp.c_str()) <<"\n";
136 vals.push_back(atof(temp.c_str()));
137 if(flag==false){
138 Sxcolscount++;
139 }
140 }
141 skip++;
142 }
143 flag=true;
144 //Sxrowcount++;
145 }
146 linecount++;
147 }
148 //cout << "csvdata header:" << names[0] << names[1] << names[2] << "";
149 //cout << "linecount:" << linecount << " " << "rowcount :" << Sxrowcount << " " << "colscount:" << Sxcolscount << "\n";
150 csvData data={linecount,Sxrowcount,Sxcolscount,xdata,vals,names,rx};
151 return data;
152}
153/*
154 * function which returns the index pos
155 * of input variables
156 */
157int getVariableIndex(vector<string> headers, string name, ofstream & logfile)
158{
159 int pos=-1;
160 for(unsigned int i=0; i<headers.size(); i++)
161 {
162 //logfile << "founded headers " << headers[i] << i << "\n";
163 if(strcmp(headers[i].c_str(),name.c_str())==0)
164 {
165 pos = i;
166 break;
167 }
168 }
169 //logfile << "founded pos " << name << ": " << pos << "\n";
170 if(pos==-1)
171 {
172 //logfile << "Variable Name not Matched :" << name;
173 logfile << "| error | " << "CoRelation-Coefficient Variable Name not Matched: " << name << " ,getVariableIndex() failed!"<< "\n";
174 logfile.close();
175 exit(1);
176 }
177 return pos;
178}
179
180/*
181 * Function which reads the csv file
182 * and stores the initial measured value X and HalfWidth confidence
183 * interval Wx and also the input variable names
184 */
185csvData readcsvInputs(const char * filename, ofstream & logfile)
186{
187 ifstream ip(filename);
188 string line;
189 vector<double> xdata;
190 vector<double> sxdata;
191 vector<string> names;
192 //vector<double> rx_ik;
193 vector< vector<string> > rx;
194 int Sxrowcount=0;
195 int linecount=1;
196 int Sxcolscount=0;
197 bool flag=false,rx_ik=false,skip0=false,skip1=false,skip2=false;
198 int myarraycount=0;
199 if(!ip.good())
200 {
201 //errorStreamPrint(LOG_STDOUT, 0, "file name not found %s.",filename);
202 logfile << "| error | " << "file name not found " << filename << "\n";
203 logfile.close();
204 exit(1);
205 }
206 while(ip.good())
207 {
208 getline(ip,line);
209 vector<string> t1;
210 if(linecount>1 && !line.empty())
211 {
212 //logfile << "array info:" << line << "\n";
213 std::replace(line.begin(), line.end(), ';', ' ');
214 std::replace(line.begin(), line.end(), ',', ' ');
215 stringstream ss(line);
216 string temp;
217 int skip=0;
218 while(ss >> temp){
219 if(skip==0)
220 {
221 skip0=true;
222 names.push_back(temp.c_str());
223 Sxrowcount++;
224 if(flag==false){
225 Sxcolscount++;
226 }
227 }
228 if(skip==1)
229 {
230 skip1=true;
231 //logfile << "xdata" << temp << " double" << atof(temp.c_str()) <<"\n";
232 xdata.push_back(atof(temp.c_str()));
233 if(flag==false){
234 Sxcolscount++;
235 }
236 }
237 if(skip==2){
238 //logfile << "sxdata" << temp << " double" << atof(temp.c_str()) <<"\n";
239 skip2=true;
240 sxdata.push_back(atof(temp.c_str()));
241 if(flag==false){
242 Sxcolscount++;
243 }
244 }
245
246 if(skip>2)
247 {
248 //logfile << "found xi " << line << "\n";
249 t1.push_back(temp.c_str());
250 rx_ik=true;
251 }
252 skip++;
253 }
254 flag=true;
255 //Sxrowcount++;
256 if(skip0==false || skip1==false || skip2==false)
257 {
258 logfile << "| error | " << filename << "| csvdata Empty, " << "DataReconciliation cannot be computed ! \n";
259 logfile.close();
260 exit(1);
261 }
262 }
263 if(rx_ik==true)
264 {
265 rx.push_back(t1);
266 }
267 linecount++;
268 }
269 //logfile << "csvdata header:" << "header length: " << names.size() << " " << names[0] << names[1] << names[2] << "" << "\n";
270 //logfile << "linecount:" << linecount << " " << "rowcount :" << Sxrowcount << " " << "colscount:" << Sxcolscount << "\n";
271 csvData data={linecount,Sxrowcount,Sxcolscount,xdata,sxdata,names,rx};
272 return data;
273}
274
275/*
276 * Function which arranges the elements in column major
277 */
278void initColumnMatrix(vector<double> data, int rows, int cols, double * tempSx)
279{
280 for (int i=0; i<rows; i++)
281 {
282 for (int j=0; j<cols;j++)
283 {
284 // store the matrix in column order
285 tempSx[j+i*rows]=data[i+j*rows];
286 }
287 }
288}
289
290/*
291 * Function to print and debug whether the matrices are stored in column major
292 */
293void printColumnAlginment(double * matrix, int rows, int cols, string name)
294{
295 cout << "\n" << "************ "<< name << " **********" << "\n";
296 for (int i=0; i < rows*cols ; i++)
297 {
298 cout << matrix[i] << " ";
299 }
300 cout << "\n";
301}
302
303/*
304 * Function to Print the matrix in row based format
305 */
306void printMatrix(double* matrix, int rows, int cols, string name, ofstream& logfile)
307{
308 logfile << "\n" << "************ "<< name << " **********" <<"\n";
309 for (int i=0;i<rows; i++)
310 {
311 for (int j=0;j<cols;j++)
312 {
313 //cout << setprecision(5);
314 logfile << std::right << setw(15) << matrix[i+j*rows];
315 logfile.flush();
316 }
317 logfile << "\n";
318 }
319 logfile << "\n";
320}
321
322/*
323 *
324 Function to Print the matrix in row based format with headers
325 */
326void printMatrixWithHeaders(double* matrix, int rows, int cols, vector<string> headers, string name, ofstream& logfile)
327{
328 logfile << "\n" << "************ "<< name << " **********" <<"\n";
329 for (int i=0;i<rows; i++)
330 {
331 logfile << std::right << setw(10) << headers[i];
332 for (int j=0;j<cols;j++)
333 {
334 //cout << setprecision(5);
335 logfile << std::right << setw(15) << matrix[i+j*rows];
336 logfile.flush();
337 //printf("% .5e ", matrix[i+j*rows]);
338 }
339 logfile << "\n";
340 }
341 logfile << "\n";
342}
343
344/*
345 *Function to Print the vecomatrix in row based format with headers
346 *based on vector arrays
347 */
348void printVectorMatrixWithHeaders(vector<double> matrix, int rows, int cols, vector<string> headers, string name, ofstream& logfile)
349{
350 logfile << "\n" << "************ "<< name << " **********" <<"\n";
351 for (int i=0;i<rows; i++)
352 {
353 logfile << std::right << setw(10) << headers[i];
354 for (int j=0;j<cols;j++)
355 {
356 //cout << setprecision(5);
357 logfile << std::right << setw(15) << matrix[i+j*rows];
358 logfile.flush();
359 //printf("% .5e ", matrix[i+j*rows]);
360 }
361 logfile << "\n";
362 }
363 logfile << "\n";
364}
365
366/*
367 *
368 Function Which gets the diagonal elements of the matrix
369 */
370void getDiagonalElements(double *matrix, int rows, int cols, double* result)
371{
372 int k=0;
373 for (int i=0;i<rows; i++)
374 {
375 for (int j=0;j<cols;j++)
376 {
377 if(i==j)
378 {
379 result[k++] = matrix[i+j*rows];
380 }
381 }
382 }
383}
384
385/*
386 * Function to transpose the Matrix
387 */
388void transposeMatrix(double * jacF, double * jacFT, int rows, int cols)
389{
390 for (int i=0;i<rows; i++)
391 {
392 for (int j=0;j<cols;j++)
393 {
394 // Perform matrix transpose store the elements in column major
395 jacFT[i*cols+j]= jacF[i+j*rows] ;
396 }
397 }
398}
399
400
401/*
402 * Matrix Multiplication using dgemm LaPack routine
403 */
404void solveMatrixMultiplication(double *matrixA, double *matrixB, int rowsa, int colsa, int rowsb, int colsb , double *matrixC, ofstream & logfile)
405{
406 char trans = 'N';
407 double one = 1.0, zero = 0.0;
408 int rowsA = rowsa;
409 int colsA = colsa;
410 int rowsB = rowsb;
411 int colsB = colsb;
412 int common = colsa;
413
414 if(colsA!=rowsB)
415 {
416 //cout << "\n Error: Column of First Matrix not equal to Rows of Second Matrix \n ";
417 //errorStreamPrint(LOG_STDOUT, 0, "solveMatrixMultiplication() Failed!, Column of First Matrix not equal to Rows of Second Matrix %i != %i.",colsA,rowsB);
418 logfile << "| error | " << "solveMatrixMultiplication() Failed!, Column of First Matrix not equal to Rows of Second Matrix " << colsA << " != "<< rowsB << "\n";
419 logfile.close();
420 exit(1);
421 }
422 // solve matrix multiplication using dgemm_ LAPACK routine
423 dgemm_(&trans, &trans, &rowsA, &colsB, &common, &one, matrixA, &rowsA, matrixB, &common, &zero, matrixC, &rowsA);
424}
425
426/*
427 * Solve the Linear System A*x=b using LAPACK Solver routine dgesv_
428 */
429void solveSystemFstar(int n, int nhrs, double * tmpMatrixD, double * tmpMatrixC, ofstream & logfile)
430{
431 int N=n; // number of rows of Matrix A
8
'N' initialized to 0
432 int NRHS=nhrs; // number of columns of Matrix B
433 int LDA=N;
434 int LDB=N;
435 int ipiv[N];
9
Declared variable-length array (VLA) has zero size
436 int info;
437 // call the external function
438 dgesv_( &N, &NRHS, tmpMatrixD, &LDA, ipiv, tmpMatrixC, &LDB, &info);
439 if( info > 0 ) {
440 //cout << "The solution could not be computed, The info satus is : " << info;
441 //errorStreamPrint(LOG_STDOUT, 0, "solveSystemFstar() Failed !, The solution could not be computed, The info satus is %i.", info);
442 logfile << "| error | " << "solveSystemFstar() Failed !, The solution could not be computed, The info satus is" << info << "\n";
443 logfile.close();
444 exit(1);
445 }
446}
447
448/*
449 * Solve the matrix Subtraction of two matrices
450 */
451void solveMatrixSubtraction(matrixData A, matrixData B, double * result, ofstream & logfile)
452{
453 if(A.rows!=B.rows && A.column!=B.column)
454 {
455 //cout << "The Matrix Dimensions are not equal to Compute ! \n";
456 //errorStreamPrint(LOG_STDOUT, 0, "solveMatrixSubtraction() Failed !, The Matrix Dimensions are not equal to Compute ! %i != %i.", A.rows,B.rows);
457 logfile << "| error | " << "solveMatrixSubtraction() Failed !, The Matrix Dimensions are not equal to Compute" << A.rows << " != " << B.rows << "\n";
458 logfile.close();
459 exit(1);
460 }
461 //printColumnAlginment(A.data,A.rows,A.column,"A-Matrix");
462 //printColumnAlginment(B.data,B.rows,B.column,"B-Matrix");
463
464 // subtract elements in cloumn major
465 for(int i=0; i < A.rows*A.column; i++)
466 {
467 result[i]=A.data[i]-B.data[i];
468 }
469}
470
471/*
472 * Solve the matrix addition of two matrices
473 */
474matrixData solveMatrixAddition(matrixData A, matrixData B, ofstream & logfile)
475{
476 double* result = (double*)calloc(A.rows*A.column,sizeof(double));
477 if(A.rows!=B.rows && A.column!=B.column)
478 {
479 //cout << "The Matrix Dimensions are not equal to Compute ! \n";
480 //errorStreamPrint(LOG_STDOUT, 0, "solveMatrixAddition() Failed !, The Matrix Dimensions are not equal to Compute ! %i != %i.", A.rows,B.rows);
481 logfile << "| error | " << "solveMatrixAddition() Failed !, The Matrix Dimensions are not equal to Compute" << A.rows << " != " << B.rows << "\n";
482 logfile.close();
483 exit(1);
484 }
485 //printColumnAlginment(A.data,A.rows,A.column,"A-Matrix");
486 //printColumnAlginment(B.data,B.rows,B.column,"B-Matrix");
487 // Add the elements in cloumn major
488 for(int i=0; i < A.rows*A.column; i++)
489 {
490 result[i]=A.data[i]+B.data[i];
491 }
492 matrixData tmpadd_a_b = {A.rows,A.column,result};
493 return tmpadd_a_b;
494}
495
496/*
497 * Function which Calculates the Matrix Multiplication
498 * of (Sx*Ft)*Fstar
499 */
500matrixData Calculate_Sx_Ft_Fstar(matrixData Sx, matrixData Ft, matrixData Fstar, ofstream & logfile)
501{
502 // Sx*Ft
503 double* tmpMatrixA = (double*)calloc(Sx.rows*Ft.column,sizeof(double));
504 solveMatrixMultiplication(Sx.data,Ft.data,Sx.rows,Sx.column,Ft.rows,Ft.column,tmpMatrixA,logfile);
505 //printMatrix1(tmpMatrixA,Sx.rows,Ft.column,"Reconciled-(Sx*Ft)");
506 //printMatrix1(Fstar.data,Fstar.rows,Fstar.column,"REconciled-FStar");
507
508 //(Sx*Ft)*Fstar
509 double* tmpMatrixB = (double*)calloc(Sx.rows*Fstar.column,sizeof(double));
510 solveMatrixMultiplication(tmpMatrixA,Fstar.data,Sx.rows,Ft.column,Fstar.rows,Fstar.column,tmpMatrixB,logfile);
511 matrixData rhsdata= {Sx.rows,Fstar.column,tmpMatrixB};
512
513 free(tmpMatrixA);
514 free(tmpMatrixB);
515 return rhsdata;
516}
517
518/*
519 * Solves the system
520 * recon_x = x - (Sx*Ft*fstar)
521 */
522matrixData solveReconciledX(matrixData x, matrixData Sx, matrixData Ft, matrixData Fstar, ofstream& logfile)
523{
524 // Sx*Ft
525 double* tmpMatrixAf = (double*)calloc(Sx.rows*Ft.column,sizeof(double));
526 solveMatrixMultiplication(Sx.data,Ft.data,Sx.rows,Sx.column,Ft.rows,Ft.column,tmpMatrixAf,logfile);
527 //printMatrix(tmpMatrixAf,Sx.rows,Ft.column,"Sx*Ft");
528 //(Sx*Ft)*fstar
529 double* tmpMatrixBf = (double*)calloc(Sx.rows*Fstar.column,sizeof(double));
530 solveMatrixMultiplication(tmpMatrixAf,Fstar.data,Sx.rows,Ft.column,Fstar.rows,Fstar.column,tmpMatrixBf,logfile);
531 //printMatrix(tmpMatrixBf,Sx.rows,Fstar.column,"(Sx*Ft*fstar)");
532
533 matrixData rhs= {Sx.rows,Fstar.column,tmpMatrixBf};
534 //matrixData rhs = Calculate_Sx_Ft_Fstar(Sx,Ft,Fstar);
535
536 double* reconciledX = (double*)calloc(x.rows*x.column,sizeof(double));
537 solveMatrixSubtraction(x,rhs,reconciledX,logfile);
538 //printMatrix(reconciledX,x.rows,x.column,"reconciled X^cap ===> (x - (Sx*Ft*fstar))");
539 if(ACTIVE_STREAM(LOG_JAC)(useStream[LOG_JAC]))
540 {
541 logfile << "Calculations of Reconciled_x ==> (x - (Sx*Ft*f*))" << "\n";
542 logfile << "====================================================";
543 printMatrix(tmpMatrixAf,Sx.rows,Ft.column,"Sx*Ft",logfile);
544 printMatrix(tmpMatrixBf,Sx.rows,Fstar.column,"(Sx*Ft*f*)",logfile);
545 printMatrix(reconciledX,x.rows,x.column,"x - (Sx*Ft*f*))",logfile);
546 logfile << "***** Completed ****** \n\n";
547 }
548 matrixData recon_x = {x.rows,x.column,reconciledX};
549 //free(reconciledX);
550 free(tmpMatrixAf);
551 free(tmpMatrixBf);
552 return recon_x;
553}
554
555/*
556 * Solves the system
557 * recon_Sx = Sx - (Sx*Ft*Fstar)
558 */
559matrixData solveReconciledSx(matrixData Sx, matrixData Ft, matrixData Fstar, ofstream& logfile)
560{
561 // Sx*Ft
562 double* tmpMatrixA = (double*)calloc(Sx.rows*Ft.column,sizeof(double));
563 solveMatrixMultiplication(Sx.data,Ft.data,Sx.rows,Sx.column,Ft.rows,Ft.column,tmpMatrixA, logfile);
564 //printMatrix(tmpMatrixA,Sx.rows,Ft.column,"Reconciled-(Sx*Ft)");
565 //printMatrix(Fstar.data,Fstar.rows,Fstar.column,"REconciled-FStar");
566
567 //(Sx*Ft)*Fstar
568 double* tmpMatrixB = (double*)calloc(Sx.rows*Fstar.column,sizeof(double));
569 solveMatrixMultiplication(tmpMatrixA,Fstar.data,Sx.rows,Ft.column,Fstar.rows,Fstar.column,tmpMatrixB, logfile);
570 //printMatrix(tmpMatrixB,Sx.rows,Fstar.column,"Reconciled-(Sx*Ft*Fstar)");
571
572 matrixData rhs= {Sx.rows,Fstar.column,tmpMatrixB};
573
574 //matrixData rhs = Calculate_Sx_Ft_Fstar(Sx,Ft,Fstar);
575 double* reconciledSx = (double*)calloc(Sx.rows*Sx.column,sizeof(double));
576 solveMatrixSubtraction(Sx,rhs,reconciledSx, logfile);
577 //printMatrix(reconciledSx,Sx.rows,Sx.column,"reconciled Sx ===> (Sx - (Sx*Ft*Fstar))");
578 if(ACTIVE_STREAM(LOG_JAC)(useStream[LOG_JAC]))
579 {
580 logfile << "Calculations of Reconciled_Sx ===> (Sx - (Sx*Ft*F*))" << "\n";
581 logfile << "============================================";
582 printMatrix(tmpMatrixA,Sx.rows,Ft.column,"(Sx*Ft)",logfile);
583 printMatrix(tmpMatrixB,Sx.rows,Fstar.column,"(Sx*Ft*F*)",logfile);
584 printMatrix(reconciledSx,Sx.rows,Sx.column,"Sx - (Sx*Ft*F*))",logfile);
585 logfile << "***** Completed ****** \n\n";
586 }
587 matrixData recon_sx ={Sx.rows,Sx.column,reconciledSx};
588 //free(reconciledSx);
589 free(tmpMatrixA);
590 free(tmpMatrixB);
591 return recon_sx;
592}
593
594/*
595 * Function Which Computes the
596 * Jacobian Matrix F
597 */
598matrixData getJacobianMatrixF(DATA* data, threadData_t *threadData, ofstream & logfile)
599{
600 // initialize the jacobian call
601 const int index = data->callback->INDEX_JAC_F;
602 ANALYTIC_JACOBIAN* jacobian = &(data->simulationInfo->analyticJacobians[index]);
603 data->callback->initialAnalyticJacobianF(data, threadData, jacobian);
604 int cols = jacobian->sizeCols;
605 int rows = jacobian->sizeRows;
606 if(cols == 0) {
607 //errorStreamPrint(LOG_STDOUT, 0, "Cannot Compute Jacobian Matrix F");
608 logfile << "| error | " << "Cannot Compute Jacobian Matrix F" << "\n";
609 logfile.close();
610 exit(1);
611 }
612 double* jacF = (double*)calloc(rows*cols,sizeof(double)); // allocate for Matrix F
613 int k=0;
614 for (int x=0; x < cols ; x++)
615 {
616 jacobian->seedVars[x] = 1.0;
617 data->callback->functionJacF_column(data, threadData, jacobian, NULL__null);
618 //cout << "Calculate one column\n:";
619 for (int y=0; y < rows ; y++)
620 {
621 jacF[k++]=jacobian->resultVars[y];
622 }
623 jacobian->seedVars[x] = 0.0;
624 }
625 matrixData Fdata ={rows,cols,jacF};
626 return Fdata;
627}
628
629/*
630 * Function Which Computes the
631 * Transpose of Jacobian Matrix FT
632 */
633matrixData getTransposeMatrix(matrixData jacF)
634{
635 int rows=jacF.column;
636 int cols=jacF.rows;
637 double* jacFT = (double*)calloc(rows*cols,sizeof(double)); // allocate for Matrix F-transpose
638 int k=0;
639 for (int i=0;i<jacF.rows; i++)
640 {
641 for (int j=0;j<jacF.column;j++)
642 {
643 // Perform matrix transpose store the elements in column major
644 //cout << (i1*jacF.rows+j1) << " index :" << (i1+j1*jacF.rows) << " value is: " << jacF.data[i1+j1*jacF.rows] << "\n";
645 jacFT[k++]= jacF.data[i+j*jacF.rows];
646
647 }
648 }
649 matrixData Ft_data ={rows,cols,jacFT};
650 return Ft_data;
651}
652
653/*
654 * function which checks and reads
655 * covariance matrix Sx from csv files
656 * and stores the data in vector format
657 */
658csvData readCovarianceMatrixSx(DATA* data, threadData_t *threadData, ofstream & logfile)
659{
660 char * Sxfile = NULL__null;
661 Sxfile = (char*)omc_flagValue[FLAG_DATA_RECONCILE_Sx];
662 if(Sxfile==NULL__null)
663 {
664 //errorStreamPrint(LOG_STDOUT, 0, "Sx file not given (eg:-sx=filename.csv), DataReconciliation cannot be computed!.");
665 logfile << "| error | " << "Sx file not given (eg:-sx=filename.csv), DataReconciliation cannot be computed!.\n";
666 logfile.close();
667 exit(1);
668 }
669 //csvData Sx_result=readcsvfiles(Sxfile,logfile);
670 csvData Sx_result=readcsvInputs(Sxfile,logfile);
671 return Sx_result;
672}
673
674/*
675 * Function which reads the vector
676 * and assign to c pointer arrays
677 */
678matrixData getCovarianceMatrixSx(csvData Sx_result, DATA* data, threadData_t *threadData)
679{
680 double* tempSx = (double*)calloc(Sx_result.rowcount*Sx_result.columncount,sizeof(double));
681 initColumnMatrix(Sx_result.sxdata , Sx_result.rowcount, Sx_result.columncount, tempSx);
682 matrixData Sx_data = {Sx_result.rowcount,Sx_result.columncount,tempSx};
683 return Sx_data;
684}
685
686/*
687 * Function which Computes
688 * covariance matrix Sx based on
689 * Half width confidence interval provided by user
690 * Sx=(Wxi/1.96)^2
691 */
692
693matrixData computeCovarianceMatrixSx(csvData Sx_result, DATA* data, threadData_t *threadData, ofstream & logfile)
694{
695 double* tempSx = (double*)calloc(Sx_result.sxdata.size()*Sx_result.sxdata.size(),sizeof(double));
696 vector<double> tmpdata;
697 int k=0;
698 for (unsigned int i=0;i<Sx_result.sxdata.size(); i++)
699 {
700 double data = pow(Sx_result.sxdata[k]/1.96,2);
701 for (unsigned int j=0;j<Sx_result.sxdata.size();j++)
702 {
703 if(i==j)
704 {
705 //tmpdata.push_back(pow(Sx_result.sxdata[k]/1.96,2));
706 //k++;
707 tmpdata.push_back(data);
708 }
709 else
710 {
711 tmpdata.push_back(0);
712 }
713 // logfile << " data " << count << "=="<< tmpdata[count++] << "\n";
714 }
715 k++;
716 }
717 //logfile << "tmpdatasize" << tmpdata.size() << "\n";
718 //logfile << "Size of vector :" << Sx_result.rx.size() << "\n";
719
720 /* check for corelation coefficient matrix and insert the elements in correct position*/
721 for (unsigned int l=0; l < Sx_result.rx.size();l++)
722 {
723 int pos1;
724 int pos2;
725 double xi;
726 double xk ;
727 for(unsigned int m=0; m<Sx_result.rx[l].size();m++)
728 {
729 if(m==0)
730 {
731 pos1 = getVariableIndex(Sx_result.headers,Sx_result.rx[l][m],logfile);
732 xi = tmpdata[(Sx_result.rowcount*pos1)+pos1];
733 //logfile << "xi =>"<< pos1 << "= "<< xi << "\n";
734 }
735
736 if(m==1)
737 {
738 pos2 = getVariableIndex(Sx_result.headers,Sx_result.rx[l][m],logfile);
739 xk = tmpdata[(Sx_result.rowcount*pos2)+pos2];
740 //logfile << "xk =>"<< pos2 << "= "<< xk << "\n";
741 }
742 if(m==2)
743 {
744 //logfile << "position:" << pos1 << ": " << pos2 << "\n";
745 //logfile << "rx_ik" << Sx_result.rx[l][m] << "*" << xi << "*" << xk << "\n";
746 //logfile << atof((Sx_result.rx[l][m]).c_str())*sqrt(xi)*sqrt(xk) << "\n";
747 double tmprx = atof((Sx_result.rx[l][m]).c_str())*sqrt(xi)*sqrt(xk);
748 // find the symmetric position and insert the elements
749 //logfile << "final position :" << (Sx_result.rowcount*pos1)+pos2 << "value is: "<< tmprx << "\n";
750 //logfile << "final position :" << (Sx_result.rowcount*pos2)+pos1 << "value is: "<< tmprx << "\n";
751 tmpdata[(Sx_result.rowcount*pos1)+pos2]=tmprx;
752 tmpdata[(Sx_result.rowcount*pos2)+pos1]=tmprx;
753 }
754 }
755 //logfile << "\n";
756 }
757 initColumnMatrix(tmpdata , Sx_result.rowcount, Sx_result.rowcount, tempSx);
758 matrixData Sx_data = {Sx_result.rowcount,Sx_result.rowcount,tempSx};
759 return Sx_data;
760}
761
762/*
763 * Function which reads the input data X from start Attribute
764 * and also stores the index of input variables which are the
765 * variables to be reconciled for Data Reconciliation
766 */
767inputData getInputDataFromStartAttribute(csvData Sx_result , DATA* data, threadData_t *threadData, ofstream & logfile)
768{
769 double *tempx = (double*)calloc(Sx_result.rowcount,sizeof(double));
770 char ** knowns = (char**)malloc(data->modelData->nInputVars * sizeof(char*));
771 vector<int> index;
772 data->callback->inputNames(data, knowns);
773 int headercount = Sx_result.headers.size();
774 /* Read data from input vars which has start attribute value set as input */
775
776 for (int h=0; h < headercount; h++)
777 {
778 tempx[h]=Sx_result.xdata[h];
779 /*
780 bool flag=false;
781 for (int in=0; in < data->modelData->nInputVars; in++)
782 {
783 if(strcmp(knowns[in], Sx_result.headers[h].c_str()) == 0)
784 {
785 //tempx[h] = data->simulationInfo->inputVars[in];
786 index.push_back(in);
787 flag=true;
788 //logfile << knowns[in] << " start value :" << data->simulationInfo->inputVars[in] << "\n";
789 //logfile << "fetch index :" << in << "\n";
790 }
791 }
792 if(flag==false)
793 {
794 logfile << "| error | " << "Input Variable Not matched or not generated: "<< Sx_result.headers[h] << " , getInputDataFromStartAttribute failed()! \n";
795 logfile.close();
796 exit(1);
797 } */
798 }
799 inputData x_data ={Sx_result.rowcount,1,tempx,index};
800 free(knowns);
801 return x_data;
802}
803
804/*
805 * Function which Copy Matrix
806 * using dcopy_ LAPACK routine
807 * this is mostly used when LAPACK routines override arrays
808 */
809matrixData copyMatrix(matrixData matdata)
810{
811 double * tmpcopymatrix = (double*)calloc(matdata.rows*matdata.column,sizeof(double));
812 int n = matdata.rows*matdata.column;
813 int inc = 1;
814 dcopy_(&n,matdata.data,&inc,tmpcopymatrix,&inc);
815
816 // for (int i=0; i < matdata.rows*matdata.column; i++)
817 // {
818 // tmpcopymatrix[i]=matdata.data[i];
819 // }
820 matrixData tmpcopymatrixdata = {matdata.rows, matdata.column, tmpcopymatrix};
821 return tmpcopymatrixdata;
822}
823
824/*
825 * Function which scales the MAtrix with constant
826 * dscal_ LAPACK_routine and result is updated in data
827 * eg: alpha = 2, data=[2,4,5]
828 * result = [4,8,10]
829 */
830void scaleVector(int rows, int cols, double alpha, double * data)
831{
832 int n=rows*cols;
833 int inc=1;
834 dscal_(&n, &alpha, data, &inc);
835}
836
837/*
838 * Function which calculates the square root of elements
839 * eg : a=[1,2,3,4]
840 * result a = [srt(1),sqrt(2).......]
841 */
842void calculateSquareRoot(double * data, int length)
843{
844 for(int i=0; i<length; i++)
845 {
846 data[i]=sqrt(data[i]);
847 }
848}
849
850/*
851 * Function which calculates
852 * J*=(recon_x-x)T*(Sx^-1)*(recon_x-x)+2.[f+F*(recon_x-x)]T*fstar
853 * where T= transpose of matrix
854 * and returns the converged value
855 */
856
857double 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)
858{
859
860 //printMatrix(conv_vector_c.data,conv_vector_c.rows,conv_vector_c.column,"Convergence_C(x,y)");
861 //printMatrix(conv_fstar.data,conv_fstar.rows,conv_fstar.column,"Convergence_f*");
862 //printMatrix(conv_recon_x.data,conv_recon_x.rows,conv_recon_x.column,"check_recon_x*");
863
864 // calculate(recon_x-x)
865 double* conv_data1 = (double*)calloc(conv_x.rows*conv_x.column,sizeof(double));
866 matrixData conv_inputs = {conv_x.rows,conv_x.column,conv_x.data};
867 solveMatrixSubtraction(conv_recon_x,conv_inputs,conv_data1,logfile);
868 matrixData conv_data1result={conv_x.rows,conv_x.column,conv_data1};
869 matrixData copy_reconx_x = copyMatrix(conv_data1result);
870
871 //printMatrix(conv_inputs.data,conv_inputs.rows,conv_inputs.column,"check_inputs");
872 //printMatrix(conv_data1result.data,conv_data1result.rows,conv_data1result.column,"(recon_X - X)");
873
874 // calculate Transpose_(recon_x-x)
875 matrixData conv_data1Transpose = getTransposeMatrix(conv_data1result);
876 //printMatrix(conv_data1Transpose.data,conv_data1Transpose.rows,conv_data1Transpose.column,"Transpose(recon_X - X)");
877
878 /* solves (Sx^-1)*(recon_x-x)
879 * Solve the inverse of matrix Sx using linear form
880 * Ax=b
881 * where A=Sx and b= (recon_x-x) to avoid inversion of Sx which is
882 * expensive
883 */
884 solveSystemFstar(conv_sx.rows,1,conv_sx.data,conv_data1result.data,logfile);
885 //printMatrix(conv_data1result.data,conv_sx.rows,conv_data1result.column,"inverse multiplication_without inverse");
886
887 double *conv_tmpmatrixlhs = (double*)calloc(conv_data1Transpose.rows*conv_data1result.column,sizeof(double));
888 /*
889 * Solve (recon_x-x)T*(Sx^-1)*(recon_x-x)
890 */
891 solveMatrixMultiplication(conv_data1Transpose.data,conv_data1result.data,conv_data1Transpose.rows,conv_data1Transpose.column,conv_data1result.rows,conv_data1result.column,conv_tmpmatrixlhs,logfile);
892 //printMatrix(conv_tmpmatrixlhs,conv_data1Transpose.rows,conv_data1result.column,"(recon_x-x)T*(Sx^-1)*(recon_x-x)");
893 matrixData struct_conv_tmpmatrixlhs = {conv_data1Transpose.rows,conv_data1result.column,conv_tmpmatrixlhs};
894
895 /*
896 * Solve rhs = 2.[f+F*(recon_x-x)]T*fstar
897 *
898 */
899 // Calculate F*(recon_x-x)
900 double * tmp_F_recon_x_x = (double*)calloc(conv_jacF.rows*copy_reconx_x.column,sizeof(double));
901 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);
902 //printMatrix(tmp_F_recon_x_x,conv_jacF.rows,copy_reconx_x.column,"F*(recon_x-x)");
903 matrixData mult_F_recon_x_x = {conv_jacF.rows,copy_reconx_x.column,tmp_F_recon_x_x};
904
905 // Calculate f + F*(recon_x-x)
906 matrixData add_f_F_recon_x_x = solveMatrixAddition(conv_vector_c, mult_F_recon_x_x, logfile);
907 //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)");
908
909 matrixData transpose_add_f_F_recon_x_x = getTransposeMatrix(add_f_F_recon_x_x);
910 //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)]");
911
912 // calculate [f + F*(recon_x-x)]T*fstar
913 double *conv_tmpmatrixrhs = (double*)calloc(transpose_add_f_F_recon_x_x.rows*conv_fstar.column,sizeof(double));
914 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);
915 //printMatrix(conv_tmpmatrixrhs, transpose_add_f_F_recon_x_x.rows, conv_fstar.column,"[f + F*(recon_x-x)]*fstar");
916
917 // scale the matrix with 2*[f + F*(recon_x-x)]T*fstar
918 scaleVector(transpose_add_f_F_recon_x_x.rows , conv_fstar.column, 2.0, conv_tmpmatrixrhs);
919 //printMatrix(conv_tmpmatrixrhs, transpose_add_f_F_recon_x_x.rows, conv_fstar.column,"2*[f + F*(recon_x-x)]*fstar");
920 matrixData struct_conv_tmpmatrixrhs= {transpose_add_f_F_recon_x_x.rows, conv_fstar.column,conv_tmpmatrixrhs};
921
922 /*
923 * solve the final J*=J*=(recon_x-x)T*(Sx^-1)*(recon_x-x)+2.[f+F*(recon_x-x)]T*fstar
924 * J*=_struct_conv_tmpmatrixlhs + struct_conv_tmpmatrixrhs
925 */
926 matrixData struct_Jstar= solveMatrixAddition(struct_conv_tmpmatrixlhs,struct_conv_tmpmatrixrhs, logfile);
927 //printMatrix(struct_Jstar.data,struct_Jstar.rows,struct_Jstar.column,"J*",logfile);
928
929 int r=data->modelData->nSetcVars; // number of setc equations
930 double val=1.0/r;
931
932 /*
933 * calculate J/r < epselon
934 */
935 scaleVector(struct_Jstar.rows,struct_Jstar.column,val,struct_Jstar.data);
936 //printMatrix(struct_Jstar.data,struct_Jstar.rows,struct_Jstar.column,"J*/r ");
937 double convergedvalue=struct_Jstar.data[0];
938
939 // free(struct_Jstar.data);
940 // free(struct_conv_tmpmatrixrhs.data);
941 // free(conv_tmpmatrixrhs);
942 // free(transpose_add_f_F_recon_x_x.data);
943 // free(add_f_F_recon_x_x.data);
944 // free(transpose_add_f_F_recon_x_x.data);
945 // free(mult_F_recon_x_x.data);
946 // free(tmp_F_recon_x_x);
947 // free(struct_conv_tmpmatrixlhs.data);
948 // free(conv_tmpmatrixlhs);
949 // free(conv_data1Transpose.data);
950 // free(copy_reconx_x.data);
951 // free(conv_data1result.data);
952 // free(conv_data1);
953
954 return convergedvalue;
955}
956
957/*
958 * Example Function which performs matrix inverse
959 * using dgetri_ and dgetrf_ LAPACK routine
960 * which is expensive one and not recommended
961 * use it when no other way to compute it
962 */
963void checkExpensiveMatrixInverse()
964{
965 double newval[3*3]={3,2,0,
966 0,0,1,
967 2,-2,1};
968
969 int N=3;
970 int LDA=N;
971 int LDB=N;
972 int ipiv[N];
973 int info=1;
974 int LWORK =N;
975 double * WORK = (double*)calloc(LWORK,sizeof(double));
976
977 dgetrf_(&N,&N,newval,&N,ipiv,&info);
978 dgetri_(&N,newval,&N,ipiv,WORK,&LWORK,&info);
979 //printMatrix(newval,3,3,"Expensive_Matrix_Inverse");
980}
981
982/*
983 * Function which performs matrix inverse without performing
984 * actual matrix inverse, Instead use the dgesv to get result
985 * Ax=b where matrix mutiplication of x=bA gives the inversed
986 * mutiplication result b with A inverse
987 */
988void checkInExpensiveMatrixInverse(ofstream & logfile)
989{
990 double newchecksx[3*3]={1,1,1,
991 0,0.95,0,
992 0,0,0.95};
993 double checksx[3*1]={-0.028,0.026,-0.004};
994 solveSystemFstar(3,1,newchecksx,checksx,logfile);
995 //printMatrix(checksx,3,1,"InExpensive_Matrix_Inverse");
996}
997
998int 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)
999{
1000 // set the inputs first
1001 /*
1002 for (int i=0; i< x.rows*x.column; i++)
1003 {
1004 data->simulationInfo->inputVars[x.index[i]]=x.data[i];
1005 //logfile << "input data:" << x.data[i]<<"\n";
1006 }*/
1007
1008 for (int i=0; i< x.rows*x.column; i++)
1
Loop condition is false. Execution continues on line 1017
1009 {
1010 data->simulationInfo->datainputVars[i]=x.data[i];
1011 //logfile << "input data:" << x.data[i]<<"\n";
1012 }
1013
1014 /* set the inputs via this special function generated for dataReconciliation
1015 * which also sets inputs for models not involving top level inputs
1016 */
1017 data->callback->data_function(data, threadData);
1018 //data->callback->input_function(data, threadData);
1019 data->callback->functionDAE(data,threadData);
1020 //data->callback->functionODE(data,threadData);
1021 data->callback->setc_function(data, threadData);
1022
1023 matrixData jacF = getJacobianMatrixF(data,threadData,logfile);
2
The value 0 is assigned to 'jacF.rows'
1024 matrixData jacFt = getTransposeMatrix(jacF);
1025
1026 printMatrix(jacF.data,jacF.rows,jacF.column,"F",logfile);
1027 printMatrix(jacFt.data,jacFt.rows,jacFt.column,"Ft",logfile);
1028
1029
1030 double* setc = (double*)calloc(data->modelData->nSetcVars,sizeof(double)); // allocate data for setc array
1031 // store the setc data to compute for convergence as setc will be overriddeen with new values
1032 double* tmpsetc = (double*)calloc(data->modelData->nSetcVars,sizeof(double));
1033
1034 /* loop to store the data C(x,y) rhs side, get the elements in reverse order */
1035 int t=0;
1036 for (int i=data->modelData->nSetcVars; i > 0; i--)
3
Assuming 'i' is <= 0
4
Loop condition is false. Execution continues on line 1044
1037 {
1038 setc[t]=data->simulationInfo->setcVars[i-1];
1039 tmpsetc[t]=data->simulationInfo->setcVars[i-1];
1040 t++;
1041 //cout << "array_setc_vars:=>" << t << ":" << data->simulationInfo->setcVars[i-1] << "\n";
1042 }
1043
1044 int nsetcvars= data->modelData->nSetcVars;
1045 matrixData vector_c = {nsetcvars,1,tmpsetc};
1046 //allocate data for matrix multiplication F*Sx
1047 double * tmpmatrixC = (double*)calloc(jacF.rows*Sx.column,sizeof(double));
1048 solveMatrixMultiplication(jacF.data,Sx.data,jacF.rows,jacF.column,Sx.rows,Sx.column,tmpmatrixC,logfile);
1049 //printMatrix(tmpmatrixC,jacF.rows,Sx.column,"F*Sx");
1050
1051 //allocate data for matrix multiplication (F*Sx)*Ftranspose
1052 double * tmpmatrixD = (double*)calloc(jacF.rows*jacFt.column,sizeof(double));
1053 solveMatrixMultiplication(tmpmatrixC,jacFt.data,jacF.rows,Sx.column,jacFt.rows,jacFt.column,tmpmatrixD,logfile);
1054 //printMatrix(tmpmatrixD,jacF.rows,jacFt.column,"F*Sx*Ft");
1055 //printMatrix(setc,nsetcvars,1,"c(x,y)");
1056 /*
1057 * Copy tmpmatrixC and tmpmatrixD to avoid loss of data
1058 * when calculating F*
1059 */
1060 matrixData cpytmpmatrixC = {jacF.rows,Sx.column,tmpmatrixC};
1061 matrixData cpytmpmatrixD = {jacF.rows,jacFt.column,tmpmatrixD};
1062 matrixData tmpmatrixC1 = copyMatrix(cpytmpmatrixC);
1063 matrixData tmpmatrixD1 = copyMatrix(cpytmpmatrixD);
1064
1065 if(ACTIVE_STREAM(LOG_JAC)(useStream[LOG_JAC]))
5
Taking false branch
1066 {
1067 logfile << "Calculations of Matrix (F*Sx*Ft) f* = c(x,y) " << "\n";
1068 logfile << "============================================\n";
1069 printMatrix(tmpmatrixC,jacF.rows,Sx.column,"F*Sx",logfile);
1070 printMatrix(tmpmatrixD,jacF.rows,jacFt.column,"F*Sx*Ft",logfile);
1071 printMatrix(setc,nsetcvars,1,"c(x,y)",logfile);
1072 }
1073 /*
1074 * calculate f* for covariance matrix (F*Sx*Ftranspose).F*= c(x,y)
1075 * matrix setc will be overridden with new values which is the output
1076 * for the calculation A *x =B
1077 * A = tmpmatrixD
1078 * B = setc
1079 */
1080 solveSystemFstar(jacF.rows,1,tmpmatrixD,setc,logfile);
6
Passing the value 0 via 1st parameter 'n'
7
Calling 'solveSystemFstar'
1081 //printMatrix(setc,jacF.rows,1,"f*");
1082 if(ACTIVE_STREAM(LOG_JAC)(useStream[LOG_JAC]))
1083 {
1084 printMatrix(setc,jacF.rows,1,"f*",logfile);
1085 logfile << "***** Completed ****** \n\n";
1086 }
1087 matrixData tmpxcap ={x.rows,1,x.data};
1088 matrixData tmpfstar = {jacF.rows,1,setc};
1089 matrixData reconciled_X = solveReconciledX(tmpxcap,Sx,jacFt,tmpfstar,logfile);
1090 //printMatrix(reconciled_X.data,reconciled_X.rows,reconciled_X.column,"reconciled_X ===> (x - (Sx*Ft*fstar))");
1091 if(ACTIVE_STREAM(LOG_JAC)(useStream[LOG_JAC]))
1092 {
1093 logfile << "Calculations of Matrix (F*Sx*Ft) F* = F*Sx " << "\n";
1094 logfile << "===============================================\n";
1095 //printMatrix(tmpmatrixC,jacF.rows,Sx.column,"F*Sx");
1096 //printMatrix(tmpmatrixD,jacF.rows,jacFt.column,"F*Sx*Ft");
1097 printMatrix(tmpmatrixC1.data,tmpmatrixC1.rows,tmpmatrixC1.column,"F*Sx",logfile);
1098 printMatrix(tmpmatrixD1.data,tmpmatrixD1.rows,tmpmatrixD1.column,"F*Sx*Ft",logfile);
1099 }
1100 /*
1101 * calculate F* for covariance matrix (F*Sx*Ftranspose).F*= (F*Sx)
1102 * tmpmatrixC1 will be overridden with new values which is the output
1103 * for the calculation A *x =B
1104 * A = tmpmatrixD
1105 * B = tmpmatrixC
1106 */
1107 solveSystemFstar(jacF.rows,Sx.column,tmpmatrixD1.data,tmpmatrixC1.data,logfile);
1108 //printMatrix(tmpmatrixC,jacF.rows,Sx.column,"Sx_F*");
1109 if(ACTIVE_STREAM(LOG_JAC)(useStream[LOG_JAC]))
1110 {
1111 printMatrix(tmpmatrixC1.data,jacF.rows,Sx.column,"F*",logfile);
1112 logfile << "***** Completed ****** \n\n";
1113 }
1114
1115 matrixData tmpFstar = {jacF.rows,Sx.column,tmpmatrixC1.data};
1116 matrixData reconciled_Sx = solveReconciledSx(Sx,jacFt,tmpFstar,logfile);
1117 //printMatrix(reconciled_Sx.data,reconciled_Sx.rows,reconciled_Sx.column,"reconciled Sx ===> (Sx - (Sx*Ft*Fstar))");
1118
1119 double value = solveConvergence(data,reconciled_X,reconciled_Sx,x,Sx,jacF,vector_c,tmpfstar,logfile);
1120 if(value > eps )
1121 {
1122 logfile << "J*/r" << "(" << value << ")" << " > " << eps << ", Value not Converged \n";
1123 logfile << "==========================================\n\n";
1124 }
1125 if(value > eps)
1126 {
1127 logfile << "Running Convergence iteration: " << iterationcount << " with the following reconciled values:" << "\n";
1128 logfile << "========================================================================" << "\n";
1129 //cout << "J*/r :=" << value << "\n";
1130 //printMatrix(jacF.data,jacF.rows,jacF.column,"F");
1131 //printMatrix(jacFt.data,jacFt.rows,jacFt.column,"Ft");
1132 //printMatrix(setc,jacF.rows,1,"f*");
1133 //printMatrix(tmpmatrixC,jacF.rows,Sx.column,"F*");
1134 printMatrixWithHeaders(reconciled_X.data,reconciled_X.rows,reconciled_X.column,csvinputs.headers,"reconciled_X ===> (x - (Sx*Ft*fstar))",logfile);
1135 printMatrixWithHeaders(reconciled_Sx.data,reconciled_Sx.rows,reconciled_Sx.column,csvinputs.headers,"reconciled_Sx ===> (Sx - (Sx*Ft*Fstar))",logfile);
1136 //free(x.data);
1137 //free(Sx.data);
1138 x.data=reconciled_X.data;
1139 //Sx.data=reconciled_Sx.data;
1140 iterationcount++;
1141 return RunReconciliation(data, threadData, x, Sx, jacF, jacFt,eps,iterationcount,csvinputs,xdiag,sxdiag,logfile);
1142 }
1143 if(value < eps && iterationcount==1)
1144 {
1145 logfile << "J*/r" << "(" << value << ")" << " > " << eps << ", Convergence iteration not required \n\n";
1146 }
1147 else
1148 {
1149 logfile << "***** Value Converged, Convergence Completed******* \n\n";
1150 }
1151 logfile << "Final Results:\n";
1152 logfile << "=============\n";
1153 logfile << "Total Iteration to Converge : " << iterationcount << "\n";
1154 logfile << "Final Converged Value(J*/r) : " << value << "\n";
1155 logfile << "Epsilon : " << eps << "\n";
1156 printMatrixWithHeaders(reconciled_X.data,reconciled_X.rows,reconciled_X.column,csvinputs.headers,"reconciled_X ===> (x - (Sx*Ft*fstar))",logfile);
1157 printMatrixWithHeaders(reconciled_Sx.data,reconciled_Sx.rows,reconciled_Sx.column,csvinputs.headers,"reconciled_Sx ===> (Sx - (Sx*Ft*Fstar))",logfile);
1158
1159 /*
1160 * Calculate half width Confidence interval
1161 * W=lambda*sqrt(Sx)
1162 * where lamba = 1.96 and
1163 * Sx - diagonal elements of reconciled_Sx
1164 */
1165 double* reconSx_diag = (double*)calloc(reconciled_Sx.rows*1,sizeof(double));
1166 getDiagonalElements(reconciled_Sx.data,reconciled_Sx.rows,reconciled_Sx.column,reconSx_diag);
1167
1168 matrixData copyreconSx_diag = {reconciled_Sx.rows,1,reconSx_diag};
1169 matrixData tmpcopyreconSx_diag = copyMatrix(copyreconSx_diag);
1170 if(ACTIVE_STREAM(LOG_JAC)(useStream[LOG_JAC]))
1171 {
1172 logfile << "Calculations of HalfWidth Confidence Interval " << "\n";
1173 logfile << "===============================================\n";
1174 printMatrix(copyreconSx_diag.data,reconciled_Sx.rows,1,"reconciled-Sx_Diagonal",logfile);
1175 }
1176 calculateSquareRoot(copyreconSx_diag.data,reconciled_Sx.rows);
1177 if(ACTIVE_STREAM(LOG_JAC)(useStream[LOG_JAC]))
1178 {
1179 printMatrix(copyreconSx_diag.data,reconciled_Sx.rows,1,"reconciled-Sx_SquareRoot",logfile);
1180 logfile << "*****Completed***********\n";
1181 }
1182 scaleVector(reconciled_Sx.rows,1,1.96,copyreconSx_diag.data);
1183 printMatrixWithHeaders(copyreconSx_diag.data,reconciled_Sx.rows,1,csvinputs.headers,"Wx-HalfWidth-Interval-(1.96)*sqrt(Sx_diagonal)",logfile);
1184
1185
1186 /*
1187 * Calculate individual tests
1188 * (recon_x - x)/sqrt(Sx-recon_Sx)
1189 */
1190 double* newSx_diag = (double*)calloc(reconciled_Sx.rows*1,sizeof(double));
1191 solveMatrixSubtraction(sxdiag,tmpcopyreconSx_diag,newSx_diag,logfile);
1192 if(ACTIVE_STREAM(LOG_JAC)(useStream[LOG_JAC]))
1193 {
1194 logfile << "Calculations of Individual Tests " << "\n";
1195 logfile << "===============================================\n";
1196 printMatrix(newSx_diag,sxdiag.rows,sxdiag.column,"Sx-recon_Sx",logfile);
1197 }
1198 calculateSquareRoot(newSx_diag,reconciled_Sx.rows);
1199 if(ACTIVE_STREAM(LOG_JAC)(useStream[LOG_JAC]))
1200 {
1201 printMatrix(newSx_diag,sxdiag.rows,sxdiag.column,"squareroot-newSx",logfile);
1202 }
1203
1204 double *newX = (double*)calloc(xdiag.rows*1,sizeof(double));
1205 solveMatrixSubtraction(reconciled_X,xdiag,newX,logfile);
1206 // calculate absolute value for this numeric analysis
1207 for (unsigned int a=0; a < sizeof(newX); a++)
1208 {
1209 newX[a] = fabs(newX[a]);
1210 }
1211 if(ACTIVE_STREAM(LOG_JAC)(useStream[LOG_JAC]))
1212 {
1213 printMatrix(newX,xdiag.rows,xdiag.column,"recon_X - X",logfile);
1214 logfile << "*********Completed***********\n";
1215 }
1216
1217 for (int val=0; val < xdiag.rows; val++)
1218 {
1219 newX[val]=newX[val]/max(newSx_diag[val],sqrt(sxdiag.data[val]/10));
1220 }
1221
1222 printMatrixWithHeaders(newX,xdiag.rows,xdiag.column,csvinputs.headers,"IndividualTests_Value- (recon_x-x)/sqrt(Sx_diag)",logfile);
1223
1224 /*
1225 * create HTML Report
1226 */
1227 ofstream myfile;
1228 time_t now = time(0);
1229 std::stringstream htmlfile;
1230 if (omc_flag[FLAG_OUTPUT_PATH])
1231 {
1232 htmlfile << string(omc_flagValue[FLAG_OUTPUT_PATH]) << "/" << data->modelData->modelName << ".html";
1233 }
1234 else
1235 {
1236 htmlfile << data->modelData->modelName << ".html";
1237 }
1238 string html= htmlfile.str();
1239 myfile.open(html.c_str());
1240
1241 /* create a csv file */
1242 ofstream csvfile;
1243 std::stringstream csv_file;
1244 if (omc_flag[FLAG_OUTPUT_PATH])
1245 {
1246 csv_file << string(omc_flagValue[FLAG_OUTPUT_PATH]) << "/" << data->modelData->modelName << "_Outputs.csv";
1247 }
1248 else
1249 {
1250 csv_file << data->modelData->modelName <<"_Outputs.csv";
1251 }
1252 string tmpcsv= csv_file.str();
1253 csvfile.open(tmpcsv.c_str());
1254
1255 /* Add Overview Data */
1256 myfile << "<!DOCTYPE html><html>\n <head> <h1> DataReconciliation Report</h1></head> \n <body> \n ";
1257 myfile << "<h2> Overview: </h2>\n";
1258 myfile << "<table> \n";
1259 myfile << "<tr> \n" << "<th align=right> ModelFile: </th> \n" << "<td>" << data->modelData->modelFilePrefix << ".mo" << "</td> </tr>\n";
1260 myfile << "<tr> \n" << "<th align=right> ModelName: </th> \n" << "<td>" << data->modelData->modelName << "</td> </tr>\n";
1261 myfile << "<tr> \n" << "<th align=right> ModelDirectory: </th> \n" << "<td>" << data->modelData->modelDir << "</td> </tr>\n";
1262 myfile << "<tr> \n" << "<th align=right> Measurement Files: </th> \n" << "<td>" << omc_flagValue[FLAG_DATA_RECONCILE_Sx] << "</td> </tr>\n";
1263 myfile << "<tr> \n" << "<th align=right> Generated: </th> \n" << "<td>" << ctime(&now) << " by "<< "<b>" << CONFIG_VERSION"OMCompiler 14867775f.14867775f.14867775f-14867775f.14867775f+14867775f" << "</b>" << "</td> </tr>\n";
1264 myfile << "</table>\n";
1265
1266 /* Add Analysis data */
1267 myfile << "<h2> Analysis: </h2>\n";
1268 myfile << "<table> \n";
1269 myfile << "<tr> \n" << "<th align=right> Number of Extracted equations: </th> \n" << "<td>" << data->modelData->nSetcVars << "</td> </tr>\n";
1270 myfile << "<tr> \n" << "<th align=right> Number of Variables to be Reconciled: </th> \n" << "<td>" << csvinputs.headers.size() << "</td> </tr>\n";
1271 myfile << "<tr> \n" << "<th align=right> Number of Iteration to Converge: </th> \n" << "<td>" << iterationcount << "</td> </tr>\n";
1272 myfile << "<tr> \n" << "<th align=right> Final Converged Value(J*/r) : </th> \n" << "<td>" << value << "</td> </tr>\n";
1273 myfile << "<tr> \n" << "<th align=right> Epsilon : </th> \n" << "<td>" << eps << "</td> </tr>\n";
1274 myfile << "<tr> \n" << "<th align=right> Final Value of the objective Function (J*) : </th> \n" << "<td>" << (value*data->modelData->nSetcVars) << "</td> </tr>\n";
1275 //myfile << "<tr> \n" << "<th align=right> Chi-square value : </th> \n" << "<td>" << quantile(complement(chi_squared(data->modelData->nSetcVars), 0.05)) << "</td> </tr>\n";
1276 if(data->modelData->nSetcVars > 200){
1277 myfile << "<tr> \n" << "<th align=right> Chi-square value : </th> \n" << "<td>" << "NOT Available for equations > 200 in setC" << "</td> </tr>\n";
1278 }
1279 else
1280 {
1281 myfile << "<tr> \n" << "<th align=right> Chi-square value : </th> \n" << "<td>" << chisquaredvalue[data->modelData->nSetcVars-1] << "</td> </tr>\n";
1282 }
1283 myfile << "<tr> \n" << "<th align=right> Result of Global Test : </th> \n" << "<td>" << "TRUE" << "</td> </tr>\n";
1284 myfile << "</table>\n";
1285
1286 /* Add Results data */
1287 myfile << "<h2> Results: </h2>\n";
1288 myfile << "<table border=2>\n";
1289 myfile << "<tr>\n" << "<th> Variables to be Reconciled </th>\n" << "<th> Initial Measured Values </th>\n" << "<th> Reconciled Values </th>\n" << "<th> Initial Uncertainty Values </th>\n" <<"<th> Reconciled Uncertainty Values </th>\n";
1290 csvfile << "Variables to be Reconciled ," << "Initial Measured Values ," << "Reconciled Values ," << "Initial Uncertainty Values ," << "Reconciled Uncertainty Values,";
1291 myfile << "<th> Results of Local Tests </th>\n" << "<th> Values of Local Tests </th>\n" << "<th> Margin to Correctness(distance from 1.96) </th>\n" << "</tr>\n";
1292 csvfile << "Results of Local Tests ," << "Values of Local Tests ," << "Margin to Correctness(distance from 1.96) ," << "\n";
1293
1294 for (unsigned int r=0; r < csvinputs.headers.size(); r++)
1295 {
1296 myfile << "<tr>\n";
1297 myfile << "<td>" << csvinputs.headers[r] << "</td>\n";
1298 csvfile << csvinputs.headers[r] << ",";
1299 myfile << "<td>" << xdiag.data[r] << "</td>\n";
1300 csvfile << xdiag.data[r] << ",";
1301 myfile << "<td>" << reconciled_X.data[r] << "</td>\n";
1302 csvfile << reconciled_X.data[r] << ",";
1303
1304 myfile << "<td>" << csvinputs.sxdata[r] << "</td>\n";
1305 csvfile << csvinputs.sxdata[r] << ",";
1306
1307 myfile << "<td>" << copyreconSx_diag.data[r] << "</td>\n";
1308 csvfile << copyreconSx_diag.data[r] << ",";
1309
1310 if(newX[r] < 1.96)
1311 {
1312 myfile << "<td>" << "TRUE" << "</td>\n";
1313 csvfile << "TRUE" << ",";
1314 }
1315 else
1316 {
1317 myfile << "<td>" << "FALSE" << "</td>\n";
1318 csvfile << "FALSE" << ",";
1319 }
1320 myfile << "<td>" << newX[r] << "</td>\n";
1321 csvfile << newX[r] << ",";
1322
1323 myfile << "<td>" << (1.96-newX[r]) << "</td>\n";
1324 csvfile << (1.96-newX[r]) << ",\n";
1325 csvfile.flush();
1326 myfile << "</tr>\n";
1327 myfile.flush();
1328 }
1329 csvfile.close();
1330 myfile << "</table>\n";
1331 myfile << "</body>\n</html>";
1332 myfile.close();
1333
1334 free(tmpFstar.data);
1335 free(tmpfstar.data);
1336 //free(tmpmatrixC);
1337 //free(tmpmatrixD);
1338 //free(setc);
1339 free(reconciled_Sx.data);
1340 free(reconciled_X.data);
1341 free(copyreconSx_diag.data);
1342 free(tmpcopyreconSx_diag.data);
1343 free(newSx_diag);
1344 free(newX);
1345 //free(jacF.data);
1346 //free(jacFt.data);
1347 //free(x.data);
1348 //free(Sx.data);
1349 return 0;
1350}
1351
1352
1353int dataReconciliation(DATA* data, threadData_t *threadData)
1354{
1355 TRACE_PUSH
1356 const char* epselon = NULL__null;
1357 epselon = (char*)omc_flagValue[FLAG_DATA_RECONCILE_Eps];
1358
1359 /*
1360 * Create a Debug Log file
1361 */
1362 ofstream logfile;
1363 std::stringstream logfilename;
1364 if (omc_flag[FLAG_OUTPUT_PATH])
1365 {
1366 logfilename << omc_flagValue[FLAG_OUTPUT_PATH] << "/" << data->modelData->modelName << "_debug.log";
1367 }
1368 else
1369 {
1370 logfilename << data->modelData->modelName << "_debug.log";
1371 }
1372 string tmplogfilename= logfilename.str();
1373 logfile.open(tmplogfilename.c_str());
1374 logfile << "| info | " << "DataReconciliation Starting!\n";
1375 logfile << "| info | " << data->modelData->modelName << "\n";
1376
1377 if(epselon==NULL__null)
1378 {
1379 //errorStreamPrint(LOG_STDOUT, 0, "Epsilon Value not given, Please specify a convergence value (eg: -eps=0.0002), DataReconciliation cannot be computed!.");
1380 logfile << "| error | " << "Epsilon Value not given, Please specify a convergence value (eg: -eps=0.0002), DataReconciliation cannot be computed!.\n";
1381 logfile.close();
1382 exit(1);
1383 }
1384 csvData Sx_data = readCovarianceMatrixSx(data, threadData,logfile);
1385 //matrixData Sx = getCovarianceMatrixSx(Sx_data, data, threadData); // Prepare the data from csv file *
1386 matrixData Sx = computeCovarianceMatrixSx(Sx_data,data,threadData,logfile); // Compute the covariance matrix from csv inputs
1387 inputData x = getInputDataFromStartAttribute(Sx_data, data, threadData, logfile); // Read the inputs from the start attribute of the modelica model
1388 matrixData jacF = getJacobianMatrixF(data, threadData, logfile); // Compute the Jacobian Matrix F
1389 matrixData jacFt = getTransposeMatrix(jacF); // Compute the Transpose of jacobian Matrix F
1390
1391 double* Sx_diag = (double*)calloc(Sx.rows*1,sizeof(double));
1392 getDiagonalElements(Sx.data,Sx.rows,Sx.column,Sx_diag);
1393 matrixData tmpSx_diag={Sx.rows,1,Sx_diag};
1394
1395 matrixData tmp_x={x.rows,x.column,x.data};
1396 matrixData x_diag=copyMatrix(tmp_x);
1397
1398 // Print the initial information
1399 logfile << "\n\nInitial Data \n" << "=============\n";
1400 printMatrixWithHeaders(x.data,x.rows,x.column,Sx_data.headers,"X",logfile);
1401 printVectorMatrixWithHeaders(Sx_data.sxdata,Sx_data.rowcount,1,Sx_data.headers,"Half-WidthConfidenceInterval",logfile);
1402 printMatrixWithHeaders(Sx.data,Sx.rows,Sx.column,Sx_data.headers,"Sx",logfile);
1403 //printMatrix(Sx_diag,Sx.rows,1,"Sx-Diagonal elements",logfile);
1404
1405 //printMatrix(jacF.data,jacF.rows,jacF.column,"F",logfile);
1406 //printMatrix(jacFt.data,jacFt.rows,jacFt.column,"Ft",logfile);
1407 // Start the Algorithm
1408 RunReconciliation(data,threadData,x,Sx,jacF,jacFt,atof(epselon),1,Sx_data,x_diag,tmpSx_diag,logfile);
1409 logfile << "| info | " << "DataReconciliation Completed! \n";
1410 logfile.flush();
1411 logfile.close();
1412 free(Sx.data);
1413 free(x.data);
1414 free(jacF.data);
1415 free(jacFt.data);
1416 free(tmpSx_diag.data);
1417 free(x_diag.data);
1418 TRACE_POP
1419 return 0;
1420}
1421
1422}