DGtal  0.9.3beta
testMPolynomial.cpp
1 
30 #include <iostream>
32 #include <iomanip>
33 #include <sstream>
34 #include "DGtal/base/Common.h"
35 #include "DGtal/math/MPolynomial.h"
36 #include "DGtal/io/readers/MPolynomialReader.h"
38 
39 using namespace std;
40 using namespace DGtal;
41 
43 // Functions for testing class Signal.
45 
49 template <typename Ring>
51 durchblick()
52 {
54  = mmonomial<Ring>( 3, 1, 0 )
55  + mmonomial<Ring>( 1, 0, 3 )
56  + mmonomial<Ring>( 0, 3, 1 )
57  + mmonomial<Ring>( 0, 0, 3 )
58  + ((Ring)5) * mmonomial<Ring>( 0, 0, 1 );
59  return p;
60 }
61 
62 double
63 durchblickC( const double & x, const double & y, const double z )
64 {
65  return x*x*x*y + x*z*z*z + y*y*y*z + z*z*z + 5.0 * z;
66 }
67 
72 bool testMPolynomialSpeed( double step = 0.01 )
73 {
74  unsigned int nbok = 0;
75  unsigned int nb = 0;
76 
77  trace.beginBlock ( "Testing block ... Evaluation speed of mpolynomials (naive)" );
78  trace.info() << setprecision( 15 ) << "step is " << step << std::endl;
79  trace.info() << "approximately " << 8.0/(step*step*step) << " computations." << std::endl;
80  MPolynomial<3, double> P = durchblick<double>();
81  double total = 0.0;
82  for ( double x = -1.0; x < 1.0; x += step )
83  {
84  for ( double y = -1.0; y < 1.0; y += step )
85  {
86  for ( double z = -1.0; z < 1.0; z += step )
87  total += P(x)(y)(z);
88  }
89  }
90  trace.info() << "Total = " << total << std::endl;
91  trace.endBlock();
92 
93  trace.beginBlock ( "Testing block ... Evaluation speed of mpolynomials" );
94  //MPolynomial<3, double> P = durchblick<double>();
95  double total1 = 0.0;
96  for ( double x = -1.0; x < 1.0; x += step )
97  {
98  MPolynomial<2, double> PX = P( x );
99  for ( double y = -1.0; y < 1.0; y += step )
100  {
101  MPolynomial<1, double> PXY = PX( y );
102  for ( double z = -1.0; z < 1.0; z += step )
103  total1 += PXY( z );
104  }
105  }
106  trace.info() << "Total1 = " << total1 << std::endl;
107  trace.endBlock();
108 
109  trace.beginBlock ( "Testing block ... Same computation in C." );
110  double total2 = 0.0;
111  for ( double x = -1.0; x < 1.0; x += step )
112  {
113  for ( double y = -1.0; y < 1.0; y += step )
114  {
115  for ( double z = -1.0; z < 1.0; z += step )
116  total2 += durchblickC( x, y, z );
117  }
118  }
119  trace.info() << "Total2 = " << total2 << std::endl;
120  trace.endBlock();
121  nbok += fabs( total1 - total ) < 1e-8 ? 1 : 0;
122  nb++;
123  trace.info() << "(" << nbok << "/" << nb << ") "
124  << "fabs( total1 - total ) < 1e-8" << std::endl;
125  nbok += fabs( total2 - total ) < 1e-8 ? 1 : 0;
126  nb++;
127  trace.info() << "(" << nbok << "/" << nb << ") "
128  << "fabs( total2 - total ) < 1e-8" << std::endl;
129 
130  trace.info() << "For information, ImaGene::Polynomial3 takes 164ms for step=0.01 and 1604ms for step = 0.005." << std::endl;
131  return nbok == nb;
132 }
133 
138 bool testMPolynomial()
139 {
140  unsigned int nbok = 0;
141  unsigned int nb = 0;
142 
143  trace.beginBlock ( "Testing block ..." );
144 
145  MPolynomial<2, int> f = mmonomial<int>(1, 2) + 3 * mmonomial<int>(4, 5);
146  trace.info() << f << std::endl;
147  int e = f(4)(2);
148  trace.info() << e << std::endl;
149  nbok += e == 24592 ? 1 : 0;
150  nb++;
151  trace.info() << derivative<1>(f) << std::endl;
152  MPolynomial<1, int> g = f(2), h = f(3);
153  trace.info() << g << " and " << h << std::endl;
154  trace.info() << gcd<double>(g, h) << std::endl; // cast polynomials to MPolynomial<1, double> first; otherwise,
155  // result will be incorrect since int is not a field
156  trace.info() << gcd(g, h) << std::endl; // to prove our point, check for yourself that the result is
157  // X_0^5 instead of X_0^2
158  nbok += gcd<double>(g,h) == mmonomial<double>(2) ? 1 : 0;
159  nb++;
160  MPolynomial<1, double> l = f(mmonomial<double>(1) - 1)(mmonomial<double>(1) - 2);
161  trace.info() << l << std::endl;
162  trace.info() << derivative<0>(l) << std::endl;
163  trace.info() << gcd(l, derivative<0>(l)) << " (gcd of two previous polys is 1)" << std::endl;
164  nbok += gcd(l, derivative<0>(l)) == 1.0 * mmonomial<double>(0) ? 1 : 0;
165  nb++;
166  trace.info() << "Durchblick (x3y+xz3+y3z+z3+5z)= " << durchblick<double>() << std::endl;
167  nbok += true ? 1 : 0;
168  nb++;
169  trace.info() << "(" << nbok << "/" << nb << ") "
170  << "true == true" << std::endl;
171  trace.endBlock();
172  MPolynomial<3, double> Q = mmonomial<double>( 0, 0, 0 )
173  + mmonomial<double>( 1, 2, 0 ) + mmonomial<double>( 4, 1, 1 );
174  std::cout << "Q(x,y,z)=1+xy^2+x^4yz = " << Q << std::endl;
175  std::cout << " degree = " << Q.degree() << std::endl;
176  std::cout << " leading = " << Q.leading() << std::endl;
177  std::cout << " Q[0] = " << Q[ 0 ] << std::endl;
178  std::cout << " Q[1] = " << Q[ 1 ] << std::endl;
179  std::cout << " Q[2] = " << Q[ 2 ] << std::endl;
180  std::cout << " Q[3] = " << Q[ 3 ] << std::endl;
181  std::cout << " Q[4] = " << Q[ 4 ] << std::endl;
182  std::cout << " dQ/dx = " << derivative<0>(Q) << std::endl;
183  std::cout << " dQ/dy = " << derivative<1>(Q) << std::endl;
184  std::cout << " dQ/dz = " << derivative<2>(Q) << std::endl;
186  P = Xe_k<3,double>( 2, 7 ) + Xe_k<3,double>( 1, 3 );
187  trace.info() << "P=" << P << std::endl;
188  return nbok == nb;
189 }
190 
191 
193 // Standard services - public :
194 
195 int main( int /*argc*/, char** /*argv*/ )
196 {
197  trace.beginBlock ( "Testing class MPolynomial" );
198 
199  bool res = testMPolynomial()
200  && testMPolynomialSpeed( 0.05 );
201  trace.emphase() << ( res ? "Passed." : "Error." ) << endl;
202  trace.endBlock();
203  return res ? 0 : 1;
204 }
205 // //
void beginBlock(const std::string &keyword="")
Trace trace
Definition: Common.h:137
const MPolyNM1 & leading() const
Definition: MPolynomial.h:1119
STL namespace.
double endBlock()
Aim: Represents a multivariate polynomial, i.e. an element of , where K is some ring or field...
Definition: MPolynomial.h:58
std::ostream & emphase()
MPolynomial< 1, Ring, Alloc > gcd(const MPolynomial< 1, Ring, Alloc > &f, const MPolynomial< 1, Ring, Alloc > &g)
Definition: MPolynomial.h:2047
DGtal is the top-level namespace which contains all DGtal functions and types.
std::ostream & info()
int degree() const
Definition: MPolynomial.h:1109