synaps/topology/TopSweep2d.h

/*********************************************************************
*      This file is part of the source code of SYNAPS library        *
*   Author(s): F.Carreras, C.Tanasescu, I. Necula, Univ. Santander   *
*              G. Gatellier, A. Labrouzy,                            *
*              O.Ruatta, B. Mourrain, GALAAD, INRIA                  *
**********************************************************************/
#ifndef synaps_sweep2d_topology_H
#define synaps_sweep2d_topology_H
//--------------------------------------------------------------------
#include <vector>
#include <map>
#include <synaps/init.h>
#include <synaps/base/Seq.h>
#include <synaps/arithm/gmp.h>
#include <synaps/linalg/VectDse.h>
#include <synaps/linalg/MatrDse.h>
#include <synaps/upol/UPolDse.h>
#include <synaps/upol/SturmSeq.h>
#include <synaps/upol/square_free.h>
#include <synaps/mpol/MPol.h>
#include <synaps/mpol/let.h>
#include <synaps/mpol/subs.h>
#include <synaps/usolve/Aberth.h>
#include <synaps/topology/point.h>
#include <synaps/topology/point_graph.h>


//--------------------------------------------------------------------
__BEGIN_NAMESPACE_SYNAPS
//--------------------------------------------------------------------
namespace projection
{
  template<class C>
  struct list_X_Y
  {
    Seq<C>             X;
    Seq<Seq<C> >       Y;
    Seq<int>           I; 
    void clear()
    {
      X.resize(0);
      Y.resize(0);
      I.resize(0);
    }

  };
  //--------------------------------------------------------------------
  template<class C>
  void print(const list_X_Y<C>& L)
  {
    for(unsigned i=0;i<L.X.size();i++)
      std::cout<<L.X[i]<<" ["<<L.Y[i]<<"] "<<L.I[i]<<std::endl;
  }
  //--------------------------------------------------------------------
  template<class C, class POL>
  POL easyquot (const POL& p, C beta, unsigned int n=1)
  {
    //    std::cout << "  easyquot "<<p<<" "<<beta<<" "<<n;
    POL r(p), l(C(1),1);
    l-= POL(beta,0);
    for(unsigned i=0;i<n;++i)
      r/=l;
    //    std::cout << " = "<<r<<std::endl;
    return r;
  }
  //--------------------------------------------------------------------
  template<class C>
  unsigned insert(Seq<C>& L, const C& x)
  {
    unsigned i=0;
    typename Seq<C>::iterator it=L.begin();
    while(it !=L.end() && x> (*it) ) 
      {it++; i++;}
    L.rep().insert(it,x);
    //    std::cout<<"insert "<<x<<" : "<<L<<" : "<<i<<std::endl;
    return i;
  }
  //--------------------------------------------------------------------
  template<class C>
  unsigned insert(list_X_Y<C>& L, const C& x, const Seq<C>& LY, int mu)
  {
    if(LY.size())
      {
        // std::cout <<"insert: "<<x<<" "<<LY<<" "<<mu<<std::endl;
        unsigned i=0;
        typename Seq<C>::iterator itX=L.X.begin();
        typename Seq<Seq<C> >::iterator itY=L.Y.begin();
        typename Seq<int>::iterator itI=L.I.begin();
        while(itX !=L.X.end() && x> (*itX) ) 
          {itX++; itY++; itI++, i++;}
        L.X.rep().insert(itX,x);
        L.Y.rep().insert(itY,LY);
        L.I.rep().insert(itI,mu);
        //      print(L);
      }
  }
  //--------------------------------------------------------------------
  template<class C,class O, class R>
  UPolDse<double> eval_to_double(const MPol<C,O,R> & p, double x)
  {
    UPolDse<double> tmp(degree(p,1)+1, AsSize());
    VECTOR::init(tmp,0);
    double c;
    for(typename MPol<C,O,R>::const_iterator it=p.begin(); it !=p.end();it++)
      {
        unsigned d0=(*it)[0], d1=(*it)[1];
        c =let::convert(it->coeff(),type::As<double>());
        for(unsigned i=0;i<d0;i++) c*=x;
        tmp[d1]+=c;
      }
    return tmp;
  }

  template<class C>
  UPolDse<double> convert_to_double(const UPolDse<C> & P)
  {
    UPolDse<double> tmp(degree(P)+1, AsSize());
    for( unsigned  k=0; k< tmp.size();++k)
      let::assign(tmp[k],P[k]);
    return tmp;
  }

  //--------------------------------------------------------------------
  template<class POL>
  POL change_x_minus_y(const POL& p)
  {
    return subs(0,POL("x0+x1"),p);
  }

  //--------------------------------------------------------------------
  template<class GRAPH>
  void change_x_plus_y(GRAPH& g, int k)
  {
    if(k!=0)
      for(unsigned i=0;i< rep(g.Vertices).size(); i++)
        {
          rep(g.Vertices)[i][0]-=rep(g.Vertices)[i][1]*k;
        }
  }
  //--------------------------------------------------------------------
  template<class GRAPH>
  void swap_x_y(GRAPH& g)
  {
    for(unsigned i=0;i< rep(g.Vertices).size(); i++)
      {
        std::swap(rep(g.Vertices)[i][0],rep(g.Vertices)[i][1]);
      }
  }

  //--------------------------------------------------------------------
}
//====================================================================
template<class O,  class SLV=SlvAberth<RR>, class X=double>
struct TopSweep2d
{
  typedef O solution_t;
  X eps;
  unsigned nbpt;
  TopSweep2d():nbpt(1),eps(1.e-4) {}
  TopSweep2d(const X& e):nbpt(1),eps(e) {}
  TopSweep2d(int n):nbpt(n),eps(1.e-4) {}
  TopSweep2d(const TopSweep2d & mth):nbpt(mth.nbpt),eps(mth.eps) {}

  typedef topology::point<double> point_t;
  typedef MatrDse<double>         m_solution_t;
  typedef Seq<VectDse<double> >   b_solution_t;
  typedef MPol<QQ>                equation_t;
  
  typedef Seq<RR>                 u_solution_t;
      
  template<class POL>
  u_solution_t u_solve(const POL &p)
  {
    return solve(p,SLV());
  }

  //--------------------------------------------------------------------
  template<class POL, class C>
  bool IsGeneric1(const POL& p, const C& alpha, const C& beta)
  {
    std::vector<C> x(2); x[0]=alpha;x[1]=beta;
    C v= MPOLDST::eval<C>(p,x);
    return (std::abs(v)< this->eps);
  }

  //--------------------------------------------------------------------
  template<class POL, class C>
  bool IsGeneric2(const POL& p, const C& alpha, const C& beta, int  k)
  {
    return true;
  }

  //--------------------------------------------------------------------
  template<class C, class POL>
  void add_regular_y(Seq<C>& S,const POL& p,const C& alpha)
  {
    UPolDse<C> pol= projection::eval_to_double(p,alpha);
    u_solution_t y_sol = this->u_solve(pol);
    C beta;
    for(unsigned j=0;j<y_sol.size();j++)
      {
        let::assign(beta,y_sol[j]);
        S.push_back(beta);
      }
    std::sort(S.begin(),S.end());
  }
  //--------------------------------------------------------------------
  template<class C, class POL>
  void add_regular_point(projection::list_X_Y<C> & L, 
                         const projection::list_X_Y<C> & L0, 
                         const POL& p,
                         unsigned N=1)
  {
    typedef C coeff_t;
    using let::convert;
    if(L0.X.size() && N)
      {
        coeff_t alpha,beta, dx;
        Seq<double>  SY;

        alpha=L0.X[0]-1;

        this->add_regular_y(SY,p,alpha);
        if(SY.size())
          {
            L.X.push_back(alpha);
            L.Y.push_back(SY);
            L.I.push_back(-1);
          }

        dx=coeff_t(1)/(N+1);
        for(unsigned l=1;l<N+1; l++)
        {
          alpha+=dx;
          SY.resize(0);
          this->add_regular_y(SY,p,alpha);
          if(SY.size())
          {
            L.X.push_back(alpha);
            L.Y.push_back(SY);
            L.I.push_back(-1);
          }
        }
        for(unsigned i=0;i<L0.X.size();i++)
          {
            L.X.push_back(L0.X[i]);
            L.I.push_back(L0.I[i]);
            L.Y.push_back(L0.Y[i]);
            if(i<int(L0.X.size()-1))
              dx=(L0.X[i+1]-L0.X[i])/(N+1);
            else 
              dx = C(1)/(N+1);
            alpha= L0.X[i]+dx;
            for(unsigned l=1;l<N+1; l++)
              {
                SY.resize(0);
                this->add_regular_y(SY,p,alpha);
                if(SY.size())
                  {
                    L.X.push_back(alpha);
                    L.Y.push_back(SY);
                    L.I.push_back(-1);
                  }
                alpha+=dx;
              }
          }
      }
    else
      {
        C alpha=-1;
        Seq<C> SY;
        L.X.push_back(alpha);
        L.I.push_back(-1);
        this->add_regular_y(SY,p,alpha);
        L.Y.push_back(SY);

        alpha=1;
        L.X.push_back(alpha);
        L.I.push_back(-1);
        SY.resize(0);
        this->add_regular_y(SY,p,alpha);
        L.Y.push_back(SY);
      }
  }


  
template<class GAMMAS, class STURMSEQ>
bool isgeneric(  GAMMAS & Gammas,  STURMSEQ & SH )
{
  typedef typename GAMMAS::value_type   upol_t;
  typedef typename STURMSEQ::value_type bpol_t;
  //  DBG_PRINT(SH);
  for ( int i = 0; i < Gammas.size(); i ++ )
    {
      bpol_t SHi =  SH[SH.size()-i-2];
      upol_t betadenum =  SHi[i+1]*(i+1);
      upol_t betanum   = -SHi[i];
      if ( betadenum != 0 )
        {
          Rational<upol_t> f(betanum,betadenum);
          Rational<upol_t> r = SHi(f);
          upol_t nm = numerator(r);
          if ( nm % Gammas[i] != 0 )  return false;
          for ( int k = 0; k < i; k++ )
            {
              SHi = diff(SHi);
              r  = SHi(f);
              nm = numerator(r);
              if ( nm % Gammas[i] != 0 ) return false;
            };
        };
    };
  return true;
};


  //--------------------------------------------------------------------
  template<class C, class POL>
  bool insert_char_points(projection::list_X_Y<C> & L, 
                          const POL & p1, const POL& p2,
                          bool both=false)
  {
    typedef POL mpol_t;
    typedef ZZ coeff_t;
    typedef UPolDse<coeff_t> upol_t;
    typedef UPolDse<UPolDse<coeff_t> > bpol_t;
    typedef std::vector<upol_t> vector_t;

    using let::assign;
    using let::convert;

    bpol_t p; assign(p,p1,1);
    bpol_t q; assign(q,p2,1);
    
    std::cout<<"p1: "<<p<<std::endl; 
    std::cout<<"p2: "<<q<<std::endl; 
    bpol_t::var_= "y";
    upol_t qq=p[degree(p)];
    if( degree(qq) !=0) return false;
    SturmSeq<upol_t> s(p,q,HABICHT());
    //    DBG_PRINT("Sturm Habicht Sequence"<<s);
    vector_t r=CoeffStHa<vector_t>(s);
    upol_t   disc=square_free(r[r.size()-1]);
    //    DBG_VAL(disc);
    vector_t Gammas;
    upol_t   Fi=disc,Fin;

    for (unsigned j=0; j<r.size()-1; ++j)
      {
        Fin= gcd(Fi,r[r.size()-j-2]);
        if(Fin !=0) Gammas.push_back(Fi/Fin);
        Fi=Fin;
      }

    if ( ! isgeneric( Gammas, s ) ) return false;

    //    DBG_PRINT("Decomposition");

    Seq<Seq<C> > Xmu, Ycrit, Sol;
    Seq<int> NbXmu;
    Seq<C> solx, soly;

    std::map<double,int> mu;
    for (unsigned j=0; j<Gammas.size(); ++j)
      {
        if (degree(Gammas[j])>=1)
          {
            UPolDse<coeff_t> tmp =Gammas[j];
            std::cout<<"Solving polynomial "<<tmp <<std::endl;
            u_solution_t SolX = u_solve(tmp);
            double r;
            for(unsigned k=0; k< SolX.size();++k)
              {
                assign(r,SolX[k]);
                mu[r]=(j+1);
                solx.push_back(r);
                // std::cout<<"x_sol: "<<r<<" mu: "<<mu[r]<<std::endl;
              }
          }
        else
          {
            Seq<double> tmpEmpty; NbXmu.push_back(0);
            Xmu.push_back(tmpEmpty);
            Ycrit.push_back(tmpEmpty);
          }
      }
    sort(solx.begin(),solx.end());
    // std::cout<<"x-crit: "<<solx<<std::endl;

    //    DBG_PRINT("Nb of critical values: ");

    for (unsigned l=0; l<solx.size(); ++l)
    {
      Seq<double> Yreg;
      double alpha=solx[l], beta;
      L.X.push_back(alpha);
      // std::cout<<"alpha: "<<alpha<<std::endl;
      int m=mu[alpha];
      if(m !=0)
        {
          // std::cout << "  x-coordinate: " << alpha<<" \n  mu: "<<m<< std::endl;

          bpol_t tp = s[s.size()-m-1];
          beta = -projection::convert_to_double(tp[m-1])(alpha)/
            projection::convert_to_double(tp[m])(alpha)/(m);
          
          unsigned l=m+1;
          if(both) l=m;
          UPolDse<double> tmp1 = 
            projection::easyquot( projection::eval_to_double(p1,alpha),beta,l);
          //   Yreg.push_back(Xmu[j][i]);
          if (degree(tmp1)>0)
            {
              std::cout<<"Solving polynomial(1) "<<tmp1<<std::endl;
              u_solution_t y_sol = u_solve(tmp1);
              //  std::cout<<"  y_sol: "<<y_sol<<std::endl;
              for(unsigned k=0; k<y_sol.size();++k)
                Yreg.push_back(convert(y_sol[k],type::As<double>()));
            }
          
          if(both)
            {
              std::cout<<"Solving polynomial(2) "<<tmp1<<std::endl;
              tmp1 = projection::easyquot(projection::eval_to_double(p2,alpha),
                                          beta,m);
              if (degree(tmp1)>0)
                {
                  u_solution_t y_sol = u_solve(tmp1);
                  for(unsigned k=0; k<y_sol.size();++k)
                    Yreg.push_back(convert(y_sol[k],type::As<double>()));
                }
              
            }
          sort(Yreg.begin(),Yreg.end());
          L.I.push_back(projection::insert(Yreg,beta));
        }
      else
        {
          // std::cout << "  x-reg: " << alpha << std::endl;
          UPolDse<double> tmp1= projection::eval_to_double(p1,alpha);
          //      std::cout <<"pol: "<<tmp1<<std::endl;
          u_solution_t y_sol = u_solve(tmp1);
          //      std::cout <<"sol: "<<y_sol<<std::endl;
          for(unsigned k=0; k<y_sol.size();++k)
            Yreg.push_back(convert(y_sol[k],type::As<double>()));

          if(both)
            {
              tmp1= projection::eval_to_double(p1,alpha);
              //  std::cout <<"pol: "<<tmp1<<std::endl;
              y_sol = u_solve(tmp1);
              //  std::cout <<"sol: "<<y_sol<<std::endl;
              for(unsigned k=0; k<y_sol.size();++k)
                Yreg.push_back(convert(y_sol[k],type::As<double>()));
            }
          sort(Yreg.begin(),Yreg.end());
          //      std::cout << "  y-reg: " << Yreg << std::endl;
          L.I.push_back(-1);
        }
      //      sort(Yreg.begin(),Yreg.end()); // Could be removed.
      L.Y.push_back(Yreg);
      // cout <<"  Y Coordinates of Points above "<<alpha<<" : " << Yreg<< endl;
    }
    return true;
  }

  //--------------------------------------------------------------------
  template<class GRAPH, class C>
  void connect(GRAPH & g, const projection::list_X_Y<C>& L)
  {
    int n0,n1,s0,s1;
    n0=g.Vertices.size();
    //    std::cout<<"n0= "<<n0<<std::endl;
    topology::point<double> p(2);
    s0=L.Y[0].size();
    //    std::cout<<"L.Y[0].size()= "<<L.Y[0].size()<<std::endl;
    p[0]=L.X[0];
    for (int j=0;j<s0;j++)
      {
        p[1]=L.Y[0][j];
        g.push_back(p);
      }

    //    std::cout<<"L.X.size()= "<<L.X.size()<<" "<<s0<<std::endl;
    for (unsigned i=1; i<L.X.size();i++)
    {
      //      std::cout<<">> "<<L.X[i-1]<<":::"<<L.X[i]<<std::endl;
      p[0]=L.X[i];
      n1=n0+s0;
      s1=L.Y[i].size();
      for (int j=0;j<s1;j++)
      {
        //      std::cout<<"L.Y[i].size()= "<<L.Y[i].size()<<std::endl;
        p[1]=L.Y[i][j];
        g.push_back(p);
        //        std::cout<<"p[1]= "<<p[1]<<std::endl;
      }
      //      std::cout<<n0<<" "<<s0<<" "<<n1<<" "<<s1<<" : "
      //               <<L.Y[i-1]<<" ; "<<L.Y[i]<<std::endl;
      int ku0,ku1,kd0,kd1;
      if (L.I[i-1]>-1)
        {
          ku0=L.I[i-1];
          kd0=ku0+1;
        }
      else
        {
          ku0=s0;
          kd0=0;
        }
      if (L.I[i]>-1)
        {
          ku1=L.I[i];
          kd1=ku1+1;
        }
      else
        {
          ku1=s1;
          kd1=0;
        }

      //   std::cout<<p[0]<<" : "<<ku0<<" "<<kd0<<" "<<ku1<<" "<<kd1<<std::endl;
      int k0,k1;

     for (k0=0;k0<ku0 && k0<ku1;k0++)
     {
       g.insert(n0+k0,n1+k0);
       std::cout<<"ru> "<<n0+k0<<" - "<<n1+k0<<" "
                <<g.Vertices[n0+k0]<<" "
                <<g.Vertices[n1+k0]
                << std::endl;
     }
     for (k1=0;k1<s0-kd0 && k1<s1-kd1;++k1)
     {
       g.insert(n0+s0-k1-1,n1+s1-k1-1);
       std::cout<<"rd> "<<n0+s0-k1-1<<" - "<<n1+s1-k1-1<<" "
                <<g.Vertices[n0+s0-k1-1]<<" "
                <<g.Vertices[n1+s1-k1-1]
                << std::endl;
     }
     // std::cout<<"n0: "<<n0<<" n1: "<<n1<<" k0: "<<k0<<" k1:"<<k1<<std::endl;
     for (int j0=k0;j0<s0-k1;j0++)
     {
       for (int j1=k0;j1<s1-k1;j1++)
       {
         g.insert(n0+j0,n1+j1);
         //      std::cout<<"Smile_inside"<<" > "<<n0+j0<<" - "<<n1+j1<<"> "
         //               <<g.Vertices[n0+j0]<<" "
         //               <<g.Vertices[n1+j1]
         //               << std::endl;
       }

     }
      n0=n1; s0=s1;
    }
  }

};

namespace topology 
{
  template<class C,class U>
  void assign(point_graph<C>& g, 
              const MPol<U>& eq, 
              TopSweep2d<C> mth)
  {
    typedef C        coeff_t;
    typedef MPol<U>  mpol_t; 

    MPol<QQ> p; let::assign(p,eq);

    projection::list_X_Y<coeff_t> L,Lr;

    int nvc=0, d= degree(p);
    bool sw=false;
    if(!mth.insert_char_points(L,p,diff(p,1)))
      {
        p = MPOLDST::swap(p,0,1);
        sw=true; 
        L.clear();
        std::cout<<"Not generic "<<nvc<<std::endl;
        while (!mth.insert_char_points(L,p,diff(p,1)) && nvc < (d*(d-1)))
          {
            p=projection::change_x_minus_y(p);
            nvc++;
            L.clear();
            std::cout<<"Not generic "<<nvc<<std::endl;
          }
      }

    print(L);
    mth.add_regular_point(Lr,L,p,mth.nbpt);
    //    std::cout<<"----------------------------------------------"<<std::endl;
    print(Lr);
    mth.connect(g,Lr);

    std::cout<<nvc<<std::endl;
    if(nvc >0)
      projection::change_x_plus_y(g,-nvc);
    if(sw)
      projection::swap_x_y(g);
    std::cout<<g.Vertices.size()<<std::endl;
  }

  template<class C,class U>
  void assign(point_graph<C>& g, 
              const MPol<U>& P1, 
              const MPol<U>& P2, 
              TopSweep2d<C> mth)
  {
    typedef  TopSweep2d<C> mth_t;
    typedef typename TopSweep2d<C>::u_solution_t u_solution_t;
    
    MPol<QQ> p1; let::assign(p1,P1);
    MPol<QQ> p2; let::assign(p2,P2);

    projection::list_X_Y<C> L1,L2,LI,LIr;
    mth.insert_char_points(L1,p1,diff(p1,1));
    mth.insert_char_points(L2,p2,diff(p2,1));

    std::cout<<"====================================="<<std::endl;
    std::cout<<"C1"<<std::endl;
    print(L1);
    std::cout<<"====================================="<<std::endl;
    std::cout<<"C2"<<std::endl;
    print(L2);


    double alpha,beta;
    int mu;
    for(unsigned i=0; i< L2.X.size();i++)
      {
        alpha=L2.X[i];
        UPolDse<double> tmp1= projection::eval_to_double(p1,alpha);
        //      std::cout <<tmp1<<std::endl;
        u_solution_t y_sol = mth.u_solve(tmp1);
        //      std::cout <<"y_sol: "<<y_sol<<std::endl;
        for(unsigned j=0;j<y_sol.size();j++)
          {
            let::assign(beta,y_sol[j]);
            mu=projection::insert(L2.Y[i],beta);
            if(mu<=L2.I[i]) L2.I[i]++;
            //std::cout <<"insert "<<y_sol[j]<<" in |"<<L2.Y[i]<< " |"<<L2.I[i]<<std::endl;
          }
      }
    for(unsigned i=0; i< L1.X.size();i++)
      {
        alpha=L1.X[i];
        UPolDse<double> tmp1= projection::eval_to_double(p2,alpha);
        //      std::cout <<"x: "<< alpha<<"\npol: "<<tmp1<<std::endl;
        u_solution_t y_sol = mth.u_solve(tmp1);
        //      std::cout <<"y_sol: "<<y_sol<<std::endl;
        for(unsigned j=0;j<y_sol.size();j++)
          {
            let::assign(beta,y_sol[j]);
            mu=projection::insert(L1.Y[i],beta);
            if(mu<=L1.I[i]) L1.I[i]++;
          }
      }

    std::cout<<"====================================="<<std::endl;
    std::cout<<"C1"<<std::endl;
    print(L1);
    std::cout<<"====================================="<<std::endl;
    std::cout<<"C2"<<std::endl;
    print(L2);
    std::cout<<"====================================="<<std::endl;
    for(unsigned i=0; i< L2.X.size();i++)
      projection::insert(L1,L2.X[i],L2.Y[i],L2.I[i]);
    print(L1);
    std::cout<<"====================================="<<std::endl;
    mth.insert_char_points(LI,p1,p2,true); 
    print(LI);
    std::cout<<">>==================================="<<std::endl;
    for(unsigned i=0; i< LI.X.size();i++)
      {
        projection::insert(L1,LI.X[i],LI.Y[i],LI.I[i]);
      }

    for(int i=0;i<L1.X.size()-1;i++)
      {
        if(L1.I[i]>=0 && L1.I[i+1]>=0)
          {
            alpha= (L1.X[i]+L1.X[i+1])/2;
            Seq<double> ly;
            UPolDse<double> tmp1= projection::eval_to_double(p1,alpha);
            u_solution_t y_sol = mth.u_solve(tmp1);
            
            for(unsigned j=0; j<y_sol.size();j++)
              {
                let::assign(beta,y_sol[j]);
                ly.push_back(beta);
              }
            tmp1= projection::eval_to_double(p2,alpha);
            y_sol = mth.u_solve(tmp1);
            for(unsigned j=0; j<y_sol.size();j++)
              {
                let::assign(beta,y_sol[j]);
                ly.push_back(beta);
              }
            sort(ly.begin(),ly.end());
            projection::insert(LIr,alpha,ly,-1);
          }
      }
    std::cout<<"<<==================================="<<std::endl;
    for(unsigned i=0; i< LIr.X.size();i++)
      {
        projection::insert(L1,LIr.X[i],LIr.Y[i],LIr.I[i]);
      }
    std::cout<<"<<==================================="<<std::endl;
    print(L1);
    std::cout<<">>==================================="<<std::endl;
    mth.connect(g,L1);
  }
}

//--------------------------------------------------------------------
__END_NAMESPACE_SYNAPS
/********************************************************************/
#endif //synaps_sweep2d_topology_H