Publicité
+ Répondre à la discussion
Affichage des résultats 1 à 7 sur 7
  1. #1
    Membre éclairé Avatar de Flo.
    Homme Profil pro Florian
    Inscrit en
    mai 2002
    Messages
    379
    Détails du profil
    Informations personnelles :
    Nom : Homme Florian
    Localisation : France, Haute Garonne (Midi Pyrénées)

    Informations forums :
    Inscription : mai 2002
    Messages : 379
    Points : 350
    Points
    350

    Par défaut [c++] Régression circulaire

    Salut,

    Pour ma première contribution (c++), je vous propose une régression circulaire sur un ensemble de points (x; y).

    Voilà la déclaration de la classe Matrix (matrix.h) :

    Code c++ :
    1
    2
    3
    4
    5
    6
    7
    8
    9
    10
    11
    12
    13
    14
    15
    16
    17
    18
    19
    20
    21
    22
    23
    24
    25
    26
    #ifndef __MATRIX__
    #define __MATRIX__
     
    #include <vector>
     
    struct Matrix
    {
    	/// @brief Solve the system AX=Y where A and Y are known and X is unknown.
    	/// @param [in] A A n*n matrix where n = A.size() / 2 = Y.size().
    	/// @param [in] X The matrix to find. X is resized if needed.
    	/// @param [in] Y A 1*n matrix where n = A.size() / 2 = Y.size().
    	/// @throw std::runtime_error If an error occured.
    	static void gaussjordan(std::vector<double> const & A, std::vector<double> & X, std::vector<double> const & Y);
     
    	/// @brief Process to the regression of the points (xlist; ylist) by a circle.
    	/// @param [in] xlist List of the x coordinates of the points that should fit a circle.
    	/// @param [in] ylist List of the y coordinates of the points that should fit a circle.
    	/// @param [in] xcenter The x coordinate of the center of the circle computed from the points (xlist; ylist).
    	/// @param [in] ycenter The y coordinate of the center of the circle computed from the points (xlist; ylist).
    	/// @param [in] radius The radius of the circle computed from the points (xlist; ylist).
    	/// @param [in] rmse The root mean square error of the regression computed as the root squared sum of the difference between the computed radius and the distances between the computed center and the points (xlist; ylist)
    	/// @throw std::runtime_error If an error occured.
    	static void circularRegression(std::vector<double> const & xlist, std::vector<double> const & ylist, double & xcenter, double & ycenter, double & radius, double & rmse);
    }; // Matrix
     
    #endif
    et voilà l'implémentation de la classe Matrix (matrix.cpp) :

    Code c++ :
    1
    2
    3
    4
    5
    6
    7
    8
    9
    10
    11
    12
    13
    14
    15
    16
    17
    18
    19
    20
    21
    22
    23
    24
    25
    26
    27
    28
    29
    30
    31
    32
    33
    34
    35
    36
    37
    38
    39
    40
    41
    42
    43
    44
    45
    46
    47
    48
    49
    50
    51
    52
    53
    54
    55
    56
    57
    58
    59
    60
    61
    62
    63
    64
    65
    66
    67
    68
    69
    70
    71
    72
    73
    74
    75
    76
    77
    78
    79
    80
    81
    82
    83
    84
    85
    86
    87
    88
    89
    90
    91
    92
    93
    94
    95
    96
    97
    98
    99
    100
    101
    102
    103
    104
    105
    106
    107
    108
    109
    110
    111
    112
    113
    114
    115
    116
    117
    118
    119
    120
    121
    122
    123
    124
    125
    126
    127
    128
    129
    130
    131
    132
    133
    134
    135
    136
    137
    138
    139
    140
    141
    142
    143
    144
    145
    146
    147
    148
    149
    150
    151
    152
    153
    154
    155
    156
    157
    158
    159
    160
    161
    162
    163
    164
    165
    166
    167
    168
    169
    #include <limits>
    #include <stdexcept>
    #include <utility>
    #include <map>
    #include <cmath>
    #include "matrix.h"
     
    void Matrix::gaussjordan(std::vector<double> const & A, std::vector<double> & X, std::vector<double> const & Y)
    {
      // get the size of the matrix and check system input
      unsigned int const n = Y.size();
      unsigned int const m = n + 1;
      if(A.size() != n * n)
        {
          throw std::runtime_error("Bad size with A and Y arguments");
        }
     
      // resize the ouptut 
      if(X.size() != n)
        {
          X.resize(n);
        }
     
      std::map<unsigned int, unsigned int> swapped;
      for(unsigned int row = 0; row < n; ++row)
        {
          swapped[row] = row;
        }
     
      // here is the system to solve
      std::vector<double> system(n * n + n);
     
      // array of pointers on each row of the system
      std::vector<double *> rows(n);
     
      // setup the system
      std::vector<double>::const_iterator a = A.begin();
      std::vector<double>::const_iterator y = Y.begin();
      std::vector<double>::iterator s = system.begin();
      for(unsigned int row = 0; row < n; ++row)
        {
          // keep the address of the first item of the current row
          rows[row] = &(*s);
     
          // set A
          for(unsigned int col = 0; col < n; ++col, ++a, ++s)
    	{
    	  *s = *a;	  
    	}
     
          // set Y
          *s++ = *y++;
        }
     
      // find X to solve  the system AX=Y
      for(unsigned int ipivot = 0; ipivot < n; ++ipivot)
        {
          double pivot = rows[ipivot][ipivot];
     
          // if the current pivot is null, search for another line with a valid pivot (for the column ipivot)
          if(fabs(pivot) <= std::numeric_limits<double>::epsilon())
    	{
    	  // browse next row, searching for a valid pivot at column ipivot
    	  for(unsigned int row = ipivot + 1; row < n; ++row)
    	    {
    	      // if found, swap the rows pointer (not the memory)
    	      pivot = rows[row][ipivot];
    	      if(fabs(pivot) > std::numeric_limits<double>::epsilon())
    		{
    		  std::swap(rows[ipivot], rows[row]);
    		  std::swap(swapped[ipivot], swapped[row]);
    		  break;
    		}	      
    	    }	  
    	}
     
          // if a valid pivot was not found, then the matrix can not be inversed
          if(fabs(pivot) <= std::numeric_limits<double>::epsilon())
    	{
    	  throw std::runtime_error("The system can not be solved");
    	}
     
          // normalize the row of the pivot so as the pivot is set to 1.0
          for(unsigned int col = 0; col < m; ++col)
    	{
    	  rows[ipivot][col] /= pivot;
    	}
     
          // for all row, except the one of the pivot, set the item of the column "ipivot" ot 0
          for(unsigned int row = 0; row < n; ++row)
    	{
    	  if(row == ipivot)
    	    {
    	      continue;
    	    }
     
    	  double const factor = rows[row][ipivot];
    	  for(unsigned int col = 0; col < m; ++col)
    	    {
    	      rows[row][col] -= factor * rows[ipivot][col];
    	    }
    	}
        }
     
      // now the system is solved, X is the last column of the system
      for(unsigned int row = 0; row < n; ++row)
        {      
          X[row] = system[swapped[row] * m + m - 1];
        }
    }
     
    void Matrix::circularRegression(std::vector<double> const & xlist, std::vector<double> const & ylist, double & xcenter, double & ycenter, double & radius, double & rmse)
    {
      xcenter = 0.0;
      ycenter = 0.0;
      radius = 0.0;
      rmse = 0.0;
     
      // compute the A (3x3 matrix) and B (3x1 matrix)
      std::vector<double> A(9);
      std::vector<double> B(3);
      std::fill(A.begin(), A.end(), 0);
      std::fill(B.begin(), B.end(), 0);
      for(unsigned int index = 0; index < xlist.size(); ++index)
        {
          double const & x = xlist[index];
          double const & y = ylist[index];
          double const xx = x * x;
          double const yy = y * y;
          double const xy = x * y;
          double const xxx = xx * x;
          double const xxy = xx * y;
          double const yyy = yy * y;
          double const xyy = x * yy;
          A[0] += 2 * xx;
          A[1] += 2 * xy;
          A[2] += x;
          A[3] += 2 * xy;
          A[4] += 2 * yy;
          A[5] += y;
          A[6] += 2 * x;
          A[7] += 2 * y;
          A[8] += 1;
          B[0] += xxx + xyy; 
          B[1] += xxy + yyy;
          B[2] += xx + yy;
        }
     
      // solve the system AX=B using a the gauss-jordan method
      std::vector<double> X(3);
      gaussjordan(A, X, B);
     
      // circle details extraction
      xcenter = X[0];
      ycenter = X[1];
      radius = sqrt(X[2] + xcenter * xcenter + ycenter * ycenter);
     
      // compute the root mean square
      for(unsigned int index = 0; index < xlist.size(); ++index)
        {
          double const & x = xlist[index];
          double const & y = ylist[index];
          double const distance = sqrt((x - xcenter) * (x - xcenter) + (y - ycenter) * (y - ycenter));
          rmse += (radius - distance) * (radius - distance);
        }
     
      rmse /= (double)xlist.size();
      rmse = sqrt(rmse);
    }
    Le tout s'utilise très facilement (en c++ je rappelle) via le code :

    Code c++ :
    1
    2
    3
    4
    5
    6
    7
    8
    9
    10
    11
    12
    13
    14
    15
    16
    17
    18
    19
    20
    21
    // variables
    std::vector<double> xlist; // ... The list of the x coordinates of the points that match a circle.
    std::vector<double> ylist; // ... The list of the y coordinates of the points that match a circle.
    double xcenter = 0.0; // To store the x coordinate of the center of the circle computed by the function.
    double ycenter = 0.0; // To store the y coordinate of the center of the circle computed by the function.
    double radius = 0.0; // To store the radius of the circle computed by the function.
    double rmse = 0.0; // To store the root mean square error resulted from the regression.
     
    try
    {
       // regression of the points by a circle
       Matrix::circularRegression(xlist, ylist, xcenter, ycenter, radius, rmse);
    }
    catch(std::exception const & error)
    {
       std::cout << error.what() << std::endl;
    }
    catch(...)
    {
       std::cout << "An undefined error has occured" << std::endl;
    }
    Et un main (main.cpp) pour tester tout ça :

    Code c++ :
    1
    2
    3
    4
    5
    6
    7
    8
    9
    10
    11
    12
    13
    14
    15
    16
    17
    18
    19
    20
    21
    22
    23
    24
    25
    26
    27
    28
    29
    30
    31
    32
    33
    34
    35
    36
    37
    38
    39
    40
    41
    42
    43
    44
    45
    46
    47
    48
    49
    50
    51
    52
    53
    54
    55
    56
    57
    58
    59
    60
    61
    62
    63
    64
    65
    66
    67
    68
    69
    70
    71
    72
    73
    74
    75
    76
    77
    78
    79
    80
    81
    82
    83
    84
    85
    86
    87
    88
    89
    90
    91
    92
    93
    94
    95
    96
    97
    98
    99
    100
    101
    102
    103
    104
    #include <iostream>
    #include <iomanip>
    #include <cmath>
    #include <stdio.h>
    #include <stdlib.h>
    #include "matrix.h"
     
    struct Circle
    {
      double xcenter;
      double ycenter;
      double radius;
    }; // Circle
     
    int main()
    {
      // for std::stw
      unsigned int const setwv = 5;
     
      // definition of the circle
      Circle const expected = {20.5, -3.7, 8.0};
      Circle computed = {0.0, 0.0, 0.0};
     
      // lists of the points on the circle used to make the regression
      std::vector<double> xlist;
      std::vector<double> ylist;
     
      // display the informations about the expected circle 
      std::cout << "  ---------------------------------------------------------------------";
      std::cout << std::endl;
      std::cout << "||" << std::right << std::setw(2 * setwv) << "Expected circle : " << std::endl;
      std::cout << "||" << std::right << std::setw(4 * setwv) <<  "x center = " << expected.xcenter << std::endl;
      std::cout << "||" << std::right << std::setw(4 * setwv) <<  "y center = " << expected.ycenter << std::endl;
      std::cout << "||" << std::right << std::setw(4 * setwv) <<  "radius = " << expected.radius << std::endl;
     
      // display the titles of the informations about the points on the circle
      std::cout << "  ---------------------------------------------------------------------";
      std::cout << std::endl;
      std::cout << "||"
    	    << std::right << std::setw(2 * setwv) << "x original" << "|" 
    	    << std::right << std::setw(2 * setwv) << "x noised" << "|"  
    	    << std::right << std::setw(2 * setwv) << "noise"
    	    << std::right << std::setw(5) << "||"
    	    << std::right << std::setw(2 * setwv) << "y original" << "|"
    	    << std::right << std::setw(2 * setwv) << "y noised" << "|" 
    	    << std::right << std::setw(2 * setwv) << "noise" << "||";
      std::cout << std::endl;
      std::cout << "  ---------------------------------------------------------------------";
      std::cout << std::endl;
     
      // compute the points used to make the regression
      for(unsigned int angle = 0; angle < 360; angle += 12)
        {
          double const radian = 2 * 3.1415 * (double)angle / 360.0;
          double const x = expected.xcenter + expected.radius * cos(radian);
          double const y = expected.ycenter + expected.radius * sin(radian);
          double const xnoise = 2.0 * (double)(rand() % 200 - 100) / 100.0;
          double const ynoise = 2.0 * (double)(rand() % 200 - 100) / 100.0;
          xlist.push_back(x + xnoise);
          ylist.push_back(y + ynoise);
     
          // values
          std::cout << "||"
    		<< std::right << std::setw(2 * setwv) << x<< "|" 
    		<< std::right << std::setw(2 * setwv) << xlist.back() << "|" 
    		<< std::right << std::setw(2 * setwv) << xnoise
    		<< std::right << std::setw(5) << "||"
    		<< std::right << std::setw(2 * setwv) << y << "|" 
    		<< std::right << std::setw(2 * setwv) << ylist.back() << "|" 
    		<< std::right << std::setw(2 * setwv) << ynoise << "||";
          std::cout << std::endl;
        }
     
      double rmse = 0.0;
     
      try
        {
          // regression of the points by a circle
          Matrix::circularRegression(xlist, ylist, computed.xcenter, computed.ycenter, computed.radius, rmse);
        }
      catch(std::exception const & error)
        {
          std::cout << error.what() << std::endl;
          return 0;
        }
      catch(...)
        {
          std::cout << "An undefined error has occured" << std::endl;
          return 0;
        }
     
      // display the informations about the computed circle 
      std::cout << "  ---------------------------------------------------------------------";
      std::cout << std::endl;
      std::cout << "||" << std::right << std::setw(2 * setwv) << "Computed circle : " << std::endl;
      std::cout << "||" << std::right << std::setw(4 * setwv) <<  "x center = " << computed.xcenter << std::endl;
      std::cout << "||" << std::right << std::setw(4 * setwv) <<  "y center = " << computed.ycenter << std::endl;
      std::cout << "||" << std::right << std::setw(4 * setwv) <<  "radius = " << computed.radius << std::endl;
      std::cout << "||" << std::right << std::setw(4 * setwv) <<  "rmse = " << rmse << std::endl;
      std::cout << "  ---------------------------------------------------------------------";
      std::cout << std::endl;
     
      return 0;
    }
    Le code ne requiert aucune autre dépendance que la Standart Template Library (native sur tous les compilateurs), pas de boost, opencv ou autres. Vous compilez sous gcc par la commande suivante :

    g++ main.cpp matrix.cpp -o test

    La sortie est du style :

    Code :
    1
    2
    3
    4
    5
    6
    7
    8
    9
    10
    11
    12
    13
    14
    15
    16
    17
    18
    19
    20
    21
    22
    23
    24
    25
    26
    27
    28
    29
    30
    31
    32
    33
    34
    35
    36
    37
    38
    39
    40
    41
    42
    43
    44
    45
      ---------------------------------------------------------------------
    ||Expected circle :
    ||         x center = 20.5
    ||         y center = -3.7
    ||           radius = 8
      ---------------------------------------------------------------------
    ||x original|  x noised|     noise   ||y original|  y noised|     noise||
      ---------------------------------------------------------------------
    ||      28.5|     27.32|     -1.18   ||      -3.7|     -4.36|     -0.66||
    ||   28.3252|   29.0052|      0.68   ||  -2.03675|  -2.03675|         0||
    ||   27.8084|   29.1884|      1.38   || -0.446197| 0.0338029|      0.48||
    ||   26.9722|   26.5322|     -0.44   ||   1.00216|   2.16216|      1.16||
    ||   25.8532|   27.0932|      1.24   ||   2.24503|   1.52503|     -0.72||
    ||   24.5002|   24.6002|       0.1   ||   3.22808|   4.12808|       0.9||
    ||   22.9724|   22.5924|     -0.38   ||   3.90836|   2.44836|     -1.46||
    ||   21.3366|   22.5566|      1.22   ||   4.25614|   4.07614|     -0.18||
    ||   19.6642|   21.5642|       1.9   ||   4.25622|   5.09622|      0.84||
    ||   18.0283|   16.5683|     -1.46   ||   3.90859|   2.62859|     -1.28||
    ||   16.5004|   18.3204|      1.82   ||   3.22845|   1.30845|     -1.92||
    ||   15.1474|   15.1874|      0.04   ||   2.24552|   3.30552|      1.06||
    ||   14.0282|   13.8682|     -0.16   ||   1.00276|   2.64276|      1.64||
    ||   13.1919|   11.6119|     -1.58   ||  -0.44552|  -0.12552|      0.32||
    ||    12.675|    13.035|      0.36   ||  -2.03603|  -2.13603|      -0.1||
    ||      12.5|     11.44|     -1.06   ||  -3.69926|  -3.17926|      0.52||
    ||   12.6747|   14.0947|      1.42   ||  -5.36252|  -4.60252|      0.76||
    ||   13.1913|   12.5713|     -0.62   ||  -6.95313|  -6.71313|      0.24||
    ||   14.0273|   13.3673|     -0.66   ||  -8.40156|  -8.42156|     -0.02||
    ||   15.1463|   13.8463|      -1.3   ||  -9.64453|  -9.76453|     -0.12||
    ||   16.4991|   16.5591|      0.06   ||  -10.6277|  -12.4077|     -1.78||
    ||   18.0269|   18.4669|      0.44   ||  -11.3081|  -10.6481|      0.66||
    ||   19.6627|   19.1227|     -0.54   ||  -11.6561|  -12.3761|     -0.72||
    ||   21.3351|   22.1551|      0.82   ||  -11.6563|  -11.4363|      0.22||
    ||    22.971|    22.031|     -0.94   ||  -11.3088|  -11.9488|     -0.64||
    ||   24.4989|   25.4389|      0.94   ||  -10.6288|  -11.7488|     -1.12||
    ||   25.8521|   25.0921|     -0.76   ||  -9.64602|  -8.50602|      1.14||
    ||   26.9714|   25.7114|     -1.26   ||  -8.40336|  -9.22336|     -0.82||
    ||   27.8078|   28.2678|      0.46   ||  -6.95516|  -6.13516|      0.82||
    ||   28.3249|   28.9049|      0.58   ||   -5.3647|   -3.8047|      1.56||
      ---------------------------------------------------------------------
    ||Computed circle :
    ||         x center = 20.3688
    ||         y center = -3.76401
    ||           radius = 8.25942
    ||             rmse = 1.01834
      ---------------------------------------------------------------------
    Bon je veux pas faire de la pub pour mon blog, mais je veux pas non plus ré-expliquer le fonctionnement du code. Le lien vers des explications détaillées est le suivant :

    http://floriansella.free.fr/index.ph...y=ti&billid=92

    Histoire de pas faire un truc creux, la régression d'un ensemble de points par un cercle (xc; yc; R) fonctionne comme la régression d'un ensemble de points par une droite. C'est une régression au sens des moindres carrés et on parvient à un système matriciel qu'il suffit de résoudre par la méthode basique de Gauss-Jordan (pivot de Gauss).

    Flo.

  2. #2
    Membre Expert
    Homme Profil pro
    Chercheur
    Inscrit en
    mars 2010
    Messages
    1 175
    Détails du profil
    Informations personnelles :
    Sexe : Homme
    Localisation : France, Paris (Île de France)

    Informations professionnelles :
    Activité : Chercheur

    Informations forums :
    Inscription : mars 2010
    Messages : 1 175
    Points : 1 684
    Points
    1 684

    Par défaut

    Bonjour,

    merci beaucoup pour cette implémentation.
    Voici quelques améliorations possibles :

    1. ne pas inverser explicitement A pour calculer y=inv(A)*x mais résoudre directement Ax=y par élimination de Gauss;
    2. faire en sorte que Matrix soit bien une matrice : en lisant rapidement, je n'ai pas vu de stockage des coefficients;
    3. encapsuler tes données;
    4. passer par class plutôt que par struct;
    5. ne pas mélanger codes d'erreur et exceptions; faire soit l'un soit l'autre.
    6. paramétrer avec des templates pour laisser à l'utilisateur le choix de la précision des calculs et éventuellement du type de données;
    7. éviter de faire du static quand ce n'est pas la peine.

    Bonne continuation

  3. #3
    Membre éclairé Avatar de Flo.
    Homme Profil pro Florian
    Inscrit en
    mai 2002
    Messages
    379
    Détails du profil
    Informations personnelles :
    Nom : Homme Florian
    Localisation : France, Haute Garonne (Midi Pyrénées)

    Informations forums :
    Inscription : mai 2002
    Messages : 379
    Points : 350
    Points
    350

    Par défaut

    Merci pour les commentaires, Aleph69.

    Ce code est plutôt destiné à mon éducation (et celle de ceux que mon code intéresse) ... j'ai pas pour objectif de le proposer chez boost (ils doivent déjà avoir des trucs pour faire ça je suppose) .

    Ceci étant dit, je vais corriger le point 5 qui est effectivement maladroit même pour un code de ce type, et proposer une autre version en considérant votre point 1, qui m'intéresse évidemment.

    Toutefois, si je n'ai rien à redire à vos remarques, qui sont sûrement pertinentes, je laisserai à tout un chacun le choix de l'implémentation d'une version propre à ses besoins, si tant est que quelques uns s'entichent de ces quelques lignes de codes.

    Flo.

  4. #4
    Membre éclairé Avatar de Flo.
    Homme Profil pro Florian
    Inscrit en
    mai 2002
    Messages
    379
    Détails du profil
    Informations personnelles :
    Nom : Homme Florian
    Localisation : France, Haute Garonne (Midi Pyrénées)

    Informations forums :
    Inscription : mai 2002
    Messages : 379
    Points : 350
    Points
    350

    Par défaut Corrections du code initial

    J'ai corrigé les points 1 et 5 décrits par Aleph69 en lieu et place du code initial dans mon premier message.

    Encore merci Aleph69.

    Flo.

  5. #5
    Membre Expert
    Homme Profil pro
    Chercheur
    Inscrit en
    mars 2010
    Messages
    1 175
    Détails du profil
    Informations personnelles :
    Sexe : Homme
    Localisation : France, Paris (Île de France)

    Informations professionnelles :
    Activité : Chercheur

    Informations forums :
    Inscription : mars 2010
    Messages : 1 175
    Points : 1 684
    Points
    1 684

    Par défaut

    J'ai encore une remarque, j'espère que ça ne pose de problème. Pour s'assurer qu'on ne provoque pas un overflow, on ne compare pas le pivot avec epsilon. La fonction dlamch de la bibliothèque lapack calcule un nombre safemin correspondant au plus petit nombre par lequel on peut diviser sans provoquer d'overflow. Je pense que c'est à cette quantité que tu devrais comparer les pivots. Moralement, il faut que c=b/a ne dépasse le plus grand nombre flottant représentable en machine MAX ( il est donné par std::numeric_limits<T>::max() en C++ ). On doit donc avoir |b/a|=|b|/|a|<=MAX, soit encore |b| <= MAX*|a|. Je pense que cela revient à passer par dlamch mais à vérifier quand même, il y a peut-être une subtilité. En tous les cas, la routine dlamch est ce qui se fait de plus robuste dans le domaine. En testant avec epsilon, tu arrêtes parfois prématurément l'algorithme alors que la matrice était bien factorisable à la précision machine.

    EDIT : ajout du lien vers la routine dlamch.

    EDIT2 : encore une chose, en C++, on préfère inclure cmath plutôt que math.h

  6. #6
    Membre éclairé Avatar de Flo.
    Homme Profil pro Florian
    Inscrit en
    mai 2002
    Messages
    379
    Détails du profil
    Informations personnelles :
    Nom : Homme Florian
    Localisation : France, Haute Garonne (Midi Pyrénées)

    Informations forums :
    Inscription : mai 2002
    Messages : 379
    Points : 350
    Points
    350

    Par défaut

    Salut Aleph69,

    J'ai remplacé l'include de math.h pour celui de cmath. Quant à ton histoire de l'epsilon, utiliser LAPACK ça ne m'intéresse pas. En revanche, je veux bien essayer de comprendre ce que tu me racontes et de l'intégrer dans mon code. J'ai donc fait une fonction qui teste un pivot sur une ligne et renvoie true si le résultat de la multiplication du dénominateur (de chaque division) par std::numeric_limits<double>::max() est inférieur au pivot, false sinon.

    Code c++ :
    1
    2
    3
    4
    5
    6
    7
    8
    9
    10
    11
    12
    13
    14
    15
    16
    17
    18
    19
    20
    21
    #define INSPECT_PIVOT 0
     
    bool isGoodPivot(double const * const row, unsigned int const size, double const pivot)
    {
    #if INSPECT_PIVOT
    	return (fabs(pivot) > std::numeric_limits<double>::epsilon());	
    #else 
    	if(fabs(pivot) < std::numeric_limits<double>::epsilon())
    	{
    		for(unsigned int col = 0; col < size; ++col)
    		{	
    			if(fabs(row[col]) > fabs(pivot * std::numeric_limits<double>::max()))
    			{			
    				return false;
    			}
    		}
    	}
     
    	return true;
    #endif	
    }

    J'ai fait un gauss-jordan sur le système AX=B suivant (sur ma machine std::numeric_limits<double>::epsilon() vaut 2.22045e-016) :

    Code :
    1
    2
    3
    4
    A[0] = +1e-17; A[1] = +2; A[2] = -9;
    A[3] = +1    ; A[4] = +6; A[5] = -8;
    A[6] = +3    ; A[7] = -5; A[8] = +3;
    B[0] = +2    ; B[1] = +1; B[2] = -3;
    Avec INSPECT_PIVOT à 0, il y a un changement de pivots sur la première ligne et les résultats sont corrects. Les valeurs de X sont :

    Code :
    1
    2
    3
    X[0] = -0.777778
    X[1] = 5.55112e-017
    X[2] = -0.222222
    Avec INSPECT_PIVOT à 1, il n'y a aucun changement de pivots mais les résultats sont faux. Les valeurs de X sont :

    Code :
    1
    2
    3
    X[0] = 0
    X[1] = 1
    X[2] = 0
    Faut dire qu'avec INSPECT_PIVOT à 1, les valeurs du système en cours de résolution montent à des + ou - 1e+17 et les opérations suivantes doivent en pâtir.

    Donc je sais pas si j'ai bien tout compris.

    Flo.

  7. #7
    Membre Expert
    Homme Profil pro
    Chercheur
    Inscrit en
    mars 2010
    Messages
    1 175
    Détails du profil
    Informations personnelles :
    Sexe : Homme
    Localisation : France, Paris (Île de France)

    Informations professionnelles :
    Activité : Chercheur

    Informations forums :
    Inscription : mars 2010
    Messages : 1 175
    Points : 1 684
    Points
    1 684

    Par défaut

    Bonjour,

    il y a plusieurs phénomènes en jeu dans ce que tu décris.

    Une première question est de savoir si une matrice est factorisable à la précision machine, c'est-à-dire si le processus d'élimination gaussienne peut être mené à son terme sans division par zéro, sans overflow et sans underflow, toute considération sur la stabilité numérique de l'algorithme mis à part. C'est une question d'informatique pur : l'algorithme arrive-t-il à son terme? Ainsi, lorsque tu divises par un pivot, il faut que cette division soit bien représentable en machine. Je te faisais juste remarquer que tu pouvais diviser par un nombre inférieur à epsilon sans nécessairement provoquer d'overflow.

    Pour tester cela, je t'ai indiqué la fonction dlamch, non pas pour l'utiliser mais pour regarder ce qui y est fait. Il y est défini un nombre sfmin correspondant au plus petit nombre par lequel on peut diviser un autre nombre sans provoquer d'overflow. Ce nombre est défini de la manière suivante :
    Code :
    1
    2
    3
    T sfmin = std::numeric_limits<T>::min();
    T small = 1./std::numeric_limits<T>::max();
    if (small>=sfmin) sfmin=small*(1+std::numeric_limits<T>::epsilon());
    Ensuite, dans ton code, tu dois uniquement comparer la valeur absolue du pivot à sfmin au lieu de epsilon; il n'est pas utile de tester toutes les divisions possibles par ton pivot. Si tu ne trouves aucun pivot admissible, alors la matrice n'est pas factorisable à la précision machine; on la dit quasi-singulière.

    La question de la stabilité numérique de l'algorithme est une autre chose. Quelle que soit le choix de ton pivot, l'erreur relative sur la solution sera au mieux bornée supérieurement par un terme supérieur à 2^(n-1) où n désigne la taille de la matrice. On connaît d'ailleurs les matrices pour lesquelles cette borne est atteinte. De ce fait, trouver une solution plus précise avec telle ou telle méthode de pivot ne la rend pas meilleure : on saura toujours trouver des matrices pour lesquelles l'erreur sur la solution explose. A noter que l'approche standard appelée "pivotage partiel par ligne" consiste à choisir pour pivot celui dont la valeur absolue est maximale. Le test de quasi-singularité n'intervient qu'après pour tester si cette valeur maximale est bien supérieure à sfmin.

    Par ailleurs, si tu veux améliorer la stabilité de ton algorithme, il faut pré-traiter ton système linéaire à l'aide d'un préconditionneur. Une méthode simple consiste par exemple à diviser chaque ligne par sa norme. Sa variante par colonne est également très utilisée. Tu peux également raffiner itérativement la solution calculée en mélangeant plusieurs degrés de précision (float/double ou double/long double); ceci te permettra également de fournir des estimations des erreurs de calcul à tes utilisateurs.

    Enfin, si le sujet t'intéresse, tu peux tenter de percer les derniers mystères de la méthode de Gauss.

Liens sociaux

Règles de messages

  • Vous ne pouvez pas créer de nouvelles discussions
  • Vous ne pouvez pas envoyer des réponses
  • Vous ne pouvez pas envoyer des pièces jointes
  • Vous ne pouvez pas modifier vos messages
  •