diff --git a/bb_test.jl b/bb_test.jl
new file mode 100644
index 00000000..11191a0d
--- /dev/null
+++ b/bb_test.jl
@@ -0,0 +1,69 @@
+using Revise
+using Test
+using Polymake
+
+function BB_iterator(c)
+    N = polytope.dim(c)
+    @info "running bb using iterator over $N-dimensional polytope"
+    bb = Polymake.BeneathBeyond(c.POINTS, c.LINEALITY_SPACE, redundant=true, triangulation=true, iscone=true)
+    tr = GC.@preserve bb begin
+        for x in bb
+            @show length(Polymake.getTriangulation(bb.algo))
+        end
+        tr = Polymake.getTriangulation(bb.algo)
+        @show last(tr)
+        tr
+    end
+    @info "bb is done!"
+    return tr
+end
+
+
+function BB_direct_call(c)
+    N = polytope.dim(c)
+    @info "running bb using direct calls over $N-dimensional polytope"
+    bb = Polymake.BeneathBeyondAlgo{Polymake.Rational}()
+    P, L = c.POINTS, c.LINEALITY_SPACE
+    tr = GC.@preserve bb P L begin    
+        Polymake.bb_initialize!(bb, P, L)
+        
+        for p in 0:6
+            @time Polymake.bb_add_point!(bb, p)
+        end
+        
+        @show Polymake.getTriangulation(bb)
+        
+        @time for p in 2:size(c.POINTS, 1)-1
+            Polymake.bb_add_point!(bb, p)
+        end
+        tr = Polymake.getTriangulation(bb)
+        @show last(tr)
+        @info "bb is done!"
+        tr
+    end
+    return tr
+end;
+
+
+
+const rs = polytope.rand_sphere(4,30);
+
+@info "running placing_triangulation"
+t1 = polytope.placing_triangulation(rs.POINTS)
+t2 = let rs = rs
+    N = polytope.dim(rs)
+    @info "running bb using bb_compute! over $N-dimensional polytope"
+    bbalgo = Polymake.BeneathBeyondAlgo{Polymake.Rational}()
+    P, L = rs.POINTS, rs.LINEALITY_SPACE
+    GC.@preserve bbalgo P L begin
+        Polymake.bb_compute!(bbalgo, P, L)
+        Polymake.getTriangulation(bbalgo)
+    end
+end
+
+@test t1 == t2
+
+t3 = BB_direct_call(rs)
+@test t1 == t3
+t4 = BB_iterator(rs)
+@test t1 == t4
diff --git a/deps/src/beneath_beyond.cpp b/deps/src/beneath_beyond.cpp
new file mode 100644
index 00000000..d25a638f
--- /dev/null
+++ b/deps/src/beneath_beyond.cpp
@@ -0,0 +1,264 @@
+#include "polymake_includes.h"
+#include <polymake/polytope/beneath_beyond_impl.h>
+#include "polymake_type_modules.h"
+
+namespace polymake { namespace polytope {
+
+template <typename E>
+class beneath_beyond_algo_for_ml: public beneath_beyond_algo<E>{
+    public:
+        typedef E value_type;
+        typedef const beneath_beyond_algo<E> Base;
+
+        beneath_beyond_algo_for_ml(): Base()
+        {
+            initialized = false;
+        };
+
+        template <typename Iterator>
+        void initialize(const Matrix<E>& rays, const Matrix<E>& lins, Iterator perm);
+        void initialize(const Matrix<E>& rays, const Matrix<E>& lins)
+        {
+        #if POLYMAKE_DEBUG
+            enable_debug_output();
+        #endif
+            initialize(rays, lins, entire(sequence(0, rays.rows())));
+        };
+
+        void process_point(Int p);
+
+        void clear();
+
+
+        // TODO: bundle all results in a structure, move all numbers into it
+        template <typename Iterator>
+        void compute(const Matrix<E>& rays, const Matrix<E>& lins, Iterator perm);
+        void compute(const Matrix<E>& rays, const Matrix<E>& lins)
+        {
+        #if POLYMAKE_DEBUG
+            enable_debug_output();
+        #endif
+            compute(rays, lins, entire(sequence(0, rays.rows())));
+        };
+
+    protected:
+        void stop_cleanup();
+
+        using Base::source_points;
+        using Base::source_linealities;
+        using Base::linealities_so_far;
+        using Base::expect_redundant;
+        using Base::source_lineality_basis;
+        using Base::linealities;
+        using Base::transform_points;
+        using Base::points;
+        using Base::generic_position;
+        using Base::triang_size;
+        using Base::AH;
+        using Base::interior_points;
+        using Base::vertices_this_step;
+        using Base::interior_points_this_step;
+        using Base::facet_normals_valid;
+        using Base::facet_normals_low_dim;
+        using Base::dual_graph;
+        using Base::vertices_so_far;
+        using Base::make_triangulation;
+        using Base::triangulation;
+        using Base::is_cone;
+        using Base::facets;
+        using compute_state = typename Base::compute_state;
+        using Base::state;
+        
+        class stop_calculation {};
+
+    private:
+        bool initialized;
+        Bitset points_added;
+};
+
+
+template <typename E>
+template <typename Iterator>
+void beneath_beyond_algo_for_ml<E>::initialize(const Matrix<E>& rays, const Matrix<E>& lins, Iterator perm)
+{
+    source_points = &rays;
+    source_linealities = &lins;
+
+    linealities_so_far.resize(0,rays.cols());
+
+    try {
+        if (lins.rows() != 0) {
+            if (expect_redundant) {
+                source_lineality_basis = basis_rows(lins);
+                linealities_so_far = lins.minor(source_lineality_basis, All);
+                linealities = &linealities_so_far;
+            } else {
+                linealities = source_linealities;
+            }
+            transform_points(); // the only place where stop_calculation could be thrown
+        } else {
+            points = source_points;
+            linealities = expect_redundant ? &linealities_so_far : source_linealities;
+        }
+
+        generic_position = !expect_redundant;
+        triang_size = 0;
+        AH = unit_matrix<E>(points->cols());
+        if (expect_redundant) {
+            interior_points.resize(points->rows());
+            vertices_this_step.resize(points->rows());
+            interior_points_this_step.resize(points->rows());
+        }
+
+        state = compute_state::zero; // moved from the main compute loop
+
+        points_added = Bitset();
+        initialized = true;
+    }
+    catch (const stop_calculation&) { 
+#if POLYMAKE_DEBUG
+        if (debug >= do_dump) cout << "stop: failed to initialize beneath_beyond_algo" << endl;
+#endif
+        // TODO: some cleanup??
+    }
+};
+
+template <typename E>
+void beneath_beyond_algo_for_ml<E>::process_point(Int p){
+    if ( !points_added.contains(p) ){
+        Base::process_point(p);
+        points_added += p;
+#if POLYMAKE_DEBUG
+        std::cout << "processed point p = " << p << std::endl;
+#endif
+    };
+};
+
+template <typename E>
+template <typename Iterator>
+void beneath_beyond_algo_for_ml<E>::compute(const Matrix<E>& rays, const Matrix<E>& lins, Iterator perm){
+    
+    initialize(rays, lins);
+
+    try
+    {
+        for (; !perm.at_end(); ++perm)
+            process_point(*perm);
+    }
+    catch (const stop_calculation&){
+#if POLYMAKE_DEBUG
+        if (debug >= do_dump) cout << "stop: degenerated to full linear space" << endl;
+#endif
+        stop_cleanup();
+    }
+
+    clear();
+
+#if POLYMAKE_DEBUG
+    if (debug >= do_dump) {
+        cout << "final ";
+        dump();
+    }
+#endif
+
+};
+
+template <typename E>
+void beneath_beyond_algo_for_ml<E>::stop_cleanup(){
+    state = compute_state::zero;
+    dual_graph.clear();
+    vertices_so_far.clear();
+    points = source_points;
+    interior_points = sequence(0, source_points->rows());
+    if (make_triangulation) {
+        triangulation.clear();
+        triang_size = 0;
+    }
+}
+
+template <typename E>
+void beneath_beyond_algo_for_ml<E>::clear(){
+
+    switch (state) {
+    case compute_state::zero:
+        if (!is_cone) {
+            // empty polyhedron
+            AH.resize(0, source_points->cols());
+            linealities_so_far.resize(0, source_points->cols());
+        }
+        break;
+    case compute_state::one:
+        // There is one empty facet in this case and the point is also a facet normal
+        facets[dual_graph.add_node()].normal = points->row(vertices_so_far.front());
+        if (make_triangulation) {
+            triang_size=1;
+            triangulation.push_back(vertices_so_far);
+        }
+        break;
+    case compute_state::low_dim:
+        if ( !facet_normals_valid )
+        {
+            try
+            {
+                facet_normals_low_dim();
+            }
+            catch(const stop_calculation& )
+            {
+                stop_cleanup();
+            }
+        }
+        break;
+    case compute_state::full_dim:
+        dual_graph.squeeze();
+        break;
+    }
+}
+
+}
+}
+
+
+template<> struct jlcxx::IsMirroredType<
+    polymake::polytope::beneath_beyond_algo<pm::Rational>> : std::false_type { };
+
+void polymake_module_add_beneath_beyond(jlcxx::Module& polymake)
+{
+    polymake
+        .add_type<jlcxx::Parametric<jlcxx::TypeVar<1>>>("_BeneathBeyondAlgo")
+        .apply<polymake::polytope::beneath_beyond_algo<pm::Rational>>([](auto wrapped) {});
+    
+    polymake
+        .add_type<jlcxx::Parametric<jlcxx::TypeVar<1>>>("BeneathBeyondAlgo")
+        .apply<polymake::polytope::beneath_beyond_algo_for_ml<pm::Rational>>([](auto wrapped) {
+            typedef typename decltype(wrapped)::type             WrappedT;
+            typedef typename decltype(wrapped)::type::value_type E;
+            wrapped.template constructor();
+
+            wrapped.method("bb_expecting_redundant", &WrappedT::expecting_redundant);
+            wrapped.method("bb_for_cone", &WrappedT::for_cone);
+            wrapped.method("bb_making_triangulation", &WrappedT::making_triangulation);
+            wrapped.method("bb_computing_vertices", &WrappedT::computing_vertices);
+
+            wrapped.method("bb_compute!", static_cast<
+                void (polymake::polytope::beneath_beyond_algo_for_ml<E>::*)(const pm::Matrix<E>&, const pm::Matrix<E>&)
+            >(&WrappedT::compute));
+
+            wrapped.method("bb_initialize!", static_cast<
+                void (polymake::polytope::beneath_beyond_algo_for_ml<E>::*)(const pm::Matrix<E>&, const pm::Matrix<E>&)
+            >(&WrappedT::initialize));
+
+            // wrapped.method("initialize", &WrappedT::initialize);
+            wrapped.method("bb_add_point!", &WrappedT::process_point);
+            wrapped.method("bb_clear!", &WrappedT::clear);
+
+            wrapped.method("getFacets", &WrappedT::getFacets);
+            wrapped.method("getVertexFacetIncidence", &WrappedT::getVertexFacetIncidence);
+            wrapped.method("getAffineHull", &WrappedT::getAffineHull);
+            wrapped.method("getVertices", &WrappedT::getVertices);
+            // wrapped.method("getNonRedundantPoints", &WrappedT::getNonRedundantPoints);
+            wrapped.method("getNonRedundantLinealities", &WrappedT::getNonRedundantLinealities);
+            wrapped.method("getLinealities", &WrappedT::getLinealities);
+            // wrapped.method("getDualGraph", &WrappedT::getDualGraph);
+            wrapped.method("getTriangulation", &WrappedT::getTriangulation);
+        });
+}
diff --git a/deps/src/polymake.cpp b/deps/src/polymake.cpp
index 06ad39a9..223c0bea 100644
--- a/deps/src/polymake.cpp
+++ b/deps/src/polymake.cpp
@@ -43,6 +43,7 @@ JLCXX_MODULE define_module_polymake(jlcxx::Module& polymake)
     polymake_module_add_polynomial(polymake);
 
     polymake_module_add_direct_calls(polymake);
+    polymake_module_add_beneath_beyond(polymake);
 
     polymake_module_add_array_polynomial(polymake, array_type);
 
diff --git a/deps/src/polymake_direct_calls.cpp b/deps/src/polymake_direct_calls.cpp
index 9b7319d1..bceeaabe 100644
--- a/deps/src/polymake_direct_calls.cpp
+++ b/deps/src/polymake_direct_calls.cpp
@@ -1,5 +1,5 @@
 #include "polymake_includes.h"
-
+#include <polymake/polytope/solve_LP.h>
 #include "polymake_type_modules.h"
 
 template<typename Scalar>
diff --git a/deps/src/polymake_includes.h b/deps/src/polymake_includes.h
index 86880ac1..45916021 100644
--- a/deps/src/polymake_includes.h
+++ b/deps/src/polymake_includes.h
@@ -21,7 +21,6 @@
 #include <polymake/TropicalNumber.h>
 #include <polymake/IncidenceMatrix.h>
 #include <polymake/Polynomial.h>
-#include <polymake/polytope/solve_LP.h>
 #include <polymake/SparseVector.h>
 
 #include <polymake/perl/calls.h>
diff --git a/deps/src/polymake_type_modules.h b/deps/src/polymake_type_modules.h
index 4fd88274..684713b2 100644
--- a/deps/src/polymake_type_modules.h
+++ b/deps/src/polymake_type_modules.h
@@ -9,7 +9,8 @@ tparametric1 polymake_module_add_array(jlcxx::Module& polymake);
 void polymake_module_add_array_polynomial(jlcxx::Module& polymake, tparametric1 arrayt);
 
 void polymake_module_add_bigobject(jlcxx::Module& polymake);
-void polymake_module_add_direct_calls(jlcxx::Module&);
+void polymake_module_add_direct_calls(jlcxx::Module& polymake);
+void polymake_module_add_beneath_beyond(jlcxx::Module& polymake);
 void polymake_module_add_incidencematrix(jlcxx::Module& polymake);
 void polymake_module_add_integer(jlcxx::Module& polymake);
 void polymake_module_add_matrix(jlcxx::Module& polymake);
diff --git a/src/Polymake.jl b/src/Polymake.jl
index 70be81cc..ac9c9982 100644
--- a/src/Polymake.jl
+++ b/src/Polymake.jl
@@ -144,6 +144,13 @@ include("tropicalnumber.jl")
 include("polynomial.jl")
 
 include("polymake_direct_calls.jl")
+include("beneath_beyond.jl")
+
+using Base.Docs
+using Markdown
+include("meta.jl")
+
+using Polymake.Meta
 
 include("generate_applications.jl")
 
diff --git a/src/beneath_beyond.jl b/src/beneath_beyond.jl
new file mode 100644
index 00000000..f881e1df
--- /dev/null
+++ b/src/beneath_beyond.jl
@@ -0,0 +1,46 @@
+mutable struct BeneathBeyond{T}
+    algo::Polymake.BeneathBeyondAlgoAllocated{T}
+    rays::Matrix{T}
+    lineality::Matrix{T}
+    perm::Vector{Int}
+
+    function BeneathBeyond(
+        rays::AbstractMatrix{T},
+        lineality::AbstractMatrix{T}=similar(rays, (0,size(rays,2)));
+        perm = collect(0:size(rays, 1)-1),
+        redundant = true,
+        triangulation = true,
+        iscone = true,
+        vertices = false,
+    ) where {T}
+
+        @assert !isempty(perm)
+
+        algo = BeneathBeyondAlgo{T}()
+        R = convert(Matrix{T}, rays)
+        L = convert(Matrix{T}, lineality)
+        p = convert(Vector{T}, perm)
+        bb = new{T}(algo, R, L, p)
+
+        bb_expecting_redundant(bb.algo, redundant)
+        bb_making_triangulation(bb.algo, triangulation)
+        bb_for_cone(bb.algo, iscone)
+        bb_computing_vertices(bb.algo, vertices)
+
+        bb_initialize!(bb.algo, R, L)
+
+        # finalizer(x->bb_clear!(x.algo), bb)
+
+        return bb
+    end
+end
+
+Base.length(bb::BeneathBeyond) = length(bb.perm)
+function Base.iterate(bb::BeneathBeyond, s = 1)
+    if s > length(bb)
+        return nothing
+    end
+    _, time, _ = @timed bb_add_point!(bb.algo, bb.perm[s])
+    @info "Iteration $s; adding point $(bb.perm[s]+1), time = $time (s)"
+    return nothing, s + 1
+end
diff --git a/src/generate_applications.jl b/src/generate_applications.jl
index 79fd9704..43c75f92 100644
--- a/src/generate_applications.jl
+++ b/src/generate_applications.jl
@@ -1,9 +1,3 @@
-using Base.Docs
-using Markdown
-include("meta.jl")
-
-using Polymake.Meta
-
 json_dir = joinpath(@__DIR__, "..", "deps", "json")
 generated_dir = joinpath(@__DIR__, "generated")
 isdir(generated_dir) || mkpath(generated_dir)