init

.gitignore

@@ -0,0 +1,18 @@
+*.log
+_*_.tex
+*.pyc
+publi.pdf
+exp.pdf
+*.aux
+toto*
+*out
+*.o
+*.so.*
+*.nav
+$.toc
+*.bbl
+*.blg
+Build*
+*~
+.ipynb_checkpoints
+Out*

CMakeLists.txt

@@ -0,0 +1,21 @@
+cmake_minimum_required(VERSION 3.0)
+project(RBF)
+enable_language(CXX)
+#
+set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS}  -O3 -Wall  -march=native \
+ -ftree-vectorize  -fstrict-aliasing -fopenmp")
+
+add_executable(
+  run
+  ../main.cpp
+  )
+
+  target_link_libraries(
+    run
+    lapack
+    blas
+    )
+
+
+
+

RBFInterp.hpp

@@ -0,0 +1,195 @@
+#include <vector>
+#include <iostream>
+#include <map>
+#include <cmath>
+#include "protos_lapack.hpp"
+#include "RBF_functions.hpp"
+using namespace std;
+using namespace lapack_c;
+#ifndef RBFInterp__hpp
+#define RBFInterp__hpp
+
+//-------------------------------------------------------------------------
+// Interpolation rbf.
+//-------------------------------------------------------------------------
+template<typename real,class Radial>  class RBFInterp{
+  double * matrix , *b;
+  real * image;
+  real eps;
+  int imgsize_x,imgsize_y;
+  double rcond,innorm; double *work; bool rcondok;
+public:
+  // (x,y) -> valeur
+  // Il ne faut pas d'ex-aequeo pour (x,y), d'où l'utilisation d'une map
+  // pour décrire les noeuds de l'interpolation.
+  typedef std::map<pair<int,int>,real> PointValue;
+
+  // consructeur:
+  RBFInterp(real epsilon):eps(epsilon){}
+  //changer epsilon
+  void set_epsilon(real epsilon){eps=epsilon;}
+  // Calcule les RBF.
+  void build(PointValue& values,bool comp_cond= false)
+  // comp_cond : on calcule ou non le conditionnement L\infinity de
+  // la matrice.
+  {
+    Radial  R(eps);
+    int np= values.size();
+   
+    matrix = new double[np*np];
+    b=  new double[np];
+
+    // construire le système linéaire (on n'utilise pas la symétrie de la
+    // matrice).
+
+    // 1) matrice:
+    auto  indix= [np](int i, int j){return i*np+j;};
+    int i=0;
+    for(typename PointValue::const_iterator I=values.begin();
+	I!=values.end();I++)
+      {
+	auto v= I->first; //la fonction est centrée en ce point.
+	int j=i;
+	for(typename PointValue::const_iterator J=I;J!=values.end();++J)
+	  {
+	    auto s = R(r(v,J->first));
+	    matrix[indix(i,j)]=s;
+	    if(i!=j)
+	      matrix[indix(j,i)] =s;
+	    ++j;
+	  }
+	++i;
+      }
+   
+    //2) second membre.
+    i=0;
+    for(typename PointValue::iterator I= values.begin();I!=values.end();I++)
+      b[i++]=I->second;
+
+    // resolution (lapack + blas).
+    int* ipiv= new int[np];
+    int nrhs=1,lda=np,ldb=np,info;
+    char trans='N'; char norm='I';
+
+    //
+    rcondok =  comp_cond;
+    if(comp_cond)
+      {
+	// norme oo de la matrice:
+	
+	work= new double[4*np];
+	// norme infinie de la matrice:
+	innorm= dlange_(&norm,&np,&np,matrix,&lda,work);
+      }
+    // factorisation PLU:
+    dgetrf_(&np,&np,matrix,&lda,ipiv,&info);
+    if(info != 0) throw  LapackException(info);
+    if(comp_cond)
+      {
+	int *iwork=new int[np];
+	dgecon_(&norm,&np,matrix,&lda,&innorm,&rcond,work,iwork,&info);
+	rcond= 1.0/rcond;
+	if(info != 0) throw  LapackException(info);
+	delete[] work; delete[] iwork;
+      }
+    // résolution du système factorisé:
+    dgetrs_(&trans,&np,&nrhs,matrix,&lda,ipiv,b,&ldb,&info);
+    if(info != 0) throw  LapackException(info);
+    //
+    delete[] matrix;
+    delete[] ipiv;
+  }
+ 
+  real test(const PointValue& interpvalues,const PointValue& testvalues)
+  {
+    // Estimer l'erreur.
+    
+    // interpvalues: les noeuds d'interpolation
+    // testvalues  : les noeuds test.
+     Radial  R(eps);
+    
+     real ret=0.0;
+     // boucle sur les noeuds test :
+     for(typename PointValue::const_iterator K=testvalues.begin();
+	 K!=testvalues.end();K++)
+       {
+	 int k=0;
+	 real u=0.0;
+	 for(typename PointValue::const_iterator L=interpvalues.begin();
+	     L!=interpvalues.end();L++)
+	   {
+	     auto dist= r(K->first,L->first);
+	     u+=b[k++]*R(dist);
+	   }
+	 ret=max(ret,abs(K->second-u));
+       }
+     return ret;
+  }
+  real Interpolator(PointValue& values,int nx,int i,int j,Radial R)
+  {
+    //auto  indix= [nx](int i, int j){return i*nx+j;};
+    real s=0.0;
+    int k=0;
+    for(typename PointValue::const_iterator K=values.begin();
+	K!=values.end();K++)
+      s+=b[k++]* R(r(K->first,make_pair(i,j)));
+    return s;
+  }
+  // interpolation sur une image nx par ny
+  void Interp(PointValue& values,real *image,int nx,int ny,int increment=1)
+  {
+    Radial  R(eps);
+    auto  indix= [nx](int i, int j){return i*nx+j;};
+
+    #pragma omp parallel for
+      for(int i=0;i<ny;i+=increment)
+	for(int j=0;j<nx;j+=increment)
+	{
+	  // real s=0.0;
+	  // int k=0;
+	  // for(typename PointValue::const_iterator K=values.begin();
+	  //     K!=values.end();K++)
+	  //   s+=b[k++]* R(r(K->first,make_pair(i,j)));
+	  //image[indix(i,j)] = s;
+	  image[indix(i,j)] =  Interpolator(values,nx,i,j,R);
+	}
+      if(increment == 2)
+	{
+	  for(int i=ny-3;i<ny;i++)
+	    for(int j=0;j<nx;j++)
+	      image[indix(i,j)] =  Interpolator(values,nx,i,j,R);
+	  for(int j=ny-3;j<nx;j++)
+	    for(int i=0;i<ny;i++)
+	      image[indix(i,j)] =  Interpolator(values,nx,i,j,R);
+	}
+	
+  }
+  ~RBFInterp()
+  {
+    delete[] b;
+   
+  }
+  double cond() const
+  {
+    // accesseur pour le conditionnement
+    if(rcondok)
+      return rcond;
+    else
+      throw "rcond pas calculé";
+  }
+ double norme_matrice() const
+  {
+    if(rcondok)
+      return innorm;
+    else
+      throw "rcond pas calculé";
+  } 
+private:
+  // distance euclidienne de 2 points .
+  inline float  r(pair<int,int> x1,pair<int,int> x2){return sqrt(r2(x1,x2));}
+  // distance euclidienne de 2 points au carré.
+  inline float  r2(pair<int,int> x1,pair<int,int> x2)
+  {return pow(get<0>(x1)-get<0>(x2),2) + pow(get<1>(x1)-get<1>(x2),2);}
+};
+  
+#endif

RBF_functions.hpp

@@ -0,0 +1,51 @@
+#ifndef RBF_functions_hpp
+#define RBF_functions_hpp
+// Les RBF
+template<typename real> class RBF{
+  // Radial Basis Functions: base.
+  protected:
+  real epsilon;
+  RBF(real eps):epsilon(eps){};
+  void set_epsilon(real eps) {epsilon = eps;}
+};
+
+
+template<typename real> class MQI: public RBF<real> {
+  // Multiquadratique inverse.
+ 
+public:
+  using RBF<real>::set_epsilon;
+  MQI(real eps): RBF<real>(eps){}
+  inline real operator()(real x){return 1./sqrt(1.+
+						RBF<real>::epsilon*x*x);}
+  };
+template<typename real> class TPS: public RBF<real> {
+  // plaque mince.
+public:
+  TPS(real eps): RBF<real>(eps){}
+  inline real operator()(real x){return exp(-RBF<real>::epsilon*x)*
+      (1.+RBF<real>::epsilon*x);}
+  };
+template<typename real> class Gauss: public RBF<real> {
+  // Gaussienne
+public:
+  Gauss(real eps): RBF<real>(eps){}
+  inline real operator()(real x){
+    if(x<1.e-10)
+      return 1.;
+    else
+      return exp(-RBF<real>::epsilon*x*x);}
+  };
+template<typename real> class TPSD: public RBF<real> {
+  // plaque mince à la Duchon..
+public:
+  TPSD(real eps): RBF<real>(eps){}
+  inline real operator()(real x){
+    if(x<1.e-9)
+      return 1.;
+    else
+      return 1/(1.+
+	      pow(RBF<real>::epsilon*x,2)*log(RBF<real>::epsilon*x));}
+  
+  };
+#endif

gnuplot.hpp

@@ -0,0 +1,30 @@
+#ifndef GNUPLOT__HPP
+#define GNUPLOT__HPP
+#include <fstream>
+// montrer les points d'interpolation.
+template<class M> void gnuplotfile(const M& Points,string fich)
+{
+  ofstream f;
+  f.open(fich);
+  for(typename M::const_iterator I=Points.begin();I!=Points.end();
+      I++)
+    {
+      auto p=I->first;
+      f<<p.first<<" "<<p.second<<endl;
+    }
+  f.close();
+}
+// montrer une "image"
+template<class real> void gnuplotfile(int nx,int ny,real* a,string fich)
+{
+  ofstream f;
+  f.open(fich);
+  for(int j=0;j<ny;j++)
+    {
+      for(int i=0;i<nx;i++)
+	f<<i<<" "<<j<<" "<<a[j*nx+i]<<endl;
+      f<<" "<<endl;
+    }
+  f.close();
+}
+#endif

interp.hpp

@@ -0,0 +1,87 @@
+#ifndef INTERP_HPP
+#define INTERP_HPP
+template<class real> inline void quadrati(real *vnodes,
+					  real& first, real& second)
+{
+
+  first  =  0.375*vnodes[0] + 0.75*vnodes[1] - 0.125*vnodes[2];
+  second = -0.125*vnodes[0] + 0.75*vnodes[1] + 0.375*vnodes[2];
+
+
+}
+template<class real> inline void interp1(real *temp, real *image,int nx,
+					 int i,int j)
+{
+  auto indix= [nx](int i, int j){return i*nx+j;};
+  int I=i-2,J=j-2;
+  //1) ligne horizontale en  haut :
+  temp[0]=image[indix(I-2,J-2)];
+  temp[1]=image[indix(I-2,J)];
+  temp[2]=image[indix(I-2,J+2)];
+  quadrati(temp,image[indix(I-2,J-1)],image[indix(I-2,J+1)]);
+	
+  //2) colonne verticale gauche :
+  temp[0]=image[indix(I-2,J-2)];
+  temp[1]=image[indix(I,J-2)];
+  temp[2]=image[indix(I+2,J-2)];
+  quadrati(temp,image[indix(I-1,J-2)],image[indix(I+1,J-2)]);
+
+  //3) colonne verticale centrale :
+  temp[0]=image[indix(I-2,J)];
+  temp[1]=image[indix(I,J)];
+  temp[2]=image[indix(I+2,J)];
+  quadrati(temp,image[indix(I-1,J)],image[indix(I+1,J)]);
+
+  //4)  colonne verticale à droite :
+  temp[0]=image[indix(I-2,J+2)];
+  temp[1]=image[indix(I,J+2)];
+  temp[2]=image[indix(I+2,J+2)];
+  quadrati(temp,image[indix(I-1,J+2)],image[indix(I+1,J+2)]);
+
+  //5) première ligne:
+  temp[0]=image[indix(I-2,J-2)];
+  temp[1]=image[indix(I-2,J)];
+  temp[2]=image[indix(I-2,J+2)]; //
+  quadrati(temp,image[indix(I-2,J-1)],image[indix(I-2,J+1)]);
+	
+  //6) deuxième ligne:
+  temp[0]=image[indix(I-1,J-2)];
+  temp[1]=image[indix(I-1,J)];
+  temp[2]=image[indix(I-1,J+2)]; //
+  quadrati(temp,image[indix(I-1,J-1)],image[indix(I-1,J+1)]);
+
+  //7) troisième ligne
+  temp[0]=image[indix(I,J-2)];
+  temp[1]=image[indix(I,J)];
+  temp[2]=image[indix(I,J+2)]; //
+  quadrati(temp,image[indix(I,J-1)],image[indix(I,J+1)]);
+	
+  //8) quatrième  ligne:
+  temp[0]=image[indix(I+1,J-2)];
+  temp[1]=image[indix(I+1,J)];
+  temp[2]=image[indix(I+1,J+2)];
+  quadrati(temp,image[indix(I+1,J-1)],image[indix(I+1,J+1)]);
+}
+// Un point sur deux étant calculé  dans image, interpoler qudratiquement
+// les autres points.
+// image : l'image
+// nx, ny: les dimensions de l'image
+template<class real> void interpquad(real* image,int nx,int ny)
+{
+  
+  // tester si nx et ny sont divisibles par 4
+  // à faire!
+  real temp[3];
+  
+  for(int i=4;i<nx;i+=4)
+    for(int j=4;j<ny;j+=4)
+      interp1(temp,image,nx,i,j);
+
+  for(int j=4;j<ny;j+=4)
+      interp1(temp,image,nx,nx-1,j);
+  for(int i=4;i<nx;i+=4)
+    interp1(temp,image,nx,i,ny-1);
+ 
+       
+}
+#endif

main.cpp

@@ -0,0 +1,189 @@
+#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;
+  //
+  typedef typename RBFInterp<real,Gauss<real>>::PointValue PointValue;
+
+  // 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++)
+    {
+      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);
+      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:
+  typedef MQI<real> typRad;
+  //  typedef TPS<real> typRad;
+  //typedef TPSD<real> typRad;
+  //typedef Gauss<real> typRad;
+
+  real eps= 0.01;
+  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");
+  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)
+	{
+	  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;
+  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;
+}
+
+    

protos_lapack.hpp

@@ -0,0 +1,90 @@
+#ifndef protos_lapack__h
+#define  protos_lapack__h
+// prototypes pour lapack et classe d'exception associée.
+namespace lapack_c
+{
+  //! prototypes for lapack routines.
+extern "C"{
+  void dgetrf_(int *n,int *m,double* a,int *lda,int *ipiv,int *info);
+}
+extern "C"{
+  void zgetrf_(int *n,int *m,double* a,int *lda,int *ipiv,int *info);
+}
+extern "C"{
+  void dgetrs_(const char *s,int *N,int *NRHS,double *A,int *LDA,int *IPIV,
+	       double *B,int *LDB,int *INFO );
+
+}
+extern "C"{
+  void zgetrs_(const char *s,int *N,int *NRHS,double *A,int *LDA,int *IPIV,
+	       double *B,int *LDB,int *INFO );
+}
+extern "C"{
+  void dgbtrf_(int *n,int *m,int *k1,int *k2,
+	       double* a,int *lda,int *ipiv,int *info);
+}
+extern "C"{
+  void zgbtrf_(int *n,int *m,int *k1,int *k2,
+	       double* a,int *lda,int *ipiv,int *info);
+}
+extern "C"{
+  void dgbtrs_(const char *s,int *N,int *k1,int *k2,
+	       int *NRHS,double *A,int *LDA,int *IPIV,
+	       double *B,int *LDB,int *INFO );
+}
+extern "C"{
+  void zgbtrs_(const char *s,int *N,int *k1,int *k2,
+	       int *NRHS,double *A,int *LDA,int *IPIV,
+	       double *B,int *LDB,int *INFO );
+}
+extern "C"{
+  void dlarnv_(int *idist,int iseed[],int *n,double *x);
+}
+extern "C"{
+  void dgehrd_(int *n,int *ilo,int *ihi,double *a,int *lda,double tau[],
+	       double work[],int *lwork,int *info);
+}
+extern "C"{
+  void dorghr_(int *n,int *ilo,int *ihi,double *a,int *lda,double tau[],
+	       double work[],int *lwork,int *info);
+}
+  extern "C"{
+    void dgeev_(const char *jobvl,const char *jobvr,int * n,double *a,
+		int *lda,double *WR, double *WI,double *VL,int *LDVL, 
+		double *VR,int  *LDVR,
+		double *WORK,int  *LWORK,int *INFO );
+
+  }
+  extern "C"{
+    void dgesv_(int *n,int *nrhs,double *A,int *lda,int *ipiv, double *B,
+		int *ldb,int *info);
+
+  }
+    extern "C"{
+      void dgbmv_(const char *trans,int *m,int *n,int *kl,int *ku,double *alpha,
+		  double *A,int *lda,double *X,int *incx,double *beta,
+		  double *y, int *incy);
+  }
+  extern "C"{
+    void dgemv_(const char *trans,int *m,int *n,double *alpha, double *A,
+		int *lda,double *X,int *incx,double *beta,
+		double *y, int *incy);
+  }
+   extern "C"{
+     double dlange_(const char *norm,int *m,int *n, double *A,
+		  int *lda,double *work);
+  }
+   extern "C"{
+     void dgecon_(const char *norm,int *n,double *a,int *lda,double *lanorme,
+		  double *rcond,double *work, int *iwork,int *info);
+		    
+  }
+  struct LapackException
+  {
+    int Info;
+    LapackException(int info):Info(info){}
+    int  what(){return Info;}
+  };
+ 
+}
+#endif

time_mes.hpp

@@ -0,0 +1,10 @@
+#ifndef TIME_MES_HPP
+#define TIME_MES_HPP
+#include <sys/time.h>
+// renvoit l'heure en secondes (pécision : micro seconde)
+double get_time() {
+  struct timeval tv;
+  gettimeofday(&tv,0);
+  return static_cast<double>(tv.tv_sec+tv.tv_usec*1e-6);
+}
+#endif