1 // This file is part of Eigen, a lightweight C++ template library 2 // for linear algebra. Eigen itself is part of the KDE project. 3 // 4 // Copyright (C) 2008 Benoit Jacob <jacob.benoit.1 (at) gmail.com> 5 // 6 // This Source Code Form is subject to the terms of the Mozilla 7 // Public License v. 2.0. If a copy of the MPL was not distributed 8 // with this file, You can obtain one at http://mozilla.org/MPL/2.0/. 9 10 #include "main.h" 11 #include <Eigen/LeastSquares> 12 13 template<typename VectorType, 14 typename HyperplaneType> 15 void makeNoisyCohyperplanarPoints(int numPoints, 16 VectorType **points, 17 HyperplaneType *hyperplane, 18 typename VectorType::Scalar noiseAmplitude) 19 { 20 typedef typename VectorType::Scalar Scalar; 21 const int size = points[0]->size(); 22 // pick a random hyperplane, store the coefficients of its equation 23 hyperplane->coeffs().resize(size + 1); 24 for(int j = 0; j < size + 1; j++) 25 { 26 do { 27 hyperplane->coeffs().coeffRef(j) = ei_random<Scalar>(); 28 } while(ei_abs(hyperplane->coeffs().coeff(j)) < 0.5); 29 } 30 31 // now pick numPoints random points on this hyperplane 32 for(int i = 0; i < numPoints; i++) 33 { 34 VectorType& cur_point = *(points[i]); 35 do 36 { 37 cur_point = VectorType::Random(size)/*.normalized()*/; 38 // project cur_point onto the hyperplane 39 Scalar x = - (hyperplane->coeffs().start(size).cwise()*cur_point).sum(); 40 cur_point *= hyperplane->coeffs().coeff(size) / x; 41 } while( cur_point.norm() < 0.5 42 || cur_point.norm() > 2.0 ); 43 } 44 45 // add some noise to these points 46 for(int i = 0; i < numPoints; i++ ) 47 *(points[i]) += noiseAmplitude * VectorType::Random(size); 48 } 49 50 template<typename VectorType> 51 void check_linearRegression(int numPoints, 52 VectorType **points, 53 const VectorType& original, 54 typename VectorType::Scalar tolerance) 55 { 56 int size = points[0]->size(); 57 assert(size==2); 58 VectorType result(size); 59 linearRegression(numPoints, points, &result, 1); 60 typename VectorType::Scalar error = (result - original).norm() / original.norm(); 61 VERIFY(ei_abs(error) < ei_abs(tolerance)); 62 } 63 64 template<typename VectorType, 65 typename HyperplaneType> 66 void check_fitHyperplane(int numPoints, 67 VectorType **points, 68 const HyperplaneType& original, 69 typename VectorType::Scalar tolerance) 70 { 71 int size = points[0]->size(); 72 HyperplaneType result(size); 73 fitHyperplane(numPoints, points, &result); 74 result.coeffs() *= original.coeffs().coeff(size)/result.coeffs().coeff(size); 75 typename VectorType::Scalar error = (result.coeffs() - original.coeffs()).norm() / original.coeffs().norm(); 76 std::cout << ei_abs(error) << " xxx " << ei_abs(tolerance) << std::endl; 77 VERIFY(ei_abs(error) < ei_abs(tolerance)); 78 } 79 80 void test_eigen2_regression() 81 { 82 for(int i = 0; i < g_repeat; i++) 83 { 84 #ifdef EIGEN_TEST_PART_1 85 { 86 Vector2f points2f [1000]; 87 Vector2f *points2f_ptrs [1000]; 88 for(int i = 0; i < 1000; i++) points2f_ptrs[i] = &(points2f[i]); 89 Vector2f coeffs2f; 90 Hyperplane<float,2> coeffs3f; 91 makeNoisyCohyperplanarPoints(1000, points2f_ptrs, &coeffs3f, 0.01f); 92 coeffs2f[0] = -coeffs3f.coeffs()[0]/coeffs3f.coeffs()[1]; 93 coeffs2f[1] = -coeffs3f.coeffs()[2]/coeffs3f.coeffs()[1]; 94 CALL_SUBTEST(check_linearRegression(10, points2f_ptrs, coeffs2f, 0.05f)); 95 CALL_SUBTEST(check_linearRegression(100, points2f_ptrs, coeffs2f, 0.01f)); 96 CALL_SUBTEST(check_linearRegression(1000, points2f_ptrs, coeffs2f, 0.002f)); 97 } 98 #endif 99 #ifdef EIGEN_TEST_PART_2 100 { 101 Vector2f points2f [1000]; 102 Vector2f *points2f_ptrs [1000]; 103 for(int i = 0; i < 1000; i++) points2f_ptrs[i] = &(points2f[i]); 104 Hyperplane<float,2> coeffs3f; 105 makeNoisyCohyperplanarPoints(1000, points2f_ptrs, &coeffs3f, 0.01f); 106 CALL_SUBTEST(check_fitHyperplane(10, points2f_ptrs, coeffs3f, 0.05f)); 107 CALL_SUBTEST(check_fitHyperplane(100, points2f_ptrs, coeffs3f, 0.01f)); 108 CALL_SUBTEST(check_fitHyperplane(1000, points2f_ptrs, coeffs3f, 0.002f)); 109 } 110 #endif 111 #ifdef EIGEN_TEST_PART_3 112 { 113 Vector4d points4d [1000]; 114 Vector4d *points4d_ptrs [1000]; 115 for(int i = 0; i < 1000; i++) points4d_ptrs[i] = &(points4d[i]); 116 Hyperplane<double,4> coeffs5d; 117 makeNoisyCohyperplanarPoints(1000, points4d_ptrs, &coeffs5d, 0.01); 118 CALL_SUBTEST(check_fitHyperplane(10, points4d_ptrs, coeffs5d, 0.05)); 119 CALL_SUBTEST(check_fitHyperplane(100, points4d_ptrs, coeffs5d, 0.01)); 120 CALL_SUBTEST(check_fitHyperplane(1000, points4d_ptrs, coeffs5d, 0.002)); 121 } 122 #endif 123 #ifdef EIGEN_TEST_PART_4 124 { 125 VectorXcd *points11cd_ptrs[1000]; 126 for(int i = 0; i < 1000; i++) points11cd_ptrs[i] = new VectorXcd(11); 127 Hyperplane<std::complex<double>,Dynamic> *coeffs12cd = new Hyperplane<std::complex<double>,Dynamic>(11); 128 makeNoisyCohyperplanarPoints(1000, points11cd_ptrs, coeffs12cd, 0.01); 129 CALL_SUBTEST(check_fitHyperplane(100, points11cd_ptrs, *coeffs12cd, 0.025)); 130 CALL_SUBTEST(check_fitHyperplane(1000, points11cd_ptrs, *coeffs12cd, 0.006)); 131 delete coeffs12cd; 132 for(int i = 0; i < 1000; i++) delete points11cd_ptrs[i]; 133 } 134 #endif 135 } 136 } 137
