dacelib / dace Goto Github PK

View Code? Open in Web Editor NEW

31.0 9.0 14.0 4.49 MB

Differential Algebra Computational Toolbox

License: Apache License 2.0

C++ 43.22% Gnuplot 0.36% CMake 2.06% C 41.44% Shell 0.32% Rich Text Format 0.09% Python 0.02% Fortran 12.49%

numerical-methods differential-algebra automatic-differentiation dace uncertainty-propagation sensitivity-analysis

dace's People

Contributors

Stargazers

Watchers

Forkers

abgandar mmassari1974 zhjsun tylerjackoliver ccoming afossa thomasclb burcakyesilirmak raexyl a-ev petercalifano wongquinn sajjikazemi

dace's Issues

Error When Using Version 2.0

Hi, I'm trying to learn dace so I downloaded the precompiled windows package and installed it at D:/DACE . My OS version is win7 sp1 64 bit and the IDE is CodeLite V12.0.0. TDM-GCC has already been installed on my computer, which includes gcc 5.1.0. When I tried to compile and run the very basic test titled "1. DA Basics" at tutorial 2, I met the following errors:

C:\Windows\system32\cmd.exe /C D:/TDM-GCC-32/bin/mingw32-make.exe -j4 SHELL=cmd.exe -e -f "TEST_01.mk" MakeIntermediateDirs && D:/TDM-GCC-32/bin/mingw32-make.exe -j4 SHELL=cmd.exe -e -f "TEST_01.mk" all
----------Building project:[ TEST_01 - Debug ]----------
D:/TDM-GCC-32/bin/g++.exe -c "E:/CodeLite_Workspace/DACE_TEST/TEST_01/main.cpp" -std=c++11 -g -O0 -Wall -o ./Debug/main.cpp.o -ID:/DACE/include -I. -I.
D:/TDM-GCC-32/bin/g++.exe -o ./Debug/TEST_01 @"TEST_01.txt" -LD:/DACE/lib -L. -ldace
./Debug/main.cpp.o: In function main': E:/CodeLite_Workspace/DACE_TEST/TEST_01/main.cpp:9: undefined reference to DACE::DA::init(unsigned int, unsigned int)'
E:/CodeLite_Workspace/DACE_TEST/TEST_01/main.cpp:10: undefined reference to DACE::DA::DA(int, double)' E:/CodeLite_Workspace/DACE_TEST/TEST_01/main.cpp:11: undefined reference to DACE::operator<<(std::ostream&, DACE::DA const&)'
E:/CodeLite_Workspace/DACE_TEST/TEST_01/main.cpp:12: undefined reference to DACE::sin(DACE::DA const&)' E:/CodeLite_Workspace/DACE_TEST/TEST_01/main.cpp:12: undefined reference to DACE::operator<<(std::ostream&, DACE::DA const&)'
E:/CodeLite_Workspace/DACE_TEST/TEST_01/main.cpp:12: undefined reference to DACE::DA::~DA()' E:/CodeLite_Workspace/DACE_TEST/TEST_01/main.cpp:10: undefined reference to DACE::DA::~DA()'
E:/CodeLite_Workspace/DACE_TEST/TEST_01/main.cpp:12: undefined reference to DACE::DA::~DA()' E:/CodeLite_Workspace/DACE_TEST/TEST_01/main.cpp:10: undefined reference to DACE::DA::~DA()'
collect2.exe: error: ld returned 1 exit status
mingw32-make.exe: *** [Debug/TEST_01] Error 1
TEST_01.mk:78: recipe for target 'Debug/TEST_01' failed
====9 errors, 0 warnings====

It seems there are errors when linking to the compiled library and I don't know how to solve it. Please tell me what I should do to fix the problems. Thank you very much!

Tutorial1/Example7 crash

Example7 crashes at line 37.

libc++abi.dylib: terminating with uncaught exception of type DACE::DACEException: daceRoot: Zero-th root does not exist (ID: 645)
Abort trap: 6

I'm running on OSX. libc++abi.dylib is version 902.0.0 (according to otool).

It runs w/o the crash if I add 1e-9 to x[0] at line 28, but the coefficients are almost certainly wrong.

Evaluate sombrero function

z
     I  COEFFICIENT              ORDER EXPONENTS
     1    9.9999999999999989e-01   0   0  0
     2    5.9104644775390618e+01   2   2  0
     3    1.3743895347200000e+11   3   1  2
     4    1.4167642029942726e+10   3   3  0
     5    7.3786976294838206e+19   4   2  2
     6   -3.4437640265295577e+19   4   4  0
     7    7.9228162514264338e+28   5   3  2
     8   -3.9349336029542597e+28   5   5  0
     9    5.3169119831396635e+37   6   0  6
    10   -2.0416942015256308e+39   6   2  4
    11    1.5312706511442231e+39   6   4  2
    12   -2.1536402429816355e+38   6   6  0
    13   -3.1056909793281687e+48   7   3  4
    14    6.3940696633227003e+47   7   5  2
    15   -2.1264739034457924e+47   7   7  0
    16    6.1299821634635554e+54   8   0  8
    17    1.9615942923083377e+56   8   2  6
    18    3.9231885846166755e+56   8   4  4
    19    9.8079714615416887e+55   8   6  2
    20    3.3681793960999757e+56   8   8  0
    21    6.8452989584562171e+65   9   1  8
    22   -1.0109980000181490e+67   9   3  6
    23    1.3479973333575320e+67   9   5  4
    24   -3.6859302083995015e+66   9   7  2
    25    3.5999976382425530e+65   9   9  0
    26    3.5336941295567687e+73  10   0 10
    27   -4.5231284858326639e+75  10   2  8
    28    5.4277541829991967e+75  10   4  6
    29    1.0855508365998393e+76  10   6  4
    30    2.0919469246976070e+75  10   8  2
    31   -3.7425681696996680e+74  10  10  0

How to load a DA object from file?

I know that we can print a DA object into a file by using the output operator "<<".

I'm wondering is there any way that we can load the DA object from the output file, e.g. with the input operator ">>".

If this is possible, I think it will be very useful in many cases. For example, we can compute the DA object and print it into a file first, next time we will load it directly and evaluate the polynomials without computing the DA object again.

Different Results for the Same Expression

Dear Developers,
I met with a strange problem. When I tried to compute the acceleration of an artificial satellite due to the Earth's central attraction, the results varied with different computation orders.
Here are my codes:

#include
#include
#include <DA/dace.h>

using namespace std;
using namespace DACE;

int main(void)
{
DA::init(6, 6);

double MU = 398600.4418e9;

AlgebraicVector<DA> s0(6);
s0[0] = 4.180620257618603e6 + DA(1);
s0[1] = -1.511228187120802e6 + DA(2);
s0[2] = 5.121306826896370e6 + DA(3);
s0[3] = 0.714518846769259e3 + DA(4);
s0[4] = 7.460015410736627e3 + DA(5);
s0[5] = 1.620242405120838e3 + DA(6);

DA x = s0[0];
DA y = s0[1];
DA z = s0[2];

DA r = sqrt(x * x + y * y + z * z);

DA ax1 = -MU/(r*r*r) * x;
DA ax2 = -(MU/(r*r)) * (x/r);

DA ay1 = -MU/(r*r*r) * y;
DA ay2 = -(MU/(r*r)) * (y/r);

DA az1 = -MU/(r*r*r) * z;
DA az2 = -(MU/(r*r)) * (z/r);


cout << ax1 << endl;
cout << ay1 << endl;
cout << az1 << endl;

cout << ax2 << endl;
cout << ay2 << endl;
cout << az2 << endl;

return 0;

}

From a mathematical viewpoint, expressions for ax1 and ax2 are definitely same, however, for ax1, it will display "ALL COEFFICIENTS ZERO" while for ax2 all things looks normal. I guess there is significant precision loss when computing the term MU/(rrr), because the order of magnitude of this term is 10^(-6) while both the numerators and denominators for MU/(r*r) and (x/r) shares the same orders of magnitude. If there do exist precision loss, this behavior seems unreasonable because 10^(-6) is not that small for scientific computing.

Wrong Results Under Debug x64 Mode

There are different behaviors between debug x64 and release x64 modes for the precompiled windows package version 2.0. For example, suppose you want to evaluate sin(1.0) using DA method and the codes are as follows:

#include
#include "stdafx.h"
#include <DA/dace.h>

using namespace std;
using namespace DACE;

int main(void)
{
DA::init(24, 1);

DA x = DA(1);
DA f = sin(x);

double y = f.evalScalar(1.0);

cout << "DA result is: " << y << endl;
cout << "Exact value is: " << sin(1.0) << endl;

return 0;

}

Suppose the Visual Studio 2017 is used, after compiling and linking, the variable y will have the wrong value 0 under debug x64 mode while it will have the right value of about 0.841471 under the release x64 mode.

Build Mac OS X Framework

On Mac OS X, build a Framework instead of / in addition to a library. This is currently tricky because of limitations in cmake preventing easy generation of non-flat header file directories.

Floating point exception handling

Problem

Currently, the DACE does not check for floating point exceptions in any of its calculations, nor are non-normal values taken into account when DA calculations are performed. This has some strange consequences, for example:

Anything multiplied by a zero DA becomes zero, even if the other DA contains infinities (overflows)
Tanh becomes zero due to this for large arguments causing overflow during sinh/cosh (see #37)
most certainly other functions may suffer from similar effects

Open Questions

This is a very difficult problem for which we first have to answer several very basic question before a solution can be implemented.

What constitutes an error?

To really solve this, we have to decide what FP exceptions actually mean in DA context. Is an overflow in order 20 reason to make the entire DA invalid? Do we care about underflows?

How to signal errors

Should these throw visible exceptions/errors the user must handle right there (like sqrt(-1)) or should DAs have a silent error flag? Should there be several error flags indicating different errors (inf, nan, ...)? How is this exposed to the different DACE interfaces?

How to propagate errors

How does each and every built-in operation propagate errors / non-normal DA input? E.g. 1/inf = 0 is perfectly valid in FP, but inf/inf yields a nan.

How do we detect these problems in the first place?

Range checks for each function to make sure it can't overflow? How would this make sure the high order coefficients don't overflow either?
Checking each and every coefficient in a result for non-normal? Surely this will be incredibly slow...
Checking the FP environment flags after each DA operation for errors? Is this portable/reliable/fast? How does this interact with other libraries manipulating the FP environment (see here)

Estimation of order sorted norms returns `nan` when zero is expected

As per documentation, DA.estimNorm() should return zero for each order sorted norm if estimation is not possible (e.g. there are less than two non-constant coefficients which are not zero).

However, calling this method on a polynomial with a single non-zero coefficient of order >0 returns a vector of nan.

This is due to the fact that that linear regression in daceEstimate is still attempted, even if a warning is thrown to inform that there is not enough data to perform it.

Wrong output allocation in matrix/vector product

When computing the product between a (row) AlgebraicVector and an AlgebraicMatrix or an AlgebraicMatrix and a (column) AlgebraicVector the size allocated for the output vector is wrong, leading to unexpected results or errors if the matrix is not square.

The allocation occurs at lines 57 and 79 of AlgebraicMatrix_t.h.

line 57: obj1.size() should be obj2.ncols()
line 79: obj2.size() should be obj1.nrows()

Strange result for the eval() function

I'm trying to evaluate a weighted-order function. For example:

#include <dace/dace.h>
#include <cmath>
#include <iostream>
#include <iomanip>

using namespace std;
using namespace DACE;

int main( void )
{
    
    DA::init( 20, 2 );

    // the first case
    DA x1 = DA(1) * DA(1) * DA(1);
    DA y1 = DA(2);
    DA z1 = exp(x1) * exp(y1);

    vector<DA> p1({DA(1), 0});
    DA f1 = z1.eval(p1);
    
    // the second case
    DA x2 = DA(1);
    DA y2 = DA(2) * DA(2) * DA(2);
    DA z2 = exp(x2) * exp(y2);

    vector<DA> p2({0, DA(2)});
    DA f2 = z2.eval(p2);

    cout << "z1 function" << endl << z1 << endl;
    cout << "z2 function" << endl << z2 << endl;
    cout << "f1 function" << endl << f1 << endl;
    cout << "f2 function" << endl << f2 << endl;

    return 0;
}

The result of f1 and f2 are supposed the same in some sense. However, it seems that some higher-order terms are missed in f2. The results are as follows.

z1 function
     I  COEFFICIENT              ORDER EXPONENTS
     1    1.0000000000000000e+00   0   0  0
     2    1.0000000000000000e+00   1   0  1
     3    5.0000000000000000e-01   2   0  2
     4    1.6666666666666666e-01   3   0  3
     5    1.0000000000000000e+00   3   3  0
     6    4.1666666666666664e-02   4   0  4
     7    1.0000000000000000e+00   4   3  1
     8    8.3333333333333332e-03   5   0  5
     9    5.0000000000000000e-01   5   3  2
    10    1.3888888888888889e-03   6   0  6
    11    1.6666666666666666e-01   6   3  3
    12    5.0000000000000000e-01   6   6  0
    13    1.9841269841269841e-04   7   0  7
    14    4.1666666666666664e-02   7   3  4
    15    5.0000000000000000e-01   7   6  1
    16    2.4801587301587302e-05   8   0  8
    17    8.3333333333333332e-03   8   3  5
    18    2.5000000000000000e-01   8   6  2
    19    2.7557319223985893e-06   9   0  9
    20    1.3888888888888889e-03   9   3  6
    21    8.3333333333333329e-02   9   6  3
    22    1.6666666666666666e-01   9   9  0
    23    2.7557319223985894e-07  10   0 10
    24    1.9841269841269841e-04  10   3  7
    25    2.0833333333333332e-02  10   6  4
    26    1.6666666666666666e-01  10   9  1
    27    2.5052108385441720e-08  11   0 11
    28    2.4801587301587302e-05  11   3  8
    29    4.1666666666666666e-03  11   6  5
    30    8.3333333333333329e-02  11   9  2
    31    2.0876756987868100e-09  12   0 12
    32    2.7557319223985893e-06  12   3  9
    33    6.9444444444444447e-04  12   6  6
    34    2.7777777777777776e-02  12   9  3
    35    4.1666666666666664e-02  12  12  0
    36    1.6059043836821616e-10  13   0 13
    37    2.7557319223985894e-07  13   3 10
    38    9.9206349206349206e-05  13   6  7
    39    6.9444444444444441e-03  13   9  4
    40    4.1666666666666664e-02  13  12  1
    41    1.1470745597729726e-11  14   0 14
    42    2.5052108385441720e-08  14   3 11
    43    1.2400793650793651e-05  14   6  8
    44    1.3888888888888887e-03  14   9  5
    45    2.0833333333333332e-02  14  12  2
    46    7.6471637318198174e-13  15   0 15
    47    2.0876756987868100e-09  15   3 12
    48    1.3778659611992946e-06  15   6  9
    49    2.3148148148148149e-04  15   9  6
    50    6.9444444444444441e-03  15  12  3
    51    8.3333333333333332e-03  15  15  0
    52    4.7794773323873859e-14  16   0 16
    53    1.6059043836821616e-10  16   3 13
    54    1.3778659611992947e-07  16   6 10
    55    3.3068783068783064e-05  16   9  7
    56    1.7361111111111110e-03  16  12  4
    57    8.3333333333333332e-03  16  15  1
    58    2.8114572543455210e-15  17   0 17
    59    1.1470745597729726e-11  17   3 14
    60    1.2526054192720860e-08  17   6 11
    61    4.1335978835978830e-06  17   9  8
    62    3.4722222222222218e-04  17  12  5
    63    4.1666666666666666e-03  17  15  2
    64    1.5619206968586228e-16  18   0 18
    65    7.6471637318198174e-13  18   3 15
    66    1.0438378493934050e-09  18   6 12
    67    4.5928865373309819e-07  18   9  9
    68    5.7870370370370373e-05  18  12  6
    69    1.3888888888888887e-03  18  15  3
    70    1.3888888888888889e-03  18  18  0
    71    8.2206352466243310e-18  19   0 19
    72    4.7794773323873859e-14  19   3 16
    73    8.0295219184108080e-11  19   6 13
    74    4.5928865373309820e-08  19   9 10
    75    8.2671957671957661e-06  19  12  7
    76    3.4722222222222218e-04  19  15  4
    77    1.3888888888888889e-03  19  18  1
    78    4.1103176233121653e-19  20   0 20
    79    2.8114572543455210e-15  20   3 17
    80    5.7353727988648631e-12  20   6 14
    81    4.1753513975736200e-09  20   9 11
    82    1.0333994708994708e-06  20  12  8
    83    6.9444444444444444e-05  20  15  5
    84    6.9444444444444447e-04  20  18  2
------------------------------------------------

z2 function
     I  COEFFICIENT              ORDER EXPONENTS
     1    1.0000000000000000e+00   0   0  0
     2    1.0000000000000000e+00   1   1  0
     3    5.0000000000000000e-01   2   2  0
     4    1.0000000000000000e+00   3   0  3
     5    1.6666666666666666e-01   3   3  0
     6    1.0000000000000000e+00   4   1  3
     7    4.1666666666666664e-02   4   4  0
     8    5.0000000000000000e-01   5   2  3
     9    8.3333333333333332e-03   5   5  0
    10    5.0000000000000000e-01   6   0  6
    11    1.6666666666666666e-01   6   3  3
    12    1.3888888888888889e-03   6   6  0
    13    5.0000000000000000e-01   7   1  6
    14    4.1666666666666664e-02   7   4  3
    15    1.9841269841269841e-04   7   7  0
    16    2.5000000000000000e-01   8   2  6
    17    8.3333333333333332e-03   8   5  3
    18    2.4801587301587302e-05   8   8  0
    19    1.6666666666666666e-01   9   0  9
    20    8.3333333333333329e-02   9   3  6
    21    1.3888888888888889e-03   9   6  3
    22    2.7557319223985893e-06   9   9  0
    23    1.6666666666666666e-01  10   1  9
    24    2.0833333333333332e-02  10   4  6
    25    1.9841269841269841e-04  10   7  3
    26    2.7557319223985894e-07  10  10  0
    27    8.3333333333333329e-02  11   2  9
    28    4.1666666666666666e-03  11   5  6
    29    2.4801587301587302e-05  11   8  3
    30    2.5052108385441720e-08  11  11  0
    31    4.1666666666666664e-02  12   0 12
    32    2.7777777777777776e-02  12   3  9
    33    6.9444444444444447e-04  12   6  6
    34    2.7557319223985893e-06  12   9  3
    35    2.0876756987868100e-09  12  12  0
    36    4.1666666666666664e-02  13   1 12
    37    6.9444444444444441e-03  13   4  9
    38    9.9206349206349206e-05  13   7  6
    39    2.7557319223985894e-07  13  10  3
    40    1.6059043836821616e-10  13  13  0
    41    2.0833333333333332e-02  14   2 12
    42    1.3888888888888887e-03  14   5  9
    43    1.2400793650793651e-05  14   8  6
    44    2.5052108385441720e-08  14  11  3
    45    1.1470745597729726e-11  14  14  0
    46    8.3333333333333332e-03  15   0 15
    47    6.9444444444444441e-03  15   3 12
    48    2.3148148148148149e-04  15   6  9
    49    1.3778659611992946e-06  15   9  6
    50    2.0876756987868100e-09  15  12  3
    51    7.6471637318198174e-13  15  15  0
    52    8.3333333333333332e-03  16   1 15
    53    1.7361111111111110e-03  16   4 12
    54    3.3068783068783064e-05  16   7  9
    55    1.3778659611992947e-07  16  10  6
    56    1.6059043836821616e-10  16  13  3
    57    4.7794773323873859e-14  16  16  0
    58    4.1666666666666666e-03  17   2 15
    59    3.4722222222222218e-04  17   5 12
    60    4.1335978835978830e-06  17   8  9
    61    1.2526054192720860e-08  17  11  6
    62    1.1470745597729726e-11  17  14  3
    63    2.8114572543455210e-15  17  17  0
    64    1.3888888888888889e-03  18   0 18
    65    1.3888888888888887e-03  18   3 15
    66    5.7870370370370373e-05  18   6 12
    67    4.5928865373309819e-07  18   9  9
    68    1.0438378493934050e-09  18  12  6
    69    7.6471637318198174e-13  18  15  3
    70    1.5619206968586228e-16  18  18  0
    71    1.3888888888888889e-03  19   1 18
    72    3.4722222222222218e-04  19   4 15
    73    8.2671957671957661e-06  19   7 12
    74    4.5928865373309820e-08  19  10  9
    75    8.0295219184108080e-11  19  13  6
    76    4.7794773323873859e-14  19  16  3
    77    8.2206352466243310e-18  19  19  0
    78    6.9444444444444447e-04  20   2 18
    79    6.9444444444444444e-05  20   5 15
    80    1.0333994708994708e-06  20   8 12
    81    4.1753513975736200e-09  20  11  9
    82    5.7353727988648631e-12  20  14  6
    83    2.8114572543455210e-15  20  17  3
    84    4.1103176233121653e-19  20  20  0
------------------------------------------------

f1 function
     I  COEFFICIENT              ORDER EXPONENTS
     1    1.0000000000000000e+00   0   0  0
     2    1.0000000000000000e+00   3   3  0
     3    5.0000000000000000e-01   6   6  0
     4    1.6666666666666666e-01   9   9  0
     5    4.1666666666666664e-02  12  12  0
     6    8.3333333333333332e-03  15  15  0
     7    1.3888888888888889e-03  18  18  0
------------------------------------------------

f2 function
     I  COEFFICIENT              ORDER EXPONENTS
     1    1.0000000000000000e+00   0   0  0
     2    1.0000000000000000e+00   3   0  3
------------------------------------------------

Bug append operator

dace/interfaces/cxx/include/dace/AlgebraicVector_t.h

Line 266 in d12a941

const size_t size = this->size();

This line of code appears to be incorrect, it should be changed to:

const size_t size = obj.size();

Michele Maestrini

Problem running Tutorial1/Example7

Hello,

I have downloaded your DACE package yesterday and successfully compiled it on Ubuntu 22.04.
The tutorial examples compiled as well (except Example11DARK78.cpp which is ignored), but when running Tutorial1/Example7, I get the following message:

`Evaluate sombrero function

terminate called after throwing an instance of 'DACE::DACEException'
what(): daceRoot: Zero-th root does not exist (ID: 645)
Aborted (core dumped)
`
I am wondering if this could be linked to the warning message appearing when compiling the whole package:

[ 96%] Built target c_core_contrib .../dace/interfaces/cxx/DA.cpp: In static member function ‘static void DACE::DA::init(unsigned int, unsigned int)’: .../dace/interfaces/cxx/DA.cpp:69:27: warning: catching polymorphic type ‘class DACE::DACEException’ by value [-Wcatch-value=] 69 | } catch(DACEException ex) { | ^~ .../dace/interfaces/cxx/DA.cpp: In static member function ‘static void DACE::DA::init(unsigned int, unsigned int)’: .../dace/interfaces/cxx/DA.cpp:69:27: warning: catching polymorphic type ‘class DACE::DACEException’ by value [-Wcatch-value=] 69 | } catch(DACEException ex) { | ^~ [ 96%] Built target c_core
All other executables run w/o error messages.

Thanks a lot for any help,

Gilles

Error while reading `AlgebraicVector`s from text file

I am trying to write/read AlgebraicVector<DA> objects to/from text files, but the operator >> is failing to correctly load the vector from file, issuing the following error:

terminate called after throwing an instance of 'DACE::DACEException'
  what():  daceRead: Invalid data (ID: 632)

This seems to be caused by the blank line between the end of one vector component (i.e. --------------) and the beginning of the new one (i.e. COEFFICIENTS ...). (Manually) removing these blank lines fixes the issue.

I'll attach a minimum example in which the ouput file test_vector_io.txt is written using the << operator, while test_vector_io_no_blanks.txt is the same file but without blank lines.

test_vector_io.txt
test_vector_io_no_blanks.txt
testIO.txt

Add automatic unit test suite

We should have automatic unit tests to run in addition to the simple CI builds to ensure builds perform as expected.

Compilation Error

Dear Developers,
The current master-branch source code (until 09/14/2018) failed to pass compilation using VS2017 under Win10 X64. Most of the errors are C2491: 'identifier' : definition of dllimport function not allowed. Please help checking the codes, thanks!

Wrong hyperbolic tangent for argument >> 0

Computing the hyperbolic tangent of a DA object with constant part |c0|>>0 returns always zero. The correct behaviour would be +1 for c0 that goes to plus infinity and -1 for c0 that goes to minus infinity.

Examples:
DACE::DA x = 100 + DACE::DA(1)
x.tanh() --> 0 but should be 1

DACE::DA x = -100 + DACE::DA(1)
x.tanh() --> 0 but should be -1

Add Tutorials to the binary packages

Build documentation in CMake

Add an option (default to off) to automatically build doxygen documentation from CMake