diff --git a/python/src/broad_phase/broad_phase.cpp b/python/src/broad_phase/broad_phase.cpp
index f62d5562c..aa8f3f0ab 100644
--- a/python/src/broad_phase/broad_phase.cpp
+++ b/python/src/broad_phase/broad_phase.cpp
@@ -6,9 +6,145 @@
 namespace py = pybind11;
 using namespace ipc;
 
+class PyBroadPhase : public BroadPhase {
+public:
+    using BroadPhase::BroadPhase; // Inherit constructors
+
+    std::string name() const override
+    {
+        PYBIND11_OVERRIDE_PURE(std::string, BroadPhase, name);
+    }
+
+    void build(
+        const Eigen::MatrixXd& vertices,
+        const Eigen::MatrixXi& edges,
+        const Eigen::MatrixXi& faces,
+        const double inflation_radius = 0) override
+    {
+        PYBIND11_OVERRIDE(
+            void, BroadPhase, build, vertices, edges, faces, inflation_radius);
+    }
+
+    void build(
+        const Eigen::MatrixXd& vertices_t0,
+        const Eigen::MatrixXd& vertices_t1,
+        const Eigen::MatrixXi& edges,
+        const Eigen::MatrixXi& faces,
+        const double inflation_radius = 0) override
+    {
+        PYBIND11_OVERRIDE(
+            void, BroadPhase, build, vertices_t0, vertices_t1, edges, faces,
+            inflation_radius);
+    }
+
+    void clear() override { PYBIND11_OVERRIDE(void, BroadPhase, clear); }
+
+    void detect_vertex_vertex_candidates(
+        std::vector<VertexVertexCandidate>& candidates) const override
+    {
+        py::gil_scoped_acquire gil; // Acquire GIL before calling Python code
+        py::function override =
+            py::get_override(this, "detect_vertex_vertex_candidates");
+        if (override) {
+            candidates = override().cast<std::vector<VertexVertexCandidate>>();
+            return;
+        }
+        throw std::runtime_error(
+            "Tried to call pure virtual function \"BroadPhase::detect_vertex_vertex_candidates\"");
+    }
+
+    void detect_edge_vertex_candidates(
+        std::vector<EdgeVertexCandidate>& candidates) const override
+    {
+        py::gil_scoped_acquire gil; // Acquire GIL before calling Python code
+        py::function override =
+            py::get_override(this, "detect_edge_vertex_candidates");
+        if (override) {
+            candidates = override().cast<std::vector<EdgeVertexCandidate>>();
+            return;
+        }
+        throw std::runtime_error(
+            "Tried to call pure virtual function \"BroadPhase::detect_edge_vertex_candidates\"");
+    }
+
+    void detect_edge_edge_candidates(
+        std::vector<EdgeEdgeCandidate>& candidates) const override
+    {
+        py::gil_scoped_acquire gil; // Acquire GIL before calling Python code
+        py::function override =
+            py::get_override(this, "detect_edge_edge_candidates");
+        if (override) {
+            candidates = override().cast<std::vector<EdgeEdgeCandidate>>();
+            return;
+        }
+        throw std::runtime_error(
+            "Tried to call pure virtual function \"BroadPhase::detect_edge_edge_candidates\"");
+    }
+
+    void detect_face_vertex_candidates(
+        std::vector<FaceVertexCandidate>& candidates) const override
+    {
+        py::gil_scoped_acquire gil; // Acquire GIL before calling Python code
+        py::function override =
+            py::get_override(this, "detect_face_vertex_candidates");
+        if (override) {
+            candidates = override().cast<std::vector<FaceVertexCandidate>>();
+            return;
+        }
+        throw std::runtime_error(
+            "Tried to call pure virtual function \"BroadPhase::detect_face_vertex_candidates\"");
+    }
+
+    void detect_edge_face_candidates(
+        std::vector<EdgeFaceCandidate>& candidates) const override
+    {
+        py::gil_scoped_acquire gil; // Acquire GIL before calling Python code
+        py::function override =
+            py::get_override(this, "detect_edge_face_candidates");
+        if (override) {
+            candidates = override().cast<std::vector<EdgeFaceCandidate>>();
+            return;
+        }
+        throw std::runtime_error(
+            "Tried to call pure virtual function \"BroadPhase::detect_edge_face_candidates\"");
+    }
+
+    void detect_face_face_candidates(
+        std::vector<FaceFaceCandidate>& candidates) const override
+    {
+        py::gil_scoped_acquire gil; // Acquire GIL before calling Python code
+        py::function override =
+            py::get_override(this, "detect_face_face_candidates");
+        if (override) {
+            candidates = override().cast<std::vector<FaceFaceCandidate>>();
+            return;
+        }
+        throw std::runtime_error(
+            "Tried to call pure virtual function \"BroadPhase::detect_face_face_candidates\"");
+    }
+
+    void
+    detect_collision_candidates(int dim, Candidates& candidates) const override
+    {
+        {
+            py::gil_scoped_acquire
+                gil; // Acquire GIL before calling Python code
+            py::function override =
+                py::get_override(this, "detect_collision_candidates");
+            if (override) {
+                candidates = override(dim).cast<Candidates>();
+                return;
+            }
+        }
+        BroadPhase::detect_collision_candidates(dim, candidates);
+    }
+};
+
 void define_broad_phase(py::module_& m)
 {
-    py::class_<BroadPhase, std::shared_ptr<BroadPhase>>(m, "BroadPhase")
+    py::class_<BroadPhase, PyBroadPhase, std::shared_ptr<BroadPhase>>(
+        m, "BroadPhase")
+        .def(py::init<>())
         .def("name", &BroadPhase::name, "Get the name of the broad phase.")
         .def(
             "build",
diff --git a/python/src/candidates/edge_edge.cpp b/python/src/candidates/edge_edge.cpp
index 6a57280ad..d4fb6ce6f 100644
--- a/python/src/candidates/edge_edge.cpp
+++ b/python/src/candidates/edge_edge.cpp
@@ -11,6 +11,12 @@ void define_edge_edge_candidate(py::module_& m)
         EdgeEdgeCandidate, CollisionStencil, ContinuousCollisionCandidate>(
         m, "EdgeEdgeCandidate")
         .def(py::init<long, long>(), py::arg("edge0_id"), py::arg("edge1_id"))
+        .def(
+            py::init([](std::tuple<long, long> edge_ids) {
+                return std::make_unique<EdgeEdgeCandidate>(
+                    std::get<0>(edge_ids), std::get<1>(edge_ids));
+            }),
+            py::arg("edge_ids"))
         .def("known_dtype", &EdgeEdgeCandidate::known_dtype)
         .def(
             "__str__",
@@ -32,4 +38,6 @@ void define_edge_edge_candidate(py::module_& m)
             "edge0_id", &EdgeEdgeCandidate::edge0_id, "ID of the first edge.")
         .def_readwrite(
             "edge1_id", &EdgeEdgeCandidate::edge1_id, "ID of the second edge.");
+
+    py::implicitly_convertible<std::tuple<long, long>, EdgeEdgeCandidate>();
 }
diff --git a/python/src/candidates/edge_face.cpp b/python/src/candidates/edge_face.cpp
index 98b606434..86f0a6851 100644
--- a/python/src/candidates/edge_face.cpp
+++ b/python/src/candidates/edge_face.cpp
@@ -10,6 +10,13 @@ void define_edge_face_candidate(py::module_& m)
 {
     py::class_<EdgeFaceCandidate>(m, "EdgeFaceCandidate")
         .def(py::init<long, long>(), py::arg("edge_id"), py::arg("face_id"))
+        .def(
+            py::init([](std::tuple<long, long> edge_and_face_id) {
+                return std::make_unique<EdgeFaceCandidate>(
+                    std::get<0>(edge_and_face_id),
+                    std::get<1>(edge_and_face_id));
+            }),
+            py::arg("edge_and_face_id"))
         .def(
             "__str__",
             [](const EdgeFaceCandidate& ev) {
@@ -29,4 +36,6 @@ void define_edge_face_candidate(py::module_& m)
         .def_readwrite("edge_id", &EdgeFaceCandidate::edge_id, "ID of the edge")
         .def_readwrite(
             "face_id", &EdgeFaceCandidate::face_id, "ID of the face");
+
+    py::implicitly_convertible<std::tuple<long, long>, EdgeFaceCandidate>();
 }
diff --git a/python/src/candidates/edge_vertex.cpp b/python/src/candidates/edge_vertex.cpp
index ac1d60fa5..3d5e82296 100644
--- a/python/src/candidates/edge_vertex.cpp
+++ b/python/src/candidates/edge_vertex.cpp
@@ -11,6 +11,13 @@ void define_edge_vertex_candidate(py::module_& m)
         EdgeVertexCandidate, CollisionStencil, ContinuousCollisionCandidate>(
         m, "EdgeVertexCandidate")
         .def(py::init<long, long>(), py::arg("edge_id"), py::arg("vertex_id"))
+        .def(
+            py::init([](std::tuple<long, long> edge_and_vertex_id) {
+                return std::make_unique<EdgeVertexCandidate>(
+                    std::get<0>(edge_and_vertex_id),
+                    std::get<1>(edge_and_vertex_id));
+            }),
+            py::arg("edge_and_vertex_id"))
         .def("known_dtype", &EdgeVertexCandidate::known_dtype)
         .def(
             "__str__",
@@ -33,4 +40,6 @@ void define_edge_vertex_candidate(py::module_& m)
             "edge_id", &EdgeVertexCandidate::edge_id, "ID of the edge")
         .def_readwrite(
             "vertex_id", &EdgeVertexCandidate::vertex_id, "ID of the vertex");
+
+    py::implicitly_convertible<std::tuple<long, long>, EdgeVertexCandidate>();
 }
diff --git a/python/src/candidates/face_face.cpp b/python/src/candidates/face_face.cpp
index dd0f80e42..c2ce5d348 100644
--- a/python/src/candidates/face_face.cpp
+++ b/python/src/candidates/face_face.cpp
@@ -10,6 +10,12 @@ void define_face_face_candidate(py::module_& m)
 {
     py::class_<FaceFaceCandidate>(m, "FaceFaceCandidate")
         .def(py::init<long, long>(), py::arg("face0_id"), py::arg("face1_id"))
+        .def(
+            py::init([](std::tuple<long, long> face_ids) {
+                return std::make_unique<FaceFaceCandidate>(
+                    std::get<0>(face_ids), std::get<1>(face_ids));
+            }),
+            py::arg("face_ids"))
         .def(
             "__str__",
             [](const FaceFaceCandidate& ff) {
@@ -30,4 +36,6 @@ void define_face_face_candidate(py::module_& m)
             "face0_id", &FaceFaceCandidate::face0_id, "ID of the first face.")
         .def_readwrite(
             "face1_id", &FaceFaceCandidate::face1_id, "ID of the second face.");
+
+    py::implicitly_convertible<std::tuple<long, long>, FaceFaceCandidate>();
 }
diff --git a/python/src/candidates/face_vertex.cpp b/python/src/candidates/face_vertex.cpp
index 89ba20e6f..3577729f2 100644
--- a/python/src/candidates/face_vertex.cpp
+++ b/python/src/candidates/face_vertex.cpp
@@ -11,6 +11,13 @@ void define_face_vertex_candidate(py::module_& m)
         FaceVertexCandidate, CollisionStencil, ContinuousCollisionCandidate>(
         m, "FaceVertexCandidate")
         .def(py::init<long, long>(), py::arg("face_id"), py::arg("vertex_id"))
+        .def(
+            py::init([](std::tuple<long, long> face_and_vertex_id) {
+                return std::make_unique<FaceVertexCandidate>(
+                    std::get<0>(face_and_vertex_id),
+                    std::get<1>(face_and_vertex_id));
+            }),
+            py::arg("face_and_vertex_id"))
         .def("known_dtype", &FaceVertexCandidate::known_dtype)
         .def(
             "__str__",
@@ -33,4 +40,6 @@ void define_face_vertex_candidate(py::module_& m)
             "face_id", &FaceVertexCandidate::face_id, "ID of the face")
         .def_readwrite(
             "vertex_id", &FaceVertexCandidate::vertex_id, "ID of the vertex");
+
+    py::implicitly_convertible<std::tuple<long, long>, FaceVertexCandidate>();
 }
diff --git a/python/src/candidates/vertex_vertex.cpp b/python/src/candidates/vertex_vertex.cpp
index f27ba9541..ba94e0eea 100644
--- a/python/src/candidates/vertex_vertex.cpp
+++ b/python/src/candidates/vertex_vertex.cpp
@@ -13,6 +13,12 @@ void define_vertex_vertex_candidate(py::module_& m)
         .def(
             py::init<long, long>(), py::arg("vertex0_id"),
             py::arg("vertex1_id"))
+        .def(
+            py::init([](std::tuple<long, long> vertex_ids) {
+                return std::make_unique<VertexVertexCandidate>(
+                    std::get<0>(vertex_ids), std::get<1>(vertex_ids));
+            }),
+            py::arg("vertex_ids"))
         .def(
             "__str__",
             [](const VertexVertexCandidate& ev) {
@@ -37,4 +43,6 @@ void define_vertex_vertex_candidate(py::module_& m)
         .def_readwrite(
             "vertex1_id", &VertexVertexCandidate::vertex1_id,
             "ID of the second vertex");
+
+    py::implicitly_convertible<std::tuple<long, long>, VertexVertexCandidate>();
 }
diff --git a/python/src/ccd/narrow_phase_ccd.cpp b/python/src/ccd/narrow_phase_ccd.cpp
index 930e1b796..8d876bb7c 100644
--- a/python/src/ccd/narrow_phase_ccd.cpp
+++ b/python/src/ccd/narrow_phase_ccd.cpp
@@ -5,9 +5,103 @@
 namespace py = pybind11;
 using namespace ipc;
 
+class PyNarrowPhaseCCD : public NarrowPhaseCCD {
+public:
+    using NarrowPhaseCCD::NarrowPhaseCCD; // Inherit constructors
+    bool point_point_ccd(
+        const VectorMax3d& p0_t0,
+        const VectorMax3d& p1_t0,
+        const VectorMax3d& p0_t1,
+        const VectorMax3d& p1_t1,
+        double& toi,
+        const double min_distance = 0.0,
+        const double tmax = 1.0) const override
+    {
+        py::gil_scoped_acquire gil; // Acquire GIL before calling Python code
+        py::function override = py::get_override(this, "point_point_ccd");
+        if (override) {
+            const py::tuple obj =
+                override(p0_t0, p1_t0, p0_t1, p1_t1, min_distance, tmax);
+            toi = obj[1].cast<double>();
+            return obj[0].cast<bool>();
+        }
+        throw std::runtime_error("pure virtual function called");
+    }
+    bool point_edge_ccd(
+        const VectorMax3d& p_t0,
+        const VectorMax3d& e0_t0,
+        const VectorMax3d& e1_t0,
+        const VectorMax3d& p_t1,
+        const VectorMax3d& e0_t1,
+        const VectorMax3d& e1_t1,
+        double& toi,
+        const double min_distance = 0.0,
+        const double tmax = 1.0) const override
+    {
+        py::gil_scoped_acquire gil; // Acquire GIL before calling Python code
+        py::function override = py::get_override(this, "point_edge_ccd");
+        if (override) {
+            const py::tuple obj = override(
+                p_t0, e0_t0, e1_t0, p_t1, e0_t1, e1_t1, min_distance, tmax);
+            toi = obj[1].cast<double>();
+            return obj[0].cast<bool>();
+        }
+        throw std::runtime_error("pure virtual function called");
+    }
+    bool point_triangle_ccd(
+        const Eigen::Vector3d& p_t0,
+        const Eigen::Vector3d& t0_t0,
+        const Eigen::Vector3d& t1_t0,
+        const Eigen::Vector3d& t2_t0,
+        const Eigen::Vector3d& p_t1,
+        const Eigen::Vector3d& t0_t1,
+        const Eigen::Vector3d& t1_t1,
+        const Eigen::Vector3d& t2_t1,
+        double& toi,
+        const double min_distance = 0.0,
+        const double tmax = 1.0) const override
+    {
+        py::gil_scoped_acquire gil; // Acquire GIL before calling Python code
+        py::function override = py::get_override(this, "point_triangle_ccd");
+        if (override) {
+            const py::tuple obj = override(
+                p_t0, t0_t0, t1_t0, t2_t0, p_t1, t0_t1, t1_t1, t2_t1,
+                min_distance, tmax);
+            toi = obj[1].cast<double>();
+            return obj[0].cast<bool>();
+        }
+        throw std::runtime_error("pure virtual function called");
+    }
+    bool edge_edge_ccd(
+        const Eigen::Vector3d& ea0_t0,
+        const Eigen::Vector3d& ea1_t0,
+        const Eigen::Vector3d& eb0_t0,
+        const Eigen::Vector3d& eb1_t0,
+        const Eigen::Vector3d& ea0_t1,
+        const Eigen::Vector3d& ea1_t1,
+        const Eigen::Vector3d& eb0_t1,
+        const Eigen::Vector3d& eb1_t1,
+        double& toi,
+        const double min_distance = 0.0,
+        const double tmax = 1.0) const override
+    {
+        py::gil_scoped_acquire gil; // Acquire GIL before calling Python code
+        py::function override = py::get_override(this, "edge_edge_ccd");
+        if (override) {
+            const py::tuple obj = override(
+                ea0_t0, ea1_t0, eb0_t0, eb1_t0, ea0_t1, ea1_t1, eb0_t1, eb1_t1,
+                min_distance, tmax);
+            toi = obj[1].cast<double>();
+            return obj[0].cast<bool>();
+        }
+        throw std::runtime_error("pure virtual function called");
+    }
+};
+
 void define_narrow_phase_ccd(py::module_& m)
 {
-    py::class_<NarrowPhaseCCD>(m, "NarrowPhaseCCD")
+    py::class_<NarrowPhaseCCD, PyNarrowPhaseCCD>(m, "NarrowPhaseCCD")
+        .def(py::init<>())
         .def(
             "point_point_ccd",
             [](const NarrowPhaseCCD& self, const VectorMax3d& p0_t0,
diff --git a/python/src/utils/thread_limiter.cpp b/python/src/utils/thread_limiter.cpp
index e3dd4190c..f142087e5 100644
--- a/python/src/utils/thread_limiter.cpp
+++ b/python/src/utils/thread_limiter.cpp
@@ -1,6 +1,7 @@
 #include <common.hpp>
 
 #include <memory>
+#include <cstdlib>
 
 #include <tbb/info.h>
 #include <tbb/global_control.h>
@@ -40,4 +41,6 @@ void define_thread_limiter(py::module_& m)
     m.def(
         "set_num_threads", &set_num_threads,
         "set maximum number of threads to use", py::arg("nthreads"));
+
+    std::atexit([]() { thread_limiter.reset(); });
 }
diff --git a/python/tests/test_candidates.py b/python/tests/test_candidates.py
new file mode 100644
index 000000000..974a4adc0
--- /dev/null
+++ b/python/tests/test_candidates.py
@@ -0,0 +1,18 @@
+import numpy as np
+
+from find_ipctk import ipctk
+from ipctk.filib import Interval
+
+from utils import load_mesh
+
+
+def test_candidates():
+    V0, E, F = load_mesh("two-cubes-close.ply")
+    V1, E, F = load_mesh("two-cubes-intersecting.ply")
+
+    mesh = ipctk.CollisionMesh(V0, E, F)
+
+    candidates = ipctk.Candidates()
+    candidates.build(mesh, V0, V1)
+
+    assert len(candidates) > 0, "No candidates generated."
\ No newline at end of file
diff --git a/python/tests/test_ccd.py b/python/tests/test_ccd.py
index 8b8f33607..aa0b0d2cb 100644
--- a/python/tests/test_ccd.py
+++ b/python/tests/test_ccd.py
@@ -18,11 +18,85 @@ def test_ccd():
     assert not ipctk.is_step_collision_free(mesh, V0, V1)
 
     toi = ipctk.compute_collision_free_stepsize(mesh, V0, V1)
-    assert 0 < toi < 1
+    assert 0 <= toi <= 1
 
     assert ipctk.is_step_collision_free(mesh, V0, (V1 - V0) * toi + V0)
 
 
+def test_custom_ccd():
+    class DumbNarrowPhaseCCD(ipctk.NarrowPhaseCCD):
+        def __init__(self):
+            ipctk.NarrowPhaseCCD.__init__(self)
+
+        def point_point_ccd(self, p0_t0, p1_t0, p0_t1, p1_t1, min_distance=0.0, tmax=1.0):
+            return True, 0.0
+
+        def point_edge_ccd(self, p_t0, e0_t0, e1_t0, p_t1, e0_t1, e1_t1, min_distance=0.0, tmax=1.0):
+            return True, 0.0
+
+        def point_triangle_ccd(self, p_t0, t0_t0, t1_t0, t2_t0, p_t1, t0_t1, t1_t1, t2_t1, min_distance=0.0, tmax=1.0):
+            return True, 0.0
+
+        def edge_edge_ccd(self, ea0_t0, ea1_t0, eb0_t0, eb1_t0, ea0_t1, ea1_t1, eb0_t1, eb1_t1, min_distance=0.0, tmax=1.0):
+            return True, 0.0
+
+    V0, E, F = load_mesh("two-cubes-close.ply")
+    V1, E, F = load_mesh("two-cubes-intersecting.ply")
+
+    mesh = ipctk.CollisionMesh(V0, E, F)
+
+    ipctk.set_num_threads(1)
+
+    assert not ipctk.is_step_collision_free(
+        mesh, V0, V1, narrow_phase_ccd=DumbNarrowPhaseCCD())
+
+    toi = ipctk.compute_collision_free_stepsize(
+        mesh, V0, V1, narrow_phase_ccd=DumbNarrowPhaseCCD())
+    assert 0 <= toi <= 1
+
+
+def test_custom_broad_phase():
+    V0, E, F = load_mesh("two-cubes-close.ply")
+    V1, E, F = load_mesh("two-cubes-intersecting.ply")
+
+    class DumbBroadPhase(ipctk.BroadPhase):
+        def __init__(self):
+            ipctk.BroadPhase.__init__(self)
+
+        def name(self):
+            return "DumbBroadPhase"
+
+        def detect_vertex_vertex_candidates(self):
+            return []
+
+        def detect_edge_vertex_candidates(self):
+            return []
+
+        def detect_edge_edge_candidates(self):
+            return [(0, 232)]
+
+        def detect_face_vertex_candidates(self):
+            return [(16, 205)]
+
+        def detect_edge_face_candidates(self):
+            return []
+
+        def detect_face_face_candidates(self):
+            return []
+
+    mesh = ipctk.CollisionMesh(V0, E, F)
+
+    broad_phase = DumbBroadPhase()
+
+    assert ipctk.is_step_collision_free(mesh, V0, V1, broad_phase=broad_phase)
+
+    toi = ipctk.compute_collision_free_stepsize(
+        mesh, V0, V1, broad_phase=broad_phase)
+    assert toi == 1
+
+    assert not ipctk.has_intersections(mesh, V0, broad_phase=broad_phase)
+
+
 def test_nonlinear_ccd():
     class LinearTrajectory(ipctk.NonlinearTrajectory):
         def __init__(self, p0, p1):