vcs.maemo.org Git - opencv/blob - tests/cxcore/src/asolvepoly.cpp

   1 // 2008-05-14, Xavier Delacour <xavier.delacour@gmail.com>
   2
   3 #include "cxcoretest.h"
   4
   5 #include <algorithm>
   6 #include <complex>
   7 #include <vector>
   8 #include <iostream>
   9
  10 #if 0
  11
  12 typedef std::complex<double> complex_type;
  13
  14 struct pred_complex {
  15     bool operator() (const complex_type& lhs, const complex_type& rhs) const
  16     {
  17         return lhs.real() != rhs.real() ? lhs.real() < rhs.real() : lhs.imag() < rhs.imag();
  18     }
  19 };
  20
  21 class CV_SolvePolyTest : public CvTest
  22 {
  23 public:
  24     CV_SolvePolyTest();
  25     ~CV_SolvePolyTest();
  26 protected:
  27     virtual void run( int start_from );
  28 };
  29
  30 CV_SolvePolyTest::CV_SolvePolyTest() : CvTest( "solve-poly", "cvSolvePoly" ) {}
  31
  32 CV_SolvePolyTest::~CV_SolvePolyTest() {}
  33
  34 void CV_SolvePolyTest::run( int )
  35 {
  36     CvRNG rng = cvRNG();
  37     int fig = 100;
  38     double range = 50;
  39
  40     for (int idx = 0, max_idx = 1000, progress = 0; idx < max_idx; ++idx)
  41     {
  42         int n = cvRandInt(&rng) % 13 + 1;
  43         std::vector<complex_type> r(n), ar(n), c(n + 1, 0);
  44         std::vector<double> a(n + 1), u(n * 2);
  45
  46         int rr_odds = 3; // odds that we get a real root
  47         for (int j = 0; j < n;)
  48         {
  49             if (cvRandInt(&rng) % rr_odds == 0 || j == n - 1)
  50                     r[j++] = cvRandReal(&rng) * range;
  51             else
  52             {
  53                     r[j] = complex_type(cvRandReal(&rng) * range,
  54                             cvRandReal(&rng) * range + 1);
  55                     r[j + 1] = std::conj(r[j]);
  56                     j += 2;
  57             }
  58         }
  59
  60         for (int j = 0, k = 1 << n, jj, kk; j < k; ++j)
  61         {
  62             int p = 0;
  63             complex_type v(1);
  64             for (jj = 0, kk = 1; jj < n && !(j & kk); ++jj, ++p, kk <<= 1)
  65                 ;
  66             for (; jj < n; ++jj, kk <<= 1)
  67             {
  68                     if (j & kk)
  69                         v *= -r[jj];
  70                     else
  71                         ++p;
  72             }
  73             c[p] += v;
  74         }
  75
  76         bool pass = false;
  77         double div;
  78         for (int maxiter = 10; !pass && maxiter < 10000; maxiter *= 2)
  79         {
  80             for (int j = 0; j < n + 1; ++j)
  81                     a[j] = c[j].real();
  82
  83             CvMat amat, umat;
  84             cvInitMatHeader(&amat, n + 1, 1, CV_64FC1, &a[0]);
  85             cvInitMatHeader(&umat, n, 1, CV_64FC2, &u[0]);
  86             cvSolvePoly(&amat, &umat, maxiter, fig);
  87
  88             for (int j = 0; j < n; ++j)
  89                     ar[j] = complex_type(u[j * 2], u[j * 2 + 1]);
  90
  91             sort(r.begin(), r.end(), pred_complex());
  92             sort(ar.begin(), ar.end(), pred_complex());
  93
  94             div = 0;
  95             double s = 0;
  96             for (int j = 0; j < n; ++j)
  97             {
  98                     s += r[j].real() + fabs(r[j].imag());
  99                     div += pow(r[j].real() - ar[j].real(), 2) + pow(r[j].imag() - ar[j].imag(), 2);
 100             }
 101             div /= s;
 102             pass = div < 1e-2;
 103         }
 104
 105         if (!pass)
 106         {
 107             ts->set_failed_test_info(CvTS::FAIL_INVALID_OUTPUT);
 108             ts->printf( CvTS::LOG, "\n" );
 109
 110             for (size_t j=0;j<r.size();++j)
 111                     ts->printf( CvTS::LOG, "r[%d]=(%g, %g)\n", j, r[j].real(), r[j].imag());
 112             ts->printf( CvTS::LOG, "\n" );
 113             for (size_t j=0;j<ar.size();++j)
 114                     ts->printf( CvTS::LOG, "ar[%d]=(%g, %g)\n", j, ar[j].real(), ar[j].imag());
 115         }
 116
 117         progress = update_progress(progress, idx-1, max_idx, 0);
 118     }
 119 }
 120
 121 CV_SolvePolyTest solve_poly_test;
 122
 123 #endif