diff --git a/Project.toml b/Project.toml
index 9eae059..c439f7a 100644
--- a/Project.toml
+++ b/Project.toml
@@ -9,6 +9,7 @@ DataStructures = "864edb3b-99cc-5e75-8d2d-829cb0a9cfe8"
 HDF5 = "f67ccb44-e63f-5c2f-98bd-6dc0ccc4ba2f"
 HiGHS = "87dc4568-4c63-4d18-b0c0-bb2238e4078b"
 JLD2 = "033835bb-8acc-5ee8-8aae-3f567f8a3819"
+JSON = "682c06a0-de6a-54ab-a142-c8b1cf79cde6"
 JuMP = "4076af6c-e467-56ae-b986-b466b2749572"
 KLU = "ef3ab10e-7fda-4108-b977-705223b18434"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
@@ -23,17 +24,17 @@ Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
 TimerOutputs = "a759f4b9-e2f1-59dc-863e-4aeb61b1ea8f"
 
 [compat]
+Conda = "1"
+DataStructures = "0.18"
+HDF5 = "0.16"
+HiGHS = "1"
+JLD2 = "0.4"
+JuMP = "1"
+KLU = "0.4"
+MathOptInterface = "1"
+OrderedCollections = "1"
+PyCall = "1"
+Requires = "1"
+Statistics = "1"
+TimerOutputs = "0.5"
 julia = "1"
-Conda="1"
-DataStructures="0.18"
-HDF5="0.16"
-HiGHS="1"
-JLD2="0.4"
-JuMP="1"
-KLU="0.4"
-MathOptInterface="1"
-OrderedCollections="1"
-PyCall="1"
-Requires="1"
-Statistics="1"
-TimerOutputs="0.5"
\ No newline at end of file
diff --git a/src/MIPLearn.jl b/src/MIPLearn.jl
index 110ad27..ed670ac 100644
--- a/src/MIPLearn.jl
+++ b/src/MIPLearn.jl
@@ -12,6 +12,7 @@ include("components.jl")
 include("extractors.jl")
 include("io.jl")
 include("problems/setcover.jl")
+include("problems/stab.jl")
 include("solvers/jump.jl")
 include("solvers/learning.jl")
 
@@ -21,6 +22,7 @@ function __init__()
     __init_extractors__()
     __init_io__()
     __init_problems_setcover__()
+    __init_problems_stab__()
     __init_solvers_jump__()
     __init_solvers_learning__()
 end
diff --git a/src/components.jl b/src/components.jl
index 73172bf..9dc2442 100644
--- a/src/components.jl
+++ b/src/components.jl
@@ -2,19 +2,21 @@
 #  Copyright (C) 2020-2023, UChicago Argonne, LLC. All rights reserved.
 #  Released under the modified BSD license. See COPYING.md for more details.
 
-global MinProbabilityClassifier = PyNULL()
-global SingleClassFix = PyNULL()
-global PrimalComponentAction = PyNULL()
-global SetWarmStart = PyNULL()
-global FixVariables = PyNULL()
 global EnforceProximity = PyNULL()
 global ExpertPrimalComponent = PyNULL()
+global FixVariables = PyNULL()
 global IndependentVarsPrimalComponent = PyNULL()
 global JointVarsPrimalComponent = PyNULL()
-global SolutionConstructor = PyNULL()
+global MemorizingCutsComponent = PyNULL()
+global MemorizingLazyComponent = PyNULL()
 global MemorizingPrimalComponent = PyNULL()
-global SelectTopSolutions = PyNULL()
 global MergeTopSolutions = PyNULL()
+global MinProbabilityClassifier = PyNULL()
+global PrimalComponentAction = PyNULL()
+global SelectTopSolutions = PyNULL()
+global SetWarmStart = PyNULL()
+global SingleClassFix = PyNULL()
+global SolutionConstructor = PyNULL()
 
 function __init_components__()
     copy!(
@@ -51,6 +53,8 @@ function __init_components__()
     )
     copy!(SelectTopSolutions, pyimport("miplearn.components.primal.mem").SelectTopSolutions)
     copy!(MergeTopSolutions, pyimport("miplearn.components.primal.mem").MergeTopSolutions)
+    copy!(MemorizingCutsComponent, pyimport("miplearn.components.cuts.mem").MemorizingCutsComponent)
+    copy!(MemorizingLazyComponent, pyimport("miplearn.components.lazy.mem").MemorizingLazyComponent)
 end
 
 export MinProbabilityClassifier,
@@ -65,4 +69,6 @@ export MinProbabilityClassifier,
     SolutionConstructor,
     MemorizingPrimalComponent,
     SelectTopSolutions,
-    MergeTopSolutions
+    MergeTopSolutions,
+    MemorizingCutsComponent,
+    MemorizingLazyComponent
diff --git a/src/problems/setcover.jl b/src/problems/setcover.jl
index fb44b63..9d03b6a 100644
--- a/src/problems/setcover.jl
+++ b/src/problems/setcover.jl
@@ -13,12 +13,11 @@ function __init_problems_setcover__()
     copy!(SetCoverGenerator, pyimport("miplearn.problems.setcover").SetCoverGenerator)
 end
 
-function build_setcover_model(data::Any; optimizer = HiGHS.Optimizer)
+function build_setcover_model_jump(data::Any; optimizer = HiGHS.Optimizer)
     if data isa String
         data = read_pkl_gz(data)
     end
     model = Model(optimizer)
-    set_silent(model)
     n_elements, n_sets = size(data.incidence_matrix)
     E = 0:n_elements-1
     S = 0:n_sets-1
@@ -32,4 +31,4 @@ function build_setcover_model(data::Any; optimizer = HiGHS.Optimizer)
     return JumpModel(model)
 end
 
-export SetCoverData, SetCoverGenerator, build_setcover_model
+export SetCoverData, SetCoverGenerator, build_setcover_model_jump
diff --git a/src/problems/stab.jl b/src/problems/stab.jl
new file mode 100644
index 0000000..831b227
--- /dev/null
+++ b/src/problems/stab.jl
@@ -0,0 +1,60 @@
+#  MIPLearn: Extensible Framework for Learning-Enhanced Mixed-Integer Optimization
+#  Copyright (C) 2020-2024, UChicago Argonne, LLC. All rights reserved.
+#  Released under the modified BSD license. See COPYING.md for more details.
+
+using JuMP
+using HiGHS
+
+global MaxWeightStableSetData = PyNULL()
+global MaxWeightStableSetGenerator = PyNULL()
+
+function __init_problems_stab__()
+    copy!(MaxWeightStableSetData, pyimport("miplearn.problems.stab").MaxWeightStableSetData)
+    copy!(MaxWeightStableSetGenerator, pyimport("miplearn.problems.stab").MaxWeightStableSetGenerator)
+end
+
+function build_stab_model_jump(data::Any; optimizer=HiGHS.Optimizer)
+    nx = pyimport("networkx")
+
+    if data isa String
+        data = read_pkl_gz(data)
+    end
+    model = Model(optimizer)
+
+    # Variables and objective function
+    nodes = data.graph.nodes
+    x = @variable(model, x[nodes], Bin)
+    @objective(model, Min, sum(-data.weights[i+1] * x[i] for i in nodes))
+
+    # Edge inequalities
+    for (i1, i2) in data.graph.edges
+        @constraint(model, x[i1] + x[i2] <= 1, base_name = "eq_edge[$i1,$i2]")
+    end
+
+    function cuts_separate(cb_data)
+        x_val = callback_value.(Ref(cb_data), x)
+        violations = []
+        for clique in nx.find_cliques(data.graph)
+            if sum(x_val[i] for i in clique) > 1.0001
+                push!(violations, sort(clique))
+            end
+        end
+        return violations
+    end
+
+    function cuts_enforce(violations)
+        @info "Adding $(length(violations)) clique cuts..."
+        for clique in violations
+            constr = @build_constraint(sum(x[i] for i in clique) <= 1)
+            submit(model, constr)
+        end
+    end
+
+    return JumpModel(
+        model,
+        cuts_separate=cuts_separate,
+        cuts_enforce=cuts_enforce,
+    )
+end
+
+export MaxWeightStableSetData, MaxWeightStableSetGenerator, build_stab_model_jump
diff --git a/src/solvers/jump.jl b/src/solvers/jump.jl
index 726c449..b7bb2b2 100644
--- a/src/solvers/jump.jl
+++ b/src/solvers/jump.jl
@@ -4,9 +4,19 @@
 
 using JuMP
 using HiGHS
+using JSON
 
 global JumpModel = PyNULL()
 
+Base.@kwdef mutable struct _JumpModelExtData
+    aot_cuts = nothing
+    cb_data = nothing
+    cuts = []
+    where::Symbol = :WHERE_DEFAULT
+    cuts_enforce::Union{Function,Nothing} = nothing
+    cuts_separate::Union{Function,Nothing} = nothing
+end
+
 # -----------------------------------------------------------------------------
 
 function _add_constrs(
@@ -35,6 +45,15 @@ function _add_constrs(
     end
 end
 
+function submit(model::JuMP.Model, constr)
+    ext = model.ext[:miplearn]
+    if ext.where == :WHERE_CUTS
+        MOI.submit(model, MOI.UserCut(ext.cb_data), constr)
+    else
+        error("not implemented")
+    end
+end
+
 function _extract_after_load(model::JuMP.Model, h5)
     if JuMP.objective_sense(model) == MOI.MIN_SENSE
         h5.put_scalar("static_sense", "min")
@@ -109,6 +128,9 @@ function _extract_after_load_constrs(model::JuMP.Model, h5)
             end
         end
     end
+    if isempty(names)
+        error("no model constraints found; note that MIPLearn ignores unnamed constraints")
+    end
     lhs = sparse(lhs_rows, lhs_cols, lhs_values, length(rhs), JuMP.num_variables(model))
     h5.put_sparse("static_constr_lhs", lhs)
     h5.put_array("static_constr_rhs", rhs)
@@ -249,17 +271,50 @@ function _extract_after_mip(model::JuMP.Model, h5)
     rhs = h5.get_array("static_constr_rhs")
     slacks = abs.(lhs * x - rhs)
     h5.put_array("mip_constr_slacks", slacks)
+
+    # Cuts
+    ext = model.ext[:miplearn]
+    h5.put_scalar("mip_cuts", JSON.json(ext.cuts))
 end
 
 function _fix_variables(model::JuMP.Model, var_names, var_values, stats)
     vars = [variable_by_name(model, v) for v in var_names]
     for (i, var) in enumerate(vars)
-        fix(var, var_values[i], force = true)
+        fix(var, var_values[i], force=true)
     end
 end
 
 function _optimize(model::JuMP.Model)
+    # Set up cut callbacks
+    ext = model.ext[:miplearn]
+    ext.cuts = []
+    function cut_callback(cb_data)
+        ext.cb_data = cb_data
+        ext.where = :WHERE_CUTS
+        if ext.aot_cuts !== nothing
+            @info "Enforcing $(length(ext.aot_cuts)) cuts ahead-of-time..."
+            violations = ext.aot_cuts
+            ext.aot_cuts = nothing
+        else
+            violations = ext.cuts_separate(cb_data)
+            for v in violations
+                push!(ext.cuts, v)
+            end
+        end
+        if !isempty(violations)
+            ext.cuts_enforce(violations)
+        end
+    end
+    if ext.cuts_separate !== nothing
+        set_attribute(model, MOI.UserCutCallback(), cut_callback)
+    end
+
+    # Optimize
+    ext.where = :WHERE_DEFAULT
     optimize!(model)
+
+    # Cleanup
+    ext.cb_data = nothing
     flush(stdout)
     Libc.flush_cstdio()
 end
@@ -291,10 +346,21 @@ end
 # -----------------------------------------------------------------------------
 
 function __init_solvers_jump__()
-    @pydef mutable struct Class
+    AbstractModel = pyimport("miplearn.solvers.abstract").AbstractModel
+    @pydef mutable struct Class <: AbstractModel
 
-        function __init__(self, inner)
+        function __init__(
+            self,
+            inner;
+            cuts_enforce::Union{Function,Nothing}=nothing,
+            cuts_separate::Union{Function,Nothing}=nothing,
+        )
+            AbstractModel.__init__(self)
             self.inner = inner
+            self.inner.ext[:miplearn] = _JumpModelExtData(
+                cuts_enforce=cuts_enforce,
+                cuts_separate=cuts_separate,
+            )
         end
 
         add_constrs(
@@ -303,7 +369,7 @@ function __init_solvers_jump__()
             constrs_lhs,
             constrs_sense,
             constrs_rhs,
-            stats = nothing,
+            stats=nothing,
         ) = _add_constrs(
             self.inner,
             from_str_array(var_names),
@@ -319,17 +385,21 @@ function __init_solvers_jump__()
 
         extract_after_mip(self, h5) = _extract_after_mip(self.inner, h5)
 
-        fix_variables(self, var_names, var_values, stats = nothing) =
+        fix_variables(self, var_names, var_values, stats=nothing) =
             _fix_variables(self.inner, from_str_array(var_names), var_values, stats)
 
         optimize(self) = _optimize(self.inner)
 
         relax(self) = Class(_relax(self.inner))
 
-        set_warm_starts(self, var_names, var_values, stats = nothing) =
+        set_warm_starts(self, var_names, var_values, stats=nothing) =
             _set_warm_starts(self.inner, from_str_array(var_names), var_values, stats)
 
         write(self, filename) = _write(self.inner, filename)
+
+        function set_cuts(self, cuts)
+            self.inner.ext[:miplearn].aot_cuts = cuts
+        end
     end
     copy!(JumpModel, Class)
 end
diff --git a/test/Project.toml b/test/Project.toml
index dce63f3..0854405 100644
--- a/test/Project.toml
+++ b/test/Project.toml
@@ -15,6 +15,7 @@ Logging = "56ddb016-857b-54e1-b83d-db4d58db5568"
 MIPLearn = "2b1277c3-b477-4c49-a15e-7ba350325c68"
 PyCall = "438e738f-606a-5dbb-bf0a-cddfbfd45ab0"
 Revise = "295af30f-e4ad-537b-8983-00126c2a3abe"
+SCIP = "82193955-e24f-5292-bf16-6f2c5261a85f"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
 [compat]
diff --git a/test/fixtures/bell5.h5 b/test/fixtures/bell5.h5
index 73bc778..afff1ab 100644
Binary files a/test/fixtures/bell5.h5 and b/test/fixtures/bell5.h5 differ
diff --git a/test/fixtures/stab-n50-00000.h5 b/test/fixtures/stab-n50-00000.h5
new file mode 100644
index 0000000..f8ae28a
Binary files /dev/null and b/test/fixtures/stab-n50-00000.h5 differ
diff --git a/test/fixtures/stab-n50-00000.pkl.gz b/test/fixtures/stab-n50-00000.pkl.gz
new file mode 100644
index 0000000..19bddd8
Binary files /dev/null and b/test/fixtures/stab-n50-00000.pkl.gz differ
diff --git a/test/src/MIPLearnT.jl b/test/src/MIPLearnT.jl
index b14b883..e3524bb 100644
--- a/test/src/MIPLearnT.jl
+++ b/test/src/MIPLearnT.jl
@@ -16,10 +16,12 @@ FIXTURES = "$BASEDIR/../fixtures"
 include("fixtures.jl")
 
 include("BB/test_bb.jl")
+include("components/test_cuts.jl")
 include("Cuts/BlackBox/test_cplex.jl")
 include("problems/test_setcover.jl")
-include("test_io.jl")
+include("problems/test_stab.jl")
 include("solvers/test_jump.jl")
+include("test_io.jl")
 include("test_usage.jl")
 
 function runtests()
@@ -27,17 +29,18 @@ function runtests()
         @testset "BB" begin
             test_bb()
         end
-        # test_cuts_blackbox_cplex()
         test_io()
         test_problems_setcover()
+        test_problems_stab()
         test_solvers_jump()
         test_usage()
+        test_cuts()
     end
 end
 
 function format()
-    JuliaFormatter.format(BASEDIR, verbose = true)
-    JuliaFormatter.format("$BASEDIR/../../src", verbose = true)
+    JuliaFormatter.format(BASEDIR, verbose=true)
+    JuliaFormatter.format("$BASEDIR/../../src", verbose=true)
     return
 end
 
diff --git a/test/src/components/test_cuts.jl b/test/src/components/test_cuts.jl
new file mode 100644
index 0000000..1562648
--- /dev/null
+++ b/test/src/components/test_cuts.jl
@@ -0,0 +1,45 @@
+#  MIPLearn: Extensible Framework for Learning-Enhanced Mixed-Integer Optimization
+#  Copyright (C) 2020-2024, UChicago Argonne, LLC. All rights reserved.
+#  Released under the modified BSD license. See COPYING.md for more details.
+
+using SCIP
+
+function gen_stab()
+    np = pyimport("numpy")
+    uniform = pyimport("scipy.stats").uniform
+    randint = pyimport("scipy.stats").randint
+    np.random.seed(42)
+    gen = MaxWeightStableSetGenerator(
+        w=uniform(10.0, scale=1.0),
+        n=randint(low=50, high=51),
+        p=uniform(loc=0.5, scale=0.0),
+        fix_graph=true,
+    )
+    data = gen.generate(1)
+    data_filenames = write_pkl_gz(data, "$BASEDIR/../fixtures", prefix="stab-n50-")
+    collector = BasicCollector(write_mps=false)
+    collector.collect(
+        data_filenames,
+        data -> build_stab_model_jump(data, optimizer=SCIP.Optimizer),
+        progress=true,
+        verbose=true,
+    )
+end
+
+function test_cuts()
+    data_filenames = ["$BASEDIR/../fixtures/stab-n50-0000$i.pkl.gz" for i in 0:0]
+    clf = pyimport("sklearn.neighbors").KNeighborsClassifier(n_neighbors=1)
+    extractor = H5FieldsExtractor(
+        instance_fields=["static_var_obj_coeffs"],
+    )
+    comp = MemorizingCutsComponent(clf=clf, extractor=extractor)
+    solver = LearningSolver(components=[comp])
+    solver.fit(data_filenames)
+    @show comp.n_features_
+    @show comp.n_targets_
+    stats = solver.optimize(
+        data_filenames[1],
+        data -> build_stab_model_jump(data, optimizer=SCIP.Optimizer),
+    )
+    @test stats["Cuts: AOT"] > 0
+end
diff --git a/test/src/fixtures.jl b/test/src/fixtures.jl
index 56f4dec..b7cc084 100644
--- a/test/src/fixtures.jl
+++ b/test/src/fixtures.jl
@@ -14,5 +14,5 @@ function fixture_setcover_data()
 end
 
 function fixture_setcover_model()
-    return build_setcover_model(fixture_setcover_data())
+    return build_setcover_model_jump(fixture_setcover_data())
 end
diff --git a/test/src/problems/test_setcover.jl b/test/src/problems/test_setcover.jl
index 5f85418..9b92f19 100644
--- a/test/src/problems/test_setcover.jl
+++ b/test/src/problems/test_setcover.jl
@@ -51,7 +51,7 @@ function test_problems_setcover_model()
     )
 
     h5 = H5File(tempname(), "w")
-    model = build_setcover_model(data)
+    model = build_setcover_model_jump(data)
     model.extract_after_load(h5)
     model.optimize()
     model.extract_after_mip(h5)
diff --git a/test/src/problems/test_stab.jl b/test/src/problems/test_stab.jl
new file mode 100644
index 0000000..60f27b3
--- /dev/null
+++ b/test/src/problems/test_stab.jl
@@ -0,0 +1,27 @@
+#  MIPLearn: Extensible Framework for Learning-Enhanced Mixed-Integer Optimization
+#  Copyright (C) 2020-2024, UChicago Argonne, LLC. All rights reserved.
+#  Released under the modified BSD license. See COPYING.md for more details.
+
+using PyCall
+using SCIP
+
+function test_problems_stab()
+    test_problems_stab_1()
+    test_problems_stab_2()
+end
+
+function test_problems_stab_1()
+    nx = pyimport("networkx")
+    data = MaxWeightStableSetData(
+        graph=nx.gnp_random_graph(25, 0.5, seed=42),
+        weights=repeat([1.0], 25),
+    )
+    h5 = H5File(tempname(), "w")
+    model = build_stab_model_jump(data, optimizer=SCIP.Optimizer)
+    model.extract_after_load(h5)
+    model.optimize()
+    model.extract_after_mip(h5)
+    @test h5.get_scalar("mip_obj_value") == -6
+    @test h5.get_scalar("mip_cuts")[1:20] == "[[0,8,11,13],[0,8,13"
+    h5.close()
+end
diff --git a/test/src/test_usage.jl b/test/src/test_usage.jl
index f02ff65..9a8cb7b 100644
--- a/test/src/test_usage.jl
+++ b/test/src/test_usage.jl
@@ -29,13 +29,13 @@ function test_usage()
 
     @debug "Collecting training data..."
     bc = BasicCollector()
-    bc.collect(data_filenames, build_setcover_model)
+    bc.collect(data_filenames, build_setcover_model_jump)
 
     @debug "Training models..."
     solver.fit(data_filenames)
 
     @debug "Solving model..."
-    solver.optimize(data_filenames[1], build_setcover_model)
+    solver.optimize(data_filenames[1], build_setcover_model_jump)
 
     @debug "Checking solution..."
     h5 = H5File(h5_filenames[1])