int main(unsigned int argc, char **argv){
  vector<unsigned int> partlist;
  cdt_skeleton cdt_skel;
  vector<cdt_skeleton> cdt_e;
  vector<unsigned int> *dimer_count;
  Triangulation TestTri;
  unsigned int size = 160;
  unsigned int dimer_size = 4;

  if(argc > 1) size = atoi(argv[1]);

  if(argc > 2) dimer_size = atoi(argv[2]); 
  //printf("No. of Dimers  |   No. of Configurations\n", size, dimer_size);

  for(unsigned int j = 2; j <= size; j = j + 2){ 
    create_CDT_size_n(&cdt_e, j);
    //dimer_count = new vector<unsigned int>((unsigned int)(j/2) + 1);
    for(unsigned int i = 0; i < cdt_e.size(); i++){
      TestTri.Clear();
      //TestTri.Create(cdt_e[i]);
      //TestTri.GenerateDimerConfigs(dimer_size, dimer_count);
    }
    cdt_e.clear();
    printf("%i\n", j);
    //printdimers(dimer_count);
    //delete dimer_count;
  }
  return 0;
}

void test_case(int n) {

	// cout << "---- " << n << " ----" << endl;

	// Read all infected people
	std::vector<K::Point_2> infected;
	infected.reserve(n);

	for(int i=0; i<n; i++) {
		// cin >> x[i] >> y[i];
		K::Point_2 p;
		cin >> p;
		infected.push_back(p);
		cout << "Read in point " << p << endl;
	}

	// Construct Delauney triangulation
	Triangulation t;
	t.insert(infected.begin(), infected.end());

	// Read all healthy people
	int m;
	cin >> m;

	for(int i=0; i<m; i++) {
		K::Point_2 escaper;
		long d;
		cin >> escaper >> d;

		// --- Find an escape path for this person ---
		// Find out at which face we are
		Face_handle current_face = t.locate(escaper);

		// Check if we are already outside
		if(t.is_infinite(current_face)) {
			cout << "y";
			continue;
		}
		// Check if we are already getting infected
		/*K::Point_2 nearest_infected = t.nearest_vertex(escaper, current_face)->point();
		cout << "Nearest infected person: " << nearest_infected << endl;
		int dx = nearest_infected.x() - escaper.x();
		int dy = nearest_infected.y() - escaper.y();
		long nearest_sqd = dx * dx + dy * dy;
		if(nearest_sqd < d) {
			cout << "n";
			continue;
		}*/

		// Recurse
		vector<Face_handle> visited;
		bool result = recurse(current_face, d, visited, t);

		if(result)
			cout << "POSSIBLE TO ESCAPE" << endl;
		else
			cout << "COULDN'T ESCAPE :(" << endl;
	}

}

int main()
{
  const Point<2> p1(-1.0, -1.0), p2(1.0, 1.0);
  Triangulation<2> triangulation;
  GridGenerator::hyper_rectangle(triangulation, p1, p2);
  triangulation.refine_global(num_levels);

  const Discretization<2> discretization(triangulation, 1);

  /**
   * Initialize a FieldType object with the return value from a function; this
   * utilizes the move constructor for FieldType.
   */
  Field<2> u = gaussian(discretization);

  std::cout << norm(u) << std::endl;

  /**
   * Reassign the FieldType object with the return value from another function;
   * this uses the move assignment operator.
   */
  u = parabola(discretization);

  std::cout << norm(u) << std::endl;

  return 0;
}

void u_new(Triangulation& T, const FT dt ) {

  for(F_v_it fv=T.finite_vertices_begin();
      fv!=T.finite_vertices_end();
      fv++) {

//  for(F_v_it fv=Tp.finite_vertices_begin();
//      fv!=Tp.finite_vertices_end();
//      fv++) {

    Vector_2 Ustar = fv->Ustar.val() ;
    Vector_2 gradp = fv->gradp.val() ;
    Vector_2 U = Ustar - dt * gradp;

    // relaxation mixing .-
    FT alpha=simu.alpha();
    Vector_2 U0=fv->U() ;
    Vector_2 U_mix = alpha*U0+ (1-alpha)*U ;

    fv->U.set( U_mix );
    fv->Delta_U.set(  U_mix - fv->Uold.val()  );
    
  }

  return;

}

//-------------------------------------------------------------------------
void TriangulationTest::TestArea()
//-------------------------------------------------------------------------
{
	Triangulation<double> triangulation;

	Point3D<double> p0;
	Point3D<double> p1;
	Point3D<double> p2;

	int pointId0 = 0;
	int pointId1 = 1;
	int pointId2 = 2;

	triangulation.points.push_back(p0);
	triangulation.points.push_back(p1);
	triangulation.points.push_back(p2);

	triangulation.AddTriangle(0,1,2);

	double area = triangulation.Area();

	CPPUNIT_ASSERT(area == 0);

	area = triangulation.Area(0);

	CPPUNIT_ASSERT(area == 0);

	area = triangulation.Area(0,1,2);

	CPPUNIT_ASSERT(area == 0);
}

// Compute all the finite voronoi vertices and the circumradius of all the finite cells.
void compute_voronoi_vertex_and_cell_radius(Triangulation& triang)
{
	bool is_there_any_problem_in_VV_computation = false;
	for (FCI cit = triang.finite_cells_begin();
			cit != triang.finite_cells_end(); cit ++)
	{
		//we tell CGAL to call our function if there is a problem
		//we also tell it not to die if things go haywire
		//CGAL::Failure_function old_ff = CGAL::set_error_handler(failure_func);
		//CGAL::Failure_behaviour old_fb = CGAL::set_error_behaviour(CGAL::CONTINUE);
		// be optimistic :-)
		//this is a global
		cgal_failed = false;
		cit->set_voronoi(triang.dual(cit));
		bool is_correct_computation = !cgal_failed;
		is_there_any_problem_in_VV_computation |= !is_correct_computation;
		if (cgal_failed)
		{
			// set cc the centroid of the cell.
			Vector cc = CGAL::NULL_VECTOR;
			for (int i = 0; i < 4; i ++)
			{
				cc = cc + (cit->vertex(i)->point() - CGAL::ORIGIN);
			}
			cc = (1./4.)*cc;
			cit->set_voronoi(CGAL::ORIGIN + cc);
		}
		//put everything back the way we found it,
		//CGAL::set_error_handler(old_ff);
		//CGAL::set_error_behaviour(old_fb);
		// set the cell radius.
		cit->set_cell_radius(CGAL::to_double((cit->vertex(0)->point()-cit->voronoi()) *(cit->vertex(0)->point()-cit->voronoi())));
	}
	return;
}

void load_alpha_on_fft( const Triangulation& T , CH_FFT& fft  ) {

  int Nb = fft.Nx();

  size_t align=fft.alignment();

  c_array al( Nb , Nb , align );

  for(F_v_it vit=T.vertices_begin();
      vit != T.vertices_end();
      vit++) {

    int nx = vit->nx.val();
    int ny = vit->ny.val();

    // "right" ordering 
    int i = ( Nb - 1 ) - ny ;
    int j = nx;

    // "wrong" ordering
    // int i = nx;
    // int j = ny;
    
    FT val =  vit->alpha0.val();
    //FT val =  vit->alpha.val();

    al(i,j) = val;

  }

  fft.set_f( al );
  
  return;
}

void load_fields_from_fft(const CH_FFT& fft , Triangulation& T  ) {

  int Nb = fft.Nx();

  c_array vx = fft.field_vel_x();
  c_array vy = fft.field_vel_y();
  c_array al = fft.field_f();

  for(F_v_it vit=T.vertices_begin();
      vit != T.vertices_end();
      vit++) {

    int nx = vit->nx.val();
    int ny = vit->ny.val();

    // "right" ordering 
    int i = ( Nb - 1 ) - ny ;
    int j = nx;

    // "wrong" ordering
    // int i = nx;
    // int j = ny;

    vit->U.set( Vector_2( real(vx(i,j)) , real(vy(i,j)) ) );

    vit->alpha.set( real( al(i,j) ) );

    //  TODO: return more fields (chem pot, pressure, force, etc)
  }

  return;
}

void Foam::DelaunayMeshTools::writeFixedPoints
(
    const fileName& fName,
    const Triangulation& t
)
{
    OFstream str(fName);

    Pout<< nl
        << "Writing fixed points to " << str.name() << endl;

    for
    (
        typename Triangulation::Finite_vertices_iterator vit =
            t.finite_vertices_begin();
        vit != t.finite_vertices_end();
        ++vit
    )
    {
        if (vit->fixed())
        {
            meshTools::writeOBJ(str, topoint(vit->point()));
        }
    }
}

void update_half_velocity( Triangulation& Tp , const bool overdamped ) {

  if (overdamped) return;
  
   for(F_v_it fv=Tp.finite_vertices_begin();
       fv!=Tp.finite_vertices_end();
       fv++) {

    Vector_2  v  = fv->U();

    // if (overdamped) 
    //    fv->U.set( v );
    // else {

    Vector_2  v0 = fv->Uold();
      //    Vector_2  v_star = fv->Ustar();

    fv->U.set(  2 * v - v0 );

      //   fv->U.set(  v + v_star - v0 );

      //    }

  }
  
  return;

}

Foam::tmp<Foam::Field<Type>> filterFarPoints
(
    const Triangulation& mesh,
    const Field<Type>& field
)
{
    tmp<Field<Type>> tNewField(new Field<Type>(field.size()));
    Field<Type>& newField = tNewField.ref();

    label added = 0;
    label count = 0;

    for
    (
        typename Triangulation::Finite_vertices_iterator vit =
            mesh.finite_vertices_begin();
        vit != mesh.finite_vertices_end();
        ++vit
    )
    {
        if (vit->real())
        {
            newField[added++] = field[count];
        }

        count++;
    }

    newField.resize(added);

    return tNewField;
}

void
club_contiguous_segment(Triangulation &triang,
	                map<int, cell_cluster> &cluster_set )
{
   for(FFI fit = triang.finite_facets_begin();
      fit != triang.finite_facets_end(); fit ++)
   {
      Cell_handle c[2] = {(*fit).first, (*fit).first->neighbor((*fit).second)};
      // if the two adjacent cells belong to the same cluster, continue.
      if(cluster_set[c[0]->id].find() == 
	 cluster_set[c[1]->id].find()) 
         continue;
      // if any one of them is not in any cluster, continue.
      if( ! cluster_set[c[0]->id].in_cluster || 
          ! cluster_set[c[1]->id].in_cluster )
         continue;

      // if any of the clusters is inside, continue.
      if( ! cluster_set[c[0]->id].outside ||
          ! cluster_set[c[1]->id].outside )
         continue;

      // merge the two clusters.
      merge_cluster(cluster_set, cluster_set[c[0]->id].find(), cluster_set[c[1]->id].find());
   }
}

//-------------------------------------------------------------------------
void TriangulationTest::TestFlipMinimize()
//-------------------------------------------------------------------------
{
	Triangulation<double> triangulation;

	std::vector<Point3D<double> >		points;

	double pos[3] = {0,0,1};

	Point3D<double> p0;
	Point3D<double> p1;
	Point3D<double> p2;

	int pointId0 = 0;
	int pointId1 = 1;
	int pointId2 = 2;

	triangulation.points.push_back(p0);
	triangulation.points.push_back(p1);
	triangulation.points.push_back(p2);

    triangulation.AddTriangle(0,1,2);

    int result = triangulation.FlipMinimize(0);

	CPPUNIT_ASSERT(result == 0);	
}

int main()
{
	//cout<<fixed<<setprecision(0);
	int cnt=0;
	while(true)
	{
		int n;
		cin>>n;
		if (n==0) break;
		
		int l,b,r,t;
		cin>>l>>b>>r>>t;

		Segment rect[4];
		rect[0] = Segment(Point(l,b),Point(r,b));
		rect[1] = Segment(Point(l,b),Point(l,t));
		rect[2] = Segment(Point(l,t),Point(r,t));
		rect[3] = Segment(Point(r,t),Point(r,b));
		
		vector<Point> points;
		for (int i=0;i<n;i++)
		{
			int x,y;
			cin>>x>>y;
			points.push_back(Point(x,y));	
		}
		Triangulation DT;
		DT.insert(points.begin(),points.end());
		

		solve(DT,rect,points);

	}
	return 0;
}

void R_s_k_2::draw_edge_footpoints(const Triangulation& mesh,
    const Edge& edge,
    const float red,
    const float green,
    const float blue)
{
  const Point& a = mesh.source_vertex(edge)->point();
  const Point& b = mesh.target_vertex(edge)->point();
  const Sample_vector& samples = edge.first->samples(edge.second);

  Sample_vector::const_iterator it;
  for (it = samples.begin(); it != samples.end(); ++it)
  {
    Sample_* sample = *it;
    Point p = sample->point();
    FT m = 0.5*(1.0 - sample->mass());

    Point q;
    if (mesh.get_plan(edge) == 0)
    {
      viewer->glColor3f(0.8f + m, m, m);
      FT Da = CGAL::squared_distance(p, a);
      FT Db = CGAL::squared_distance(p, b);
      if (Da < Db) q = a;
      else         q = b;
    }
    else
    {
      viewer->glColor3f(red + m, green + m, blue + m);
      FT t = sample->coordinate();
      q = CGAL::ORIGIN + (1.0 - t)*(a - CGAL::ORIGIN) + t*(b - CGAL::ORIGIN);
    }
    draw_segment(p, q);
  }
}

	void print_triangles(Triangulation dt){ 
		Point p;
		Triangle tri;
		Triangulation::Finite_faces_iterator it;
		for (it = dt.finite_faces_begin(); it != dt.finite_faces_end(); it++) {
			tri = dt.triangle(it);
			cout << "Triangle: " << tri << endl;
		}
	}	

int
main(int argc,char* argv[])
{
  const char* filename = (argc > 1) ? argv[1] : "data/points.xy";
  std::ifstream input(filename);
  Triangulation t;
  Filter is_finite(t);
  Finite_triangulation ft(t, is_finite, is_finite);

  Point p ;
  while(input >> p){
    t.insert(p);
  }

  vertex_iterator vit, ve;
  // Associate indices to the vertices
  int index = 0;
  // boost::tie assigns the first and second element of the std::pair
  // returned by boost::vertices to the variables vit and ve
  for(boost::tie(vit,ve)=boost::vertices(ft); vit!=ve; ++vit ){
    vertex_descriptor  vd = *vit;
    vertex_id_map[vd]= index++;
    }

  // Dijkstra's shortest path needs property maps for the predecessor and distance
  // We first declare a vector
  std::vector<vertex_descriptor> predecessor(boost::num_vertices(ft));
  // and then turn it into a property map
  boost::iterator_property_map<std::vector<vertex_descriptor>::iterator,
                               VertexIdPropertyMap>
    predecessor_pmap(predecessor.begin(), vertex_index_pmap);

  std::vector<double> distance(boost::num_vertices(ft));
  boost::iterator_property_map<std::vector<double>::iterator,
                               VertexIdPropertyMap>
    distance_pmap(distance.begin(), vertex_index_pmap);

  // start at an arbitrary vertex
  vertex_descriptor source = *boost::vertices(ft).first;
  std::cout << "\nStart dijkstra_shortest_paths at " << source->point() <<"\n";

  boost::dijkstra_shortest_paths(ft, source,
				 distance_map(distance_pmap)
				 .predecessor_map(predecessor_pmap)
				 .vertex_index_map(vertex_index_pmap));

  for(boost::tie(vit,ve)=boost::vertices(ft); vit!=ve; ++vit ){
    vertex_descriptor vd = *vit;
    std::cout << vd->point() << " [" <<  vertex_id_map[vd] << "] ";
    std::cout << " has distance = "  << boost::get(distance_pmap,vd)
	      << " and predecessor ";
    vd =  boost::get(predecessor_pmap,vd);
    std::cout << vd->point() << " [" <<  vertex_id_map[vd] << "]\n ";
  }

  return 0;
}

vector<int>
compute_smax(Triangulation& triang, 
             map<int, cell_cluster> &cluster_set,
             const double& mr)
{
   for(ACI cit = triang.all_cells_begin();
      cit != triang.all_cells_end(); cit ++)
   {
      cluster_set[cit->id] = cell_cluster(cit->id);
      cit->visited = false;
   }

   int max_cnt = 0;
   for(FCI cit = triang.finite_cells_begin();
      cit != triang.finite_cells_end(); cit ++)
   {
      if( ! is_maxima(cit) ) continue;

      #ifndef __OUTSIDE__
      if( cit->outside ) continue;
      #endif

      #ifndef __INSIDE__
      if( ! cit->outside ) continue;
      #endif
      
      if(max_cnt++%1000 == 0) cerr << "+";
      grow_maximum(cit, triang, cluster_set);
   }
   cerr << ".";

   // club_segment(triang, cluster_set, mr );
   // cerr << ".";
   club_contiguous_segment(triang, cluster_set );
   cerr << ".";

   // Compute the volume of each cluster. Remember after merging the 
   // 'rep' field is more useful than cluster_id.
   vector<int> cluster_ids;
   vector<double> cluster_vol;
   cluster_ids.clear();
   cluster_vol.clear();
   calc_cluster_volume_and_store_with_cluster_rep(triang, 
	  					  cluster_set,
			                          cluster_vol, 
						  cluster_ids);
   cerr << ".";

   // Sort the clusters with respect to the volumes.
   vector<int> sorted_indices;
   sorted_indices.clear();
   sort_cluster_wrt_volume(cluster_vol, cluster_ids, sorted_indices);
   cerr << ".";
   return sorted_indices;
}

Face_handle test_point_location(const Triangulation &t,
                                const Point &query,
                                const Triangulation::Locate_type &lt_in)
{
  Triangulation::Locate_type lt, lt2;
  int li, li2;
  Face_handle fh;
  CGAL::Bounded_side bs;
  CGAL::Oriented_side os;

  fh = t.locate(query, lt, li);
  CGAL_assertion(lt == lt_in);
  if (lt_in == Triangulation::EMPTY) {
    CGAL_assertion(fh == Face_handle());
    return fh;
  }

  bs = t.side_of_face(query, fh, lt2, li2);
  os = t.oriented_side(fh, query);
  CGAL_USE(bs);
  CGAL_USE(os);
  
  CGAL_assertion(lt2 == lt_in);

  switch (lt_in)
    {
    case Triangulation::VERTEX:
    case Triangulation::EDGE:
    {
      CGAL_assertion(fh != Face_handle());
      CGAL_assertion(bs == CGAL::ON_BOUNDARY);
      CGAL_assertion(os == CGAL::ON_ORIENTED_BOUNDARY);

      CGAL_assertion(li == li2);
      break;
    }
    case Triangulation::FACE:
    {
      CGAL_assertion(fh != Face_handle());
      CGAL_assertion(bs == CGAL::ON_BOUNDED_SIDE);
      CGAL_assertion(os == CGAL::ON_POSITIVE_SIDE);
      break;
    }
    case Triangulation::EMPTY:
    {
      // Handled above
      CGAL_assertion(false);
      break;
    }
    case Triangulation::OUTSIDE_CONVEX_HULL: CGAL_error();
    case Triangulation::OUTSIDE_AFFINE_HULL: CGAL_error();
    }

  return fh;
}

	void print_circles(Triangulation dt){ 
		Point p;
		Triangle tri;
		Triangulation::Finite_faces_iterator it;
		for (it = dt.finite_faces_begin(); it != dt.finite_faces_end(); it++) {
			p=dt.circumcenter(it);
			tri = dt.triangle(it);
			cout << "Circle center located at: " << p << " and with radius: "
				  << squared_distance(p, tri.vertex(0)) << endl;
		}
	}