Fix cpplint on examples.

Bug fix on persistence_from_simple_simplex_tree.
Add persistence examples tests.



git-svn-id: svn+ssh://scm.gforge.inria.fr/svnroot/gudhi/trunk@843 636b058d-ea47-450e-bf9e-a15bfbe3eedb

Former-commit-id: 9b7b73abb4a5d6bb110376deb689247bfdae035c

3 files changed, 237 insertions, 222 deletions

diff --git a/src/Contraction/example/Garland_heckbert.cpp b/src/Contraction/example/Garland_heckbert.cpp
index a41f65aa..70f29b6a 100644
--- a/src/Contraction/example/Garland_heckbert.cpp
+++ b/src/Contraction/example/Garland_heckbert.cpp
@@ -1,7 +1,4 @@
-/*
- * Garland_heckbert.h
- *  Created on: Feb 10, 2015
- * This file is part of the Gudhi Library. The Gudhi library 
+/* This file is part of the Gudhi Library. The Gudhi library 
  *    (Geometric Understanding in Higher Dimensions) is a generic C++ 
  *    library for computational topology.
  *
@@ -28,12 +25,13 @@
 #ifndef GARLAND_HECKBERT_H_
 #define GARLAND_HECKBERT_H_
 
+#include <gudhi/Point.h>
+#include <gudhi/Edge_contraction.h>
+#include <gudhi/Skeleton_blocker.h>
+#include <gudhi/Off_reader.h>
+
 #include <boost/timer/timer.hpp>
 #include <iostream>
-#include "gudhi/Point.h"
-#include "gudhi/Edge_contraction.h"
-#include "gudhi/Skeleton_blocker.h"
-#include "gudhi/Off_reader.h"
 
 #include "Garland_heckbert/Error_quadric.h"
 
@@ -51,7 +49,6 @@ struct Geometry_trait {
  */
 struct Garland_heckbert_traits : public Skeleton_blocker_simple_geometric_traits<Geometry_trait> {
  public:
-
   struct Garland_heckbert_vertex : public Simple_geometric_vertex {
     Error_quadric<Geometry_trait::Point> quadric;
   };
@@ -68,6 +65,7 @@ typedef Skeleton_blocker_contractor<Complex> Complex_contractor;
  */
 class GH_placement : public Gudhi::contraction::Placement_policy<EdgeProfile> {
   Complex& complex_;
+
  public:
   typedef Gudhi::contraction::Placement_policy<EdgeProfile>::Placement_type Placement_type;
 
@@ -91,8 +89,8 @@ class GH_placement : public Gudhi::contraction::Placement_policy<EdgeProfile> {
  */
 class GH_cost : public Gudhi::contraction::Cost_policy<EdgeProfile> {
   Complex& complex_;
- public:
 
+ public:
   typedef Gudhi::contraction::Cost_policy<EdgeProfile>::Cost_type Cost_type;
 
   GH_cost(Complex& complex) : complex_(complex) { }
@@ -115,13 +113,13 @@ class GH_cost : public Gudhi::contraction::Cost_policy<EdgeProfile> {
  */
 class GH_visitor : public Gudhi::contraction::Contraction_visitor<EdgeProfile> {
   Complex& complex_;
- public:
 
+ public:
   GH_visitor(Complex& complex) : complex_(complex) { }
 
-  //Compute quadrics for every vertex v
-  //The quadric of v consists in the sum of quadric
-  //of every triangles passing through v weighted by its area
+  // Compute quadrics for every vertex v
+  // The quadric of v consists in the sum of quadric
+  // of every triangles passing through v weighted by its area
 
   void on_started(Complex & complex) override {
     for (auto v : complex.vertex_range()) {
@@ -147,7 +145,8 @@ class GH_visitor : public Gudhi::contraction::Contraction_visitor<EdgeProfile> {
 
 int main(int argc, char *argv[]) {
   if (argc != 4) {
-    std::cerr << "Usage " << argv[0] << " input.off output.off N to load the file input.off, contract N edges and save the result to output.off.\n";
+    std::cerr << "Usage " << argv[0] << " input.off output.off N to load the file input.off, contract N edges and save "
+        << "the result to output.off.\n";
     return EXIT_FAILURE;
   }
 
@@ -172,8 +171,7 @@ int main(int argc, char *argv[]) {
                                 new GH_cost(complex),
                                 new GH_placement(complex),
                                 contraction::make_link_valid_contraction<EdgeProfile>(),
-                                new GH_visitor(complex)
-                                );
+                                new GH_visitor(complex));
 
   std::cout << "Contract " << num_contractions << " edges" << std::endl;
   contractor.contract_edges(num_contractions);
@@ -183,7 +181,7 @@ int main(int argc, char *argv[]) {
       complex.num_edges() << " edges and" <<
       complex.num_triangles() << " triangles." << std::endl;
 
-  //write simplified complex
+  // write simplified complex
   Skeleton_blocker_off_writer<Complex> off_writer(argv[2], complex);
 
   return EXIT_SUCCESS;
@@ -191,4 +189,4 @@ int main(int argc, char *argv[]) {
 
 
 
-#endif /* GARLAND_HECKBERT_H_ */
+#endif  // GARLAND_HECKBERT_H_
diff --git a/src/Contraction/example/Garland_heckbert/Error_quadric.h b/src/Contraction/example/Garland_heckbert/Error_quadric.h
index 72134c9d..a033aa00 100644
--- a/src/Contraction/example/Garland_heckbert/Error_quadric.h
+++ b/src/Contraction/example/Garland_heckbert/Error_quadric.h
@@ -1,164 +1,182 @@
-/*

- * Error_quadric.h

+/* This file is part of the Gudhi Library. The Gudhi library 

+ *    (Geometric Understanding in Higher Dimensions) is a generic C++ 

+ *    library for computational topology.

  *

- *  Created on: 24 janv. 2014

- *      Author: dsalinas

+ *    Author(s):       David Salinas

+ *

+ *    Copyright (C) 2014  INRIA Sophia Antipolis-M�diterran�e (France)

+ *

+ *    This program is free software: you can redistribute it and/or modify

+ *    it under the terms of the GNU General Public License as published by

+ *    the Free Software Foundation, either version 3 of the License, or

+ *    (at your option) any later version.

+ *

+ *    This program is distributed in the hope that it will be useful,

+ *    but WITHOUT ANY WARRANTY; without even the implied warranty of

+ *    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the

+ *    GNU General Public License for more details.

+ *

+ *    You should have received a copy of the GNU General Public License

+ *    along with this program.  If not, see <http://www.gnu.org/licenses/>.

+ * 

  */

 

 #ifndef ERROR_QUADRIC_H_

 #define ERROR_QUADRIC_H_

 

-#include <vector>

-#include <utility>

 #include <boost/optional/optional.hpp>

 

+#include <vector>

+#include <utility>

 

-template <typename Point> class Error_quadric{

-private :

-	double coeff[10];

-

-public :

-	Error_quadric(){

-		clear();

-	}

-

-	/**

-	 * Quadric corresponding to the L2 distance to the plane.

-	 *

-	 * According to the notation of Garland Heckbert, they

-	 * denote a quadric symetric matrix as :

-	 * Q = [ q11 q12 q13 q14]

-	 *     [ q12 q22 q23 q24]

-	 *     [ q13 q23 q33 q34]

-	 *     [ q14 q24 q34 q44]

-	 *

-	 * which is represented by a vector with 10 elts that

-	 * are denoted ci for clarity with :

-	 * Q = [ c0  c1  c2  c3 ]

-	 *     [ c1  c4  c5  c6 ]

-	 *     [ c2  c5  c7  c8 ]

-	 *     [ c3  c6  c8  c9 ]

-	 *

-	 * The constructor return the quadrics that represents

-	 * the squared distance to the plane defined by triangle p0,p1,p2

-	 * times the area of triangle p0,p1,p2.

-	 */

-	Error_quadric(const Point & p0,const Point & p1,const Point & p2){

-

-		Point normal(unit_normal(p0,p1,p2));

-		double a=normal[0];

-		double b=normal[1];

-		double c=normal[2];

-		double d= -a*p0[0]-b*p0[1]-c*p0[2];

-		coeff[0] = a*a ;

-		coeff[1] = a*b ;

-		coeff[2] = a*c ;

-		coeff[3] = a*d ;

-		coeff[4] = b*b ;

-		coeff[5] = b*c ;

-		coeff[6] = b*d ;

-		coeff[7] = c*c ;

-		coeff[8] = c*d ;

-		coeff[9] = d*d ;

-

-		double area_p0p1p2 = std::sqrt(squared_area(p0,p1,p2));

-		for(auto& x : coeff)

-			x*= area_p0p1p2;

-	}

-

-

-	inline double squared_area(const Point& p0,const Point& p1,const Point& p2) {

-		//if (x1,x2,x3) = p1-p0 and (y1,y2,y3) = p2-p0

-		//then the squared area is = (u^2+v^2+w^2)/4

-		//with:		u = x2 * y3 - x3 * y2;

-		//		    v = x3 * y1 - x1 * y3;

-		//	     	w = x1 * y2 - x2 * y1;

-		Point p0p1(p1-p0);

-		Point p0p2(p2-p0);

-		double A = p0p1[1] * p0p2[2] - p0p1[2] * p0p2[1];

-		double B = p0p1[2] * p0p2[0] - p0p1[0] * p0p2[2];

-		double C = p0p1[0] * p0p2[1] - p0p1[1] * p0p2[0];

-		return 1./4. * (A*A+B*B+C*C);

-	}

-

-

-	void clear(){

-		for(auto& x:coeff)

-			x=0;

-	}

-

-	Error_quadric& operator+=(const Error_quadric& other){

-		if(this!=&other)

-			for(int i = 0 ; i < 10; ++i)

-				coeff[i] += other.coeff[i];

-		return *this;

-	}

-

-	/**

-	 * @return The quadric quost defined by the scalar product v^T Q v where Q is the quadratic form of Garland/Heckbert

-	 */

-	inline double cost(const Point& point) const{

-		double cost =

-				coeff[0]*point.x()*point.x()+coeff[4]*point.y()*point.y()+coeff[7]*point.z()*point.z()

-				+2*(coeff[1]*point.x()*point.y()+coeff[5]*point.y()*point.z()+coeff[2]*point.z()*point.x())

-				+2*(coeff[3]*point.x()+coeff[6]*point.y()+coeff[8]*point.z())

-				+coeff[9];

-		if(cost<0) return 0;

-		else {

-			return cost;

-		}

-	}

-

-	inline double grad_determinant() const{

-		return

-				coeff[0] * coeff[4] * coeff[7]

-				                            -  coeff[0] * coeff[5] * coeff[5]

-				                                                           -  coeff[1] * coeff[1] * coeff[7]

-				                                                                                          +2*coeff[1] * coeff[5] * coeff[2]

-				                                                                                                                         -  coeff[4] * coeff[2] * coeff[2];

-	}

-

-	/**

-	 * Return the point such that it minimizes the gradient of the quadric.

-	 * Det must be passed with the determinant value of the gradient (should be non zero).

-	 */

-	inline Point solve_linear_gradient(double det) const{

-		return Point({

-				(-coeff[1]*coeff[5]*coeff[8]+coeff[1]*coeff[7]*coeff[6]+coeff[2]*coeff[8]*coeff[4]-coeff[2]*coeff[5]*coeff[6]-coeff[3]*coeff[4]*coeff[7]+coeff[3]*coeff[5]*coeff[5])/ det,

-				(coeff[0]*coeff[5]*coeff[8]-coeff[0]*coeff[7]*coeff[6]-coeff[5]*coeff[2]*coeff[3]-coeff[1]*coeff[2]*coeff[8]+coeff[6]*coeff[2]*coeff[2]+coeff[1]*coeff[3]*coeff[7])/det,

-				(-coeff[8]*coeff[0]*coeff[4]+coeff[8]*coeff[1]*coeff[1]+coeff[2]*coeff[3]*coeff[4]+coeff[5]*coeff[0]*coeff[6]-coeff[5]*coeff[1]*coeff[3]-coeff[1]*coeff[2]*coeff[6])/det

-		});

-	}

-

-

-	/**

-	 * returns the point that minimizes the quadric.

-	 * It inverses the quadric if its determinant is higher that a given threshold .

-	 * If the determinant is lower than this value the returned value is uninitialized.

-	 */

-	boost::optional<Point> min_cost(double scale=1) const{

-		//		const double min_determinant = 1e-4 * scale*scale;

-		const double min_determinant = 1e-5;

-		boost::optional<Point> pt_res;

-		double det = grad_determinant();

-		if (std::abs(det)>min_determinant)

-			pt_res = solve_linear_gradient(det);

-		return pt_res;

-	}

-

-	friend std::ostream& operator<< (std::ostream& stream, const Error_quadric& quadric) {

-		stream << "\n[ "<<quadric.coeff[0]<<","<<quadric.coeff[1]<<","<<quadric.coeff[2]<<","<<quadric.coeff[3]<<";\n";

-		stream <<    " "<<quadric.coeff[1]<<","<<quadric.coeff[4]<<","<<quadric.coeff[5]<<","<<quadric.coeff[6]<<";\n";

-		stream <<    " "<<quadric.coeff[2]<<","<<quadric.coeff[5]<<","<<quadric.coeff[7]<<","<<quadric.coeff[8]<<";\n";

-		stream <<    " "<<quadric.coeff[3]<<","<<quadric.coeff[6]<<","<<quadric.coeff[8]<<","<<quadric.coeff[9]<<"]";

-		return stream;

-	}

-

-

+template <typename Point> class Error_quadric {

+ private:

+  double coeff[10];

+

+ public:

+  Error_quadric() {

+    clear();

+  }

+

+  /**

+   * Quadric corresponding to the L2 distance to the plane.

+   *

+   * According to the notation of Garland Heckbert, they

+   * denote a quadric symetric matrix as :

+   * Q = [ q11 q12 q13 q14]

+   *     [ q12 q22 q23 q24]

+   *     [ q13 q23 q33 q34]

+   *     [ q14 q24 q34 q44]

+   *

+   * which is represented by a vector with 10 elts that

+   * are denoted ci for clarity with :

+   * Q = [ c0  c1  c2  c3 ]

+   *     [ c1  c4  c5  c6 ]

+   *     [ c2  c5  c7  c8 ]

+   *     [ c3  c6  c8  c9 ]

+   *

+   * The constructor return the quadrics that represents

+   * the squared distance to the plane defined by triangle p0,p1,p2

+   * times the area of triangle p0,p1,p2.

+   */

+  Error_quadric(const Point & p0, const Point & p1, const Point & p2) {

+    Point normal(unit_normal(p0, p1, p2));

+    double a = normal[0];

+    double b = normal[1];

+    double c = normal[2];

+    double d = -a * p0[0] - b * p0[1] - c * p0[2];

+    coeff[0] = a*a;

+    coeff[1] = a*b;

+    coeff[2] = a*c;

+    coeff[3] = a*d;

+    coeff[4] = b*b;

+    coeff[5] = b*c;

+    coeff[6] = b*d;

+    coeff[7] = c*c;

+    coeff[8] = c*d;

+    coeff[9] = d*d;

+

+    double area_p0p1p2 = std::sqrt(squared_area(p0, p1, p2));

+    for (auto& x : coeff)

+      x *= area_p0p1p2;

+  }

+

+  inline double squared_area(const Point& p0, const Point& p1, const Point& p2) {

+    // if (x1,x2,x3) = p1-p0 and (y1,y2,y3) = p2-p0

+    // then the squared area is = (u^2+v^2+w^2)/4

+    // with: u = x2 * y3 - x3 * y2;

+    //       v = x3 * y1 - x1 * y3;

+    //       w = x1 * y2 - x2 * y1;

+    Point p0p1(p1 - p0);

+    Point p0p2(p2 - p0);

+    double A = p0p1[1] * p0p2[2] - p0p1[2] * p0p2[1];

+    double B = p0p1[2] * p0p2[0] - p0p1[0] * p0p2[2];

+    double C = p0p1[0] * p0p2[1] - p0p1[1] * p0p2[0];

+    return 1. / 4. * (A * A + B * B + C * C);

+  }

+

+  void clear() {

+    for (auto& x : coeff)

+      x = 0;

+  }

+

+  Error_quadric& operator+=(const Error_quadric& other) {

+    if (this != &other) {

+      for (int i = 0; i < 10; ++i)

+        coeff[i] += other.coeff[i];

+    }

+    return *this;

+  }

+

+  /**

+   * @return The quadric quost defined by the scalar product v^T Q v where Q is the quadratic form of Garland/Heckbert

+   */

+  inline double cost(const Point& point) const {

+    double cost =

+        coeff[0] * point.x() * point.x() + coeff[4] * point.y() * point.y() + coeff[7] * point.z() * point.z()

+        + 2 * (coeff[1] * point.x() * point.y() + coeff[5] * point.y() * point.z() + coeff[2] * point.z() * point.x())

+        + 2 * (coeff[3] * point.x() + coeff[6] * point.y() + coeff[8] * point.z())

+        + coeff[9];

+    if (cost < 0) {

+      return 0;

+    } else {

+      return cost;

+    }

+  }

+

+  inline double grad_determinant() const {

+    return

+    coeff[0] * coeff[4] * coeff[7]

+        - coeff[0] * coeff[5] * coeff[5]

+        - coeff[1] * coeff[1] * coeff[7]

+        + 2 * coeff[1] * coeff[5] * coeff[2]

+        - coeff[4] * coeff[2] * coeff[2];

+  }

+

+  /**

+   * Return the point such that it minimizes the gradient of the quadric.

+   * Det must be passed with the determinant value of the gradient (should be non zero).

+   */

+  inline Point solve_linear_gradient(double det) const {

+    return Point({

+                 (-coeff[1] * coeff[5] * coeff[8] + coeff[1] * coeff[7] * coeff[6] + coeff[2] * coeff[8] * coeff[4] -

+                     coeff[2] * coeff[5] * coeff[6] - coeff[3] * coeff[4] * coeff[7] + coeff[3] * coeff[5] * coeff[5])

+                     / det,

+                 (coeff[0] * coeff[5] * coeff[8] - coeff[0] * coeff[7] * coeff[6] - coeff[5] * coeff[2] * coeff[3] -

+                     coeff[1] * coeff[2] * coeff[8] + coeff[6] * coeff[2] * coeff[2] + coeff[1] * coeff[3] * coeff[7])

+                     / det,

+                 (-coeff[8] * coeff[0] * coeff[4] + coeff[8] * coeff[1] * coeff[1] + coeff[2] * coeff[3] * coeff[4] +

+                     coeff[5] * coeff[0] * coeff[6] - coeff[5] * coeff[1] * coeff[3] - coeff[1] * coeff[2] * coeff[6])

+                     / det

+    });

+  }

+

+  /**

+   * returns the point that minimizes the quadric.

+   * It inverses the quadric if its determinant is higher that a given threshold .

+   * If the determinant is lower than this value the returned value is uninitialized.

+   */

+  boost::optional<Point> min_cost(double scale = 1) const {

+    // const double min_determinant = 1e-4 * scale*scale;

+    const double min_determinant = 1e-5;

+    boost::optional<Point> pt_res;

+    double det = grad_determinant();

+    if (std::abs(det) > min_determinant)

+      pt_res = solve_linear_gradient(det);

+    return pt_res;

+  }

+

+  friend std::ostream& operator<<(std::ostream& stream, const Error_quadric& quadric) {

+    stream << "\n[ " << quadric.coeff[0] << "," << quadric.coeff[1] << "," << quadric.coeff[2] << "," <<

+        quadric.coeff[3] << ";\n";

+    stream << " " << quadric.coeff[1] << "," << quadric.coeff[4] << "," << quadric.coeff[5] << "," <<

+        quadric.coeff[6] << ";\n";

+    stream << " " << quadric.coeff[2] << "," << quadric.coeff[5] << "," << quadric.coeff[7] << "," <<

+        quadric.coeff[8] << ";\n";

+    stream << " " << quadric.coeff[3] << "," << quadric.coeff[6] << "," << quadric.coeff[8] << "," <<

+        quadric.coeff[9] << "]";

+    return stream;

+  }

 };

 

-

-

-

-#endif /* ERROR_QUADRIC_H_ */

-

+#endif  // ERROR_QUADRIC_H_

diff --git a/src/Contraction/example/Rips_contraction.cpp b/src/Contraction/example/Rips_contraction.cpp
index bd0a8b8c..d21246ed 100644
--- a/src/Contraction/example/Rips_contraction.cpp
+++ b/src/Contraction/example/Rips_contraction.cpp
@@ -19,24 +19,23 @@
  *    You should have received a copy of the GNU General Public License
  *    along with this program.  If not, see <http://www.gnu.org/licenses/>.
  */
+#include <gudhi/Edge_contraction.h>
+#include <gudhi/Skeleton_blocker.h>
+#include <gudhi/Off_reader.h>
+#include <gudhi/Point.h>
+
 #include <boost/timer/timer.hpp>
 #include <iostream>
-#include "gudhi/Edge_contraction.h"
-#include "gudhi/Skeleton_blocker.h"
-#include "gudhi/Off_reader.h"
-#include "gudhi/Point.h"
 
 using namespace std;
 using namespace Gudhi;
 using namespace skbl;
 using namespace contraction;
 
-
-struct Geometry_trait{
-	typedef Point_d Point;
+struct Geometry_trait {
+  typedef Point_d Point;
 };
 
-
 typedef Geometry_trait::Point Point;
 typedef Skeleton_blocker_simple_geometric_traits<Geometry_trait> Complex_geometric_traits;
 typedef Skeleton_blocker_geometric_complex< Complex_geometric_traits > Complex;
@@ -44,62 +43,62 @@ typedef Edge_profile<Complex> Profile;
 typedef Skeleton_blocker_contractor<Complex> Complex_contractor;
 
 template<typename ComplexType>
-void build_rips(ComplexType& complex, double offset){
-	if (offset<=0) return;
-	auto vertices = complex.vertex_range();
-	for (auto p = vertices.begin(); p != vertices.end(); ++p)
-		for (auto q = p; ++q != vertices.end(); /**/){
-			if ( squared_dist(complex.point(*p), complex.point(*q)) < 4 * offset * offset)
-				complex.add_edge(*p,*q);
-		}
+void build_rips(ComplexType& complex, double offset) {
+  if (offset <= 0) return;
+  auto vertices = complex.vertex_range();
+  for (auto p = vertices.begin(); p != vertices.end(); ++p)
+    for (auto q = p; ++q != vertices.end(); /**/) {
+      if (squared_dist(complex.point(*p), complex.point(*q)) < 4 * offset * offset)
+        complex.add_edge(*p, *q);
+    }
 }
 
-int main (int argc, char *argv[])
-{
-	if (argc!=3){
-		std::cerr << "Usage "<<argv[0]<<" ../../../data/meshes/SO3_10000.off 0.3 to load the file ../../data/SO3_10000.off and contract the Rips complex built with paremeter 0.3.\n";
-		return -1;
-	}
+int main(int argc, char *argv[]) {
+  if (argc != 3) {
+    std::cerr << "Usage " << argv[0] << " ../../../data/meshes/SO3_10000.off 0.3 to load the file " <<
+        "../../data/SO3_10000.off and contract the Rips complex built with paremeter 0.3.\n";
+    return -1;
+  }
 
-	Complex complex;
+  Complex complex;
 
-	// load only the points
-	Skeleton_blocker_off_reader<Complex> off_reader(argv[1],complex,true);
-	if(!off_reader.is_valid()){
-		std::cerr << "Unable to read file:"<<argv[1]<<std::endl;
-		return EXIT_FAILURE;
-	}
+  // load only the points
+  Skeleton_blocker_off_reader<Complex> off_reader(argv[1], complex, true);
+  if (!off_reader.is_valid()) {
+    std::cerr << "Unable to read file:" << argv[1] << std::endl;
+    return EXIT_FAILURE;
+  }
 
-	std::cout << "Build the Rips complex with "<<complex.num_vertices()<<" vertices"<<std::endl;
+  std::cout << "Build the Rips complex with " << complex.num_vertices() << " vertices" << std::endl;
 
-	build_rips(complex,atof(argv[2]));
+  build_rips(complex, atof(argv[2]));
 
-	boost::timer::auto_cpu_timer t;
+  boost::timer::auto_cpu_timer t;
 
-	std::cout << "Initial complex has "<<
-			complex.num_vertices()<<" vertices and "<<
-			complex.num_edges()<<" edges"<<std::endl;
+  std::cout << "Initial complex has " <<
+      complex.num_vertices() << " vertices and " <<
+      complex.num_edges() << " edges" << std::endl;
 
-	Complex_contractor contractor(complex,
-			new Edge_length_cost<Profile>,
-			contraction::make_first_vertex_placement<Profile>(),
-			contraction::make_link_valid_contraction<Profile>(),
-			contraction::make_remove_popable_blockers_visitor<Profile>());
-	contractor.contract_edges();
+  Complex_contractor contractor(complex,
+                                new Edge_length_cost<Profile>,
+                                contraction::make_first_vertex_placement<Profile>(),
+                                contraction::make_link_valid_contraction<Profile>(),
+                                contraction::make_remove_popable_blockers_visitor<Profile>());
+  contractor.contract_edges();
 
-	std::cout << "Counting final number of simplices \n";
-	unsigned num_simplices = std::distance(complex.simplex_range().begin(),complex.simplex_range().end());
+  std::cout << "Counting final number of simplices \n";
+  unsigned num_simplices = std::distance(complex.simplex_range().begin(), complex.simplex_range().end());
 
-	std::cout << "Final complex has "<<
-			complex.num_vertices()<<" vertices, "<<
-			complex.num_edges()<<" edges, "<<
-			complex.num_blockers()<<" blockers and "<<
-			num_simplices<<" simplices"<<std::endl;
+  std::cout << "Final complex has " <<
+      complex.num_vertices() << " vertices, " <<
+      complex.num_edges() << " edges, " <<
+      complex.num_blockers() << " blockers and " <<
+      num_simplices << " simplices" << std::endl;
 
 
-	std::cout << "Time to simplify and enumerate simplices:\n";
+  std::cout << "Time to simplify and enumerate simplices:\n";
 
-	return EXIT_SUCCESS;
+  return EXIT_SUCCESS;
 }