diff --git a/src/analysis/sampling/hit_and_run.jl b/src/analysis/sampling/hit_and_run.jl
index fa2b82c910e81ffec2536f00ef081fb39747dfdb..03609850d1934266ed403bed32192272f5ae8a2c 100644
--- a/src/analysis/sampling/hit_and_run.jl
+++ b/src/analysis/sampling/hit_and_run.jl
@@ -47,17 +47,25 @@ samples = hit_and_run(100_000, opt_model; keepevery=10, samplesize=5000)
 ```
 """
 function hit_and_run(
-    N::Int,
-    opt_model;
+    model,
+    optimizer;
+    modifications = [],
+    N = 1000,
     keepevery = _constants.sampling_keep_iters,
     samplesize = _constants.sampling_size,
-    random_objective = false,
+    num_warmup_points = length(reactions(model))
 )
 
-    lbs, ubs = get_bound_vectors(opt_model) # get actual ub and lb constraints, can't use model function because the user may have changed them in the function arguments
-
-    wpoints = _get_warmup_points(opt_model; random_objective = random_objective)
+    # generate warmup points
+    ws = flux_variability_analysis(
+        model,
+        collect(1:10),
+        Tulip.Optimizer;
+        bounds = objective_bounds(0.90),
+        ret = m -> COBREXA.JuMP.value.(m[:x])
+    )
 
+    
     nwpts = size(wpoints, 2) # number of warmup points generated
     samples = zeros(size(wpoints, 1), samplesize) # sample storage
     current_point = zeros(size(wpoints, 1))
diff --git a/src/analysis/sampling/warmup.jl b/src/analysis/sampling/warmup.jl
new file mode 100644
index 0000000000000000000000000000000000000000..4aa832ca3731fcd37964e911bb89142a1ed3bbbb
--- /dev/null
+++ b/src/analysis/sampling/warmup.jl
@@ -0,0 +1,38 @@
+"""
+warmup(
+    model::MetabolicModel,
+    optimizer;
+    warmup_points::Vector{Int},
+    modifications = [],
+    workers = [myid()],
+
+Generates warmup points for samplers by sequentially minimizing and maximizing 
+reactions at `warmup_points`. Very similar to [`flux_variability_analysis`](@ref).
+"""
+function warmup(
+    model::MetabolicModel,
+    optimizer;
+    warmup_points::Vector{Int},
+    modifications = [],
+    workers = [myid()],
+)
+    # create optimization problem, apply constraints, load on all workers
+    optmodel = make_optimization_model(model, optimizer)
+    for mod in modifications
+        mod(model, optmodel)
+    end
+    
+    map(fetch, save_at.(workers, :warmup_model, Ref(:($optmodel))))
+
+    ret = m -> value.(m[:x]) # get all the fluxes
+    fluxes = dpmap(
+        rid -> :($COBREXA._FVA_optimize_reaction(warmup_model, $rid, $ret)),
+        CachingPool(workers),
+        [-warmup_points warmup_points],
+    )
+
+    # free the data on workers
+    map(fetch, remove_from.(workers, :warmup_model))
+
+    return fluxes
+end