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