self_avoiding_polygon/approx/polygon.cpp

#include "polygon.hpp"


//---------------
// Polygon::init
//---------------

void
Polygon::init(){
  for(size_t i=0;i<grid_size;++i){
    grid[i]=empty;
    grid_histo[i]=0;
  }
  for(size_t i=0;i<grid_width;++i){
    grid[pos(i,0)]=border;
    grid[pos(i,grid_height-1)]=border;
  }
  for(size_t i=0;i<grid_height;++i){
    grid[pos(0,i)]=border;
    grid[pos(grid_width-1,i)]=border;
  }
  for(size_t i=0;i<xinit;++i){
    grid[pos(i,1)]=border;
  }
  head=pos(xinit,yinit);
  grid[head]=1;
}

//--------------------
// Polygon::operator=
//--------------------

Polygon&
Polygon::operator=(const Polygon& P){
  head=P.head;
  length=P.length;
  memcpy(grid,P.grid,grid_size*sizeof(size_t));
  memcpy(grid_histo,P.grid_histo,grid_size*sizeof(size_t));
  memcpy(step_histo,P.step_histo,(max_len+1)*sizeof(size_t));
  return *this;
}

//------------------
// Polygon::display
//------------------

void
Polygon::display() const{
    for(size_t z=0;z<grid_height;++z){
    size_t y=grid_height-z-1;
    for(size_t x=0;x<grid_width;++x){
      size_t v=grid[pos(x,y)];
      if(v==border) cout<<"\033[47m \033[0m";
      else if(v==empty or step_histo[v]!=grid_histo[pos(x,y)] or v>length) cout<<" ";
      else cout<<"\033[42m\033[30m"<<" "<<"\033[0m";
    }
    cout<<endl;
  }
}

//-------------
// Polygon::fp
//-------------

Reel
Polygon::fp() const{
  map<size_t,InfoVertex> infos;
  size_t id=0;
  int64 x=xinit;
  int64 y=yinit;
  size_t c=pos(x,y);
  for(size_t t=1;t<=length;++t){
    size_t l=left(c);
    size_t r=right(c);
    size_t u=up(c);
    size_t d=down(c);
    auto it=infos.find(c);
    if(it==infos.end()){
      infos[c]=InfoVertex(id++,x,y,l,r,u,d);
    }
    else{
      InfoVertex& info=it->second;
      info.v[0]=l;
      info.v[1]=r;
      info.v[2]=u;
      info.v[3]=d;
      info.deg=4;
    }
    it=infos.find(l);
    if(it==infos.end()) infos[l]=InfoVertex(id++,x-1,y);
    it=infos.find(r);
    if(it==infos.end()) infos[r]=InfoVertex(id++,x+1,y);
    it=infos.find(u);
    if(it==infos.end()) infos[u]=InfoVertex(id++,x,y+1);
    it=infos.find(d);
    if(it==infos.end()) infos[d]=InfoVertex(id++,x,y-1);
    size_t s=step_histo[t+1];
    if(grid[l]==t+1 and s==grid_histo[l]){
      c=l;
      --x;
    }
    else if(grid[r]==t+1 and s==grid_histo[r]){
      c=r;
      ++x;
    }
    else if(grid[u]==t+1 and s==grid_histo[u]){
      c=u;
      ++y;
    }
    else if(grid[d]==t+1 and s==grid_histo[d]){
      c=d;
      --y;
    }
  }
  size_t ne=infos.size();
  AvxMatrix AvxB;
  AvxMatrix AvxC;
  AvxB.clear();
  AvxC.clear();
  for(auto it=infos.begin();it!=infos.end();++it){
    InfoVertex& info=it->second;
    size_t i=info.id;
    for(size_t k=0;k<info.deg;++k){
      size_t j=infos[info.v[k]].id;
      AvxB.get(i,j)=1;
      AvxB.get(j,i)=1;
    }
    for(auto it2=infos.begin();it2!=infos.end();++it2){
      InfoVertex& info2=it2->second;
      size_t j=info2.id;
      int64 dx=abs(info.x-info2.x);
      int64 dy=abs(info.y-info2.y);
      AvxC.get(i,j)=get_coeff(dx,dy);
    }
  }
  AvxMatrix AvxM;
  AvxM.clear();
  AvxM.from_C_B(AvxC,AvxB,ne);
 
  Reel det=AvxM.Gauss(ne,ne+1);

  Reel res=0;
  for(size_t i=0;i<ne;++i){
    res+=AvxB.get_diag_square_sym(i,ne)*AvxM.get(i,ne);
  }
  res*=(det*0.25);
  return res;
}

//----------------------
// Polygon::next_vertex
//----------------------

size_t
Polygon::next_vertex(size_t t,size_t c) const{
  size_t l=left(c);
  size_t r=right(c);
  size_t u=up(c);
  size_t d=down(c);
  size_t s=step_histo[t+1];
  if(grid[l]==t+1 and s==grid_histo[l]) return l;
  if(grid[r]==t+1 and s==grid_histo[r]) return r;
  if(grid[u]==t+1 and s==grid_histo[u]) return u;
  if(grid[d]==t+1 and s==grid_histo[d]) return d;
  assert(false);
  return 0;
}