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@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
13template<typename VectorType,
14         typename HyperplaneType>
15void 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
50template<typename VectorType>
51void 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
64template<typename VectorType,
65         typename HyperplaneType>
66void 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
80void 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