data.cc

#include <iostream>
#include <fstream>
#include <set>
#include <functional>
#include <vector>
#include <algorithm>
#include "functions.h"
#include "data.h"




// ------- class vertices ----------------------------

int vertices::n=0;

int vertices::d=0;

vertices_lookup vertices::lookup=vertices_lookup();


void vertices::set_dimensions(int N, int D)
{
  n=N; d=D;
  lookup.generate_tables(n,d);
}


vertices::vertices()
{
}


vertices::vertices(const simplex& S)
{
  simplex_to_vertices(S);
}


const simplex vertices::vertices_to_simplex() const
{
  simplex s=1;
  vertex k=1, l;
  const_iterator l_iterator=begin();
  for (vertex i=1; i<=d; ++i) {
    l = (*l_iterator)+1;
    for (vertex j=k; j<=l-1; ++j)
      s=s+lookup.get_binomial(n-j,d-i);
    k=l+1;
    ++l_iterator;
  }
  return(s);
}


const vertices & vertices::simplex_to_vertices(const simplex & S)
{
  (*this) = lookup.simplex_to_vertices(S);
  return(*this);
}


std::ostream & operator << (std::ostream & out, const vertices & v)
{
  out << *( v.begin() );
  vertices::const_iterator i=v.begin();   i++;
  for (; i!=v.end(); ++i)
    out << "_" << *i;
  return(out);
}


bool vertices_order::operator()(const vertices & a, const vertices & b) const
{
  size_t n1=a.size(), n2=b.size();
  if (n1<n2) return(true);
  if (n1>n2) return(false);
  // a,b have same size
  for (vertices::iterator i=a.begin(), j=b.begin();
       i!=a.end() && j!=b.end(); ++i, ++j) {
    if (*i < *j) return(true);
    if (*i > *j) return(false);
  }
  // a == b
  return(false);
}


bool operator==(const vertices & a, const vertices & b)
{
  return( std::set<vertex,std::less<vertex> >(a) ==
          std::set<vertex,std::less<vertex> >(b) );
}


// ------- class vertices_lookup --------------------

const vertices & vertices_lookup::simplex_to_vertices(const simplex & S) const
{
  return( SimplexToVertices[S-1] );
}


// this one is messy, but it only reverses vertices_to_simplex()
vertices vertices_lookup::manual_vertices_to_simplex(const simplex & S) const
{
  vertices result;
  simplex s=S;
  simplex b;
  vertex i,j,l=0,k;
  for (k=1; k<d; k++) {
    l++;  i=l;   j=1;
    b=binomial(n-l,d-k);
    while (s>b && b>0) {
      j++;  l++;
      s=s-b;
      b=binomial(n-l,d-k);
    }
    result.insert(result.begin(),l-1);
  }
  result.insert(result.begin(),s+l-1);
  return(result);
}

// tweaked for special (e.g. negative) values
int vertices_lookup::get_binomial(int i, int j) const
{
  if (i>=0 && i<=n && j>=0 && j<=d && j<=i)
    return(fast_binomial[i][j]);
  else
    return(1);
}


void vertices_lookup::generate_tables(int N, int D)
{
  n=N;  d=D;
  fast_binomial.erase(fast_binomial.begin(),fast_binomial.end());
  for(int i=0; i<=n; ++i) {
    std::vector<int> row;
    for(int j=0; j<=i && j<=d; j++)
      row.push_back( binomial(i,j) );
    fast_binomial.push_back( row );
  }
  SimplexToVertices.erase(SimplexToVertices.begin(),SimplexToVertices.end());
  for(int s=1; s<=binomial(n,d); ++s)
    SimplexToVertices.push_back(manual_vertices_to_simplex(s));
}




// ------- class compact_simplices -------------------

compact_simplices::compact_simplices()
{
  reserve(50);
}

compact_simplices::~compact_simplices()
{
}

std::ostream & operator << (std::ostream & out, const compact_simplices & t)
{
  out << "[" << *( t.begin() );
  for (compact_simplices::const_iterator i=t.begin()+1; i!=t.end(); ++i)
    out << "," << *i;
  out << "]";
  return(out);
}

const bool operator==(const compact_simplices & a, const compact_simplices & b)
{
  // friendship is not inherited
  return( std::vector<simplex>(a) == std::vector<simplex>(b)  );
}

hash_value compact_simplices::hash_function() const
{
  hash_value result=0;
  int n=101;
  for (const_iterator i=begin(); i!=end(); ++i, n+=37)
    result += n* (*i) + n*n;
  return(result);
}





// ------- class simplices ---------------------------

simplices::simplices(const compact_simplices & s)
  :compact_simplices(s)
{
  decompress();    // de-compress vertex information
}

simplices::simplices(const std::set<vertices,vertices_order> & s)
{
  v.erase(v.begin(), v.end());
  copy(s.begin(), s.end(), back_inserter(v) );
  compress();
}


simplices::~simplices()
{
}

void simplices::compress()
{
  erase(begin(),end());
  for (std::vector<vertices>::const_iterator i=v.begin(); i!=v.end(); i++)
    push_back( (*i).vertices_to_simplex() );
  std::sort(begin(),end());
}


void simplices::decompress()
{
  v.erase(v.begin(),v.end());
  for (const_iterator i=begin(); i!=end(); i++)
    v.push_back( vertices().simplex_to_vertices(*i) );
}


bool simplices::starshaped(const vertex origin) const
{
  bool result = true;
  for (std::vector<vertices>::const_iterator
         vi=v.begin(); vi!=v.end(); vi++) {
    result=result && (find(vi->begin(),vi->end(),origin)!=vi->end());
  }
  return(result);
}


bool simplices::fine() const
{
  vertices support;
  for (std::vector<vertices>::const_iterator
         vi=v.begin(); vi!=v.end(); vi++) {
    support.insert(vi->begin(), vi->end());
  }
  return support.full_set();
}


std::ostream & operator << (std::ostream & out, const simplices & s)
{
  out << "[" << s.v.front();
  for (std::vector<vertices>::const_iterator
         si=s.v.begin()+1; si!=s.v.end(); si++)
    out << ", " << *si;
  out << "]";
  return(out);
}

// ------- class flip --------------------------------

flip::flip()
{
  deltaplus.reserve(10);
  deltaminus.reserve(10);
}

flip::flip(const std::vector<vertices>& pos, const std::vector<vertices>& neg)
  : deltaplus(pos), deltaminus(neg)
{}

flip::flip(const flip & copyfrom)
{
  deltaplus.reserve(10);
  deltaminus.reserve(10);
  (*this)=copyfrom;
}

flip::~flip()
{
}

flip & flip::operator =(const flip & copyfrom)
{
  deltaplus  = copyfrom.get_deltaplus();
  deltaminus = copyfrom.get_deltaminus();
  return(*this);
}

std::ostream & operator << (std::ostream & out, const flip & f)
{
  out << "< ";
  for (std::vector<vertices>::const_iterator i=f.get_deltaplus().begin();
       i!=f.get_deltaplus().end(); ++i) {
    std::cout << *i << " ";
  }
  out << "| ";
  for (std::vector<vertices>::const_iterator i=f.get_deltaminus().begin();
       i!=f.get_deltaminus().end(); ++i) {
    std::cout << *i << " ";
  }
  out << ">";
  return(out);
}

void flip::mirror()
{
  swap(deltaplus,deltaminus);
}




// ------- class flips -------------------------------

flips::flips()
{
}

flips::~flips()
{
}


// ------- class goodcircuit -------------------------


goodcircuit::goodcircuit(const simplices & s, const flip & f)
  :supporter(f)
{
  const std::vector<vertices> & deltaplus = f.get_deltaplus();
  const std::vector<vertices> & v = s.get_vertices();
  supported.reserve(deltaplus.size());
  good=true;


  // find simplices that include members of deltaplus
  std::vector<vertices> empty_result;
  for (size_t n=0; n<deltaplus.size(); n++) {
    supported.push_back(empty_result);  supported[n].reserve(10);
    bool found_match=false;
    for (std::vector<vertices>::const_iterator i=v.begin(); i!=v.end(); ++i) {
      if (includes( (*i).begin(), (*i).end(),
                    deltaplus[n].begin(), deltaplus[n].end() )) {
        found_match=true;
        supported[n].push_back(*i);
      }
    }
    good=good && found_match;
    if (!good) return;
  }

  // check that all links are equal
  for (size_t n=0; n<deltaplus.size(); n++) {
    std::set<vertices,vertices_order> l;
    for (std::vector<vertices>::const_iterator i=supported[n].begin();
             i!=supported[n].end(); ++i) {
      vertices one_simplex;
      std::insert_iterator<vertices>
        ins(one_simplex,one_simplex.begin());
      set_difference( (*i).begin(), (*i).end(),
                       deltaplus[n].begin(), deltaplus[n].end(),
                       ins );
      l.insert(l.begin(),one_simplex);
    }
    link.push_back(l);
    if (n>0 && link[0]!=link[n]) {
      good=false;
      return;
    }
  }
}


void goodcircuit::do_flip(const simplices & s, const flip & f)
{
  bistellarneighbor.erase(bistellarneighbor.begin(), bistellarneighbor.end());

  const std::set<vertices,vertices_order> & use_link = link[0];
  const std::vector<vertices> & v = s.get_vertices();
  const std::vector<vertices> & deltaplus  = f.get_deltaplus();
  const std::vector<vertices> & deltaminus = f.get_deltaminus();

  // start with all simplices
  std::set<vertices,vertices_order> tri;
  copy(v.begin(), v.end(), inserter(tri,tri.begin()) );

  // remove all suported simplices
  std::set<vertices,vertices_order> to_remove;

  for (std::set<vertices,vertices_order>::const_iterator i=use_link.begin();
       i!=use_link.end(); ++i)
    for (std::vector<vertices>::const_iterator j=deltaplus.begin();
         j!=deltaplus.end(); ++j) {
      vertices one_simplex;
      set_union( (*j).begin(), (*j).end(),
                 (*i).begin(), (*i).end(),
                 inserter(one_simplex,one_simplex.begin()) );
      to_remove.insert(to_remove.begin(),one_simplex);
    }
  set_difference( tri.begin(), tri.end(),
                  to_remove.begin(), to_remove.end(),
                  inserter(bistellarneighbor,bistellarneighbor.begin()),
                  vertices_order() );

  // now add the flip of the supported simplices
  for (std::set<vertices,vertices_order>::const_iterator i=use_link.begin();
       i!=use_link.end(); ++i)
    for (std::vector<vertices>::const_iterator j=deltaminus.begin();
         j!=deltaminus.end(); j++) {
      vertices one_simplex;
      set_union( (*j).begin(), (*j).end(),
                 (*i).begin(), (*i).end(),
                 inserter(one_simplex,one_simplex.begin()) );
      bistellarneighbor.insert( bistellarneighbor.begin(), one_simplex );
    }
}