diff --git a/Project.toml b/Project.toml
index 5caf2f0..ae3cc07 100644
--- a/Project.toml
+++ b/Project.toml
@@ -9,6 +9,7 @@ DataFrames = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0"
 DocStringExtensions = "ffbed154-4ef7-542d-bbb7-c09d3a79fcae"
 JSON = "682c06a0-de6a-54ab-a142-c8b1cf79cde6"
 Logging = "56ddb016-857b-54e1-b83d-db4d58db5568"
+Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
 Tables = "bd369af6-aec1-5ad0-b16a-f7cc5008161c"
 Unicode = "4ec0a83e-493e-50e2-b9ac-8f72acf5a8f5"
 
diff --git a/src/GardenOptim.jl b/src/GardenOptim.jl
index 7773f20..800f3d8 100644
--- a/src/GardenOptim.jl
+++ b/src/GardenOptim.jl
@@ -6,7 +6,7 @@ using DocStringExtensions
 using Tables
 
 export loadplants, loadgarden, loadclassification, loadcosts
-export update!, randomgardenevolution!, outputgarden
+export update!, fillgardenrandomly!, gardenevolution!, outputgarden
 
 @template (FUNCTIONS, METHODS, MACROS) =
     """
diff --git a/src/mcmc.jl b/src/mcmc.jl
index f235528..efa9c91 100644
--- a/src/mcmc.jl
+++ b/src/mcmc.jl
@@ -1,4 +1,5 @@
 using Logging
+using Random
 
 "Return a random index to be filled from the garden mask."
 function randomindex(mask::Matrix{Bool})::Int
@@ -66,21 +67,18 @@ function update!(
     garden
 end
 
-"Fill the garden randomly with a predefined number of plants."
-function randomfillgarden!(garden::Matrix{Int}, mask::Matrix{Bool}, plantcount::Int)
-    garden[mask] = rand(1:plantcount, sum(mask))
+"Fill the garden randomly with a predefined number for each plant."
+function fillgardenrandomly!(garden::Matrix{Int}, mask::Matrix{Bool}, plants::DataFrame)
+    cells = vcat([repeat([plant], count) for (plant, count) in eachrow(plants)]...)
+    # fill the remaining slots with random plants
+    diffcount = sum(mask) - length(cells)
+    cells = vcat(cells, rand(cells, diffcount))
+    garden[mask] = shuffle!(indexin(cells, plants.name))
     garden
 end
 
-"Update the garden for a given number of steps, starting from a random initialisation."
-function randomgardenevolution!(
-    garden::Matrix{Int},
-    mask::Matrix{Bool},
-    costs::Matrix{Float64};
-    steps::Int = 10000
-)
-    m = size(costs, 1)
-    garden = randomfillgarden!(garden, mask, m)
+"Update the garden for a given number of steps."
+function gardenevolution!(garden::Matrix{Int}, mask::Matrix{Bool}, costs::Matrix{Float64}; steps::Int = 10000)
     for i = 1:steps
         update!(garden, mask, costs, 10.0)
     end