[717] | 1 | #include "../bdm/math/square_mat.h" |
---|
| 2 | #include "../bdm/math/chmat.h" |
---|
| 3 | #include "base/user_info.h" |
---|
| 4 | #include "../square_mat_point.h" |
---|
| 5 | #include "UnitTest++.h" |
---|
| 6 | #include "TestReporterStdout.h" |
---|
| 7 | #include <iostream> |
---|
| 8 | #include <iomanip> |
---|
| 9 | #include <stdlib.h> |
---|
| 10 | #include <string.h> |
---|
| 11 | |
---|
| 12 | using std::cout; |
---|
| 13 | using std::cerr; |
---|
| 14 | using std::endl; |
---|
| 15 | |
---|
| 16 | using bdm::fsqmat; |
---|
| 17 | using bdm::chmat; |
---|
| 18 | using bdm::ldmat; |
---|
| 19 | using bdm::shared_ptr; |
---|
| 20 | using bdm::UIFile; |
---|
| 21 | using bdm::UI; |
---|
| 22 | |
---|
| 23 | const char *agenda_file_name = "agenda.cfg"; |
---|
| 24 | double epsilon = 0.00001; |
---|
| 25 | bool fast = false; |
---|
| 26 | |
---|
| 27 | namespace UnitTest { |
---|
| 28 | |
---|
| 29 | // can't include mat_checks.h because CheckClose is different in this file |
---|
| 30 | extern bool AreClose ( const itpp::vec &expected, const itpp::vec &actual, |
---|
| 31 | double tolerance ); |
---|
| 32 | |
---|
| 33 | extern bool AreClose ( const itpp::mat &expected, const itpp::mat &actual, |
---|
| 34 | double tolerance ); |
---|
| 35 | |
---|
| 36 | void CheckClose ( TestResults &results, const itpp::mat &expected, |
---|
| 37 | const itpp::mat &actual, double tolerance, |
---|
| 38 | TestDetails const& details ) { |
---|
[1064] | 39 | if ( !AreClose ( expected, actual, tolerance ) ) { |
---|
| 40 | MemoryOutStream stream; |
---|
| 41 | stream << "failed at " << expected.rows() |
---|
| 42 | << " x " << expected.cols(); |
---|
[717] | 43 | |
---|
[1064] | 44 | results.OnTestFailure ( details, stream.GetText() ); |
---|
| 45 | } |
---|
[717] | 46 | } |
---|
| 47 | |
---|
| 48 | } |
---|
| 49 | |
---|
| 50 | typedef void ( *FTestMatrix ) ( int, square_mat_point * ); |
---|
| 51 | |
---|
| 52 | template<typename TMatrix> |
---|
| 53 | void test_matrix ( int index, square_mat_point *point ) { |
---|
[1064] | 54 | Real_Timer tt; |
---|
[717] | 55 | |
---|
[1064] | 56 | cout << "agenda[" << index << "]:" << endl; |
---|
| 57 | mat A = point->get_matrix(); |
---|
| 58 | int sz = A.rows(); |
---|
| 59 | CHECK_EQUAL ( A.cols(), sz ); |
---|
[717] | 60 | |
---|
[1064] | 61 | tt.tic(); |
---|
| 62 | TMatrix sq_mat ( A ); |
---|
| 63 | double elapsed = tt.toc(); |
---|
| 64 | cout << "ctor(" << sz << " x " << sz << "): " << elapsed << " s" << endl; |
---|
[717] | 65 | |
---|
[1064] | 66 | tt.tic(); |
---|
| 67 | mat res = sq_mat.to_mat(); |
---|
| 68 | elapsed = tt.toc(); |
---|
[717] | 69 | |
---|
[1064] | 70 | if ( !fast ) { |
---|
| 71 | CHECK_CLOSE ( A, res, epsilon ); |
---|
| 72 | } |
---|
[717] | 73 | |
---|
[1064] | 74 | cout << "to_mat: " << elapsed << " s" << endl; |
---|
[717] | 75 | |
---|
[1064] | 76 | vec v = point->get_vector(); |
---|
| 77 | double w = point->get_scalar(); |
---|
| 78 | TMatrix sq_mat2 = sq_mat; |
---|
[717] | 79 | |
---|
[1064] | 80 | tt.tic(); |
---|
| 81 | sq_mat2.opupdt ( v, w ); |
---|
| 82 | elapsed = tt.toc(); |
---|
[717] | 83 | |
---|
[1064] | 84 | if ( !fast ) { |
---|
| 85 | mat expA = A + w * outer_product ( v, v ); |
---|
| 86 | CHECK_CLOSE ( expA, sq_mat2.to_mat(), epsilon ); |
---|
| 87 | } |
---|
[717] | 88 | |
---|
[1064] | 89 | cout << "opupdt: " << elapsed << " s" << endl; |
---|
[717] | 90 | |
---|
[1064] | 91 | TMatrix invmat ( sz ); |
---|
[717] | 92 | |
---|
[1064] | 93 | tt.tic(); |
---|
| 94 | sq_mat.inv ( invmat ); |
---|
| 95 | elapsed = tt.toc(); |
---|
[717] | 96 | |
---|
[1064] | 97 | mat invA; |
---|
| 98 | if ( !fast ) { |
---|
| 99 | invA = inv ( A ); |
---|
| 100 | CHECK_CLOSE ( invA, invmat.to_mat(), epsilon ); |
---|
| 101 | } |
---|
[717] | 102 | |
---|
[1064] | 103 | cout << "inv: " << elapsed << " s" << endl; |
---|
[717] | 104 | |
---|
[1064] | 105 | tt.tic(); |
---|
| 106 | double ld = sq_mat.logdet(); |
---|
| 107 | elapsed = tt.toc(); |
---|
[717] | 108 | |
---|
[1064] | 109 | if ( !fast ) { |
---|
| 110 | double d = det ( A ); |
---|
| 111 | CHECK_CLOSE ( log ( d ), ld, epsilon ); |
---|
| 112 | } |
---|
[717] | 113 | |
---|
[1064] | 114 | cout << "logdet: " << elapsed << " s" << endl; |
---|
[717] | 115 | |
---|
[1064] | 116 | tt.tic(); |
---|
| 117 | double q = sq_mat.qform ( ones ( sz ) ); |
---|
| 118 | elapsed = tt.toc(); |
---|
[717] | 119 | |
---|
[1064] | 120 | if ( !fast ) { |
---|
| 121 | CHECK_CLOSE ( sumsum ( A ), q, epsilon ); |
---|
| 122 | } |
---|
[717] | 123 | |
---|
[1064] | 124 | cout << "qform(1): " << elapsed << " s" << endl; |
---|
[717] | 125 | |
---|
[1064] | 126 | tt.tic(); |
---|
| 127 | q = sq_mat.qform ( v ); |
---|
| 128 | elapsed = tt.toc(); |
---|
[717] | 129 | |
---|
[1064] | 130 | if ( !fast ) { |
---|
| 131 | double r = ( A * v ) * v; |
---|
| 132 | CHECK_CLOSE ( r, q, epsilon ); |
---|
| 133 | } |
---|
[717] | 134 | |
---|
[1064] | 135 | cout << "qform(v): " << elapsed << " s" << endl; |
---|
[717] | 136 | |
---|
[1064] | 137 | tt.tic(); |
---|
| 138 | q = sq_mat.invqform ( v ); |
---|
| 139 | elapsed = tt.toc(); |
---|
[717] | 140 | |
---|
[1064] | 141 | if ( !fast ) { |
---|
| 142 | double r = ( invA * v ) * v; |
---|
| 143 | CHECK_CLOSE ( r, q, epsilon ); |
---|
| 144 | } |
---|
[717] | 145 | |
---|
[1064] | 146 | cout << "invqform: " << elapsed << " s" << endl; |
---|
[717] | 147 | |
---|
[1064] | 148 | TMatrix twice = sq_mat; |
---|
[717] | 149 | |
---|
[1064] | 150 | tt.tic(); |
---|
| 151 | twice += sq_mat; |
---|
| 152 | elapsed = tt.toc(); |
---|
[717] | 153 | |
---|
[1064] | 154 | if ( !fast ) { |
---|
| 155 | res = 2 * A; |
---|
| 156 | CHECK_CLOSE ( res, twice.to_mat(), epsilon ); |
---|
| 157 | } |
---|
[717] | 158 | |
---|
[1064] | 159 | cout << "+=: " << elapsed << " s" << endl; |
---|
[717] | 160 | |
---|
[1064] | 161 | sq_mat2 = sq_mat; |
---|
[717] | 162 | |
---|
[1064] | 163 | tt.tic(); |
---|
| 164 | sq_mat2.mult_sym ( A ); |
---|
| 165 | elapsed = tt.toc(); |
---|
[717] | 166 | |
---|
[1064] | 167 | if ( !fast ) { |
---|
| 168 | res = ( A * A ) * A.T(); |
---|
| 169 | CHECK_CLOSE ( res, sq_mat2.to_mat(), epsilon ); |
---|
| 170 | } |
---|
[717] | 171 | |
---|
[1064] | 172 | cout << "^2: " << elapsed << " s" << endl; |
---|
[717] | 173 | } |
---|
| 174 | |
---|
| 175 | void test_agenda ( FTestMatrix test ) { |
---|
[1064] | 176 | UIFile fag ( agenda_file_name ); |
---|
| 177 | Array<shared_ptr<square_mat_point> > mag; |
---|
| 178 | UI::get ( mag, fag, "agenda", UI::compulsory ); |
---|
| 179 | int sz = mag.size(); |
---|
| 180 | CHECK ( sz > 0 ); |
---|
| 181 | for ( int i = 0; i < sz; ++i ) { |
---|
| 182 | test ( i, mag ( i ).get() ); |
---|
| 183 | } |
---|
[717] | 184 | } |
---|
| 185 | |
---|
| 186 | SUITE ( ldmat ) { |
---|
[1064] | 187 | TEST ( agenda ) { |
---|
| 188 | test_agenda ( test_matrix<ldmat> ); |
---|
| 189 | } |
---|
[717] | 190 | } |
---|
| 191 | |
---|
| 192 | SUITE ( fsqmat ) { |
---|
[1064] | 193 | TEST ( agenda ) { |
---|
| 194 | test_agenda ( test_matrix<fsqmat> ); |
---|
| 195 | } |
---|
[717] | 196 | } |
---|
| 197 | |
---|
| 198 | SUITE ( chmat ) { |
---|
[1064] | 199 | TEST ( agenda ) { |
---|
| 200 | test_agenda ( test_matrix<chmat> ); |
---|
| 201 | } |
---|
[717] | 202 | } |
---|
| 203 | |
---|
| 204 | int main ( int argc, char const *argv[] ) { |
---|
[1064] | 205 | bool unknown = false; |
---|
| 206 | int update_next = 0; // 1 suite, 2 epsilon, 3 agenda file |
---|
| 207 | const char *suite = "ldmat"; |
---|
| 208 | const char **param = argv + 1; |
---|
| 209 | while ( *param && !unknown ) { |
---|
| 210 | if ( update_next ) { |
---|
| 211 | if ( update_next == 1 ) { |
---|
| 212 | suite = *param; |
---|
| 213 | } else if ( update_next == 2 ) { |
---|
| 214 | double eps = atof ( *param ); |
---|
| 215 | if ( eps > 0 ) { |
---|
| 216 | epsilon = eps; |
---|
| 217 | } else { |
---|
| 218 | cerr << "invalid epsilon value ignored" << endl; |
---|
| 219 | } |
---|
| 220 | } else { |
---|
| 221 | agenda_file_name = *param; |
---|
| 222 | } |
---|
[717] | 223 | |
---|
[1064] | 224 | update_next = 0; |
---|
| 225 | } else { |
---|
| 226 | if ( !strcmp ( *param, "-a" ) ) { |
---|
| 227 | update_next = 3; |
---|
| 228 | } else if ( !strcmp ( *param, "-c" ) ) { |
---|
| 229 | update_next = 1; |
---|
| 230 | } else if ( !strcmp ( *param, "-e" ) ) { |
---|
| 231 | update_next = 2; |
---|
| 232 | } else if ( !strcmp ( *param, "-f" ) ) { |
---|
| 233 | fast = true; |
---|
| 234 | } else { |
---|
| 235 | unknown = true; |
---|
| 236 | } |
---|
| 237 | } |
---|
[717] | 238 | |
---|
[1064] | 239 | ++param; |
---|
| 240 | } |
---|
[717] | 241 | |
---|
[1064] | 242 | if ( unknown || update_next ) { |
---|
| 243 | cerr << "usage: " << argv[0] << " [ -f ] [ -e epsilon ] [ -a agenda_input.cfg ] [ -c class ]" << endl; |
---|
| 244 | } else { |
---|
| 245 | UnitTest::TestReporterStdout reporter; |
---|
| 246 | UnitTest::TestRunner runner ( reporter ); |
---|
| 247 | return runner.RunTestsIf ( UnitTest::Test::GetTestList(), |
---|
| 248 | suite, |
---|
| 249 | UnitTest::True(), |
---|
| 250 | 0 ); |
---|
| 251 | } |
---|
[717] | 252 | } |
---|