test/vector_expression.cc

Code
Comments
Other
Rev Date Author Line
3605 27 Jan 17 peter 1 // $Id$
3605 27 Jan 17 peter 2
3605 27 Jan 17 peter 3 /*
3605 27 Jan 17 peter 4   Copyright (C) 2017 Peter Johansson
3605 27 Jan 17 peter 5
3605 27 Jan 17 peter 6   This file is part of the yat library, http://dev.thep.lu.se/yat
3605 27 Jan 17 peter 7
3605 27 Jan 17 peter 8   The yat library is free software; you can redistribute it and/or
3605 27 Jan 17 peter 9   modify it under the terms of the GNU General Public License as
3605 27 Jan 17 peter 10   published by the Free Software Foundation; either version 3 of the
3605 27 Jan 17 peter 11   License, or (at your option) any later version.
3605 27 Jan 17 peter 12
3605 27 Jan 17 peter 13   The yat library is distributed in the hope that it will be useful,
3605 27 Jan 17 peter 14   but WITHOUT ANY WARRANTY; without even the implied warranty of
3605 27 Jan 17 peter 15   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
3605 27 Jan 17 peter 16   General Public License for more details.
3605 27 Jan 17 peter 17
3605 27 Jan 17 peter 18   You should have received a copy of the GNU General Public License
3605 27 Jan 17 peter 19   along with yat. If not, see <http://www.gnu.org/licenses/>.
3605 27 Jan 17 peter 20 */
3605 27 Jan 17 peter 21
3605 27 Jan 17 peter 22 #include <config.h>
3605 27 Jan 17 peter 23
3605 27 Jan 17 peter 24 #include "Suite.h"
3605 27 Jan 17 peter 25
3605 27 Jan 17 peter 26 #include "yat/utility/Matrix.h"
3605 27 Jan 17 peter 27 #include "yat/utility/Vector.h"
3605 27 Jan 17 peter 28 #include "yat/utility/VectorView.h"
3605 27 Jan 17 peter 29
3607 27 Jan 17 peter 30 #include <cmath>
3607 27 Jan 17 peter 31
3605 27 Jan 17 peter 32 using namespace theplu::yat;
3605 27 Jan 17 peter 33 using namespace utility;
3605 27 Jan 17 peter 34
3605 27 Jan 17 peter 35 // create a vector that looks random very deterministically
3605 27 Jan 17 peter 36 Vector make_vector(size_t n, double seed)
3605 27 Jan 17 peter 37 {
3605 27 Jan 17 peter 38   Vector res(n);
3605 27 Jan 17 peter 39   double val = 0;
3605 27 Jan 17 peter 40   for (size_t i=0; i<n; ++i) {
3605 27 Jan 17 peter 41     val += sin(seed + i);
3605 27 Jan 17 peter 42     res(i) = val;
3605 27 Jan 17 peter 43   }
3605 27 Jan 17 peter 44   return res;
3605 27 Jan 17 peter 45 }
3605 27 Jan 17 peter 46
3605 27 Jan 17 peter 47
3605 27 Jan 17 peter 48 // create a matrix that looks random very deterministically
3605 27 Jan 17 peter 49 Matrix make_matrix(size_t row, size_t column, double seed)
3605 27 Jan 17 peter 50 {
3605 27 Jan 17 peter 51   Matrix res(row, column);
3605 27 Jan 17 peter 52   double val = 0;
3605 27 Jan 17 peter 53   for (size_t i=0; i<row; ++i)
3605 27 Jan 17 peter 54     for (size_t j=0; j<column; ++j) {
3605 27 Jan 17 peter 55       val += sin(seed);
3605 27 Jan 17 peter 56       seed += 1.0;
3605 27 Jan 17 peter 57       res(i, j) = val;
3605 27 Jan 17 peter 58     }
3605 27 Jan 17 peter 59   return res;
3605 27 Jan 17 peter 60 }
3605 27 Jan 17 peter 61
3605 27 Jan 17 peter 62
3605 27 Jan 17 peter 63 bool check(const VectorBase& lhs, const VectorBase& rhs, test::Suite& suite,
3605 27 Jan 17 peter 64            unsigned int N=1)
3605 27 Jan 17 peter 65 {
3605 27 Jan 17 peter 66   if (lhs.size() != rhs.size()) {
3605 27 Jan 17 peter 67     suite.err() << "error: wrong dimensions: comparing "
3605 27 Jan 17 peter 68                 << lhs.size() << " with " << rhs.size() << "\n";
3605 27 Jan 17 peter 69     suite.add(false);
3605 27 Jan 17 peter 70     return false;
3605 27 Jan 17 peter 71   }
3605 27 Jan 17 peter 72   if (!suite.equal_range(lhs.begin(), lhs.end(), rhs.begin(), N)) {
3605 27 Jan 17 peter 73     suite.add(false);
3605 27 Jan 17 peter 74     return false;
3605 27 Jan 17 peter 75   }
3605 27 Jan 17 peter 76   return true;
3605 27 Jan 17 peter 77 }
3605 27 Jan 17 peter 78
3605 27 Jan 17 peter 79
3605 27 Jan 17 peter 80 class Foo
3605 27 Jan 17 peter 81 {
3605 27 Jan 17 peter 82 public:
3605 27 Jan 17 peter 83   Foo(void) : raw_(3), view_(raw_) {}
3605 27 Jan 17 peter 84   VectorMutable& vec(void) { return vec_; }
3605 27 Jan 17 peter 85   VectorMutable& view(void) { return view_; }
3605 27 Jan 17 peter 86 private:
3605 27 Jan 17 peter 87   Vector vec_;
3605 27 Jan 17 peter 88   Vector raw_;
3605 27 Jan 17 peter 89   VectorView view_;
3605 27 Jan 17 peter 90 };
3605 27 Jan 17 peter 91
3605 27 Jan 17 peter 92
3677 14 Aug 17 peter 93 void vector_func(const VectorBase& vec)
3677 14 Aug 17 peter 94 {
3677 14 Aug 17 peter 95 }
3677 14 Aug 17 peter 96
3677 14 Aug 17 peter 97
3605 27 Jan 17 peter 98 int main(int argc, char* argv[])
3605 27 Jan 17 peter 99 {
3605 27 Jan 17 peter 100   test::Suite suite(argc, argv);
3605 27 Jan 17 peter 101   Vector x = make_vector(3, 1);
3605 27 Jan 17 peter 102   Vector y = make_vector(3, 2);
3605 27 Jan 17 peter 103   Vector z = make_vector(3, 3);
3605 27 Jan 17 peter 104
3605 27 Jan 17 peter 105   // Vector Vector
3605 27 Jan 17 peter 106
3605 27 Jan 17 peter 107   // addition operator+
3605 27 Jan 17 peter 108   {
3605 27 Jan 17 peter 109     suite.out() << "test: Vector + Vector\n";
3605 27 Jan 17 peter 110     Vector oldres(x);
3605 27 Jan 17 peter 111     oldres += y;
3605 27 Jan 17 peter 112     Vector res = x + y;
3605 27 Jan 17 peter 113     check(res, oldres, suite);
3605 27 Jan 17 peter 114     oldres += z;
3605 27 Jan 17 peter 115     res = x + y + z;
3605 27 Jan 17 peter 116     check(res, oldres, suite);
3677 14 Aug 17 peter 117 #ifdef YAT_TICKET897
3677 14 Aug 17 peter 118     vector_func(x + y);
3677 14 Aug 17 peter 119 #endif
3605 27 Jan 17 peter 120   }
3605 27 Jan 17 peter 121
3605 27 Jan 17 peter 122   // subtraction operator-
3605 27 Jan 17 peter 123   {
3605 27 Jan 17 peter 124     suite.out() << "test: Vector - Vector\n";
3605 27 Jan 17 peter 125     Vector oldres(x);
3605 27 Jan 17 peter 126     oldres -= y;
3605 27 Jan 17 peter 127     oldres -= z;
3605 27 Jan 17 peter 128     Vector res = x - y - z;
3605 27 Jan 17 peter 129     check(res, oldres, suite);
3677 14 Aug 17 peter 130 #ifdef YAT_TICKET897
3677 14 Aug 17 peter 131     vector_func(x - y);
3677 14 Aug 17 peter 132 #endif
3605 27 Jan 17 peter 133   }
3605 27 Jan 17 peter 134
3605 27 Jan 17 peter 135   // combining addition and subtraction
3605 27 Jan 17 peter 136   {
3605 27 Jan 17 peter 137     suite.out() << "test: Vector + Vector - Vector\n";
3605 27 Jan 17 peter 138     Vector oldres(x);
3605 27 Jan 17 peter 139     oldres += y;
3605 27 Jan 17 peter 140     oldres -= z;
3605 27 Jan 17 peter 141     Vector res = x + y - z;
3605 27 Jan 17 peter 142     check(res, oldres, suite);
3677 14 Aug 17 peter 143 #ifdef YAT_TICKET897
3677 14 Aug 17 peter 144     vector_func(x + y - z);
3677 14 Aug 17 peter 145 #endif
3605 27 Jan 17 peter 146   }
3605 27 Jan 17 peter 147
3605 27 Jan 17 peter 148   // Vector scalar operations
3605 27 Jan 17 peter 149   double c = 1.602;
3605 27 Jan 17 peter 150   double e = 2.718;
3605 27 Jan 17 peter 151
3605 27 Jan 17 peter 152   // multiply
3605 27 Jan 17 peter 153   {
3605 27 Jan 17 peter 154     suite.out() << "test: scalar * vector * scalar\n";
3605 27 Jan 17 peter 155     Vector oldres(x);
3605 27 Jan 17 peter 156     oldres *= e;
3605 27 Jan 17 peter 157     oldres *= c;
3605 27 Jan 17 peter 158     Vector res = c * x * e;
3605 27 Jan 17 peter 159     check(res, oldres, suite);
3677 14 Aug 17 peter 160 #ifdef YAT_TICKET897
3677 14 Aug 17 peter 161     vector_func(c * x);
3677 14 Aug 17 peter 162 #endif
3605 27 Jan 17 peter 163   }
3605 27 Jan 17 peter 164
3605 27 Jan 17 peter 165   // Matrix stuff
3605 27 Jan 17 peter 166
3605 27 Jan 17 peter 167   suite.out() << " === Matrix stuff ===\n";
3605 27 Jan 17 peter 168
3605 27 Jan 17 peter 169   Matrix A = make_matrix(3, 3, 1);
3605 27 Jan 17 peter 170   Matrix B = make_matrix(3, 3, 2);
3605 27 Jan 17 peter 171   Matrix C = make_matrix(4, 3, 3);
3605 27 Jan 17 peter 172   Matrix D = make_matrix(3, 4, 4);
3605 27 Jan 17 peter 173
3605 27 Jan 17 peter 174   // Matrix * Vector multiplication
3605 27 Jan 17 peter 175   {
3605 27 Jan 17 peter 176     suite.out() << "Matrix * Vector\n";
3605 27 Jan 17 peter 177     Vector oldres = A * (B * y);
3605 27 Jan 17 peter 178
3605 27 Jan 17 peter 179     suite.out() << "A * B * y\n";
3605 27 Jan 17 peter 180     Vector res1 = A * B * y;
3605 27 Jan 17 peter 181     check(res1, oldres, suite, 2);
3605 27 Jan 17 peter 182
3605 27 Jan 17 peter 183     suite.out() << "A * (B * y)\n";
3605 27 Jan 17 peter 184     Vector res2 = A * (B * y);
3605 27 Jan 17 peter 185     check(res2, oldres, suite, 2);
3605 27 Jan 17 peter 186
3605 27 Jan 17 peter 187     suite.out() << "(A * B) * y\n";
3605 27 Jan 17 peter 188     Vector res3 = (A * B) * y;
3605 27 Jan 17 peter 189     check(res3, oldres, suite, 2);
3677 14 Aug 17 peter 190
3677 14 Aug 17 peter 191 #ifdef YAT_TICKET897
3677 14 Aug 17 peter 192     vector_func(B*y);
3677 14 Aug 17 peter 193 #endif
3605 27 Jan 17 peter 194   }
3605 27 Jan 17 peter 195
3605 27 Jan 17 peter 196   // Vector * Matrix
3605 27 Jan 17 peter 197   {
3605 27 Jan 17 peter 198     suite.out() << "Vector * Matrix\n";
3605 27 Jan 17 peter 199     Vector res1 = A * x;
3605 27 Jan 17 peter 200     Matrix At(A);
3605 27 Jan 17 peter 201     At.transpose();
3605 27 Jan 17 peter 202     Vector res2 = x * At;
3605 27 Jan 17 peter 203     check(res1, res2, suite);
3677 14 Aug 17 peter 204 #ifdef YAT_TICKET897
3677 14 Aug 17 peter 205     vector_func(A * x);
3677 14 Aug 17 peter 206 #endif
3605 27 Jan 17 peter 207   }
3605 27 Jan 17 peter 208
3605 27 Jan 17 peter 209   {
3605 27 Jan 17 peter 210     suite.out() << "Vector * Matrix * Vector\n";
3605 27 Jan 17 peter 211     Vector vec = A * x;
3605 27 Jan 17 peter 212     double scalar = vec * y;
3605 27 Jan 17 peter 213     double res = y * A * x;
3605 27 Jan 17 peter 214     if (!suite.equal(scalar, res))
3605 27 Jan 17 peter 215       suite.add(false);
3605 27 Jan 17 peter 216   }
3605 27 Jan 17 peter 217
3605 27 Jan 17 peter 218   // Assignment into Vector
3605 27 Jan 17 peter 219   {
3605 27 Jan 17 peter 220     suite.out() << "assignment into vector\n";
3605 27 Jan 17 peter 221     Vector res = A * x;
3605 27 Jan 17 peter 222     Vector res2;
3605 27 Jan 17 peter 223     res2 = A * x;
3605 27 Jan 17 peter 224     check(res, res2, suite);
3605 27 Jan 17 peter 225   }
3605 27 Jan 17 peter 226
3605 27 Jan 17 peter 227   // Assignment into VectorView
3605 27 Jan 17 peter 228   {
3605 27 Jan 17 peter 229     suite.out() << "assignment into view\n";
3605 27 Jan 17 peter 230     Matrix M(A);
3605 27 Jan 17 peter 231     M.row_view(0) = B * x;
3605 27 Jan 17 peter 232     Vector res = B*x;
3605 27 Jan 17 peter 233     check(M.row_const_view(0), res, suite);
3605 27 Jan 17 peter 234   }
3605 27 Jan 17 peter 235
3605 27 Jan 17 peter 236   // Assignment into a VectorMutable
3605 27 Jan 17 peter 237   {
3605 27 Jan 17 peter 238     suite.out() << "assignment into abstract VectorMutable\n";
3605 27 Jan 17 peter 239     Vector res = 2*B*y;
3605 27 Jan 17 peter 240     Foo foo;
3605 27 Jan 17 peter 241     foo.vec() = 2*B*y;
3605 27 Jan 17 peter 242     check(foo.vec(), res, suite);
3605 27 Jan 17 peter 243     res = B*z;
3605 27 Jan 17 peter 244     foo.view() = B*z;
3605 27 Jan 17 peter 245     check(foo.view(), res, suite);
3605 27 Jan 17 peter 246   }
3605 27 Jan 17 peter 247
3605 27 Jan 17 peter 248   // combining everything
3605 27 Jan 17 peter 249   {
3605 27 Jan 17 peter 250     suite.out() << "=== Chombah dooobah ===\n";
3605 27 Jan 17 peter 251     // OK, we have matrices A, B, C, D and vectors x, y, z, scalars e and c
3605 27 Jan 17 peter 252     Vector res = A*x + 2*B*y + 3*D*C*z + e*(x+y+3*z) * (x*A*y) + A*B*x - B*x;
3605 27 Jan 17 peter 253
3605 27 Jan 17 peter 254     // let's calculate with old interface
3605 27 Jan 17 peter 255     Vector oldres = A*x;
3605 27 Jan 17 peter 256      Vector vec2 = B*y;
3605 27 Jan 17 peter 257     vec2 *= 2;
3605 27 Jan 17 peter 258     oldres += vec2;
3605 27 Jan 17 peter 259     Vector vec3 = D*(C*z);
3605 27 Jan 17 peter 260     vec3 *= 3;
3605 27 Jan 17 peter 261     oldres += vec3;
3605 27 Jan 17 peter 262     Vector vec4 = z;
3605 27 Jan 17 peter 263     vec4 *= 3;
3605 27 Jan 17 peter 264     vec4 += y;
3605 27 Jan 17 peter 265     vec4 += x;
3605 27 Jan 17 peter 266     vec4 *= e * (x*A*y);
3605 27 Jan 17 peter 267     oldres += vec4;
3605 27 Jan 17 peter 268     Vector vec5 = A*(B*x);
3605 27 Jan 17 peter 269     oldres += vec5;
3605 27 Jan 17 peter 270     Vector vec6 = B*x;
3605 27 Jan 17 peter 271     vec6 *= -1;
3605 27 Jan 17 peter 272     oldres += vec6;
3605 27 Jan 17 peter 273
3605 27 Jan 17 peter 274     check(res, oldres, suite, 2);
3605 27 Jan 17 peter 275   }
3605 27 Jan 17 peter 276
3608 27 Jan 17 peter 277   // operator+=
3608 27 Jan 17 peter 278   {
3608 27 Jan 17 peter 279     suite.out() << "test: Vector += VectorExpression\n";
3608 27 Jan 17 peter 280     Vector res1 = x + A*z;
3608 27 Jan 17 peter 281     Vector res2(x);
3608 27 Jan 17 peter 282     res2 += A*z;
3608 27 Jan 17 peter 283     check(res1, res2, suite);
3608 27 Jan 17 peter 284   }
3608 27 Jan 17 peter 285
3608 27 Jan 17 peter 286
3608 27 Jan 17 peter 287   // operator-=
3608 27 Jan 17 peter 288   {
3608 27 Jan 17 peter 289     suite.out() << "test: Vector -= VectorExpression\n";
3608 27 Jan 17 peter 290     Vector res1 = x - A*z;
3608 27 Jan 17 peter 291     Vector res2(x);
3608 27 Jan 17 peter 292     res2 -= A*z;
3608 27 Jan 17 peter 293     check(res1, res2, suite);
3608 27 Jan 17 peter 294   }
3608 27 Jan 17 peter 295
3610 27 Jan 17 peter 296   // negation
3610 27 Jan 17 peter 297   {
3610 27 Jan 17 peter 298     suite.out() << "test: negation\n";
3610 27 Jan 17 peter 299     Vector res = -x;
3610 27 Jan 17 peter 300     suite.out() << "res: " << res << "\n";
3610 27 Jan 17 peter 301     res = -res;
3610 27 Jan 17 peter 302     suite.out() << "res: " << res << "\n";
3610 27 Jan 17 peter 303     suite.out() << "x: " << x << "\n";
3610 27 Jan 17 peter 304     check(res, x, suite);
3610 27 Jan 17 peter 305   }
3608 27 Jan 17 peter 306
3610 27 Jan 17 peter 307
3605 27 Jan 17 peter 308   // testing assignment operator
3605 27 Jan 17 peter 309   return suite.return_value();
3605 27 Jan 17 peter 310 }

Generated on Sat Oct 14 02:21:45 2023 (UTC) by svndigest 0.11.2pre (dev build 1672)