diff --git a/deps/build.jl b/deps/build.jl
index f3effee..4a0a082 100644
--- a/deps/build.jl
+++ b/deps/build.jl
@@ -5,7 +5,7 @@ function install_miplearn()
     Conda.update()
     pip = joinpath(dirname(pyimport("sys").executable), "pip")
     isfile(pip) || error("$pip: invalid path")
-    run(`$pip install miplearn==0.2.0.dev8`)
+    run(`$pip install miplearn==0.2.0.dev9`)
 end
 
 install_miplearn()
diff --git a/src/modeling/jump_instance.jl b/src/modeling/jump_instance.jl
index 465750c..8519f8f 100644
--- a/src/modeling/jump_instance.jl
+++ b/src/modeling/jump_instance.jl
@@ -7,6 +7,7 @@ using JuMP
 
 @pydef mutable struct PyJuMPInstance <: miplearn.Instance
     function __init__(self, model)
+        init_miplearn_ext(model)
         self.model = model
         self.samples = []
     end
@@ -22,25 +23,25 @@ using JuMP
     function get_variable_features(self, var_name)
         model = self.model
         v = variable_by_name(model, var_name)
-        return model.ext[:miplearn][:variable_features][v]
+        return get(model.ext[:miplearn][:variable_features], v, [0.0])
     end
 
     function get_variable_category(self, var_name)
         model = self.model
         v = variable_by_name(model, var_name)
-        return model.ext[:miplearn][:variable_categories][v]
+        return get(model.ext[:miplearn][:variable_categories], v, "default")
     end
 
     function get_constraint_features(self, cname)
         model = self.model
         c = constraint_by_name(model, cname)
-        return model.ext[:miplearn][:constraint_features][c]
+        return get(model.ext[:miplearn][:constraint_features], c, [0.0])
     end
 
     function get_constraint_category(self, cname)
         model = self.model
         c = constraint_by_name(model, cname)
-        return model.ext[:miplearn][:constraint_categories][c]
+        return get(model.ext[:miplearn][:constraint_categories], c, "default")
     end
 end
 
@@ -50,7 +51,8 @@ struct JuMPInstance
 end
 
 
-function JuMPInstance(model::Model)
+function JuMPInstance(model)
+    model isa Model || error("model should be a JuMP.Model. Found $(typeof(model)) instead.")
     return JuMPInstance(PyJuMPInstance(model))
 end
 
diff --git a/src/modeling/jump_solver.jl b/src/modeling/jump_solver.jl
index 8864f3d..9dc8fef 100644
--- a/src/modeling/jump_solver.jl
+++ b/src/modeling/jump_solver.jl
@@ -316,49 +316,49 @@ function get_variables(
     values, rc = nothing, nothing
 
     # Variable names
-    names = Tuple(JuMP.name.(vars))
+    names = JuMP.name.(vars)
 
     # Primal values
     if !isempty(data.solution)
-        values = Tuple([data.solution[v] for v in vars])
+        values = [data.solution[v] for v in vars]
     end
 
     if with_static
         # Lower bounds
-        lb = Tuple(
+        lb = [
             JuMP.is_binary(v) ? 0.0 :
                 JuMP.has_lower_bound(v) ? JuMP.lower_bound(v) :
                     -Inf
             for v in vars
-        )
+        ]
 
         # Upper bounds
-        ub = Tuple(
+        ub = [
             JuMP.is_binary(v) ? 1.0 :
                 JuMP.has_upper_bound(v) ? JuMP.upper_bound(v) :
                     Inf
             for v in vars
-        )
+        ]
 
         # Variable types
-        types = Tuple(
+        types = [
             JuMP.is_binary(v) ? "B" :
                 JuMP.is_integer(v) ? "I" :
                     "C"
             for v in vars
-        )
+        ]
 
         # Objective function coefficients
         obj = objective_function(data.model)
-        obj_coeffs = Tuple(
+        obj_coeffs = [
             v ∈ keys(obj.terms) ? obj.terms[v] : 0.0
             for v in vars
-        )
+        ]
     end
 
-    rc = isempty(data.reduced_costs) ? nothing : Tuple(data.reduced_costs)
+    rc = isempty(data.reduced_costs) ? nothing : data.reduced_costs
 
-    return miplearn.features.VariableFeatures(
+    vf = miplearn.features.VariableFeatures(
         names=names,
         lower_bounds=lb,
         upper_bounds=ub,
@@ -367,6 +367,7 @@ function get_variables(
         reduced_costs=rc,
         values=values,
     )
+    return vf
 end
 
 
@@ -406,7 +407,7 @@ function get_constraints(
                 if ftype == JuMP.AffExpr
                     push!(
                         lhs,
-                        Tuple(
+                        [
                             (
                                 MOI.get(
                                     constr.model.moi_backend,
@@ -420,7 +421,7 @@ function get_constraints(
                                 MOI.ConstraintFunction(),
                                 constr.index,
                             ).terms
-                        )
+                        ]
                     )
                     if stype == MOI.EqualTo{Float64}
                         push!(senses, "=")
@@ -441,17 +442,12 @@ function get_constraints(
         end
     end
 
-    function to_tuple(x)
-        x !== nothing || return nothing
-        return Tuple(x)
-    end
-
     return miplearn.features.ConstraintFeatures(
-        names=to_tuple(names),
-        senses=to_tuple(senses),
-        lhs=to_tuple(lhs),
-        rhs=to_tuple(rhs),
-        dual_values=to_tuple(dual_values),
+        names=names,
+        senses=senses,
+        lhs=lhs,
+        rhs=rhs,
+        dual_values=dual_values,
     )
 end
 
@@ -471,23 +467,35 @@ end
         )
     end
 
-    add_constraints(self, cf) =
+    function add_constraints(self, cf)
+        lhs = cf.lhs
+        if lhs isa Matrix
+            # Undo incorrect automatic conversion performed by PyCall
+            lhs = [col[:] for col in eachcol(lhs)]
+        end
         add_constraints(
             self.data,
-            lhs=[[term for term in constr] for constr in cf.lhs],
-            rhs=[r for r in cf.rhs],
-            senses=[s for s in cf.senses],
-            names=[n for n in cf.names],
+            lhs=lhs,
+            rhs=cf.rhs,
+            senses=cf.senses,
+            names=cf.names,
         )
+    end
 
-    are_constraints_satisfied(self, cf; tol=1e-5) =
-        tuple(are_constraints_satisfied(
+    function are_constraints_satisfied(self, cf; tol=1e-5)
+        lhs = cf.lhs
+        if lhs isa Matrix
+            # Undo incorrect automatic conversion performed by PyCall
+            lhs = [col[:] for col in eachcol(lhs)]
+        end
+        return are_constraints_satisfied(
             self.data,
-            lhs=[[term for term in constr] for constr in cf.lhs],
-            rhs=[r for r in cf.rhs],
-            senses=[s for s in cf.senses],
+            lhs=lhs,
+            rhs=cf.rhs,
+            senses=cf.senses,
             tol=tol,
-        )...)
+        )
+    end
 
     build_test_instance_infeasible(self) =
         build_test_instance_infeasible()
diff --git a/src/modeling/macros.jl b/src/modeling/macros.jl
index 716746b..b4755c6 100644
--- a/src/modeling/macros.jl
+++ b/src/modeling/macros.jl
@@ -5,6 +5,7 @@
 function init_miplearn_ext(model)::Dict
     if :miplearn ∉ keys(model.ext)
         model.ext[:miplearn] = Dict{Symbol, Any}()
+        model.ext[:miplearn][:instance_features] = [0.0]
         model.ext[:miplearn][:variable_features] = Dict{VariableRef, Vector{Float64}}()
         model.ext[:miplearn][:variable_categories] = Dict{VariableRef, String}()
         model.ext[:miplearn][:constraint_features] = Dict{ConstraintRef, Vector{Float64}}()
diff --git a/test/modeling/learning_solver_test.jl b/test/modeling/learning_solver_test.jl
index 8593a97..b2f79a3 100644
--- a/test/modeling/learning_solver_test.jl
+++ b/test/modeling/learning_solver_test.jl
@@ -6,43 +6,55 @@ using JuMP
 using MIPLearn
 using Gurobi
 
-@testset "macros" begin
-    weights = [1.0, 2.0, 3.0]
-    prices = [5.0, 6.0, 7.0]
-    capacity = 3.0
+@testset "LearningSolver" begin
+    @testset "model with annotations" begin
+        # Create standard JuMP model
+        weights = [1.0, 2.0, 3.0]
+        prices = [5.0, 6.0, 7.0]
+        capacity = 3.0
+        model = Model()
 
-    # Create standard JuMP model
-    model = Model()
-    n = length(weights)
-    @variable(model, x[1:n], Bin)
-    @objective(model, Max, sum(x[i] * prices[i] for i in 1:n))
-    @constraint(model, c1, sum(x[i] * weights[i] for i in 1:n) <= capacity)
+        n = length(weights)
+        @variable(model, x[1:n], Bin)
+        @objective(model, Max, sum(x[i] * prices[i] for i in 1:n))
+        @constraint(model, c1, sum(x[i] * weights[i] for i in 1:n) <= capacity)
 
-    # Add ML information to the model
-    @feature(model, [5.0])
-    @feature(c1, [1.0, 2.0, 3.0])
-    @category(c1, "c1")
-    for i in 1:n
-        @feature(x[i], [weights[i]; prices[i]])
-        @category(x[i], "type-$i")
+        # Add ML information to the model
+        @feature(model, [5.0])
+        @feature(c1, [1.0, 2.0, 3.0])
+        @category(c1, "c1")
+        for i in 1:n
+            @feature(x[i], [weights[i]; prices[i]])
+            @category(x[i], "type-$i")
+        end
+    
+        # Should store ML information
+        @test model.ext[:miplearn][:variable_features][x[1]] == [1.0, 5.0]
+        @test model.ext[:miplearn][:variable_features][x[2]] == [2.0, 6.0]
+        @test model.ext[:miplearn][:variable_features][x[3]] == [3.0, 7.0]
+        @test model.ext[:miplearn][:variable_categories][x[1]] == "type-1"
+        @test model.ext[:miplearn][:variable_categories][x[2]] == "type-2"
+        @test model.ext[:miplearn][:variable_categories][x[3]] == "type-3"
+        @test model.ext[:miplearn][:constraint_features][c1] == [1.0, 2.0, 3.0]
+        @test model.ext[:miplearn][:constraint_categories][c1] == "c1"
+        @test model.ext[:miplearn][:instance_features] == [5.0]
+
+        solver = LearningSolver(Gurobi.Optimizer)
+        instance = JuMPInstance(model)
+        stats = solve!(solver, instance)
+        @test stats["mip_lower_bound"] == 11.0
+        @test length(instance.py.samples) == 1
+        fit!(solver, [instance])
+        solve!(solver, instance)
     end
- 
-    # Should store ML information
-    @test model.ext[:miplearn][:variable_features][x[1]] == [1.0, 5.0]
-    @test model.ext[:miplearn][:variable_features][x[2]] == [2.0, 6.0]
-    @test model.ext[:miplearn][:variable_features][x[3]] == [3.0, 7.0]
-    @test model.ext[:miplearn][:variable_categories][x[1]] == "type-1"
-    @test model.ext[:miplearn][:variable_categories][x[2]] == "type-2"
-    @test model.ext[:miplearn][:variable_categories][x[3]] == "type-3"
-    @test model.ext[:miplearn][:constraint_features][c1] == [1.0, 2.0, 3.0]
-    @test model.ext[:miplearn][:constraint_categories][c1] == "c1"
-    @test model.ext[:miplearn][:instance_features] == [5.0]
 
-    solver = LearningSolver(Gurobi.Optimizer)
-    instance = JuMPInstance(model)
-    stats = solve!(solver, instance)
-    @test stats["mip_lower_bound"] == 11.0
-    @test length(instance.py.samples) == 1
-    fit!(solver, [instance])
-    solve!(solver, instance)
+    @testset "plain model" begin
+        model = Model()
+        @variable(model, x, Bin)
+        @variable(model, y, Bin)
+        @objective(model, Max, x + y)
+        solver = LearningSolver(Gurobi.Optimizer)
+        instance = JuMPInstance(model)
+        stats = solve!(solver, instance)
+    end
 end