test/pca2.cc

Code
Comments
Other
Rev Date Author Line
4363 08 Sep 23 peter 1 // $Id$
4363 08 Sep 23 peter 2
4363 08 Sep 23 peter 3 /*
4363 08 Sep 23 peter 4   Copyright (C) 2023 Peter Johansson
4363 08 Sep 23 peter 5
4363 08 Sep 23 peter 6   This file is part of the yat library, https://dev.thep.lu.se/yat
4363 08 Sep 23 peter 7
4363 08 Sep 23 peter 8   The yat library is free software; you can redistribute it and/or
4363 08 Sep 23 peter 9   modify it under the terms of the GNU General Public License as
4363 08 Sep 23 peter 10   published by the Free Software Foundation; either version 3 of the
4363 08 Sep 23 peter 11   License, or (at your option) any later version.
4363 08 Sep 23 peter 12
4363 08 Sep 23 peter 13   The yat library is distributed in the hope that it will be useful,
4363 08 Sep 23 peter 14   but WITHOUT ANY WARRANTY; without even the implied warranty of
4363 08 Sep 23 peter 15   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
4363 08 Sep 23 peter 16   General Public License for more details.
4363 08 Sep 23 peter 17
4363 08 Sep 23 peter 18   You should have received a copy of the GNU General Public License
4363 08 Sep 23 peter 19   along with yat. If not, see <https://www.gnu.org/licenses/>.
4363 08 Sep 23 peter 20 */
4363 08 Sep 23 peter 21
4363 08 Sep 23 peter 22 #include <config.h>
4363 08 Sep 23 peter 23
4363 08 Sep 23 peter 24 #include "Suite.h"
4363 08 Sep 23 peter 25
4363 08 Sep 23 peter 26 #include "yat/utility/Matrix.h"
4363 08 Sep 23 peter 27 #include "yat/utility/PCA.h"
4363 08 Sep 23 peter 28 #include "yat/utility/Vector.h"
4369 13 Sep 23 peter 29 #include "yat/utility/VectorConstView.h"
4363 08 Sep 23 peter 30 #include "yat/utility/VectorView.h"
4363 08 Sep 23 peter 31
4363 08 Sep 23 peter 32 #include "yat/random/random.h"
4363 08 Sep 23 peter 33
4363 08 Sep 23 peter 34 #include <cmath>
4363 08 Sep 23 peter 35 #include <fstream>
4363 08 Sep 23 peter 36 #include <string>
4363 08 Sep 23 peter 37
4363 08 Sep 23 peter 38 using namespace theplu::yat;
4363 08 Sep 23 peter 39
4363 08 Sep 23 peter 40 void test1(test::Suite&);
4363 08 Sep 23 peter 41 void test2(test::Suite&);
4369 13 Sep 23 peter 42 void test3(test::Suite&);
4369 13 Sep 23 peter 43 void test4(test::Suite&);
4363 08 Sep 23 peter 44
4363 08 Sep 23 peter 45 int main(int argc, char* argv[])
4363 08 Sep 23 peter 46 {
4363 08 Sep 23 peter 47   test::Suite suite(argc, argv);
4369 13 Sep 23 peter 48   test1(suite);
4369 13 Sep 23 peter 49   test2(suite);
4369 13 Sep 23 peter 50   test3(suite);
4369 13 Sep 23 peter 51   test4(suite);
4363 08 Sep 23 peter 52   return suite.return_value();
4363 08 Sep 23 peter 53 }
4363 08 Sep 23 peter 54
4363 08 Sep 23 peter 55
4363 08 Sep 23 peter 56 utility::Vector noise(size_t n)
4363 08 Sep 23 peter 57 {
4363 08 Sep 23 peter 58   utility::Vector x(n);
4363 08 Sep 23 peter 59   random::Gaussian gaussian;
4363 08 Sep 23 peter 60   std::generate(x.begin(), x.end(), gaussian);
4363 08 Sep 23 peter 61   return x;
4363 08 Sep 23 peter 62 }
4363 08 Sep 23 peter 63
4363 08 Sep 23 peter 64
4363 08 Sep 23 peter 65 bool eye_test(test::Suite& suite, const utility::Matrix& X, size_t N=1)
4363 08 Sep 23 peter 66 {
4363 08 Sep 23 peter 67   bool ok=true;
4363 08 Sep 23 peter 68   for (size_t i=0; i<X.rows(); ++i)
4363 08 Sep 23 peter 69     for (size_t j=0; j<X.columns(); ++j) {
4363 08 Sep 23 peter 70       if (i==j && suite.equal(X(i,j), 1.0, N))
4363 08 Sep 23 peter 71         continue;
4363 08 Sep 23 peter 72       else if (i!=j && suite.equal_fix(X(i,j), 0.0, N * 1e-8))
4363 08 Sep 23 peter 73         continue;
4363 08 Sep 23 peter 74       ok = false;
4363 08 Sep 23 peter 75       suite.err() << "error: element " << i << " : " << j << "\n";
4363 08 Sep 23 peter 76     }
4363 08 Sep 23 peter 77
4363 08 Sep 23 peter 78   return suite.add(ok);
4363 08 Sep 23 peter 79 }
4363 08 Sep 23 peter 80
4363 08 Sep 23 peter 81
4369 13 Sep 23 peter 82 void pca_test(test::Suite& suite, size_t dimension, size_t n)
4363 08 Sep 23 peter 83 {
4369 13 Sep 23 peter 84   suite.out() << "=== pca test ===\n";
4363 08 Sep 23 peter 85   assert(dimension > 2);
4363 08 Sep 23 peter 86   suite.out() << "dimension: " << dimension << "\n";
4363 08 Sep 23 peter 87   utility::Vector offset(dimension);
4363 08 Sep 23 peter 88   offset(0) = 3.0;
4363 08 Sep 23 peter 89   offset(1) = 1.0;
4363 08 Sep 23 peter 90   offset(2) = 0.5;
4363 08 Sep 23 peter 91
4363 08 Sep 23 peter 92   utility::Vector pc1(dimension);
4363 08 Sep 23 peter 93   pc1(0) = 1.0;
4363 08 Sep 23 peter 94   pc1(1) = 2.0;
4363 08 Sep 23 peter 95   pc1 *= 1.0 / std::sqrt(pc1*pc1);
4369 13 Sep 23 peter 96   suite.out() << "PC1 : " << utility::VectorConstView(pc1, 0, 3) << "\n";
4369 13 Sep 23 peter 97
4363 08 Sep 23 peter 98   utility::Vector pc2(dimension);
4363 08 Sep 23 peter 99   pc2(0) = 1.0;
4363 08 Sep 23 peter 100   pc2(1) = -0.5;
4363 08 Sep 23 peter 101   pc2(2) = 0.4;
4363 08 Sep 23 peter 102   pc2 *= 1.0 / std::sqrt(pc2*pc2);
4369 13 Sep 23 peter 103   double ev1 = 2.0;
4363 08 Sep 23 peter 104   double ev2 = 1.0;
4369 13 Sep 23 peter 105   double epsilon = 1e-4;
4369 13 Sep 23 peter 106   suite.out() << "PC2 : " << utility::VectorConstView(pc2, 0, 3) << "\n";
4363 08 Sep 23 peter 107   assert(ev2 < ev1);
4363 08 Sep 23 peter 108
4363 08 Sep 23 peter 109   if (std::abs(pc1 * pc2) > 1e-8) {
4363 08 Sep 23 peter 110     throw std::runtime_error("pc1 and pc2 are not orthogonal");
4363 08 Sep 23 peter 111   }
4363 08 Sep 23 peter 112
4369 13 Sep 23 peter 113   // create data from n samples
4363 08 Sep 23 peter 114   random::Gaussian gaussian;
4369 13 Sep 23 peter 115   suite.out() << "# data points: " << n << "\n";
4363 08 Sep 23 peter 116   utility::Matrix X(dimension, n);
4363 08 Sep 23 peter 117   for (size_t i=0; i<n; ++i) {
4363 08 Sep 23 peter 118     auto x = X.column_view(i);
4369 13 Sep 23 peter 119     x = offset + gaussian() * pc1 * std::sqrt(ev1) +
4369 13 Sep 23 peter 120       gaussian() * pc2 * std::sqrt(ev2) +
4369 13 Sep 23 peter 121       std::sqrt(epsilon) * noise(dimension);
4363 08 Sep 23 peter 122   }
4363 08 Sep 23 peter 123   utility::PCA pca(X);
4363 08 Sep 23 peter 124   suite.out() << "# eigen values: " << pca.eigenvalues().size() << "\n";
4363 08 Sep 23 peter 125   if (pca.eigenvalues().size() != std::min(dimension, n)) {
4363 08 Sep 23 peter 126     suite.add(false);
4369 13 Sep 23 peter 127     suite.err() << "wong number of eigenvalues; expected "
4363 08 Sep 23 peter 128                 << std::min(dimension, n) << "\n";
4363 08 Sep 23 peter 129   }
4369 13 Sep 23 peter 130   suite.out() << "1st eigenvalue: " << pca.eigenvalues()(0) << "\n";
4369 13 Sep 23 peter 131   suite.add(suite.equal_fix(pca.eigenvalues()(0), ev1, 0.5));
4369 13 Sep 23 peter 132   suite.out() << "2nd eigenvalue: " << pca.eigenvalues()(1) << "\n";
4369 13 Sep 23 peter 133   suite.add(suite.equal_fix(pca.eigenvalues()(1), ev2, 0.5));
4363 08 Sep 23 peter 134
4363 08 Sep 23 peter 135   suite.out() << "eigenvectors: "
4363 08 Sep 23 peter 136               << pca.eigenvectors().rows() << " x "
4363 08 Sep 23 peter 137               << pca.eigenvectors().columns() << "\n";
4363 08 Sep 23 peter 138   if (pca.eigenvectors().rows() != pca.eigenvalues().size() ||
4363 08 Sep 23 peter 139       pca.eigenvectors().columns() != dimension) {
4363 08 Sep 23 peter 140     suite.add(false);
4363 08 Sep 23 peter 141     suite.err() << "wrong dimension of PCA::eigenvectors()\n";
4369 13 Sep 23 peter 142     suite.err() << "expected: " << pca.eigenvalues().size()
4369 13 Sep 23 peter 143                 << " x " << dimension << "\n";
4369 13 Sep 23 peter 144     return;
4363 08 Sep 23 peter 145   }
4363 08 Sep 23 peter 146
4363 08 Sep 23 peter 147   suite.out() << "top left corner of eigen vectors:\n";
4363 08 Sep 23 peter 148   for (size_t i=0; i<5; ++i) {
4363 08 Sep 23 peter 149     suite.out() << "pc" << (i+1);
4363 08 Sep 23 peter 150     for (size_t j=0; j<5; ++j)
4363 08 Sep 23 peter 151       suite.out() << " " << pca.eigenvectors()(i,j);
4363 08 Sep 23 peter 152     suite.out() << "\n";
4363 08 Sep 23 peter 153   }
4363 08 Sep 23 peter 154
4363 08 Sep 23 peter 155   if (!eye_test(suite,
4363 08 Sep 23 peter 156                 pca.eigenvectors() * transpose(pca.eigenvectors()),
4363 08 Sep 23 peter 157                 100))
4369 13 Sep 23 peter 158     suite.err() << "error: eigenvectors not orthonormal\n";
4363 08 Sep 23 peter 159
4363 08 Sep 23 peter 160
4363 08 Sep 23 peter 161   utility::Vector p1 = pca.eigenvectors() * pc1;
4363 08 Sep 23 peter 162   suite.out() << "pc1\n";
4363 08 Sep 23 peter 163   for (size_t i=0; i<5; ++i) {
4363 08 Sep 23 peter 164     suite.out() << i << " " << p1(i) << "\n";
4363 08 Sep 23 peter 165     if (i == 0 && std::abs(p1(i)) < 0.99) {
4363 08 Sep 23 peter 166       suite.err() << "error: p1(0)\n";
4363 08 Sep 23 peter 167       suite.add(false);
4363 08 Sep 23 peter 168     }
4369 13 Sep 23 peter 169     if (i != 0 && std::abs(p1(i)) > 0.02) {
4363 08 Sep 23 peter 170       suite.err() << "error: p1(" << i << ")\n";
4363 08 Sep 23 peter 171       suite.add(false);
4363 08 Sep 23 peter 172     }
4363 08 Sep 23 peter 173   }
4363 08 Sep 23 peter 174   // mirror if inferred pc is negated
4363 08 Sep 23 peter 175   if (p1(0) < 0)
4363 08 Sep 23 peter 176     pc1 *= -1.0;
4363 08 Sep 23 peter 177
4363 08 Sep 23 peter 178   utility::Vector p2 = pca.eigenvectors() * pc2;
4363 08 Sep 23 peter 179   suite.out() << "pc2\n";
4363 08 Sep 23 peter 180   for (size_t i=0; i<5; ++i) {
4363 08 Sep 23 peter 181     suite.out() << i << " " << p2(i) << "\n";
4363 08 Sep 23 peter 182     if (i == 1 && std::abs(p2(i)) < 0.99) {
4363 08 Sep 23 peter 183       suite.add(false);
4363 08 Sep 23 peter 184       suite.err() << "error: p2(0)\n";
4363 08 Sep 23 peter 185     }
4369 13 Sep 23 peter 186     if (i != 1 && std::abs(p2(i)) > 0.02) {
4363 08 Sep 23 peter 187       suite.err() << "error: p2(" << i << ")\n";
4363 08 Sep 23 peter 188       suite.add(false);
4363 08 Sep 23 peter 189     }
4363 08 Sep 23 peter 190   }
4363 08 Sep 23 peter 191   // mirror if inferred pc is negated
4363 08 Sep 23 peter 192   if (p2(1) < 0)
4363 08 Sep 23 peter 193     pc2 *= -1.0;
4363 08 Sep 23 peter 194
4363 08 Sep 23 peter 195
4363 08 Sep 23 peter 196   suite.out() << "test projection\n";
4363 08 Sep 23 peter 197   utility::Matrix X2(dimension, 5);
4363 08 Sep 23 peter 198   utility::Matrix alpha(5, 2);
4363 08 Sep 23 peter 199   alpha(1,0) = 1.0;
4363 08 Sep 23 peter 200   alpha(2,1) = 1.0;
4363 08 Sep 23 peter 201   alpha(3,0) = 1.0;
4363 08 Sep 23 peter 202   alpha(3,1) = 1.0;
4363 08 Sep 23 peter 203   alpha(4,0) = 10.0;
4363 08 Sep 23 peter 204   alpha(3,1) = -1.0;
4363 08 Sep 23 peter 205   for (size_t i=0; i<X2.columns(); ++i)
4363 08 Sep 23 peter 206     X2.column_view(i) = offset + alpha(i, 0) * pc1 + alpha(i, 1) * pc2;
4363 08 Sep 23 peter 207
4363 08 Sep 23 peter 208   utility::Matrix Y = pca.projection(X2);
4363 08 Sep 23 peter 209   suite.out() << "Y dim: " << Y.rows() << " x " << Y.columns() << "\n";
4363 08 Sep 23 peter 210
4363 08 Sep 23 peter 211   for (size_t i=0; i<Y.rows(); ++i) {
4363 08 Sep 23 peter 212     for (size_t j=0; j<Y.columns(); ++j) {
4363 08 Sep 23 peter 213
4363 08 Sep 23 peter 214       if (i > 1 || alpha(j,i)==0.0) {
4369 13 Sep 23 peter 215         if (!suite.add(suite.equal_fix(Y(i,j), 0.0, 0.3)))
4363 08 Sep 23 peter 216           suite.err() << "error: element: " << i << " " << j
4363 08 Sep 23 peter 217                       << " " << Y(i,j) << "\n";
4363 08 Sep 23 peter 218       }
4363 08 Sep 23 peter 219       else {
4363 08 Sep 23 peter 220         if (!suite.add(suite.equal_fix(Y(i,j), alpha(j,i), 0.2)))
4363 08 Sep 23 peter 221           suite.err() << "error: element: " << i << " " << j
4363 08 Sep 23 peter 222                       << " " << Y(i,j) << "\n";
4363 08 Sep 23 peter 223       }
4363 08 Sep 23 peter 224     }
4363 08 Sep 23 peter 225   }
4370 13 Sep 23 peter 226
4370 13 Sep 23 peter 227   utility::Vector z = pca.projection(X2.column_const_view(1));
4370 13 Sep 23 peter 228   if (!suite.equal_range_fix(z.begin(), z.end(), Y.begin_column(1), 1e-6)) {
4370 13 Sep 23 peter 229     suite.err() << "error: z not equal 2nd column of Y\n";
4370 13 Sep 23 peter 230     suite.out() << "z: " << z << "\n";
4370 13 Sep 23 peter 231     suite.add(false);
4370 13 Sep 23 peter 232   }
4363 08 Sep 23 peter 233 }
4363 08 Sep 23 peter 234
4363 08 Sep 23 peter 235
4363 08 Sep 23 peter 236 void test1(test::Suite& suite)
4363 08 Sep 23 peter 237 {
4369 13 Sep 23 peter 238   try {
4369 13 Sep 23 peter 239     pca_test(suite, 5000, 100);
4369 13 Sep 23 peter 240   }
4369 13 Sep 23 peter 241   catch (std::exception& e) {
4369 13 Sep 23 peter 242     suite.add(false);
4369 13 Sep 23 peter 243     suite.err() << "test 1 failed: " << e.what() << "\n";
4369 13 Sep 23 peter 244   }
4363 08 Sep 23 peter 245 }
4363 08 Sep 23 peter 246
4363 08 Sep 23 peter 247
4363 08 Sep 23 peter 248 void test2(test::Suite& suite)
4363 08 Sep 23 peter 249 {
4369 13 Sep 23 peter 250   try {
4369 13 Sep 23 peter 251     pca_test(suite, 6, 10000);
4369 13 Sep 23 peter 252   }
4369 13 Sep 23 peter 253   catch (std::exception& e) {
4369 13 Sep 23 peter 254     suite.add(false);
4369 13 Sep 23 peter 255     suite.err() << "test 2 failed: " << e.what() << "\n";
4369 13 Sep 23 peter 256   }
4363 08 Sep 23 peter 257 }
4369 13 Sep 23 peter 258
4369 13 Sep 23 peter 259
4369 13 Sep 23 peter 260 void test3(test::Suite& suite)
4369 13 Sep 23 peter 261 {
4369 13 Sep 23 peter 262   try {
4369 13 Sep 23 peter 263     suite.out() << "test3\n";
4369 13 Sep 23 peter 264     random::Gaussian gaussian;
4369 13 Sep 23 peter 265     // create one dimensional data, standard gaussian
4369 13 Sep 23 peter 266     size_t N = 100000;
4369 13 Sep 23 peter 267     utility::Matrix X(1, N);
4369 13 Sep 23 peter 268     std::generate(X.begin(), X.end(), gaussian);
4369 13 Sep 23 peter 269     utility::Matrix Cov = X * transpose(X);
4369 13 Sep 23 peter 270     Cov *= (1.0 / N);
4369 13 Sep 23 peter 271     suite.out() << "Cov: " << Cov << "\n";
4369 13 Sep 23 peter 272     utility::PCA pca(X);
4369 13 Sep 23 peter 273     suite.out() << "eigenvalue: " << pca.eigenvalues()(0) << "\n";
4369 13 Sep 23 peter 274     if (!suite.equal_fix(Cov(0,0), pca.eigenvalues()(0), 1e-5)) {
4369 13 Sep 23 peter 275       suite.add(false);
4369 13 Sep 23 peter 276       suite.err() << "Cov != PCA::eigenvalues()(0)\n";
4369 13 Sep 23 peter 277     }
4369 13 Sep 23 peter 278     if (!suite.equal_fix(1.0, pca.eigenvalues()(0), 0.01)) {
4369 13 Sep 23 peter 279       suite.add(false);
4369 13 Sep 23 peter 280       suite.err() << "PCA::eigenvalues()(0) != 1.0\n";
4369 13 Sep 23 peter 281     }
4369 13 Sep 23 peter 282   }
4369 13 Sep 23 peter 283   catch (std::exception& e) {
4369 13 Sep 23 peter 284     suite.add(false);
4369 13 Sep 23 peter 285     suite.err() << "test 3 failed: " << e.what() << "\n";
4369 13 Sep 23 peter 286   }
4369 13 Sep 23 peter 287 }
4369 13 Sep 23 peter 288
4369 13 Sep 23 peter 289
4369 13 Sep 23 peter 290 void test4(test::Suite& suite)
4369 13 Sep 23 peter 291 {
4369 13 Sep 23 peter 292   try {
4369 13 Sep 23 peter 293     suite.out() << "test4\n";
4369 13 Sep 23 peter 294     double s = 3.2;
4369 13 Sep 23 peter 295     random::Gaussian gaussian(s);
4369 13 Sep 23 peter 296     // create one dimensional data, standard gaussian
4369 13 Sep 23 peter 297     size_t N = 100;
4369 13 Sep 23 peter 298     utility::Matrix X(120, N);
4369 13 Sep 23 peter 299     // only generate first dimension/row
4369 13 Sep 23 peter 300     std::generate(X.begin_row(0), X.end_row(0), gaussian);
4369 13 Sep 23 peter 301
4369 13 Sep 23 peter 302     utility::PCA pca(X);
4369 13 Sep 23 peter 303     suite.out() << "1st eigenvalue: " << pca.eigenvalues()(0) << "\n";
4369 13 Sep 23 peter 304     if (!suite.equal_fix(s*s, pca.eigenvalues()(0), s*s*0.5)) {
4369 13 Sep 23 peter 305       suite.add(false);
4369 13 Sep 23 peter 306       suite.err() << "PCA::eigenvalues()(0) != 1.0\n";
4369 13 Sep 23 peter 307     }
4369 13 Sep 23 peter 308
4369 13 Sep 23 peter 309     utility::Matrix Cov = X * transpose(X);
4369 13 Sep 23 peter 310     Cov *= (1.0 / N);
4369 13 Sep 23 peter 311     suite.out() << "Cov(0,0): " << Cov(0,0) << "\n";
4369 13 Sep 23 peter 312     if (!suite.equal_fix(Cov(0,0), pca.eigenvalues()(0), 0.01)) {
4369 13 Sep 23 peter 313       suite.add(false);
4369 13 Sep 23 peter 314       suite.err() << "Cov != PCA::eigenvalues()(0)\n";
4369 13 Sep 23 peter 315     }
4369 13 Sep 23 peter 316   }
4369 13 Sep 23 peter 317   catch (std::exception& e) {
4369 13 Sep 23 peter 318     suite.add(false);
4369 13 Sep 23 peter 319     suite.err() << "test 3 failed: " << e.what() << "\n";
4369 13 Sep 23 peter 320   }
4369 13 Sep 23 peter 321 }

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