diff --git a/CMakeLists.txt b/CMakeLists.txt
index 143319a..7b732e5 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -10,7 +10,8 @@ project(simplicial CXX)
 #   * boost under /usr/include/boost
 #   * CGAL under /usr/include/CGAL
 #
-set(EIGEN3_INCLUDE_DIR "/usr/include/eigen3")
+set(EIGEN3_INCLUDE_DIR "/usr/local/include/eigen3")
+set(EIGEN3_INCLUDE_DIR "/home/crvs/eigen-3.3.2")
 include_directories(${EIGEN3_INCLUDE_DIR})
 set(CGAL_INCLUDE_DIR "/usr/local/include")
 include_directories(${CGAL_INCLUDE_DIR})
@@ -21,19 +22,61 @@ include(${CGAL_USE_FILE})
 # end_requirements
 
 
+# Using CLANG instead of gcc for reasons
+set(CMAKE_C_COMPILER             "/usr/bin/clang-3.8")
+set(CMAKE_C_FLAGS                "-Wall -std=c99")
+set(CMAKE_C_FLAGS_DEBUG          "-g")
+set(CMAKE_C_FLAGS_MINSIZEREL     "-Os -DNDEBUG")
+set(CMAKE_C_FLAGS_RELEASE        "-O4 -DNDEBUG")
+set(CMAKE_C_FLAGS_RELWITHDEBINFO "-O2 -g")
+
+set(CMAKE_CXX_COMPILER             "/usr/bin/clang++-3.8")
+set(CMAKE_CXX_FLAGS                "-Wall -std=c++11")
+set(CMAKE_CXX_FLAGS_DEBUG          "-g")
+set(CMAKE_CXX_FLAGS_MINSIZEREL     "-Os -DNDEBUG")
+set(CMAKE_CXX_FLAGS_RELEASE        "-O4 -DNDEBUG")
+set(CMAKE_CXX_FLAGS_RELWITHDEBINFO "-O2 -g")
+# Linker flags for
+# set(CMAKE_CXX_FLAGS "-lspatialindex_c -lspatialindex ${CMAKE_CXX_FLAGS}")
+
+set(CMAKE_AR      "/usr/bin/llvm-ar")
+set(CMAKE_LINKER  "/usr/bin/llvm-ld")
+set(CMAKE_NM      "/usr/bin/llvm-nm")
+set(CMAKE_OBJDUMP "/usr/bin/llvm-objdump")
+set(CMAKE_RANLIB  "/usr/bin/llvm-ranlib")
 
 set(CMAKE_BUILD_TYPE Debug)
-include_directories("./lib")
 
-set(CMAKE_CXX_FLAGS "-std=c++11 ${CMAKE_CXX_FLAGS}")
+include_directories("./lib") # where the header files are
 
-set(STUPID_SOURCES "./src/main.cpp")
-add_executable(my_stupid_test ${STUPID_SOURCES})
+add_library(scomplex SHARED "./lib/scomplex/simplicial_complex.cpp")
+add_library(coeffflow SHARED "./lib/scomplex/coeff_flow.cpp")
+add_library(pathsnap SHARED "./lib/scomplex/path_snapper.cpp")
 
-set(YOTHER_STUPID_SOURCES "./src/graphtest.cpp")
-add_executable(yother_stupid_test ${YOTHER_STUPID_SOURCES})
+add_executable(test0 "src/test0.cpp")
+target_link_libraries(test0 scomplex)
 
-# set(CMAKE_CXX_FLAGS "-lspatialindex_c -lspatialindex ${CMAKE_CXX_FLAGS}")
- set(OTHER_STUPID_SOURCES "./src/libtest.cpp")
-add_executable(my_other_stupid_test ${OTHER_STUPID_SOURCES})
+add_executable(test1 "src/test1.cpp")
+target_link_libraries(test1 scomplex pathsnap)
+
+add_executable(test2 "src/test2.cpp")
+target_link_libraries(test2 scomplex pathsnap)
+
+add_executable(test3 "src/test3.cpp")
+target_link_libraries(test3 scomplex pathsnap)
+
+add_executable(test4 "src/test4.cpp")
+target_link_libraries(test4 scomplex pathsnap)
+
+add_executable(test5 "src/test5.cpp")
+target_link_libraries(test5 scomplex pathsnap)
+
+# set(STUPID_SOURCES "./src/main.cpp")
+# add_executable(my_stupid_test ${STUPID_SOURCES})
+
+# set(YOTHER_STUPID_SOURCES "./src/graphtest.cpp")
+# add_executable(yother_stupid_test ${YOTHER_STUPID_SOURCES})
+
+# set(OTHER_STUPID_SOURCES "./src/libtest.cpp")
+# add_executable(my_other_stupid_test ${OTHER_STUPID_SOURCES})
 
diff --git a/lib/scomplex/Simplicial_Complex.h b/lib/scomplex/Simplicial_Complex.h
deleted file mode 100644
index 780608c..0000000
--- a/lib/scomplex/Simplicial_Complex.h
+++ /dev/null
@@ -1,293 +0,0 @@
-/**
- * @file Simplicial_Complex.h
- * @brief  utilities to create and manipulate geometric simplicial complexes
- * @author Joao Carvalho
- * @version 0.1
- * @date 2017-01-20
- *
- * Copyright (C)
- * 2017 - Joao Carvalho
- * 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 2
- * 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, write to the Free Software
- * Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.
- *
- */
-#ifndef QUOT_SIMPLICIAL_COMPLEX
-#define QUOT_SIMPLICIAL_COMPLEX
-
-#include <Eigen/Sparse>
-
-#include <gudhi/Simplex_tree/Simplex_tree_siblings.h>
-#include <gudhi/Simplex_tree.h>
-
-#include <scomplex/base_utils.hpp>
-// #include <scomplex/utils.h>
-
-#include <math.h>
-#include <vector>
-#include <set>
-
-#include <typeinfo>
-
-// point_t needs to have an empty constructor, and a size() function call,
-// further it needs a begin() function call. Mostly it needs to behave like a
-// std::vector
-
-namespace simplicial {
-typedef std::pair<int, std::vector<double>> chain_t;
-
-template <typename point_t>
-class SimplicialComplex {
-   private:
-    // the type of matrix to be used
-    typedef typename Eigen::SparseMatrix<double> matrix_t;
-
-    struct SimpleOptions : Gudhi::Simplex_tree_options_full_featured {
-        // simplex tree options
-        typedef int Vertex_handle;
-    };
-    typedef typename Gudhi::Simplex_tree<SimpleOptions> SimplexTree;
-    typedef typename Gudhi::Simplex_tree<SimpleOptions>::Simplex_handle
-        Simplex_handle;
-
-    typedef std::vector<Simplex_handle *> level_t;
-    typedef std::vector<level_t> levels_t;
-
-   public:
-    std::vector<std::vector<int>> get_level(int dimen) {
-        std::vector<std::vector<int>> level;
-        for (auto simp : simplices.complex_simplex_range()) {
-            if (simplices.dimension(simp) == dimen) {
-                std::vector<int> v_simp;
-                for (auto v : simplices.simplex_vertex_range(simp)) {
-                    v_simp.push_back(v);
-                }
-                level.push_back(v_simp);
-            }
-        }
-        return level;
-    }
-
-    std::vector<point_t> points;
-    SimplexTree simplices;
-
-    int get_dimension(int level) { return dimensions.at(level); }
-
-    SimplicialComplex<point_t>() { geometric_q = true; }
-
-    /**
-     * @brief  builds a simplicial complex object out of points and "triangles"
-     *
-     * @param points_a  points
-     * @param tris      the top dimensional simplicial complex structure
-     */
-    SimplicialComplex<point_t>(std::list<point_t> points_a,
-                               std::list<std::list<int>> tris) {
-        geometric_q = true;
-        points = std::vector<point_t>();
-        for (auto pt : points_a) {
-            points.push_back(pt);
-        }
-        for (auto s : tris) {
-            // cant have simplices with repeated vertices so we dedupe the lists
-            simplices.insert_simplex_and_subfaces(dedupe_list(s));
-            int d = s.size() - 1;
-            if (simplices.dimension() < d) {
-                simplices.set_dimension(d);
-            }
-            // error (something is out-of-range)
-            // calculate_matrices();
-        }
-
-        // we need to count i--simplices in order to get a key for them
-        std::vector<int> count = std::vector<int>(simplices.dimension() + 1, 0);
-        levels_t levels;
-        levels = levels_t(simplices.dimension() + 1, level_t());
-
-        auto simplex_range = simplices.complex_simplex_range();
-        for (auto s : simplex_range) {
-            int d = simplices.dimension(s);
-            simplices.assign_key(s, count.at(d)++);
-            levels.at(d).push_back(&s);
-            auto v_range = simplices.simplex_vertex_range(s);
-            auto sbd = simplices.boundary_simplex_range(s);
-        }
-        dimensions = count;
-    }
-
-    matrix_t get_boundary(int d) {
-        if (boundary_matrices.size() == 0) {
-            calculate_matrices();
-        }
-        if (d < boundary_matrices.size()) {
-            return boundary_matrices.at(d);
-        } else {
-            return matrix_t(0, 0);
-        }
-    }
-
-    SimplicialComplex<point_t> quotient(int f(point_t)) {
-        int n_points = 1;
-        auto corresp = std::vector<int>();
-
-        // making the correspondence between points in the original complex and
-        // points in the quotient. point index 0 corresponds to the "virtual"
-        // point (in case it exists).
-        std::list<point_t> q_points;
-        bool geometric = true;
-
-        for (auto p : points) {
-            if (f(p) != 0) {
-                corresp.push_back(n_points++);
-                q_points.push_back(p);
-            } else {
-                if (geometric) {
-                    // need to realize that it is not geometric, (flip geometric
-                    // and add the virtual point
-                    auto pos = q_points.begin();
-                    q_points.insert(pos, point_t());
-                    geometric = false;
-                }
-                corresp.push_back(0);
-            }
-        }
-        if (geometric) {
-            for (int i = 0; i < corresp.size(); i++) {
-                corresp.at(i) -= 1;
-            }
-        }
-
-        // producing list of simplices
-        std::list<std::list<int>> simp_list;
-
-        for (auto s : simplices.complex_simplex_range()) {
-            std::list<int> s_q;
-            for (int v : simplices.simplex_vertex_range(s)) {
-                s_q.push_back(corresp.at(v));
-            }
-            simp_list.push_back(dedupe_list(s_q));
-        }
-
-        auto quotient_sc = SimplicialComplex<point_t>(q_points, simp_list);
-        quotient_sc.geometric_q = geometric;
-
-        return quotient_sc;
-    }
-
-    /**
-     * @brief  function deciding if a simplicial complex is geometric, note
-     * that
-     * quotient simplicial complexes are not geometric.
-     */
-    bool is_geometric() { return geometric_q; }
-
-   private:
-    bool geometric_q;
-    std::vector<int> dimensions;
-
-    std::vector<matrix_t> boundary_matrices;
-
-    std::set<int> vertex_set(Simplex_handle s) {
-        std::set<int> v_set;
-        for (auto v : simplices.simplex_vertex_range(s)) {
-            v_set.insert(s);
-        }
-        return v_set;
-    }
-
-    // need to do the boundary inclusion crap (you know what I mean)
-    int boundary_index(Simplex_handle s_1, Simplex_handle s_2) {
-        // s1 has to include into s2
-
-        /*
-        int dim1 = simplices.dimension(s_1);
-        int dim2 = simplices.dimension(s_2);
-        decltype(s_1) bdry_h;
-        decltype(s_2) simp_h;
-        if (dim2 == dim1 - 1) {
-            simp_h = decltype(s_1)(s_1);
-            bdry_h = decltype(s_2)(s_2);
-        } else if (dim1 == dim2 - 1) {
-            simp_h = decltype(s_2)(s_2);
-            bdry_h = decltype(s_1)(s_1);
-        } else {
-            std::cerr << "CAN'T BE A BOUNDARY RELATION";
-            throw;
-        }
-
-        bool found = false;
-        int index = dim1;
-
-        // <--- this is not working need to find the parent!!
-        // simp_h = simp_h.parent();
-
-        auto someth = simp_h->second.children()->parent();
-        auto someth2 = someth.first;
-
-        std::cout << "simp_h: " << typeid(simp_h).name() << std::endl;
-        std::cout << "simp_h-fst: " << someth2 << std::endl;
-        std::cout << "someth: " << typeid(someth).name() << std::endl;
-        // std::cout << "someth-fst: " << someth2 << std::endl;
-
-        while (not found) { // <-- this is not working yet
-        if (simp_h == bdry_h) {found = true; index--;} else {index--;simp_h =
-        simp_h->parent;bdry_h=bdry_h->parent;}
-        }
-        */
-
-        auto it_1 = simplices.simplex_vertex_range(s_1).begin();
-        auto end_1 = simplices.simplex_vertex_range(s_1).end();
-
-        auto it_2 = simplices.simplex_vertex_range(s_2).begin();
-        auto end_2 = simplices.simplex_vertex_range(s_2).end();
-
-        int orient = 0;
-        int unmatch = 0;
-
-        for (int p = 0; p <= simplices.dimension(s_2); p++) {
-            if (*it_1 != *it_2) {
-                orient = p++;
-                unmatch++;
-
-            } else {
-                if (it_1 != end_1) {
-                    it_1++;
-                }
-                it_2++;
-            }
-        }
-        return pow(-1, orient);
-    }
-
-    void calculate_matrices() {
-        boundary_matrices = std::vector<matrix_t>();
-        for (int k = 0; k < simplices.dimension(); k++) {
-            boundary_matrices.push_back(
-                matrix_t(dimensions.at(k), dimensions.at(k + 1)));
-        }
-        for (auto s : simplices.complex_simplex_range()) {
-            int j = simplices.key(s);
-            for (auto bs : simplices.boundary_simplex_range(s)) {
-                int i = simplices.key(bs);
-                int k = simplices.dimension(bs);
-                boundary_matrices.at(k).coeffRef(i, j) = boundary_index(bs, s);
-            }
-        }
-    }
-
-    bool is_point(point_t pt) { return not(pt.size() == 0); }
-
-};  // class SimplicialComplex
-
-};      // namespace simplicial
-#endif  // SIMPLICIAL_COMPLEX
diff --git a/lib/scomplex/base_utils.hpp b/lib/scomplex/base_utils.hpp
deleted file mode 100644
index 443666d..0000000
--- a/lib/scomplex/base_utils.hpp
+++ /dev/null
@@ -1,47 +0,0 @@
-#ifndef BASE_UTILS
-#define BASE_UTILS
-
-#include <vector>
-#include <list>
-#include <initializer_list>
-
-std::list<int> dedupe_list(std::list<int>& list) {
-    // O(n log(n))
-    std::set<int> no_reps;
-    std::list<int> no_reps_list;
-    for (int el : list) {
-        no_reps.insert(el);
-    }
-    for (int el : no_reps) {
-        no_reps_list.push_back(el);
-    }
-    return no_reps_list;
-}
-
-std::list<int> dedupe_list(std::vector<int> list) {
-    // O(n log(n))
-    std::set<int> no_reps;
-    std::list<int> no_reps_list;
-    for (int el : list) {
-        no_reps.insert(el);
-    }
-    for (int el : no_reps) {
-        no_reps_list.push_back(el);
-    }
-    return no_reps_list;
-}
-
-std::list<int> dedupe_list(std::initializer_list<int> list) {
-    // O(n log(n))
-    std::set<int> no_reps;
-    std::list<int> no_reps_list;
-    for (int el : list) {
-        no_reps.insert(el);
-    }
-    for (int el : no_reps) {
-        no_reps_list.push_back(el);
-    }
-    return no_reps_list;
-}
-
-#endif
diff --git a/lib/scomplex/chain_calc.hpp b/lib/scomplex/chain_calc.hpp
new file mode 100644
index 0000000..b88f9fc
--- /dev/null
+++ b/lib/scomplex/chain_calc.hpp
@@ -0,0 +1,183 @@
+#pragma once
+
+#include <scomplex/simplicial_complex.hpp>
+#include <scomplex/types.hpp>
+
+#include <Eigen/IterativeLinearSolvers>
+#include <Eigen/Sparse>
+#include <cmath>
+#include <exception>
+#include <tuple>
+#include <vector>
+
+#include <memory>
+#include <boost/geometry.hpp>
+#include <boost/geometry/geometries/point_xy.hpp>
+#include <boost/geometry/geometries/polygon.hpp>
+#include <fstream>
+
+namespace gsimp {
+
+class non_zero_chain : public std::exception {};
+class cant_draw_dimension : public std::exception {};
+
+class bounding_chain {
+    std::shared_ptr<simplicial_complex> s_comp;
+    std::vector<std::unique_ptr<matrix_t>> boundary_matrices;
+    void populate_matrices();
+
+   public:
+    bounding_chain(simplicial_complex& sc);
+    bounding_chain(std::shared_ptr<simplicial_complex> sc);
+    bounding_chain(std::vector<point_t>& points, std::vector<cell_t>& tris);
+    ~bounding_chain();
+    chain_t get_bounding_chain(chain_t&);
+    // void draw_chain(std::string, chain_t);
+};
+
+//---------------------------------------
+// implementation
+
+bounding_chain::~bounding_chain() {}
+
+bounding_chain::bounding_chain(std::shared_ptr<simplicial_complex> sc) {
+    s_comp = sc;
+    populate_matrices();
+}
+
+bounding_chain::bounding_chain(simplicial_complex& sc) {
+    s_comp.reset(new simplicial_complex(sc));
+    populate_matrices();
+}
+
+bounding_chain::bounding_chain(std::vector<point_t>& points,
+                               std::vector<cell_t>& tris) {
+    s_comp.reset(new simplicial_complex(points, tris));
+    populate_matrices();
+}
+
+void bounding_chain::populate_matrices() {
+    for (int d = 0; d < s_comp->dimension(); ++d) {
+        auto level_matrix = s_comp->get_boundary_matrix(d);
+        boundary_matrices.push_back(
+            std::unique_ptr<matrix_t>(new matrix_t(level_matrix)));
+    }
+}
+//
+// round vectors for comparison
+vector_t round_vec(vector_t vec) {
+    vector_t rounded_vec(vec.rows());
+    for (int i = 0; i < vec.rows(); ++i) {
+        int coef = round(vec.coeffRef(i));
+        if (coef != 0) rounded_vec.coeffRef(i) = coef;
+    }
+    return rounded_vec;
+}
+
+bool equals(vector_t vec1, vector_t vec2) {
+    for (int i = 0; i < vec1.rows(); ++i) {
+        if (vec1.coeffRef(i) != vec2.coeffRef(i)) {
+            std::cout << i << ": " << vec1.coeffRef(i) << " "
+                      << vec2.coeffRef(i) << '\n';
+            return false;
+        }
+    }
+    return true;
+}
+
+// // TODO: drowing facilities
+// //  - draw chain
+// //  - draw 1-chain with arrows
+// //  - put it into its own header file
+// void bounding_chain::draw_chain(std::string filename, chain_t chain) {
+//     typedef boost::geometry::model::d2::point_xy<double> point_type;
+//     typedef boost::geometry::model::ring<point_type> polygon_type;
+//
+//     std::vector<polygon_type> polys;
+//
+//     for (auto tri : s_comp->get_level(2)) {
+//         std::vector<point_type> polygon_vec;
+//         for (auto k : tri) {
+//             point_t pt = s_comp->get_points().at(k);
+//             point_type pt_g(pt[0], pt[1]);
+//             polygon_vec.push_back(pt_g);
+//         }
+//         polygon_type poly(polygon_vec.begin(), polygon_vec.end());
+//         polys.push_back(poly);
+//     }
+//
+//     // now we want to draw this
+//     std::ofstream svg(filename.c_str());
+//
+//     boost::geometry::svg_mapper<point_type> mapper(svg, 150, 150);
+//
+//     int chain_d = std::get<0>(chain);
+//     vector_t chain_v = std::get<1>(chain);
+//     for (int i = 0; i < chain_v.rows(); ++i) {
+//         polygon_type poly = polys.at(i);
+//         mapper.add(poly);
+//         if (chain_d == 2) {
+//             if (chain_v.coeffRef(i) == 0)
+//                 mapper.map(
+//                     poly,
+//                     "fill-opacity:0.5;fill:rgb(153,204,0);stroke:rgb(153,"
+//                     "204,0);stroke-width:2");
+//             else
+//                 mapper.map(
+//                     poly,
+//                     "fill-opacity:0.5;fill:rgb(230,104,58);stroke:rgb(230,"
+//                     "104,58);stroke-width:2");
+//         } else {
+//             mapper.map(poly,
+//                        "fill-opacity:0.5;fill:rgb(153,204,0);stroke:rgb(153,"
+//                        "204,0);stroke-width:2");
+//         }
+//
+//         if (chain_d == 1) {
+//             auto lines = s_comp->get_level(1);
+//             for (int i = 0; i < lines.size(); ++i) {
+//                 if (chain_v.coeffRef(i) != 0) {
+//                     std::vector<point_type> line_vec;
+//                     for (auto k : lines.at(i)) {
+//                         point_t pt = s_comp->get_points().at(k);
+//                         point_type pt_g(pt[0], pt[1]);
+//                         line_vec.push_back(pt_g);
+//                     }
+//                     polygon_type line_g(line_vec.begin(), line_vec.end());
+//                     mapper.map(
+//                         line_g,
+//                         "fill-opacity:0.5;fill:rgb(250,0,0);stroke:rgb(250,"
+//                         "0,0);stroke-width:2");
+//                 }
+//             }
+//         }
+//     }
+// }
+
+chain_t bounding_chain::get_bounding_chain(chain_t& chain) {
+    int chain_d;
+    vector_t chain_v;
+    std::tie<int, vector_t>(chain_d, chain_v) = chain;
+
+    if (chain_d >= s_comp->dimension()) throw non_zero_chain();
+
+    Eigen::LeastSquaresConjugateGradient<matrix_t> lscg;
+    lscg.compute(*(boundary_matrices.at(chain_d)));
+
+    vector_t bound_chain(round_vec(lscg.solve(chain_v)));
+    vector_t result(round_vec(*(boundary_matrices.at(chain_d)) * bound_chain));
+
+    if (equals(result, chain_v))
+        return chain_t(chain_d - 1, bound_chain);
+    else {
+        // std::cout << "where to write the svg file?";
+        // std::string filename;
+        // std::cin >> filename;
+        // draw_chain(filename, chain_t(1, bound_chain));
+        // std::cout << '\n';
+        // std::cout << result << '\n';
+        // std::cout << chain_v << '\n';
+        throw non_zero_chain();
+    }
+}
+};
diff --git a/lib/scomplex/coeff_flow.hpp b/lib/scomplex/coeff_flow.hpp
new file mode 100644
index 0000000..e4e613f
--- /dev/null
+++ b/lib/scomplex/coeff_flow.hpp
@@ -0,0 +1,104 @@
+#include <tuple>
+#include <iostream>
+#include <queue>
+#include <scomplex/types.hpp>
+#include <scomplex/simplicial_complex.hpp>
+
+namespace gsimp {
+
+typedef std::tuple<cell_t, cell_t, double> q_elem_t;
+typedef std::queue<q_elem_t> queue_t;
+
+class out_of_context : std::exception {};
+class no_bounding_chain : std::exception {};
+
+chain_t coeff_flow(simplicial_complex& s_comp,  //
+                   chain_t p,                   //
+                   cell_t sigma_0,              //
+                   double c_0) {                //
+
+    if (chain_dim(p) != s_comp.dimension() - 1) throw out_of_context();
+    // 01
+    chain_t c_chain = s_comp.new_chain(s_comp.dimension());
+
+    std::vector<bool> seen_sigma(
+        (size_t)s_comp.get_level_size(s_comp.dimension()), false);
+
+    std::vector<bool> seen_tau(
+        (size_t)s_comp.get_level_size(s_comp.dimension() - 1), false);
+
+    // 04
+    size_t tau_i, sigma_i(s_comp.cell_to_index(sigma_0));
+    chain_val(c_chain, sigma_i) = c_0;
+    seen_sigma.at(sigma_i) = true;
+
+    // 05
+    queue_t queue;
+    for (cell_t tau : s_comp.cell_boundary(sigma_0))
+        queue.emplace(sigma_0, tau, c_0);
+
+    // 08
+    while (not queue.empty()) {
+        cell_t sigma;
+        cell_t tau;
+        double c;
+        std::tie(sigma, tau, c) = queue.front();
+        queue.pop();
+        sigma_i = s_comp.cell_to_index(sigma);
+        tau_i = s_comp.cell_to_index(tau);
+
+        // 10
+        if (seen_sigma.at(sigma_i) && chain_val(c_chain, sigma_i) != c)
+            throw no_bounding_chain();
+        else if (not(seen_tau.at(tau_i) && seen_sigma.at(sigma_i))) {
+            // 14
+            seen_tau.at(tau_i) = true;
+            seen_sigma.at(sigma_i) = true;
+            chain_val(c_chain, sigma_i) = c;
+
+            // 15
+            std::vector<std::pair<int, cell_t>> tau_cofaces;
+            for (auto coface : s_comp.get_cof_and_ind(tau)) {
+                size_t coface_i = s_comp.cell_to_index(std::get<1>(coface));
+                if (coface_i != sigma_i) tau_cofaces.push_back(coface);
+            }
+
+            // 16
+            double predicted_bdry =
+                s_comp.boundary_inclusion_index(tau, sigma) * c;
+            if (tau_cofaces.size() == 0 &&              //
+                predicted_bdry != chain_val(p, tau_i))  //
+                throw no_bounding_chain();              // 18
+            else if (tau_cofaces.size() != 0) {
+                int sigma_p_ind;
+                cell_t sigma_p;
+                std::tie(sigma_p_ind, sigma_p) = tau_cofaces[0];  // 20
+                double c_p = sigma_p_ind*(chain_val(p,tau_i)-s_comp.boundary_inclusion_index(tau, sigma)*c);
+                for (cell_t tau_p : s_comp.cell_boundary(sigma_p))
+                    if (not seen_tau.at(s_comp.cell_to_index(tau_p)))
+                        queue.emplace(sigma_p, tau_p, c_p);
+            }
+        }
+    }
+    return c_chain;
+}
+
+chain_t coeff_flow_embedded(simplicial_complex& s_comp, chain_t p) {
+    if (chain_dim(p) != s_comp.dimension() - 1) throw out_of_context();
+
+    cell_t sigma;
+    double c;
+    auto level = s_comp.get_level(s_comp.dimension() - 1);
+    for (auto tau : s_comp.get_level(s_comp.dimension() - 1)) {
+        if (s_comp.get_cofaces(tau).size() < 2) {
+            int sigma_ind;
+            size_t tau_i = s_comp.cell_to_index(tau);
+            std::tie(sigma_ind, sigma) = s_comp.get_cof_and_ind(tau)[0];
+            c = sigma_ind * chain_val(p, tau_i);
+            chain_t sol = coeff_flow(s_comp, p, sigma, c);
+            return sol;
+        }
+    }
+    throw out_of_context();
+}
+};  // namespace gsimp
diff --git a/lib/scomplex/drawing_utils.hpp b/lib/scomplex/drawing_utils.hpp
new file mode 100644
index 0000000..8bee7fc
--- /dev/null
+++ b/lib/scomplex/drawing_utils.hpp
@@ -0,0 +1,147 @@
+#pragma once
+
+#include<limits>
+#include <cstdio>
+#include <fstream>
+#include <iostream>
+#include <memory>
+
+#include <vector>
+#include <string>
+
+namespace gsimp {
+
+class drawing {
+    struct impl;
+    std::shared_ptr<impl> pimpl;
+    std::string filename, temp_filename;
+   public:
+    drawing(std::string, int, int);
+    drawing &operator=(const drawing &other);
+    void add_arrow_definition();
+    void draw_arrow(std::vector<double>, std::vector<double>, std::string, double);
+    void draw_edge(std::vector<double>, std::vector<double>, std::string, double);
+    void draw_triangle(std::vector<std::vector<double>>, std::string);
+};
+
+struct drawing::impl {
+    std::ofstream file;
+    double maxx = -std::numeric_limits<double>::infinity();
+    double maxy = -std::numeric_limits<double>::infinity();
+    double minx = std::numeric_limits<double>::infinity();
+    double miny = std::numeric_limits<double>::infinity();
+    int size_x, size_y;
+    std::string element_definitions;
+    bool has_arrows;
+    std::string filename;
+    std::string temp_filename;
+
+
+    impl(std::string filename, int sizex, int sizey)
+        : size_x(sizex), size_y(sizey) {
+        this->filename = filename + ".svg";
+        temp_filename = filename + "-temp.svg";
+        file.open(temp_filename);
+    }
+
+    ~impl() {
+        file.close();
+        std::ifstream temp_file(temp_filename);
+        std::ofstream file_finished;
+        file_finished.open(filename);
+
+        file_finished << "<svg xmlns='http://www.w3.org/2000/svg' version='1.1' viewBox='"    //
+                      << 0 << " "            //
+                      << 0 << " "            //
+                      << maxx - minx << " "  //
+                      << maxy - miny << "' height = '500px' width='500px'>\n";
+        file_finished << "<g transform='translate(" << -minx << "," << -miny << ")' >\n";
+        for (std::string file_line; std::getline(temp_file, file_line);) {
+            file_finished << file_line + "\n";
+        }
+        file_finished << "</g>\n";
+        file_finished << "</svg>\n";
+        file_finished.close();
+        std::remove(temp_filename.c_str());
+    }
+
+    void adjust_to_point(std::vector<double> point) {
+        minx = fmin(point[0] , minx);
+        maxx = fmax(point[0] , maxx);
+        miny = fmin(point[1] , miny);
+        maxy = fmax(point[1] , maxy);
+    }
+};
+
+drawing& drawing::operator=(const drawing& other) {
+    pimpl = std::move(other.pimpl);
+    return *this;
+}
+
+drawing::drawing(std::string filename, int sizex, int sizey) {
+    pimpl = std::shared_ptr<impl>(new impl(filename, sizex, sizey));
+}
+
+void drawing::draw_triangle(std::vector<std::vector<double>> coords,
+                            std::string color) {
+    if (coords.size() != 3) throw "A triangle should have 3 vertices";
+
+    std::string tri("<polyline points='");
+    for (auto coord : coords) {
+        double x(coord[0]), y(coord[1]);
+        tri += std::to_string(x) + ", " + std::to_string(y) + " ";
+        pimpl->adjust_to_point(coord);
+    }
+    tri += "' stroke='none' fill='" + color + "'/>\n";
+    pimpl->file << tri;
+}
+
+void drawing::add_arrow_definition() {
+    std::string arrow_defn = "<defs>\n";                               //
+    arrow_defn += "<marker id='head' orient= auto ";                   //
+    arrow_defn += "markerHeight='6' markerWidth='6'",                  //
+        arrow_defn += "refX='0' refY='1.5'>\n";                        //
+    arrow_defn += "<path d='M 0 0 V 3 L 4 1.5 Z' fill='#ff0000'/>\n";  //
+    arrow_defn += "</marker>\n";                                       //
+    arrow_defn += "</defs>\n";
+
+    pimpl->file << arrow_defn;
+}
+
+void drawing::draw_edge(std::vector<double> src, std::vector<double> trg,
+                        std::string color, double stroke) {
+    std::string line_s = "<path stroke-linecap='round' ";  //
+    line_s += "stroke-width='";                                      //
+    line_s += std::to_string(stroke);                               //
+    line_s += "' fill='none' stroke='";                              //
+    line_s += color;                                                //
+    line_s += "' d='M";                                          //
+
+    std::string bgn(""), end(",");
+    for (int i = 0; i < src.size(); ++i) {
+        bgn += std::to_string(src[i]) + " ";
+        end += std::to_string(trg[i]) + " ";
+    }
+    line_s += bgn + end + "'/>\n";
+    pimpl->file << line_s;
+}
+
+void drawing::draw_arrow(std::vector<double> src, std::vector<double> trg,
+                         std::string color, double stroke) {
+    std::string arrow =                                              //
+        "<path id = 'arrow' stroke-linecap='round' stroke-width='";  //
+    arrow += std::to_string(stroke);                                 //
+    arrow += "' fill='none' stroke=";                                //
+    arrow += color;                                                  //
+    arrow += " d = 'M ";                                             //
+
+    std::string bgn(" "), mid(" L "), end(" L ");
+    for (int i = 0; i < src.size(); ++i) {
+        bgn += std::to_string(src[i]) + " ";
+        mid += std::to_string(trg[i] - src[i] / 2) + " ";
+        end += std::to_string(trg[i]) + "' ";
+    }
+    arrow += bgn + mid + end + " />";
+    (pimpl->file) << arrow;
+}
+};
diff --git a/lib/scomplex/graph_utils.hpp b/lib/scomplex/graph_utils.hpp
index 3cfd5f1..5ce5c26 100644
--- a/lib/scomplex/graph_utils.hpp
+++ b/lib/scomplex/graph_utils.hpp
@@ -1,13 +1,8 @@
 #pragma once
 
 #include <vector>
-// reverse
-#include <algorithm>
-
-// to compute norms
-#include "Eigen/Eigen"
-
-#include "scomplex/Simplicial_Complex.h"
+#include <scomplex/types.hpp>
+#include <scomplex/simplicial_complex.hpp>
 
 // graph libraries
 #include "boost/graph/adjacency_list.hpp"
@@ -15,97 +10,100 @@
 #include "boost/graph/named_function_params.hpp"
 #include "boost/property_map/vector_property_map.hpp"
 
+namespace gsimp {
+
+/*
+typedefs:
+    graph_t
+functions:
+    calculate_one_skelleton_graph(simplicial_complex& s_comp )
+    complete_path(const graph_t& g, std::vector<vertex_t> vec)
+    shortest_path(const graph_t& g, vertex_t s, vertex_t t)
+*/
+
 typedef typename boost::adjacency_list<             //
     boost::vecS,                                    //
     boost::vecS,                                    //
     boost::undirectedS,                             //
-    boost::no_property,                             //
+    boost::no_property,                             // vertex property
     boost::property<boost::edge_weight_t, double>>  //
-    Graph;
+    graph_t;
 
-typedef typename boost::graph_traits<Graph>::vertex_descriptor vertex_t;
+typedef typename boost::graph_traits<graph_t>::vertex_descriptor vertex_t;
 
-typedef std::vector<double> point_t;
-
-Graph calculate_one_skelleton_graph(
-    simplicial::SimplicialComplex<point_t>& scomp  //
-    ) {
-    int num_points{scomp.get_dimension(0)};
-    int num_edges{scomp.get_dimension(1)};
+graph_t calculate_one_skelleton_graph(simplicial_complex& s_comp  //
+                                      ) {
+    size_t num_points(s_comp.get_level_size(0));
+    size_t num_edges(s_comp.get_level_size(1));
 
     std::vector<std::pair<int, int>> edges(num_edges);
     std::vector<double> weights(num_edges);
 
-    for (std::vector<int> edge : scomp.get_level(1)) {
+    for (cell_t edge : s_comp.get_level(1)) {
         int src_index = edge.at(0);
         int trg_index = edge.at(1);
 
-        point_t src_point = scomp.points.at(edge.at(0));
-        point_t trg_point = scomp.points.at(edge.at(1));
+        point_t src_point = s_comp.get_points().at(edge.at(0));
+        point_t trg_point = s_comp.get_points().at(edge.at(1));
 
         edges.push_back(std::make_pair(src_index, trg_index));
 
-        Eigen::VectorXd diff(src_point.size());
-
-        for (int i = 0; i < src_point.size(); i++)
-            diff(i) = src_point.at(i) - trg_point.at(i);
-
-        weights.push_back(diff.norm());
+        double norm = 0;
+        for (int i = 0; i < src_point.size(); ++i)
+            norm += pow(src_point.at(i) - trg_point.at(i), 2);
+        weights.push_back(sqrt(norm));
     }
 
-    const Graph one_skelleton(edges.begin(),    //
-                              edges.end(),      //
-                              weights.begin(),  //
-                              num_points);      //
-
+    const graph_t one_skelleton(edges.begin(),    //
+                                edges.end(),      //
+                                weights.begin(),  //
+                                num_points);      //
     return one_skelleton;
 }
 
-// calculate shortest paths:
-std::vector<vertex_t> shortest_path(const Graph& g, vertex_t s, vertex_t t) {
+/**
+* @brief  find the shortest path between two vertices in a graph.
+*
+* @param g: graph (passed by ref)
+* @param s: source vertex (type: decltype(g)::vertex_descriptor)
+* @param t: source vertex (type: decltype(g)::vertex_descriptor)
+*
+* @returns: a vector of vertices
+*            (type: std::vector<decltype(g)::vertex_descriptor>)
+*/
+std::vector<vertex_t> shortest_path(const graph_t& g, vertex_t s, vertex_t t) {
+    // make a predecessor map
     boost::vector_property_map<vertex_t> predecessors(boost::num_vertices(g));
-
     boost::dijkstra_shortest_paths(g, s,
                                    boost::predecessor_map(&predecessors[0]));
-
+    // construct the path (going backwards from the target to the source)
     vertex_t it = t;
     std::vector<vertex_t> s_t_path{};
-    while (it != s) {
+    do {
         s_t_path.push_back(it);
         it = predecessors[it];
-    }
-    s_t_path.push_back(it);
+    } while (it != s);
     std::reverse(s_t_path.begin(), s_t_path.end());
     return s_t_path;
-
-    // std::cout << "I'm calculating shortest paths\n";
 }
 
-std::vector<vertex_t> complete_path(const Graph& g, std::vector<vertex_t> vec) {
+std::vector<vertex_t> complete_path(const graph_t& g,
+                                    std::vector<vertex_t> vec) {
     std::vector<vertex_t>::iterator s{vec.begin()};
-    std::list<vertex_t> full_path{*s};
-    while (s != vec.end()) {
-        std::vector<vertex_t>::iterator t = std::next(s);
-        // calculate the path between s and t
+    std::vector<vertex_t>::iterator t = std::next(s);
+
+    // start calculatinng the full path
+    std::vector<vertex_t> full_path;
+    full_path.push_back(*s);
+    while (t != vec.end()) {
+        // get the next vertex calculate and add the path between s and t
+        // (contains vertices other than s) to the full path
         std::vector<vertex_t> path_portion{shortest_path(g, *s, *t)};
-
-        // add the path between s and t to the full path
-        for (vertex_t vert : path_portion) {
-            // allways skip the first element since we are calculating the paths
-            // between s,t then s',t' where s' = t and that would make use of a
-            // non-existing edge.
-            if (vert != *path_portion.begin()) {
-                full_path.push_back(vert);
-            }
-        }
-
+        // quick vector merging
+        full_path.insert(full_path.end(),path_portion.begin(),path_portion.end());
         std::advance(s, 1);
+        std::advance(t, 1);
     }
-    // initialize with size so we don't have to relocate
-    std::vector<vertex_t> full_path_vect(full_path.size());
-    for (vertex_t vert : full_path) {
-        full_path_vect.push_back(vert);
-    }
-
-    return full_path_vect;
+    return full_path;
 }
+};
diff --git a/lib/scomplex/nn_utils.hpp b/lib/scomplex/nn_utils.hpp
index fb42161..9dd8676 100644
--- a/lib/scomplex/nn_utils.hpp
+++ b/lib/scomplex/nn_utils.hpp
@@ -1,46 +1,98 @@
 #pragma once
 
+#include <scomplex/types.hpp>
+
 #include <vector>
 #include <algorithm>
 
+// TODO: go through the includes and find out which ones are needed
 #include <CGAL/Kd_tree.h>
-
-// #include <CGAL/Simple_cartesian.h>
 #include <CGAL/Search_traits_d.h>
+#include <CGAL/Search_traits_adapter.h>
 #include <CGAL/Cartesian_d.h>
 #include <CGAL/K_neighbor_search.h>
-// #include <CGAL/Euclidean_distance.h>
 #include <CGAL/Orthogonal_k_neighbor_search.h>
+#include <boost/tuple/tuple.hpp>
 
-typedef std::vector<double> point_t;
+namespace gsimp {
 
-typedef CGAL::Cartesian_d<double> cart;
-typedef cart::Point_d point_d;
-typedef CGAL::Search_traits_d<cart> traits;
-typedef CGAL::Orthogonal_k_neighbor_search<traits> neighbor_search_t;
-typedef neighbor_search_t::Tree tree_t;
+/*
+types used outside:
+    tree_t
+functions provided:
+    make_tree(tree_t& tree, std::vector<point_t>& point_list);
+    nearest_neighbor(tree_t& tree, point_t point);
+    nearest_neighbor_index(tree_t& tree, point_t point);
+    snap_points(tree_t& tree, std::vector<point_t> points);
+    snap_points_to_indexes( tree_t& tree, std::vector<point_t> points);
+*/
+
+// kd-tree point type, d-dimensional point
+typedef CGAL::Cartesian_d<double>::Point_d point_d;
 
-// typedef CGAL::Kd_tree<traits> kd_tree;
-// typedef CGAL::K_neighbor_search<traits> k_neighbor;
-// typedef CGAL::Euclidean_distance<traits> euclidean;
+typedef CGAL::Search_traits_d<CGAL::Cartesian_d<double>> cartesian_traits;
 
-// template <int Dimension>
-// typedef cart::Point_d<Dimension> point_d;
+typedef CGAL::Search_traits_adapter<     // define the search tratis
+    boost::tuple<point_d, size_t>,       // pairing d-point with its index
+    CGAL::Nth_of_tuple_property_map<     // custom kernel
+        0,                               // position of the points to compare
+        boost::tuple<point_d, size_t>>,  // type tha will be used
+    cartesian_traits>                    // traits to use
+    traits;                              // name
 
-void make_tree(tree_t& tree, std::list<point_t>& point_list) {
-    tree.reserve(point_list.size());
-    std::list<point_d> insertable_points;
-    for (auto pt : point_list) {
+// type of search tree
+typedef CGAL::Orthogonal_k_neighbor_search<traits> neighbor_search_t;
+
+typedef neighbor_search_t::Tree tree_t;
+
+void make_tree(tree_t& tree, std::vector<point_t>& point_list) {
+    // tree.reserve(point_list.size());
+    for (size_t ind = 0; ind < point_list.size(); ++ind) {
+        auto pt = point_list.at(ind);
         point_d pt_d(pt.size(), pt.begin(), pt.end());
-        tree.insert(pt_d);
+        tree.insert(boost::make_tuple(pt_d, ind));
     }
-    //    return tree;
-};
+}
 
-point_t nearest_neighbour(tree_t& tree, point_t point) {
+// given a point snap it to the tree
+point_t nearest_neighbor(tree_t& tree, point_t point) {
     point_d pt(point.size(), point.begin(), point.end());
+    // "1" refers to the number of nearest neighbors to find
     neighbor_search_t search(tree, pt, 1);
-    point_d result_pt = search.begin()->first;
+    point_d result_pt;
+    size_t ind;
+    boost::tie<point_d, size_t>(result_pt, ind) = search.begin()->first;
     point_t point_vec(result_pt.cartesian_begin(), result_pt.cartesian_end());
     return point_vec;
 }
+
+// given a vector of points snap them to the tree
+std::vector<point_t> snap_points(tree_t& tree, std::vector<point_t> points) {
+    std::vector<point_t> snapped_points;
+    for (auto point : points) {
+        snapped_points.push_back(nearest_neighbor(tree, point));
+    }
+    return snapped_points;
+}
+
+size_t nearest_neighbor_index(tree_t& tree, point_t point) {
+    point_d pt(point.size(), point.begin(), point.end());
+    // "1" refers to the number of nearest neighbors to find
+    neighbor_search_t search(tree, pt, 1);
+    point_d result_pt;
+    size_t ind;
+    boost::tie<point_d, size_t>(result_pt, ind) = search.begin()->first;
+    return ind;
+}
+
+std::vector<size_t> snap_points_to_indexes(  //
+    tree_t& tree,                            //
+    std::vector<point_t> points) {
+    std::vector<size_t> snapped_points;
+    for (auto point : points) {
+        snapped_points.push_back(  //
+            nearest_neighbor_index(tree, point));
+    }
+    return snapped_points;
+}
+};
diff --git a/lib/scomplex/path_snapper.cpp b/lib/scomplex/path_snapper.cpp
new file mode 100644
index 0000000..7687e2a
--- /dev/null
+++ b/lib/scomplex/path_snapper.cpp
@@ -0,0 +1,142 @@
+#include <scomplex/path_snapper.hpp>
+#include <iostream>
+
+#include <Eigen/Sparse>
+#include <scomplex/types.hpp>
+#include <scomplex/simplicial_complex.hpp>
+#include <scomplex/nn_utils.hpp>
+#include <scomplex/graph_utils.hpp>
+
+namespace gsimp {
+
+struct path_snapper::impl {
+    tree_t point_tree;     // defined in nn_utils.hpp
+    graph_t vertex_graph;  // defined in graph_utils.hpp
+    std::shared_ptr<simplicial_complex> s_comp;
+
+    impl(std::shared_ptr<simplicial_complex> p_sc) {
+        s_comp = p_sc;
+        auto points = s_comp->get_points();
+        make_tree(point_tree, points);
+        vertex_graph = calculate_one_skelleton_graph(*s_comp);
+    };
+
+    impl(simplicial_complex& sc) {
+        s_comp.reset(new simplicial_complex(sc));
+        auto points = s_comp->get_points();
+        make_tree(point_tree, points);
+        vertex_graph = calculate_one_skelleton_graph(*s_comp);
+    }
+
+    impl(std::vector<point_t>& pts, std::vector<cell_t>& cells) {
+        s_comp = std::shared_ptr<simplicial_complex>(
+            new simplicial_complex(pts, cells));
+        make_tree(point_tree, pts);
+        vertex_graph = calculate_one_skelleton_graph(*s_comp);
+    }
+
+    ~impl(){};
+
+    std::vector<size_t> snap_path(std::vector<point_t> path) {
+        auto way_points = snap_points_to_indexes(point_tree, path);
+        auto snapped_path = complete_path(vertex_graph, way_points);
+        return snapped_path;
+    }
+
+    std::vector<std::pair<cell_t, int>> index_pairs(std::vector<point_t> path) {
+        auto vertex_path = snap_path(path);
+        std::vector<std::pair<cell_t, int>> pair_seq;
+        for (auto it = vertex_path.begin(); std::next(it) != vertex_path.end();
+             ++it) {
+            //
+            //
+            size_t src(*it), trg(*(std::next(it)));            //
+            if (src < trg) {                                   //
+                pair_seq.push_back(                            //
+                    std::make_pair<std::vector<size_t>, int>(  //
+                        std::vector<size_t>{src, trg}, 1));    //
+            } else {                                           //
+                pair_seq.push_back(                            //
+                    std::make_pair<std::vector<size_t>, int>(  //
+                        std::vector<size_t>{trg, src}, -1));   //
+            }                                                  //
+        }
+        return pair_seq;
+    }
+};
+
+path_snapper::path_snapper(std::shared_ptr<simplicial_complex> sc) {
+    p_impl = std::shared_ptr<impl>(new impl(sc));
+};
+
+path_snapper::path_snapper(simplicial_complex& sc) {
+    p_impl = std::shared_ptr<impl>(new impl(sc));
+}
+
+path_snapper::path_snapper(std::vector<point_t>& pts,
+                           std::vector<cell_t>& cells) {
+    p_impl = std::shared_ptr<impl>(new impl(pts, cells));
+}
+
+path_snapper::~path_snapper() {}
+
+path_snapper::path_snapper(path_snapper& other) { p_impl = other.p_impl; }
+
+path_snapper& path_snapper::operator=(const path_snapper& other) {
+    p_impl = other.p_impl;
+    return *this;
+}
+
+std::vector<point_t> path_snapper::snap_path_to_points(
+    std::vector<point_t> path) {
+    auto index_path = p_impl->snap_path(path);
+    std::vector<point_t> point_path;
+    for (size_t p : index_path)
+        point_path.push_back(p_impl->s_comp->get_points().at(p));
+    return point_path;
+}
+
+std::vector<size_t> path_snapper::snap_path_to_indices(
+    std::vector<point_t> path) {
+    return p_impl->snap_path(path);
+}
+
+chain_t path_snapper::snap_path_to_chain(std::vector<point_t> path) {
+    auto pair_seq = p_impl->index_pairs(path);
+    chain_t rep = p_impl->s_comp->new_chain(1);
+    for (auto pair : pair_seq) {
+        size_t index = p_impl->s_comp->cell_to_index(std::get<0>(pair));
+        chain_val(rep, index) += std::get<1>(pair);
+    }
+    return rep;
+}
+
+std::vector<point_t> path_snapper::index_sequence_to_point(
+    std::vector<size_t> ind_path) {
+    std::vector<point_t> pt_path;
+    for (size_t ind : ind_path)
+        pt_path.push_back(p_impl->s_comp->get_points().at(ind));
+    return pt_path;
+}
+
+std::vector<size_t> path_snapper::point_sequence_to_index(
+    std::vector<point_t> pt_path) {
+    std::vector<size_t> ind_path =
+        snap_points_to_indexes(p_impl->point_tree, pt_path);
+    return ind_path;
+}
+
+chain_t path_snapper::index_sequence_to_chain(std::vector<size_t> ind_path) {
+    auto it = ind_path.begin();
+    chain_t chain = p_impl->s_comp->new_chain(1);
+    for (; std::next(it) != ind_path.end(); ++it) {
+        auto ind = p_impl->s_comp->cell_to_index({*it, *std::next(it)});
+        chain_val(chain, ind) += (*it < *std::next(it)) ? 1 : -1;
+    }
+    return chain;
+}
+
+chain_t path_snapper::point_sequence_to_chain(std::vector<point_t> pt_path) {
+    return index_sequence_to_chain(point_sequence_to_index(pt_path));
+}
+};
diff --git a/lib/scomplex/path_snapper.hpp b/lib/scomplex/path_snapper.hpp
new file mode 100644
index 0000000..43025ea
--- /dev/null
+++ b/lib/scomplex/path_snapper.hpp
@@ -0,0 +1,31 @@
+#pragma once
+
+#include <scomplex/simplicial_complex.hpp>
+#include <scomplex/types.hpp>
+
+namespace gsimp {
+
+class path_snapper {
+   private:
+    struct impl;
+    std::shared_ptr<impl> p_impl;
+
+   public:
+    // constructors and such
+    path_snapper(std::vector<point_t>&, std::vector<cell_t>&);
+    path_snapper(std::shared_ptr<simplicial_complex>);
+    path_snapper(simplicial_complex&);
+    ~path_snapper();
+    path_snapper(path_snapper&);
+    path_snapper& operator=(const path_snapper&);
+    // the things we want to do
+    std::vector<point_t> snap_path_to_points(std::vector<point_t>);
+    std::vector<size_t> snap_path_to_indices(std::vector<point_t>);
+    chain_t snap_path_to_chain(std::vector<point_t>);
+    // interconversion
+    std::vector<point_t> index_sequence_to_point(std::vector<size_t>);
+    std::vector<size_t> point_sequence_to_index(std::vector<point_t>);
+    chain_t index_sequence_to_chain(std::vector<size_t>);
+    chain_t point_sequence_to_chain(std::vector<point_t>);
+};
+};
diff --git a/lib/scomplex/qhull_parsing.hpp b/lib/scomplex/qhull_parsing.hpp
index 7a88b97..446b2db 100644
--- a/lib/scomplex/qhull_parsing.hpp
+++ b/lib/scomplex/qhull_parsing.hpp
@@ -1,4 +1,3 @@
-
 #pragma once
 
 #include <fstream>
@@ -12,40 +11,37 @@
 #include <tuple>
 
 typedef std::vector<double> point_t;
-typedef std::vector<int> cell_t;
+typedef std::vector<size_t> cell_t;
 
 template <typename data_t>
 std::vector<data_t> tokenize(std::string str) {
     std::istringstream sstream(str);
-    std::vector<data_t> tokens{std::istream_iterator<data_t>{sstream},
+    std::vector<data_t> tokens{//
+                               std::istream_iterator<data_t>{sstream},
                                std::istream_iterator<data_t>{}};
     return tokens;
 }
 
-template <typename data_t>
-void output_vector(std::vector<data_t> vec) {
-    std::copy(vec.begin(), vec.end(),
-              std::ostream_iterator<data_t>(std::cout, " "));
+template <typename Iterator>
+void cout_iterator(Iterator first, Iterator last) {
+    std::copy(first, last,
+              std::ostream_iterator<
+                  typename std::iterator_traits<Iterator>::value_type>(
+                  std::cout, " "));
     std::cout << '\n';
 }
 
 std::pair<std::vector<point_t>, std::vector<cell_t>> parse_qhull_file(
     std::string filename) {
     std::ifstream file_stream(filename);
-    std::string line;
 
-    int current_line_number{0};
-    int dimensionality;
-    //
-    int num_points;
-    int num_cells;
-    //
-    int cur_point_num;
-    int cur_cell_num;
+    if (!file_stream.is_open()) {
+        throw std::runtime_error("failed to open \"" + filename + "\"\n");
+    }
 
-    int line_number{0};
-    int total_points{0};
-    int total_cells{0};
+    std::string line;
+
+    int dimensionality, total_points{0}, total_cells{0};
 
     std::vector<point_t> points;
     std::vector<cell_t> cells;
@@ -55,17 +51,13 @@ std::pair<std::vector<point_t>, std::vector<cell_t>> parse_qhull_file(
         std::getline(file_stream, line);
         // case of the first line of the file (get dimensions)
         if (line_number == 0) {
-            auto tokens = tokenize<int>(line);
+            auto tokens = tokenize<size_t>(line);
             dimensionality = tokens.at(0);
-            // std::cout << "found first line! dimensions == " << dimensionality
-            //          << '\n';
         }
         // second line of the file (get number of points)
         if (line_number == 1) {
-            auto tokens = tokenize<int>(line);
+            auto tokens = tokenize<size_t>(line);
             total_points = tokens.at(0);
-            // std::cout << "found second line! number of points == "
-            //          << total_points << '\n';
         }
         // lines pertaining to points
         if (line_number > 1 && line_number < total_points + 2) {
@@ -75,39 +67,16 @@ std::pair<std::vector<point_t>, std::vector<cell_t>> parse_qhull_file(
         }
         // first line after points (get number of cells)
         if (line_number == total_points + 2) {
-            auto tokens = tokenize<int>(line);
+            auto tokens = tokenize<size_t>(line);
             total_cells = tokens.at(0);
-            // std::cout << "found cell count! number of cells == " <<
-            // total_cells
-            //          << '\n';
         }
         // lines pertaining to cells
         if (line_number > total_points + 2 &&
             line_number < total_points + total_cells + 3) {
-            auto cell = tokenize<int>(line);
+            auto cell = tokenize<size_t>(line);
             cells.push_back(cell);
         }
     }
+    file_stream.close();
     return std::make_pair(points, cells);
 }
-
-/*
- * usage example
- *
-int main() {
-    std::vector<point_t> points;
-    std::vector<cell_t> cells;
-    std::tie(points, cells) = parse_qhull_file("qh-test.dat");
-
-    std::cout << points.size() << '\n';
-    for (point_t point : points) {
-        output_vector<double>(point);
-    }
-    std::cout << cells.size() << '\n';
-    for (cell_t cell : cells) {
-        output_vector<int>(cell);
-    }
-
-    return 0;
-}
-*/
diff --git a/lib/scomplex/simplicial_complex.cpp b/lib/scomplex/simplicial_complex.cpp
new file mode 100644
index 0000000..19ef127
--- /dev/null
+++ b/lib/scomplex/simplicial_complex.cpp
@@ -0,0 +1,366 @@
+#include <scomplex/simplicial_complex.hpp>
+#include <scomplex/types.hpp>
+
+#include <iterator>  // for debuging purposes
+
+#include <gudhi/Simplex_tree.h>
+#include <Eigen/Sparse>
+
+#include <tuple>
+#include <memory>  // smart pointers
+#include <functional>
+#include <cmath>
+
+namespace gsimp {
+
+struct simplicial_complex::impl {
+    // auxiliary types
+    struct SimpleOptions : Gudhi::Simplex_tree_options_full_featured {
+        typedef size_t Vertex_handle;
+    };
+    typedef Gudhi::Simplex_tree<SimpleOptions> simp_tree;
+    typedef Gudhi::Simplex_tree<SimpleOptions>::Simplex_handle simp_handle;
+    typedef std::vector<simp_handle*> level_t;
+    typedef std::vector<std::unique_ptr<level_t>> levels_t;
+
+    // member variables
+    std::vector<point_t> points;
+    simp_tree simplices;
+    bool quotient_q;
+    std::vector<matrix_t> boundary_matrices;
+    levels_t levels;
+
+    impl() : quotient_q(false){};
+
+    impl(std::vector<point_t>& arg_points, std::vector<cell_t>& arg_tris)
+        : points(arg_points), quotient_q(false) {
+        // create the simplex tree
+        for (auto tri : arg_tris) {
+            simplices.insert_simplex_and_subfaces(dedupe_vec(tri));
+            int d = tri.size() - 1;
+            if (simplices.dimension() < d) simplices.set_dimension(d);
+        }
+
+        // assign a key to each simplex in each level
+        std::vector<size_t> count(simplices.dimension() + 1, 0);
+        for (int i = 0; i < simplices.dimension() + 1; ++i) {
+            level_t* level = new level_t();
+            levels.push_back(std::unique_ptr<level_t>(level));
+        }
+
+        auto simplex_range = simplices.complex_simplex_range();
+        for (auto s : simplex_range) {
+            int d = simplices.dimension(s);
+            simplices.assign_key(s, count.at(d)++);
+            levels.at(d)->push_back(new simp_handle(s));
+        }
+    }
+
+    ~impl() {
+        // get rid of the levels so they won't dangle
+        for (int i = levels.size(); i <= 0; --i) levels.at(i).reset();
+    };
+
+    size_t get_level_size(int level) { return levels.at(level)->size(); }
+
+    // calculate the index of s_1 in the boundary of s_2
+    int boundary_index(simp_handle s_1, simp_handle s_2) {
+        auto it_1 = simplices.simplex_vertex_range(s_1).begin();
+        auto end_1 = simplices.simplex_vertex_range(s_1).end();
+        //
+        auto it_2 = simplices.simplex_vertex_range(s_2).begin();
+        //
+        int orient = 0;
+        int unmatch = 0;
+        //
+        for (int p = 0; p <= simplices.dimension(s_2); p++) {
+            if (*it_1 != *it_2) {
+                orient = p++;
+                unmatch++;
+            } else {
+                if (it_1 != end_1) it_1++;
+                it_2++;
+            }
+        }
+        return pow(-1, orient);
+    }
+
+    std::vector<size_t> dedupe_vec(std::vector<size_t>& vec) {
+        std::set<size_t> no_reps;
+        std::vector<size_t> no_reps_list;
+        for (size_t el : vec) {
+            no_reps.insert(el);
+        }
+        for (size_t el : no_reps) {
+            no_reps_list.push_back(el);
+        }
+        return no_reps_list;
+    }
+
+    void calculate_matrices() {
+        boundary_matrices = std::vector<matrix_t>();
+        for (int k = 0; k < simplices.dimension(); k++) {
+            boundary_matrices.push_back(
+                matrix_t(get_level_size(k), get_level_size(k + 1)));
+        }
+        for (auto s : simplices.complex_simplex_range()) {
+            int j = simplices.key(s);
+            for (auto bs : simplices.boundary_simplex_range(s)) {
+                int i = simplices.key(bs);
+                int k = simplices.dimension(bs);
+                boundary_matrices.at(k).coeffRef(i, j) = boundary_index(bs, s);
+            }
+        }
+    }
+
+    std::vector<cell_t> get_level(int level) {
+        std::vector<cell_t> level_cells;
+        for (auto simp : *levels.at(level)) {
+            cell_t v_simp;
+            for (auto v : simplices.simplex_vertex_range(*simp)) {
+                v_simp.push_back(v);
+            }
+            level_cells.push_back(v_simp);
+        }
+        return level_cells;
+    }
+
+    simp_handle index_to_handle(int d, size_t tau) {
+        return *(levels.at(d)->at(tau));
+    }
+
+    size_t handle_to_index(simp_handle tau) { return simplices.key(tau); }
+
+    simp_handle cell_to_handle(cell_t tau) {
+        auto sh = simplices.find(tau);
+        return sh;
+    }
+
+    cell_t handle_to_cell(simp_handle tau) {
+        cell_t cell;
+        for (auto v : simplices.simplex_vertex_range(tau)) cell.push_back(v);
+        return cell;
+    }
+};  // struct impl
+
+std::vector<std::pair<int, cell_t>> simplicial_complex::get_bdry_and_ind(
+    cell_t cell) {
+    impl::simp_handle simp = p_impl->cell_to_handle(cell);
+    auto c_boundary = p_impl->simplices.boundary_simplex_range(simp);
+    std::vector<std::pair<int, cell_t>> boundary_and_indices;
+    for (auto face : c_boundary)
+        boundary_and_indices.push_back(              //
+            std::make_pair<int, cell_t>(             //
+                p_impl->boundary_index(face, simp),  //
+                p_impl->handle_to_cell(face)));      //
+    return boundary_and_indices;
+};
+
+
+
+std::vector<std::pair<int, size_t>> simplicial_complex::get_bdry_and_ind_index(
+    int d, size_t cell) {
+    impl::simp_handle simp = p_impl->index_to_handle(d, cell);
+    auto c_boundary = p_impl->simplices.boundary_simplex_range(simp);
+    std::vector<std::pair<int, size_t>> boundary_and_indices;
+    for (auto face : c_boundary)
+        boundary_and_indices.push_back(              //
+            std::make_pair<int, size_t>(             //
+                p_impl->boundary_index(face, simp),  //
+                p_impl->handle_to_index(face)));     //
+    return boundary_and_indices;
+};
+
+std::vector<size_t> simplicial_complex::cell_boundary_index(int d,
+                                                            size_t cell) {
+    impl::simp_handle simp = p_impl->index_to_handle(d, cell);
+    auto c_boundary = p_impl->simplices.boundary_simplex_range(simp);
+    std::vector<size_t> s_boundary;
+    for (auto c : c_boundary) s_boundary.push_back(p_impl->handle_to_index(c));
+    return s_boundary;
+}
+
+std::vector<cell_t> simplicial_complex::cell_boundary(cell_t cell) {
+    impl::simp_handle simp = p_impl->cell_to_handle(cell);
+    auto c_boundary = p_impl->simplices.boundary_simplex_range(simp);
+    std::vector<cell_t> s_boundary;
+    for (auto c : c_boundary) s_boundary.push_back(p_impl->handle_to_cell(c));
+    return s_boundary;
+}
+
+int simplicial_complex::boundary_inclusion_index(cell_t c1, cell_t c2) {
+    return p_impl->boundary_index(p_impl->cell_to_handle(c1),   //
+                                  p_impl->cell_to_handle(c2));  //
+};
+int simplicial_complex::boundary_inclusion_index(int d1, size_t s1,    //
+                                                 int d2, size_t s2) {  //
+    cell_t c1 = index_to_cell(d1, s1);
+    cell_t c2 = index_to_cell(d2, s2);
+    return p_impl->boundary_index(p_impl->cell_to_handle(c1),
+                                  p_impl->cell_to_handle(c2));
+};
+
+std::vector<std::pair<int, cell_t>> simplicial_complex::get_cof_and_ind(
+    cell_t cell) {
+    std::vector<std::pair<int, cell_t>> c_cofaces;
+    for (auto face : get_cofaces(cell))
+        c_cofaces.push_back(                                   //
+            std::make_pair(                                    //
+                boundary_inclusion_index(cell, face), face));  //
+    return c_cofaces;
+}
+
+std::vector<std::pair<int, size_t>> simplicial_complex::get_cof_and_ind_index(
+    int d, size_t c) {
+    std::vector<std::pair<int, size_t>> c_cofaces;
+    for (auto face : get_cofaces_index(d, c))
+        c_cofaces.push_back(                                          //
+            std::make_pair(                                           //
+                boundary_inclusion_index(d, c, d + 1, face), face));  //
+    return c_cofaces;
+}
+
+int simplicial_complex::get_level_size(int level) {
+    return p_impl->get_level_size(level);
+}
+
+simplicial_complex::simplicial_complex(std::vector<point_t>& arg_points,
+                                       std::vector<cell_t>& arg_tris) {
+    p_impl = std::shared_ptr<impl>(new impl(arg_points, arg_tris));
+}
+
+simplicial_complex::simplicial_complex(const simplicial_complex& other) {
+    p_impl = other.p_impl;
+}
+
+simplicial_complex& simplicial_complex::operator=(
+    const simplicial_complex& other) {
+    p_impl = other.p_impl;
+    return *this;
+}
+
+simplicial_complex::~simplicial_complex() {}
+
+std::vector<point_t> simplicial_complex::get_points() { return p_impl->points; }
+
+std::vector<cell_t> simplicial_complex::get_level(int level) {
+    return p_impl->get_level(level);
+}
+
+matrix_t simplicial_complex::get_boundary_matrix(int d) {
+    // uninstantiated boundary matrices
+    if (p_impl->boundary_matrices.size() == 0) p_impl->calculate_matrices();
+
+    // now they have to be instantiated, get them
+    if (0 <= d && d < p_impl->boundary_matrices.size())
+        return p_impl->boundary_matrices.at(d);
+    else
+        throw NoChain();
+}
+
+int simplicial_complex::dimension() { return p_impl->simplices.dimension(); }
+
+cell_t simplicial_complex::index_to_cell(int d, size_t ind) {
+    auto sh = p_impl->levels.at(d)->at(ind);
+    return p_impl->handle_to_cell(*sh);
+}
+
+size_t simplicial_complex::cell_to_index(cell_t simp) {
+    auto sh = p_impl->simplices.find(simp);
+    return p_impl->simplices.key(sh);
+}
+
+bool simplicial_complex::is_quotient() { return p_impl->quotient_q; }
+
+std::vector<size_t> simplicial_complex::get_cofaces_index(int d, size_t face) {
+    std::vector<size_t> s_cofaces;
+    impl::simp_handle face_h = *(p_impl->levels.at(d)->at(face));
+    cell_t my_cell = p_impl->handle_to_cell(face_h);
+    // codimension 1 faces
+    auto range = p_impl->simplices.cofaces_simplex_range(face_h, 1);
+    for (auto tau : range) s_cofaces.push_back(p_impl->simplices.key(tau));
+    return s_cofaces;
+}
+
+std::vector<cell_t> simplicial_complex::get_cofaces(cell_t face) {
+    std::vector<cell_t> s_cofaces;
+    auto face_h = p_impl->cell_to_handle(face);
+    // codimension 1 faces
+    auto range = p_impl->simplices.cofaces_simplex_range(face_h, 1);
+    for (auto tau : range) s_cofaces.push_back(p_impl->handle_to_cell(tau));
+    return s_cofaces;
+}
+
+chain_t simplicial_complex::new_chain(int d) {
+    vector_t v(get_level_size(d));
+    return chain_t(d,v);
+}
+
+
+// TODO: get top dimensional cells from the simplicial complex to add to the
+// simplex tree (instead of all the simplices).
+simplicial_complex simplicial_complex::quotient(int char_fun(point_t)) {
+    int n_points = 1;
+    auto corresp = std::vector<size_t>();
+
+    // making the correspondence between points in the original complex
+    // and points in the quotient. point index 0 corresponds to the
+    // "virtual" point (in case it exists).
+    std::vector<point_t> points;
+    bool quotient_q = p_impl->quotient_q;
+
+    for (auto p : p_impl->points) {
+        if (char_fun(p) != 0) {
+            corresp.push_back(n_points++);
+            points.push_back(p);
+        } else {
+            if (not quotient_q) {
+                // need to realize that it is not geometric, (flip
+                // quotient_q and add the virtual point
+                points.insert(points.begin(), point_t());
+                quotient_q = true;
+            }
+            corresp.push_back(0);
+        }
+    }
+
+    if (not quotient_q)
+        for (int i = 0; i < corresp.size(); i++) corresp.at(i) -= 1;
+
+    // producing list of simplices
+    std::vector<cell_t> simp_list;
+
+    for (auto s : p_impl->simplices.complex_simplex_range()) {
+        cell_t s_q;
+        for (int v : p_impl->simplices.simplex_vertex_range(s))
+            s_q.push_back(corresp.at(v));
+        simp_list.push_back(s_q);
+    }
+
+    simplicial_complex quotient_sc(points, simp_list);
+    quotient_sc.p_impl->quotient_q = quotient_q;
+
+    return quotient_sc;
+}
+
+};  // namespace gsimp
+
+// DEPRECATED
+
+/*
+ *std::set<size_t> vertex_set(simp_handle s) {
+ *    std::set<size_t> v_set;
+ *    for (auto v : simplices.simplex_vertex_range(s)) {
+ *        v_set.insert(v);
+ *    }
+ *    return v_set;
+ *}
+ */
+
+/* MAY BE INCLUDED AGAIN WHEN WE START DOING MORE WITH QUOTIENTS
+ *
+ * bool is_point(point_t pt) {
+ *     return (not quotient_q ? true : not(pt.size() == 0));
+ * }
+ *
+ */
diff --git a/lib/scomplex/simplicial_complex.hpp b/lib/scomplex/simplicial_complex.hpp
new file mode 100644
index 0000000..a292cc5
--- /dev/null
+++ b/lib/scomplex/simplicial_complex.hpp
@@ -0,0 +1,52 @@
+#pragma once
+
+#include <scomplex/types.hpp>
+#include <memory>
+
+namespace gsimp {
+
+class NoChain {};
+
+class simplicial_complex {
+    // implementation details
+    struct impl;
+    std::shared_ptr<impl> p_impl;
+
+   public:
+    // constructor (no default)
+    simplicial_complex(std::vector<point_t>&, std::vector<cell_t>&);
+    simplicial_complex(const simplicial_complex&);
+    simplicial_complex& operator=(const simplicial_complex&);
+    // destructor
+    ~simplicial_complex();
+    // basic info
+    std::vector<point_t> get_points();
+    int dimension();
+    bool is_quotient();
+    // level-wise info
+    chain_t new_chain(int d);
+    // create chains
+    std::vector<cell_t> get_level(int);
+    int get_level_size(int);
+    // inclusion index
+    int boundary_inclusion_index(cell_t, cell_t);
+    int boundary_inclusion_index(int, size_t, int, size_t);
+    // cell boundaries
+    std::vector<cell_t> cell_boundary(cell_t);
+    std::vector<size_t> cell_boundary_index(int, size_t);
+    std::vector<std::pair<int, cell_t>> get_bdry_and_ind(cell_t);
+    std::vector<std::pair<int, size_t>> get_bdry_and_ind_index(int, size_t);
+    // treating cofaces
+    std::vector<cell_t> get_cofaces(cell_t);
+    std::vector<size_t> get_cofaces_index(int, size_t);
+    std::vector<std::pair<int, cell_t>> get_cof_and_ind(cell_t);
+    std::vector<std::pair<int, size_t>> get_cof_and_ind_index(int, size_t);
+    // boundary matrices
+    matrix_t get_boundary_matrix(int);
+    // cells and indices back and forth
+    cell_t index_to_cell(int, size_t);
+    size_t cell_to_index(cell_t);
+    //
+    simplicial_complex quotient(int char_fun(point_t));
+};  // class simplicial_complex
+};  // namespace gsimp
diff --git a/lib/scomplex/types.hpp b/lib/scomplex/types.hpp
new file mode 100644
index 0000000..3db2d56
--- /dev/null
+++ b/lib/scomplex/types.hpp
@@ -0,0 +1,63 @@
+#pragma once
+
+#include <iostream>
+#include <Eigen/Sparse>
+#include <functional>
+#include <vector>
+
+namespace gsimp {
+// geometric types
+typedef std::vector<double> point_t;
+typedef std::vector<size_t> cell_t;
+
+// linear algebra types
+typedef typename Eigen::SparseMatrix<double> matrix_t;
+typedef typename Eigen::SparseVector<double> vector_t;
+
+// chains
+typedef typename std::pair<int, vector_t> chain_t;
+
+int& chain_dim(chain_t& p) { return std::get<0>(p); }
+vector_t& chain_rep(chain_t& p) { return std::get<1>(p); }
+double& chain_val(chain_t& p, size_t i) { return std::get<1>(p).coeffRef(i); }
+
+chain_t add(chain_t chain1, chain_t chain2) {
+    if (std::get<0>(chain1) == std::get<0>(chain2))
+        return chain_t(std::get<0>(chain1),
+                       std::get<1>(chain1) + std::get<1>(chain2));
+    std::cout << "chains must have the same dimension to be added (+)";
+    throw std::exception();
+}
+
+void add_to(chain_t chain1, chain_t chain2) {
+    if (std::get<0>(chain1) == std::get<0>(chain2))
+        std::get<1>(chain1) += std::get<1>(chain2);
+    else {
+        std::cout << "chains must have the same dimension to be added (+=)";
+        throw std::exception();
+    }
+}
+
+chain_t subtract(chain_t chain1, chain_t chain2) {
+    if (std::get<0>(chain1) == std::get<0>(chain2))
+        return chain_t(std::get<0>(chain1),
+                       std::get<1>(chain1) - std::get<1>(chain2));
+    std::cout << "chains must have the same dimension to be added (-)";
+    throw std::exception();
+}
+
+void subtract_to(chain_t chain1, chain_t chain2) {
+    if (std::get<0>(chain1) == std::get<0>(chain2))
+        std::get<1>(chain1) -= std::get<1>(chain2);
+    else {
+        std::cout << "chains must have the same dimension to be added (-=)";
+        throw std::exception();
+    }
+}
+
+chain_t prod(double coef, chain_t chain) {
+    return chain_t(std::get<0>(chain), std::get<1>(chain) * coef);
+}
+
+void prod_to(double coef, chain_t chain) { std::get<1>(chain) *= coef; }
+}  // namespace gsimp
diff --git a/lib/scomplex/utils.h b/lib/scomplex/utils.h
deleted file mode 100644
index d96fa9d..0000000
--- a/lib/scomplex/utils.h
+++ /dev/null
@@ -1,4 +0,0 @@
-#pragma once
-
-#include "scomplex/Simplicial_Complex.h"
-#include "scomplex/graph_utils.hpp"
diff --git a/lib/scomplex/utils.hpp b/lib/scomplex/utils.hpp
new file mode 100644
index 0000000..0b3d61b
--- /dev/null
+++ b/lib/scomplex/utils.hpp
@@ -0,0 +1,98 @@
+#include <scomplex/types.hpp>
+#include <scomplex/drawing_utils.hpp>
+#include <scomplex/simplicial_complex.hpp>
+
+namespace gsimp {
+
+class DimensionException : public std::exception {};
+
+void draw_complex(simplicial_complex& scomp, drawing& my_drawing,
+                  double line_thickness, std::string fill_color = "#00ff00",
+                  std::string line_color = "#000000") {
+    if (scomp.dimension() > 2) DimensionException();
+
+    auto points = scomp.get_points();
+    auto edges = scomp.get_level(1);
+    auto tris = scomp.get_level(2);
+
+    for (auto cell : tris) {
+        std::vector<std::vector<double>> cell_pt;
+        for (auto i : cell) cell_pt.push_back(points[i]);
+        my_drawing.draw_triangle(cell_pt, fill_color);
+    }
+
+    for (auto cell : scomp.get_level(1)) {
+        std::vector<std::vector<double>> cell_pt;
+        for (auto i : cell) cell_pt.push_back(points[i]);
+        my_drawing.draw_edge(cell_pt[0], cell_pt[1], line_color,
+                             line_thickness);
+    }
+}
+
+void draw_path(path_snapper& p_snap, drawing& my_drawing,
+               std::vector<std::vector<double>> path, double stroke,
+               std::string path_color = "#ff0000") {
+    auto snapped = p_snap.snap_path_to_points(path);
+    auto it = snapped.begin();
+    auto trail = it++;
+    for (; it != snapped.end(); ++it, ++trail) {
+        auto init = *trail;
+        auto dest = *it;
+        my_drawing.draw_edge(init, dest, path_color, stroke);
+    }
+}
+
+void draw_2_chain(simplicial_complex& s_comp, drawing& my_drawing,
+                  chain_t chain, std::string chain_color = "#ff0000") {
+    if (chain_dim(chain) != 2) DimensionException();
+
+    auto points = s_comp.get_points();
+    auto level_2 = s_comp.get_level(2);
+    Eigen::SparseVector<double>::InnerIterator it(chain_rep(chain));
+    for (; it; ++it) {
+        cell_t cell = level_2[it.index()];
+        std::vector<point_t> cell_pt;
+        for (size_t i : cell) cell_pt.push_back(points[i]);
+        // TODO: get rid of these thresholds
+        if (it.value() > 0.01 || it.value() < -0.01)
+        my_drawing.draw_triangle(cell_pt, chain_color);
+    }
+}
+void draw_1_chain(simplicial_complex& s_comp, drawing& my_drawing,
+                  chain_t chain, double stroke,
+                  std::string chain_color = "#ff0000") {
+    if (chain_dim(chain) != 1) DimensionException();
+
+    auto points = s_comp.get_points();
+    auto level_1 = s_comp.get_level(1);
+    Eigen::SparseVector<double>::InnerIterator it(chain_rep(chain));
+    for (; it; ++it) {
+        cell_t edge = level_1[it.index()];
+        my_drawing.draw_edge(points[edge[0]], points[edge[1]], chain_color,
+                             stroke);
+    }
+}
+
+void draw_path(simplicial_complex& scomp, drawing& my_drawing,
+               std::vector<std::vector<double>> path, double stroke,
+               std::string path_color = "#ff0000") {
+    path_snapper p_snap(scomp);
+    auto snapped = p_snap.snap_path_to_points(path);
+    auto it = snapped.begin();
+    auto trail = it++;
+    for (; it != snapped.end(); ++it, ++trail) {
+        auto init = *trail;
+        auto dest = *it;
+        my_drawing.draw_edge(init, dest, path_color, stroke);
+    }
+}
+
+drawing draw_complex(simplicial_complex& scomp, std::string filename,
+                     double stroke, int width = 500, int height = 500,
+                     std::string fill_color = "#00ff00",
+                     std::string line_color = "#000000") {
+    drawing my_drawing(filename, width, height);
+    draw_complex(scomp, my_drawing, stroke, fill_color, line_color);
+    return my_drawing;
+}
+};
diff --git a/src/graphtest.cpp b/src/graphtest.cpp
deleted file mode 100644
index 6329c70..0000000
--- a/src/graphtest.cpp
+++ /dev/null
@@ -1,15 +0,0 @@
-
-#include <iostream>
-
-#include "boost/graph/adjacency_list.hpp"
-#include "boost/graph/dijkstra_shortest_paths.hpp"
-
-typedef boost::adjacency_list<boost::vecS, boost::vecS, boost::undirectedS>
-    Graph;
-
-int main() {
-    Graph g(10);
-    boost::add_edge(0, 1, g);
-
-    return 0;
-}
diff --git a/src/libtest.cpp b/src/libtest.cpp
deleted file mode 100644
index 9a7a696..0000000
--- a/src/libtest.cpp
+++ /dev/null
@@ -1,160 +0,0 @@
-#include <scomplex/Simplicial_Complex.h>
-#include <scomplex/utils.h>
-#include <scomplex/nn_utils.hpp>
-#include <scomplex/trace.h>
-#include <scomplex/qhull_parsing.hpp>
-
-#include "boost/graph/graph_traits.hpp"
-
-#include <iostream>
-#include <vector>
-#include <string>
-#include <tuple>
-
-typedef std::vector<double> point_t;
-typedef std::list<int> simp_t;
-typedef std::pair<int, std::vector<int>> chain_t;
-
-int f(point_t a) { return a.at(0) < .5 ? 0 : 1; }
-
-int main() {
-    std::list<point_t> point_list;
-    std::list<simp_t> simp_list;
-    // ---
-
-    point_list.push_back(point_t({0, 0}));
-    point_list.push_back(point_t({0, 5}));
-    point_list.push_back(point_t({1, 0}));
-    point_list.push_back(point_t({1, 1}));
-
-    simp_list.push_back(simp_t({0, 1, 2}));
-    simp_list.push_back(simp_t({1, 2, 3}));
-
-    simplicial::SimplicialComplex<point_t> sc(point_list, simp_list);
-    for (auto pt : sc.points) {
-        std::cout << pt.at(0) << " " << pt.at(1) << std::endl;
-    }
-
-    auto range = sc.simplices.complex_simplex_range();
-    for (auto s : range) {
-        std::cout << sc.simplices.dimension(s) << " , " << sc.simplices.key(s)
-                  << " : ";
-        for (auto v : sc.simplices.simplex_vertex_range(s)) {
-            std::cout << v << " ";
-        }
-        // iterate over boundary
-        std::cout << std::endl
-                  << "    boundary: ";
-        for (auto sb : sc.simplices.boundary_simplex_range(s)) {
-            auto v_range = sc.simplices.simplex_vertex_range(sb);
-            for (auto v : v_range) {
-                std::cout << v << " ";
-            }
-            std::cout << " ; ";
-        }
-        std::cout << std::endl
-                  << "  coboundary: ";
-        for (auto sb : sc.simplices.cofaces_simplex_range(s, 1)) {
-            auto v_range = sc.simplices.simplex_vertex_range(sb);
-            for (auto v : v_range) {
-                std::cout << v << " ";
-            }
-            std::cout << " ; ";
-        }
-        std::cout << std::endl;
-    }
-
-    std::cout << "now quotient" << std::endl;
-    auto sq = sc.quotient(f);
-
-    /* the points are working, what's up with the simplices??
-    for (auto p : sq.points) {
-        std::cout << "point: " << p.at(0) << " " << p.at(1) << std::endl;
-    }
-    */
-
-    range = sq.simplices.complex_simplex_range();
-    std::cout << "got this far" << std::endl;
-    for (auto s : range) {
-        std::cout << sq.simplices.dimension(s) << " , " << sq.simplices.key(s)
-                  << " : ";
-        for (auto v : sq.simplices.simplex_vertex_range(s)) {
-            std::cout << v << " ";
-        }
-        // iterate over boundary
-        std::cout << std::endl
-                  << "    boundary: ";
-        for (auto sb : sq.simplices.boundary_simplex_range(s)) {
-            auto v_range = sq.simplices.simplex_vertex_range(sb);
-            for (auto v : v_range) {
-                std::cout << v << " ";
-            }
-            std::cout << " ; ";
-        }
-        std::cout << std::endl
-                  << "  coboundary: ";
-        for (auto sb : sq.simplices.cofaces_simplex_range(s, 1)) {
-            auto v_range = sq.simplices.simplex_vertex_range(sb);
-            for (auto v : v_range) {
-                std::cout << v << " ";
-            }
-            std::cout << " ; ";
-        }
-        std::cout << std::endl;
-    }
-    std::cout << "got this far" << std::endl;
-
-    // checking the matrices
-    auto g = sc.get_boundary(1);
-    std::cout << g << std::endl;
-    g = sc.get_boundary(0);
-    std::cout << g << std::endl;
-
-    g = sq.get_boundary(1);
-    std::cout << g << std::endl;
-    g = sq.get_boundary(0);
-    std::cout << g << std::endl;
-
-    std::cout << "printing level: " << std::endl;
-    for (auto s : sq.get_level(1)) {
-        std::cout << s.size() << ": ";
-        std::cout << s.at(0) << " ";
-        std::cout << s.at(1) << " ";
-        std::cout << std::endl;
-    }
-
-    auto G = calculate_one_skelleton_graph(sc);
-    auto p = shortest_path(G, 0, 3);
-
-    typedef typename decltype(G)::vertex_descriptor vertex_t;
-
-    /*
-    std::cout << '\n';
-    std::copy(p.begin(), p.end(),
-              std::ostream_iterator<vertex_t>{std::cout, " "});
-    std::cout << '\n';
-    */
-
-    tree_t rt;
-    make_tree(rt, point_list);
-
-    point_t q{0.1, 0.2};
-    point_t q_prime = nearest_neighbour(rt, q);
-    // get_r_tree_from_points<2>(point_list);
-    std::ostream_iterator<double> outstr(std::cout, " ");
-    std::copy(q_prime.begin(), q_prime.end(), outstr);
-
-    std::vector<point_t> points_v;
-    std::vector<cell_t> cells_v;
-    std::tie(points_v, cells_v) =
-        parse_qhull_file("/home/crvs/qhull-test/qh-test.dat");
-    std::list<point_t> points(points_v.begin(), points_v.end());
-    std::list<std::list<int>> cells;
-    for (auto cell_v : cells_v) {
-        std::list<int> cell(cell_v.begin(), cell_v.end());
-        cells.push_back(cell);
-    }
-
-    simplicial::SimplicialComplex<point_t> rsc(points, cells);
-    return 0;
-}
diff --git a/src/main.cpp b/src/main.cpp
deleted file mode 100644
index c296cbb..0000000
--- a/src/main.cpp
+++ /dev/null
@@ -1,68 +0,0 @@
-#include <gudhi/Simplex_tree.h>
-
-#include <iostream>
-#include <vector>
-
-int main() {
-    struct Myoptions : Gudhi::Simplex_tree_options_full_featured {
-        typedef int Vertex_handle;
-    };
-    Gudhi::Simplex_tree<Myoptions> simplexTree;
-    auto d = {0, 1, 2};
-    simplexTree.insert_simplex_and_subfaces(d);
-
-    auto range = simplexTree.complex_simplex_range();
-    int count = 0;
-
-    for (auto s : range) {
-        simplexTree.assign_key(s, count++);
-        auto v_range = simplexTree.simplex_vertex_range(s);
-        std::cout << simplexTree.filtration(s) << " , " << simplexTree.key(s)
-                  << " : ";
-        for (auto v : v_range) {
-            std::cout << v << " ";
-        }
-        std::cout << ": ";
-        auto sbd = simplexTree.boundary_simplex_range(s);
-        for (auto sb : sbd) {
-            auto v_range = simplexTree.simplex_vertex_range(sb);
-            for (auto v : v_range) {
-                std::cout << v << " ";
-            }
-            std::cout << " ; ";
-        }
-        std::cout << std::endl;
-    }
-
-    // testing the effects of reinsertion
-    simplexTree.insert_simplex_and_subfaces(d);
-    for (auto s : range) {
-        auto v_range = simplexTree.simplex_vertex_range(s);
-        std::cout << simplexTree.filtration(s) << " , " << simplexTree.key(s)
-                  << " : ";
-        for (auto v : v_range) {
-            std::cout << v << " ";
-        }
-        std::cout << ": ";
-        auto sbd = simplexTree.boundary_simplex_range(s);
-        for (auto sb : sbd) {
-            auto v_range = simplexTree.simplex_vertex_range(sb);
-            for (auto v : v_range) {
-                std::cout << v << " ";
-            }
-            std::cout << " ; ";
-        }
-        std::cout << std::endl;
-    }
-
-    auto p = simplexTree.find({0, 1, 2});
-    for (auto v : simplexTree.simplex_vertex_range(p)) {
-        std::cout << v << " ";
-    }
-    std::cout << std::endl;
-
-    std::vector<int> vec;
-    std::cout << "vector size " << vec.size() << std::endl;
-
-    return 0;
-}
diff --git a/src/test0.cpp b/src/test0.cpp
new file mode 100644
index 0000000..3dbf707
--- /dev/null
+++ b/src/test0.cpp
@@ -0,0 +1,59 @@
+/**
+ * SIMPLICIAL COMPLEX TEST:
+ *     takes as argument a path to a file containing the description of a
+ *   triangulation as output by qhull (the last column of the point coordinates
+ *   is ignored).
+ *     constantly asks the user to input a number which will be interpreted as a
+ *   1-cell index, and outputs the cofaces of the given cell, (together with the
+ *   boundary index) and the faces of the cell.
+ **/
+#include <vector>
+#include <iostream>
+#include <iterator>
+
+#include <scomplex/types.hpp>
+#include <scomplex/qhull_parsing.hpp>
+#include <scomplex/simplicial_complex.hpp>
+
+int my_char(gsimp::point_t point) { return 1; }
+
+int& do_this(int& a) { return a; }
+
+int main(int argc, char* argv[]) {
+    std::vector<gsimp::point_t> points_v;
+    std::vector<gsimp::cell_t> cells_v;
+
+    std::string filename;
+    if (argc == 1) {
+        std::cout << "insert file name: ";
+        std::cin >> filename;
+    } else
+        filename = argv[1];
+    std::tie(points_v, cells_v) = parse_qhull_file(filename);
+
+    gsimp::simplicial_complex s_comp(points_v, cells_v);
+
+    std::cout << "complex dimension: " << s_comp.dimension() << "\n";
+    for (int i = 0; i <= s_comp.dimension(); ++i) {
+        std::cout << "level " << i << " has " << s_comp.get_level_size(i)
+                  << " cells\n";
+    }
+
+    size_t i;
+    std::cout << "get 1-cell number: ";
+    while (std::cin >> i) {
+        std::vector<std::pair<int, size_t>> my_cofaces =
+            s_comp.get_cof_and_ind_index(1, i);
+        for (auto s : my_cofaces) {
+            std::cout << std::get<0>(s) << ", ";
+            std::cout << std::get<1>(s) << " ;  ";
+        }
+        std::cout << '\n';
+        std::vector<size_t> faces = s_comp.cell_boundary_index(1, i);
+        std::copy(faces.begin(), faces.end(),
+                  std::ostream_iterator<size_t>(std::cout, " "));
+        std::cout << "\nget 1-cell number: ";
+    }
+
+    return 0;
+}
diff --git a/src/test1.cpp b/src/test1.cpp
new file mode 100644
index 0000000..960ec85
--- /dev/null
+++ b/src/test1.cpp
@@ -0,0 +1,55 @@
+/**
+ * PATH SNAPPING TEST:
+ *     takes as argument a path to a file containing the description of a
+ *   triangulation as output by qhull (the last column of the point coordinates
+ *   is ignored).
+ *     comes up with a path in this complex (hard coded) and snaps it to the
+ *   complex, the only point of this is to test the libraries for errors.
+ **/
+
+#include <iostream>
+#include <iterator>
+
+#include <scomplex/types.hpp>
+#include <scomplex/simplicial_complex.hpp>
+#include <scomplex/qhull_parsing.hpp>
+
+// testing
+#include <scomplex/path_snapper.hpp>
+
+int my_char(gsimp::point_t point) { return 1; }
+
+int& do_this(int& a) { return a; }
+
+int main(int argc, char* argv[]) {
+    std::vector<point_t> points_v;
+    std::vector<cell_t> cells_v;
+
+    std::string filename;
+    if (argc == 1) {
+        std::cout << "insert file name: ";
+        std::cin >> filename;
+    } else
+        filename = argv[1];
+    std::tie(points_v, cells_v) = parse_qhull_file(filename);
+
+    auto s_comp = std::shared_ptr<gsimp::simplicial_complex>(
+        new gsimp::simplicial_complex(points_v, cells_v));
+
+    std::cout << "complex dimension: " << s_comp->dimension() << "\n";
+    for (int i = 0; i <= s_comp->dimension(); ++i) {
+        std::cout << "level " << i << " has " << s_comp->get_level_size(i)
+                  << " cells\n";
+    }
+
+    gsimp::path_snapper p_snap(s_comp);
+    std::vector<size_t> snapped =
+        p_snap.snap_path_to_indices({{-2, -2}, {-2, 2}, {2, 2}, {-2, -2}});
+    std::copy(snapped.begin(), snapped.end(),
+              std::ostream_iterator<size_t>(std::cout, " "));
+    std::cout << '\n';
+
+    auto chain = p_snap.index_sequence_to_chain(snapped);
+
+    return 0;
+}
diff --git a/src/test2.cpp b/src/test2.cpp
new file mode 100644
index 0000000..e27c8df
--- /dev/null
+++ b/src/test2.cpp
@@ -0,0 +1,60 @@
+/**
+ * PATH SNAPPING TEST:
+ *     takes as argument a path to a file containing the description of a
+ *   triangulation as output by qhull (the last column of the point coordinates
+ *   is ignored).
+ *     outputs the complex to an svg file called my_complex
+ **/
+
+#include <iostream>
+#include <iterator>
+
+#include <scomplex/types.hpp>
+#include <scomplex/qhull_parsing.hpp>
+#include <scomplex/simplicial_complex.hpp>
+#include <scomplex/path_snapper.hpp>
+#include <scomplex/drawing_utils.hpp>
+
+int my_char(gsimp::point_t point) { return 1; }
+
+int main(int argc, char* argv[]) {
+    std::vector<point_t> points_v;
+    std::vector<cell_t> cells_v;
+
+    std::string filename;
+    if (argc == 1) {
+        std::cout << "insert file name: ";
+        std::cin >> filename;
+    } else
+        filename = argv[1];
+    std::tie(points_v, cells_v) = parse_qhull_file(filename);
+
+    gsimp::simplicial_complex s_comp(points_v, cells_v);
+
+    gsimp::drawing drawing("my_complex", 500, 500);
+
+    for (auto cell : s_comp.get_level(2)) {
+        std::vector<std::vector<double>> cell_pt;
+        for (auto i : cell) cell_pt.push_back(points_v[i]);
+        drawing.draw_triangle(cell_pt, "#00ff00");
+    }
+
+    for (auto cell : s_comp.get_level(1)) {
+        std::vector<std::vector<double>> cell_pt;
+        for (auto i : cell) cell_pt.push_back(points_v[i]);
+        drawing.draw_edge(cell_pt[0], cell_pt[1], "#000000", 0.002);
+    }
+
+    gsimp::path_snapper p_snap(s_comp);
+    std::vector<size_t> snapped =
+        p_snap.snap_path_to_indices({{-2, -2}, {-2, 2}, {2, 2}, {-2, -2}});
+    auto it = snapped.begin();
+    auto trail = it++;
+    for (; it != snapped.end(); ++it, ++trail) {
+        auto init = s_comp.get_points()[*trail];
+        auto dest = s_comp.get_points()[*it];
+        drawing.draw_edge(init, dest, "#ff0000", 0.003);
+    }
+
+    return 0;
+}
diff --git a/src/test3.cpp b/src/test3.cpp
new file mode 100644
index 0000000..e92babc
--- /dev/null
+++ b/src/test3.cpp
@@ -0,0 +1,48 @@
+/**
+ * PATH SNAPPING TEST:
+ *     takes as argument a path to a file containing the description of a
+ *   triangulation as output by qhull (the last column of the point coordinates
+ *   is ignored).
+ *     outputs the complex to an svg file called my_complex
+ **/
+
+#include <iostream>
+#include <iterator>
+
+#include <scomplex/types.hpp>
+#include <scomplex/qhull_parsing.hpp>
+#include <scomplex/simplicial_complex.hpp>
+#include <scomplex/path_snapper.hpp>
+#include <scomplex/utils.hpp>
+
+// testing
+#include <scomplex/drawing_utils.hpp>
+
+int my_char(gsimp::point_t point) { return 1; }
+
+int main(int argc, char* argv[]) {
+    std::vector<point_t> points_v;
+    std::vector<cell_t> cells_v;
+
+    std::string filename;
+    if (argc == 1) {
+        std::cout << "insert file name: ";
+        std::cin >> filename;
+    } else
+        filename = argv[1];
+    std::tie(points_v, cells_v) = parse_qhull_file(filename);
+
+    gsimp::simplicial_complex s_comp(points_v, cells_v);
+
+    gsimp::drawing drawing("my_complex", 500, 500);
+
+    gsimp::draw_complex(s_comp, drawing, 0.003, "#00ff00", "#000000");
+
+    { // scope for p_snap
+        gsimp::path_snapper p_snap(s_comp);
+        gsimp::draw_path(p_snap, drawing, {{-2, -2}, {2, -2}, {2, 2}, {-2, -2}},
+                         0.003);
+    }
+
+    return 0;
+}
diff --git a/src/test4.cpp b/src/test4.cpp
new file mode 100644
index 0000000..6faebde
--- /dev/null
+++ b/src/test4.cpp
@@ -0,0 +1,61 @@
+/**
+ * PATH SNAPPING TEST:
+ *     takes as argument a path to a file containing the description of a
+ *   triangulation as output by qhull (the last column of the point coordinates
+ *   is ignored).
+ *     outputs the complex to an svg file called my_complex
+ **/
+
+#include <iostream>
+#include <iterator>
+
+#include <scomplex/types.hpp>
+#include <scomplex/qhull_parsing.hpp>
+#include <scomplex/simplicial_complex.hpp>
+#include <scomplex/path_snapper.hpp>
+#include <scomplex/drawing_utils.hpp>
+#include <scomplex/utils.hpp>
+
+// testing chain_calc
+#include <scomplex/chain_calc.hpp>
+
+int my_char(gsimp::point_t point) { return 1; }
+
+int main(int argc, char* argv[]) {
+    std::vector<point_t> points_v;
+    std::vector<cell_t> cells_v;
+
+    std::string filename;
+    if (argc == 1) {
+        std::cout << "insert file name: ";
+        std::cin >> filename;
+    } else
+        filename = argv[1];
+    std::tie(points_v, cells_v) = parse_qhull_file(filename);
+
+    gsimp::simplicial_complex s_comp(points_v, cells_v);
+    gsimp::drawing drawing("my_complex", 500, 500);
+    gsimp::draw_complex(s_comp, drawing, 0.003, "#00ff00", "#000000");
+
+    std::vector<point_t> path2 = {{-2, -2}, { 2, -2}, { 2,  2}, {-2, -2}};
+    std::vector<point_t> path  = {{-2, -2}, {.3, .3}, {.2, -2}, {-2, -2}};
+
+    // auxiliary structures
+    gsimp::path_snapper p_snap(s_comp);
+    gsimp::bounding_chain c_calc(s_comp);
+
+    gsimp::draw_path(p_snap, drawing, path, 0.003);
+    gsimp::draw_path(p_snap, drawing, path2, 0.003);
+
+    gsimp::chain_t path_chain = p_snap.snap_path_to_chain(path);
+    gsimp::add_to(p_snap.snap_path_to_chain(path2),path_chain);
+    try {
+        gsimp::chain_t bound_chain = c_calc.get_bounding_chain(path_chain);
+        draw_2_chain(s_comp, drawing, bound_chain,
+                     "#ff00ff' fill-opacity='0.8");
+    } catch (const std::exception& e) {
+        std::cout << "oh well, nothing I can do";
+    }
+
+    return 0;
+}
diff --git a/src/test5.cpp b/src/test5.cpp
new file mode 100644
index 0000000..3d1e72c
--- /dev/null
+++ b/src/test5.cpp
@@ -0,0 +1,58 @@
+/**
+ * PATH SNAPPING TEST:
+ *     takes as argument a path to a file containing the description of a
+ *   triangulation as output by qhull (the last column of the point coordinates
+ *   is ignored).
+ *     outputs the complex to an svg file called my_complex
+ **/
+
+#include <iostream>
+#include <iterator>
+
+#include <scomplex/coeff_flow.hpp>
+#include <scomplex/types.hpp>
+#include <scomplex/qhull_parsing.hpp>
+#include <scomplex/simplicial_complex.hpp>
+#include <scomplex/path_snapper.hpp>
+#include <scomplex/drawing_utils.hpp>
+#include <scomplex/utils.hpp>
+
+// testing chain_calc
+#include <scomplex/chain_calc.hpp>
+
+int my_char(gsimp::point_t point) { return 1; }
+
+int main(int argc, char* argv[]) {
+    std::vector<point_t> points_v;
+    std::vector<cell_t> cells_v;
+
+    std::string filename;
+    if (argc == 1) {
+        std::cout << "insert file name: ";
+        std::cin >> filename;
+    } else
+        filename = argv[1];
+    std::tie(points_v, cells_v) = parse_qhull_file(filename);
+
+    gsimp::simplicial_complex s_comp(points_v, cells_v);
+    gsimp::drawing drawing("my_complex", 500, 500);
+    gsimp::draw_complex(s_comp, drawing, 0.003, "#00ff00", "#000000");
+
+    std::vector<point_t> path2 = {{-2, -2}, { 2, -2}, { 2,  2}, {-2, -2}};
+    std::vector<point_t> path  = {{-2, -2}, {.3, .3}, {.2, -2}, {-2, -2}};
+
+    // auxiliary structures
+    gsimp::path_snapper p_snap(s_comp);
+    gsimp::bounding_chain c_calc(s_comp);
+
+    gsimp::draw_path(p_snap, drawing, path, 0.003);
+    gsimp::draw_path(p_snap, drawing, path2, 0.003);
+
+    gsimp::chain_t path_chain = p_snap.snap_path_to_chain(path);
+    gsimp::add_to(p_snap.snap_path_to_chain(path2),path_chain);
+
+    gsimp::chain_t bound_chain = gsimp::coeff_flow_embedded(s_comp,path_chain);
+    gsimp::draw_2_chain(s_comp, drawing, bound_chain, "#ff00ff' fill-opacity='0.8");
+
+    return 0;
+}
diff --git a/src/utils.hpp b/src/utils.hpp
new file mode 100644
index 0000000..7348ac5
--- /dev/null
+++ b/src/utils.hpp
@@ -0,0 +1,2 @@
+#include <scomplex/simplicial_complex.hpp>
+namespace gsimp