mirror of
https://gitlab.com/libeigen/eigen.git
synced 2025-04-21 09:09:36 +08:00
Add missing test files for Polynomials module.
This commit is contained in:
Manuel Yguel
parent
9a4a08da46
commit
5ef83d6a6c
263
unsupported/test/polynomialsolver.cpp
Normal file
263
unsupported/test/polynomialsolver.cpp
Normal file
@ -0,0 +1,263 @@
|
||||
// This file is part of Eigen, a lightweight C++ template library
|
||||
// for linear algebra.
|
||||
//
|
||||
// Copyright (C) 2010 Manuel Yguel <manuel.yguel@gmail.com>
|
||||
//
|
||||
// Eigen is free software; you can redistribute it and/or
|
||||
// modify it under the terms of the GNU Lesser General Public
|
||||
// License as published by the Free Software Foundation; either
|
||||
// version 3 of the License, or (at your option) any later version.
|
||||
//
|
||||
// Alternatively, you can redistribute it and/or
|
||||
// modify it under the terms of the GNU General Public License as
|
||||
// published by the Free Software Foundation; either version 2 of
|
||||
// the License, or (at your option) any later version.
|
||||
//
|
||||
// Eigen is distributed in the hope that it will be useful, but WITHOUT ANY
|
||||
// WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
|
||||
// FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License or the
|
||||
// GNU General Public License for more details.
|
||||
//
|
||||
// You should have received a copy of the GNU Lesser General Public
|
||||
// License and a copy of the GNU General Public License along with
|
||||
// Eigen. If not, see <http://www.gnu.org/licenses/>.
|
||||
|
||||
#include "main.h"
|
||||
#include <unsupported/Eigen/Polynomials>
|
||||
#include <iostream>
|
||||
#include <algorithm>
|
||||
|
||||
#ifdef HAS_GSL
|
||||
#include "gsl_helper.h"
|
||||
#endif
|
||||
|
||||
using namespace std;
|
||||
|
||||
template<int Size>
|
||||
struct ei_increment_if_fixed_size
|
||||
{
|
||||
enum {
|
||||
ret = (Size == Dynamic) ? Dynamic : Size+1
|
||||
};
|
||||
};
|
||||
|
||||
|
||||
|
||||
|
||||
template<int Deg, typename POLYNOMIAL, typename SOLVER>
|
||||
bool aux_evalSolver( const POLYNOMIAL& pols, SOLVER& psolve )
|
||||
{
|
||||
typedef typename POLYNOMIAL::Scalar Scalar;
|
||||
|
||||
typedef typename SOLVER::RootsType RootsType;
|
||||
typedef Matrix<Scalar,Deg,1> EvalRootsType;
|
||||
|
||||
const int deg = pols.size()-1;
|
||||
|
||||
psolve.compute( pols );
|
||||
const RootsType& roots( psolve.roots() );
|
||||
EvalRootsType evr( deg );
|
||||
for( int i=0; i<roots.size(); ++i ){
|
||||
evr[i] = std::abs( poly_eval( pols, roots[i] ) ); }
|
||||
|
||||
bool evalToZero = evr.isZero( test_precision<Scalar>() );
|
||||
if( !evalToZero )
|
||||
{
|
||||
cerr << "WRONG root: " << endl;
|
||||
cerr << "Polynomial: " << pols.transpose() << endl;
|
||||
cerr << "Roots found: " << roots.transpose() << endl;
|
||||
cerr << "Abs value of the polynomial at the roots: " << evr.transpose() << endl;
|
||||
cerr << endl;
|
||||
}
|
||||
|
||||
#ifdef HAS_GSL
|
||||
if (ei_is_same_type< Scalar, double>::ret)
|
||||
{
|
||||
typedef GslTraits<Scalar> Gsl;
|
||||
RootsType gslRoots(deg);
|
||||
Gsl::eigen_poly_solve( pols, gslRoots );
|
||||
EvalRootsType gslEvr( deg );
|
||||
for( int i=0; i<gslRoots.size(); ++i )
|
||||
{
|
||||
gslEvr[i] = std::abs( poly_eval( pols, gslRoots[i] ) );
|
||||
}
|
||||
bool gslEvalToZero = gslEvr.isZero( test_precision<Scalar>() );
|
||||
if( !evalToZero )
|
||||
{
|
||||
if( !gslEvalToZero ){
|
||||
cerr << "GSL also failed" << endl; }
|
||||
else{
|
||||
cerr << "GSL did NOT failed" << endl; }
|
||||
cerr << "GSL roots found: " << gslRoots.transpose() << endl;
|
||||
cerr << "Abs value of the polynomial at the GSL roots: " << gslEvr.transpose() << endl;
|
||||
cerr << endl;
|
||||
}
|
||||
}
|
||||
#endif //< HAS_GSL
|
||||
|
||||
|
||||
std::vector<Scalar> rootModuli( roots.size() );
|
||||
Map< EvalRootsType > aux( &rootModuli[0], roots.size() );
|
||||
aux = roots.array().abs();
|
||||
std::sort( rootModuli.begin(), rootModuli.end() );
|
||||
bool distinctModuli=true;
|
||||
for( size_t i=1; i<rootModuli.size() && distinctModuli; ++i )
|
||||
{
|
||||
if( ei_isApprox( rootModuli[i], rootModuli[i-1] ) ){
|
||||
distinctModuli = false; }
|
||||
}
|
||||
VERIFY( evalToZero || !distinctModuli );
|
||||
|
||||
return distinctModuli;
|
||||
}
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
template<int Deg, typename POLYNOMIAL>
|
||||
void evalSolver( const POLYNOMIAL& pols )
|
||||
{
|
||||
typedef typename POLYNOMIAL::Scalar Scalar;
|
||||
|
||||
typedef PolynomialSolver<Scalar, Deg > PolynomialSolverType;
|
||||
|
||||
PolynomialSolverType psolve;
|
||||
aux_evalSolver<Deg, POLYNOMIAL, PolynomialSolverType>( pols, psolve );
|
||||
}
|
||||
|
||||
|
||||
|
||||
|
||||
template< int Deg, typename POLYNOMIAL, typename ROOTS, typename REAL_ROOTS >
|
||||
void evalSolverSugarFunction( const POLYNOMIAL& pols, const ROOTS& roots, const REAL_ROOTS& real_roots )
|
||||
{
|
||||
typedef typename POLYNOMIAL::Scalar Scalar;
|
||||
|
||||
typedef PolynomialSolver<Scalar, Deg > PolynomialSolverType;
|
||||
|
||||
PolynomialSolverType psolve;
|
||||
if( aux_evalSolver<Deg, POLYNOMIAL, PolynomialSolverType>( pols, psolve ) )
|
||||
{
|
||||
//It is supposed that
|
||||
// 1) the roots found are correct
|
||||
// 2) the roots have distinct moduli
|
||||
|
||||
typedef typename POLYNOMIAL::Scalar Scalar;
|
||||
typedef typename REAL_ROOTS::Scalar Real;
|
||||
|
||||
typedef PolynomialSolver<Scalar, Deg > PolynomialSolverType;
|
||||
typedef typename PolynomialSolverType::RootsType RootsType;
|
||||
typedef Matrix<Scalar,Deg,1> EvalRootsType;
|
||||
|
||||
//Test realRoots
|
||||
std::vector< Real > calc_realRoots;
|
||||
psolve.realRoots( calc_realRoots );
|
||||
VERIFY( calc_realRoots.size() == (size_t)real_roots.size() );
|
||||
|
||||
const Scalar psPrec = ei_sqrt( test_precision<Scalar>() );
|
||||
|
||||
for( size_t i=0; i<calc_realRoots.size(); ++i )
|
||||
{
|
||||
bool found = false;
|
||||
for( size_t j=0; j<calc_realRoots.size()&& !found; ++j )
|
||||
{
|
||||
if( ei_isApprox( calc_realRoots[i], real_roots[j] ), psPrec ){
|
||||
found = true; }
|
||||
}
|
||||
VERIFY( found );
|
||||
}
|
||||
|
||||
//Test greatestRoot
|
||||
VERIFY( ei_isApprox( roots.array().abs().maxCoeff(),
|
||||
ei_abs( psolve.greatestRoot() ), psPrec ) );
|
||||
|
||||
//Test smallestRoot
|
||||
VERIFY( ei_isApprox( roots.array().abs().minCoeff(),
|
||||
ei_abs( psolve.smallestRoot() ), psPrec ) );
|
||||
|
||||
bool hasRealRoot;
|
||||
//Test absGreatestRealRoot
|
||||
Real r = psolve.absGreatestRealRoot( hasRealRoot );
|
||||
VERIFY( hasRealRoot == (real_roots.size() > 0 ) );
|
||||
if( hasRealRoot ){
|
||||
VERIFY( ei_isApprox( real_roots.array().abs().maxCoeff(), ei_abs(r), psPrec ) ); }
|
||||
|
||||
//Test absSmallestRealRoot
|
||||
r = psolve.absSmallestRealRoot( hasRealRoot );
|
||||
VERIFY( hasRealRoot == (real_roots.size() > 0 ) );
|
||||
if( hasRealRoot ){
|
||||
VERIFY( ei_isApprox( real_roots.array().abs().minCoeff(), ei_abs( r ), psPrec ) ); }
|
||||
|
||||
//Test greatestRealRoot
|
||||
r = psolve.greatestRealRoot( hasRealRoot );
|
||||
VERIFY( hasRealRoot == (real_roots.size() > 0 ) );
|
||||
if( hasRealRoot ){
|
||||
VERIFY( ei_isApprox( real_roots.array().maxCoeff(), r, psPrec ) ); }
|
||||
|
||||
//Test smallestRealRoot
|
||||
r = psolve.smallestRealRoot( hasRealRoot );
|
||||
VERIFY( hasRealRoot == (real_roots.size() > 0 ) );
|
||||
if( hasRealRoot ){
|
||||
VERIFY( ei_isApprox( real_roots.array().minCoeff(), r, psPrec ) ); }
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
template<typename _Scalar, int _Deg>
|
||||
void polynomialsolver(int deg)
|
||||
{
|
||||
typedef ei_increment_if_fixed_size<_Deg> Dim;
|
||||
typedef Matrix<_Scalar,Dim::ret,1> PolynomialType;
|
||||
typedef Matrix<_Scalar,_Deg,1> EvalRootsType;
|
||||
|
||||
cout << "Standard cases" << endl;
|
||||
PolynomialType pols = PolynomialType::Random(deg+1);
|
||||
evalSolver<_Deg,PolynomialType>( pols );
|
||||
|
||||
cout << "Hard cases" << endl;
|
||||
_Scalar multipleRoot = ei_random<_Scalar>();
|
||||
EvalRootsType allRoots = EvalRootsType::Constant(deg,multipleRoot);
|
||||
roots_to_monicPolynomial( allRoots, pols );
|
||||
evalSolver<_Deg,PolynomialType>( pols );
|
||||
|
||||
cout << "Test sugar" << endl;
|
||||
EvalRootsType realRoots = EvalRootsType::Random(deg);
|
||||
roots_to_monicPolynomial( realRoots, pols );
|
||||
evalSolverSugarFunction<_Deg>(
|
||||
pols,
|
||||
realRoots.template cast <
|
||||
std::complex<
|
||||
typename NumTraits<_Scalar>::Real
|
||||
>
|
||||
>(),
|
||||
realRoots );
|
||||
}
|
||||
|
||||
|
||||
template<typename _Scalar> void polynomialsolver_scalar()
|
||||
{
|
||||
CALL_SUBTEST_1( (polynomialsolver<_Scalar,1>(1)) );
|
||||
CALL_SUBTEST_2( (polynomialsolver<_Scalar,2>(2)) );
|
||||
CALL_SUBTEST_3( (polynomialsolver<_Scalar,3>(3)) );
|
||||
CALL_SUBTEST_4( (polynomialsolver<_Scalar,4>(4)) );
|
||||
CALL_SUBTEST_5( (polynomialsolver<_Scalar,5>(5)) );
|
||||
CALL_SUBTEST_6( (polynomialsolver<_Scalar,6>(6)) );
|
||||
CALL_SUBTEST_7( (polynomialsolver<_Scalar,7>(7)) );
|
||||
CALL_SUBTEST_8( (polynomialsolver<_Scalar,8>(8)) );
|
||||
|
||||
CALL_SUBTEST_9( (polynomialsolver<_Scalar,Dynamic>(
|
||||
ei_random<int>(9,45)
|
||||
)) );
|
||||
}
|
||||
|
||||
void test_polynomialsolver()
|
||||
{
|
||||
for(int i = 0; i < g_repeat; i++)
|
||||
{
|
||||
polynomialsolver_scalar<double>();
|
||||
polynomialsolver_scalar<float>();
|
||||
}
|
||||
}
|
||||
124
unsupported/test/polynomialutils.cpp
Normal file
124
unsupported/test/polynomialutils.cpp
Normal file
@ -0,0 +1,124 @@
|
||||
// This file is part of Eigen, a lightweight C++ template library
|
||||
// for linear algebra.
|
||||
//
|
||||
// Copyright (C) 2010 Manuel Yguel <manuel.yguel@gmail.com>
|
||||
//
|
||||
// Eigen is free software; you can redistribute it and/or
|
||||
// modify it under the terms of the GNU Lesser General Public
|
||||
// License as published by the Free Software Foundation; either
|
||||
// version 3 of the License, or (at your option) any later version.
|
||||
//
|
||||
// Alternatively, you can redistribute it and/or
|
||||
// modify it under the terms of the GNU General Public License as
|
||||
// published by the Free Software Foundation; either version 2 of
|
||||
// the License, or (at your option) any later version.
|
||||
//
|
||||
// Eigen is distributed in the hope that it will be useful, but WITHOUT ANY
|
||||
// WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
|
||||
// FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License or the
|
||||
// GNU General Public License for more details.
|
||||
//
|
||||
// You should have received a copy of the GNU Lesser General Public
|
||||
// License and a copy of the GNU General Public License along with
|
||||
// Eigen. If not, see <http://www.gnu.org/licenses/>.
|
||||
|
||||
#include "main.h"
|
||||
#include <Eigen/Polynomials>
|
||||
#include <iostream>
|
||||
|
||||
using namespace std;
|
||||
|
||||
template<int Size>
|
||||
struct ei_increment_if_fixed_size
|
||||
{
|
||||
enum {
|
||||
ret = (Size == Dynamic) ? Dynamic : Size+1
|
||||
};
|
||||
};
|
||||
|
||||
template<typename _Scalar, int _Deg>
|
||||
void realRoots_to_monicPolynomial_test(int deg)
|
||||
{
|
||||
typedef ei_increment_if_fixed_size<_Deg> Dim;
|
||||
typedef Matrix<_Scalar,Dim::ret,1> PolynomialType;
|
||||
typedef Matrix<_Scalar,_Deg,1> EvalRootsType;
|
||||
|
||||
PolynomialType pols(deg+1);
|
||||
EvalRootsType roots = EvalRootsType::Random(deg);
|
||||
realRoots_to_monicPolynomial( roots, pols );
|
||||
|
||||
EvalRootsType evr( deg );
|
||||
for( int i=0; i<roots.size(); ++i ){
|
||||
evr[i] = std::abs( poly_eval( pols, roots[i] ) ); }
|
||||
|
||||
bool evalToZero = evr.isZero( test_precision<_Scalar>() );
|
||||
if( !evalToZero ){
|
||||
cerr << evr.transpose() << endl; }
|
||||
VERIFY( evalToZero );
|
||||
}
|
||||
|
||||
template<typename _Scalar> void realRoots_to_monicPolynomial_scalar()
|
||||
{
|
||||
CALL_SUBTEST_2( (realRoots_to_monicPolynomial_test<_Scalar,2>(2)) );
|
||||
CALL_SUBTEST_3( (realRoots_to_monicPolynomial_test<_Scalar,3>(3)) );
|
||||
CALL_SUBTEST_4( (realRoots_to_monicPolynomial_test<_Scalar,4>(4)) );
|
||||
CALL_SUBTEST_5( (realRoots_to_monicPolynomial_test<_Scalar,5>(5)) );
|
||||
CALL_SUBTEST_6( (realRoots_to_monicPolynomial_test<_Scalar,6>(6)) );
|
||||
CALL_SUBTEST_7( (realRoots_to_monicPolynomial_test<_Scalar,7>(7)) );
|
||||
CALL_SUBTEST_8( (realRoots_to_monicPolynomial_test<_Scalar,17>(17)) );
|
||||
|
||||
int deg = ei_random<int>(18,26);
|
||||
CALL_SUBTEST_9( (realRoots_to_monicPolynomial_test<_Scalar,Dynamic>(deg)) );
|
||||
}
|
||||
|
||||
|
||||
|
||||
|
||||
template<typename _Scalar, int _Deg>
|
||||
void CauchyBounds(int deg)
|
||||
{
|
||||
typedef ei_increment_if_fixed_size<_Deg> Dim;
|
||||
typedef Matrix<_Scalar,Dim::ret,1> PolynomialType;
|
||||
typedef Matrix<_Scalar,_Deg,1> EvalRootsType;
|
||||
|
||||
PolynomialType pols(deg+1);
|
||||
EvalRootsType roots = EvalRootsType::Random(deg);
|
||||
realRoots_to_monicPolynomial( roots, pols );
|
||||
_Scalar M = cauchy_max_bound( pols );
|
||||
_Scalar m = cauchy_min_bound( pols );
|
||||
_Scalar Max = roots.array().abs().maxCoeff();
|
||||
_Scalar min = roots.array().abs().minCoeff();
|
||||
bool eval = (M >= Max) && (m <= min);
|
||||
if( !eval )
|
||||
{
|
||||
cerr << "Roots: " << roots << endl;
|
||||
cerr << "Bounds: (" << m << ", " << M << ")" << endl;
|
||||
cerr << "Min,Max: (" << min << ", " << Max << ")" << endl;
|
||||
}
|
||||
VERIFY( eval );
|
||||
}
|
||||
|
||||
template<typename _Scalar> void CauchyBounds_scalar()
|
||||
{
|
||||
CALL_SUBTEST_2( (CauchyBounds<_Scalar,2>(2)) );
|
||||
CALL_SUBTEST_3( (CauchyBounds<_Scalar,3>(3)) );
|
||||
CALL_SUBTEST_4( (CauchyBounds<_Scalar,4>(4)) );
|
||||
CALL_SUBTEST_5( (CauchyBounds<_Scalar,5>(5)) );
|
||||
CALL_SUBTEST_6( (CauchyBounds<_Scalar,6>(6)) );
|
||||
CALL_SUBTEST_7( (CauchyBounds<_Scalar,7>(7)) );
|
||||
CALL_SUBTEST_8( (CauchyBounds<_Scalar,17>(17)) );
|
||||
|
||||
int deg = ei_random<int>(18,26);
|
||||
CALL_SUBTEST_9( (CauchyBounds<_Scalar,Dynamic>(deg)) );
|
||||
}
|
||||
|
||||
void test_polynomialutils()
|
||||
{
|
||||
for(int i = 0; i < g_repeat; i++)
|
||||
{
|
||||
realRoots_to_monicPolynomial_scalar<double>();
|
||||
realRoots_to_monicPolynomial_scalar<float>();
|
||||
CauchyBounds_scalar<double>();
|
||||
CauchyBounds_scalar<float>();
|
||||
}
|
||||
}
|
||||
Loading…
x
Reference in New Issue
Block a user