Radial-Functions-Interpolation/main.cpp

#include "RBFInterp.hpp"
#include "interp.hpp"
#include "gnuplot.hpp"
#include <iostream>
#include <fstream>
#include <cstdlib>
#include <ctime>
#include <algorithm>
#include <string>
#include "time_mes.hpp"
using namespace std;
// Fonction à interpoler z= f(x,y)
float finit(int x,int y)
{
  return 1.0+float(x*sin(float(x)/1000.)+2*y)/1000.;
}

int main()
{

  typedef float real;
  //  voir RBFInterp qui définit ceci:
  typedef typename RBFInterp<real,Gauss<real>>::PointValue PointValue;
  // qui est une map (i,j) -> valeur flottante.

  // taille de l'"image":
  int nx=1500,ny=1500;
  cout<<"nx : "<<nx<<" ny : "<<ny<<endl;
  real* image= new real[nx*ny];
  //
  
  std::srand(std::time(nullptr)); // seed with current time.

  // Choisir npts points d'interpolation au hasard  dans l'image:
  PointValue Values;
  int npts= 50;
  for(int i=0;i<npts;i++)
    {
      // la position des points est choisie au hasard:
      auto vrand=float(std::rand())/RAND_MAX;
      int xx= static_cast<int>(nx*vrand);
      vrand=float(std::rand())/RAND_MAX;
      int yy= static_cast<int>(ny*vrand);
      auto k=make_pair(xx,yy);
      // valeur de l'image en ce point.
      Values[k]=finit(xx,yy);
      if(xx<0 || xx>=nx||yy<0 || yy>=ny)
	// vérifier que les points sont dans l'image.
	throw "pas bien";
     }
  cout<<"nb. values : "<<Values.size()<<endl;
  //
  double start,stop;
 
  // Choisir au hasard une partie des points d'interpolation:
  PointValue Test,NoTest;
  float portion_in_Notest = 0.8;
  for(typename PointValue::const_iterator I =Values.begin();I!=Values.end();I++)
    {
      float srand=float(std::rand())/RAND_MAX;
      if(srand<portion_in_Notest)
	// les points qui vont servir a construire l'interpolation.
	NoTest[I->first]= I->second;
      else
	// les points de test.
	Test[I->first]= I->second;
    }
  cout<<"estimation d'epsilon. nb. noeuds : "<<NoTest.size()<<
    " nb. points de test : "<<Test.size()<<endl;
  
  // interpolant:
  // définition de la fontion radiale utilisée (MQI semble mailleure)
  typedef MQI<real> typRad;
  //typedef TPS<real> typRad;
  //typedef TPSD<real> typRad;
  //typedef Gauss<real> typRad;

  real eps= 0.01; //valeur de départ de eps
  RBFInterp<real,typRad> A(eps);

  //
  
  // Les points d'interpolation :
  gnuplotfile(Values,"points");


  //--1--- Détermination d'un epsilon "optimal" pour les RBF.
  ofstream f; f.open("eps-err");//pour faire un graphique precision(epsilon)
  start=get_time();
  real epsm= eps;
  real testm=1.e10;
  for(int iter=1;iter<=25;iter++)
    {
      eps/= 10.;
      A.set_epsilon(eps);
      bool calcond=false; // on calcule ou non le conditionnement (pas
                          // nécessaire, mais peut donner un renseignement
                          //utile).
      try{
	A.build(NoTest,calcond);
      }
      catch(const LapackException& e )
	{
	  // normalement quand epsilon est  trop petit, le système
	  // linéaire devient singulier. On va donc à priori venir ici
	  // plutôt que de sortir normalement de la boucle.
	  eps*=10;
	  break;
	}
      auto test= A.test(NoTest,Test);
      cout<<"eps : "<<eps<<" test : "<<test<<endl;;
      if(calcond)
	{
	  // si on a calculé le conditionnement de la matrice:
	  if(calcond)
	    {
	      cout<<"norme matrice : "<<A.norme_matrice();
	      cout<<" conditionnement : "<<A.cond();
	      cout<<endl<<endl;
	    }
	}
      f<<eps<<" "<<test<<endl;
      if(test<testm)
	{
	  // la valeur epsm choisie pour eps.
	  testm=test;
	  epsm=eps;
	}
    }
  stop=get_time();
  cout <<endl<<"--- (durée)  preparation :"<<stop-start<<endl;
  f.close();
  cout<<"estimation ok"<<endl;
  
  // On introduit le meilleur epsilon trouvé pour les RBF.
  A.set_epsilon(epsm);
  cout<<"epsilon utilisé : "<<epsm<<endl;
  start=get_time();
  try{
    A.build(Values);
  }
  catch(const LapackException& e )
   {
    cout<<"Lapack exception: "<<e.Info<<endl;
    throw;
  }
  stop=get_time();
  cout <<"--- (durée)  interpolants : "<<stop-start<<endl;
  cout<<"interpolants ok"<<endl;
  
  start=get_time();
  int increment=2; //pour utiliser l'interp. quad.
  // (un point sur 2 dans caque direction)
  A.Interp(Values,image,nx,ny,increment);
  stop=get_time();
  cout <<"--- (durée)  interpolation RBF : "<<stop-start<<endl;
  cout<<"interpolation ok"<<endl;
  start=get_time();
  interpquad(image,nx,ny);
  stop=get_time();
  cout <<"--- (durée)  interpolation quad. : "<<stop-start<<endl;
  string aa;cout<<"g pour sortir les fichiers gnuplot ou tout autre caractère sinon :";
  cin>>aa;
  if(aa=="g")
    {
      cout<<"gnuplot files :"<<endl;
      // interpolation sur la grille :
      gnuplotfile(nx,ny,image,"image");
    }
  // erreur :
  real* diff= new real[nx*ny];

  auto indix= [nx](int i, int j){return i*nx+j;};
#pragma omp parallel for
  for(int i=0;i<nx;i++)
    for(int j=0;j<ny;j++)
      diff[indix(i,j)]=abs(image[indix(i,j)]-finit(i,j));
  if(aa=="g")
    gnuplotfile(nx,ny,diff,"diff");

  // erreur relative
#pragma omp parallel for
  for(int i=0;i<nx;i++)
    for(int j=0;j<ny;j++)
      diff[indix(i,j)]/=abs(finit(i,j));
  auto vmax=max_element(diff,diff+nx*ny);
  cout<<"erreur relative max: "<<*vmax<<endl;
  if(aa=="g")
    gnuplotfile(nx,ny,diff,"diffrel");
 
  delete[] image;
  delete[] diff;
}