Update SpatialRust and run scripts

Project.toml

@@ -1,19 +1,18 @@
 name = "SpatialRust"
 uuid = "8b437046-e5dd-42a9-af37-0c38301ea563"
 authors = ["Manuela Vanegas Ferro"]
-version = "0.4.0"
+version = "1.0.0"
 
 [deps]
 Agents = "46ada45e-f475-11e8-01d0-f70cc89e6671"
-BenchmarkTools = "6e4b80f9-dd63-53aa-95a3-0cdb28fa8baf"
 CSV = "336ed68f-0bac-5ca0-87d4-7b16caf5d00b"
 DataFrames = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0"
 Distributed = "8ba89e20-285c-5b6f-9357-94700520ee1b"
 Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f"
 DrWatson = "634d3b9d-ee7a-5ddf-bec9-22491ea816e1"
 Profile = "9abbd945-dff8-562f-b5e8-e1ebf5ef1b79"
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
-StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
+Revise = "295af30f-e4ad-537b-8983-00126c2a3abe"
 Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
 StatsBase = "2913bbd2-ae8a-5f71-8c99-4fb6c76f3a91"
 

scripts/ParameterRuns.jl

@@ -1,44 +1,31 @@
 # Load Packages
 import Pkg
 Pkg.activate(".") 
-using Agents, CSV, DataFrames, Distributed, DrWatson, Statistics
-using SpatialRust
+using CSV, DataFrames, DrWatson, Statistics, SpatialRust
 
-# Define parameter options
-mean_temp = collect(20.0:0.5:25.0) # mean temperature values: min:step:max
-rain_prob = collect(0.4:0.1:0.9) # rain probability values
-wind_prob = collect(0.1:0.1:0.9) # wind probability values
+# Define parameter options and # of steps
+mean_temp = collect(19.0:0.5:24.0) # mean temperature values: min:step:max
+rain_prob = collect(0.3:0.2:0.9) # rain probability values
+wind_prob = collect(0.3:0.2:0.9) # wind probability values
 reps = 5 # number of replicates per parameter combination
+years = 2
+
+# Define which parameter columns are included in output dataframe and their order
+parsorder = [:rep, :wind_prob, :rain_prob, :mean_temp]
 
 # Create dictionary with parameter conditions
 conds = Dict(
     :mean_temp => mean_temp,
     :rain_prob => rain_prob,
     :wind_prob => wind_prob,
-    :reps => collect(1:reps),
-
-    # other settings
-    :shade_d => [6],
-    :barrier_arr => (1,1,0,0),
-    :target_shade => 0.5,
-    :prune_period => 365,
-    :fungicide_period => 365,
-    :barrier_rows => 2,
-    :shade_g_rate => 0.05,
-    :steps => 1095,
-    :rust_gr => 1.63738,
-    :cof_gr => 0.393961,
-    :spore_pct => 0.821479,
-    :fruit_load => 0.597133,
-    :uv_inact => 0.166768,
-    :rain_washoff => 0.23116,
-    :rain_distance => 0.80621,
-    :wind_distance => 3.29275,
-    :exhaustion => 0.17458
+    :rep => collect(1:reps),
+    :steps => years * 365,
 )
 
-# Run simulations (in parallel)
-results = parameters_experiment(conds)
+# Generate dataframe with all the parameter permutations
+pars = DataFrame(dict_list(conds))
+# Run simulations
+results = parameters_experiment(pars, parsorder)
 
 # Write results
 CSV.write("/srv/results/parameterexp.csv", results)

src/ABM/CGrowerSteps.jl

@@ -1,60 +1,83 @@
-function harvest!(model::ABM)
-    harvest = 0.0
-    # ids = model.current.coffee_ids
-    # for c in (model[id] for id in ids)
-    for c in Iterators.filter(c -> c isa Coffee, allagents(model))
-        harvest += c.production / model.pars.harvest_cycle
-        c.production = 1.0
-        # if plant.fung_this_cycle
-        #     plant.fung_this_cycle = false
-        #     plant.productivity = plant.productivity / 0.8
-        # end
-        # if plant.pruned_this_cycle
-        #     plant.pruned_this_cycle = false
-        #     plant.productivity = plant.productivity / 0.9
-        # end
+function harvest!(model::SpatialRustABM)
+    yprod = sum(map(c -> c.production, model.agents))
+    model.current.prod += yprod
+    cost = model.current.costs += model.mngpars.fixed_costs +
+        yprod * (model.mngpars.other_costs * (1.0 -  (model.current.shadeacc / 365.0) * mean(model.shade_map)) + 0.012)
+    model.current.shadeacc = 0.0
+
+    model.current.fung_count = 0
+    map(a -> new_harvest_cycle!(a, model.mngpars.lesion_survive), model.agents)
+    return nothing
+end
+
+function new_harvest_cycle!(c::Coffee, surv_p::Float64)
+    c.production = 0.0
+    c.deposited *= surv_p
+    surv_n = trunc(Int, c.n_lesions * surv_p)
+    if surv_n == 0
+        c.n_lesions = 0
+        empty!(c.ages)
+        empty!(c.areas)
+        empty!(c.spores)
+        if c.deposited < 0.05
+            c.deposited = 0.0
+            c.rusted = false
+        end
+    else
+        lost = c.n_lesions - surv_n
+        c.n_lesions = surv_n
+        deleteat!(c.ages, 1:lost)
+        deleteat!(c.areas, 1:lost)
+        deleteat!(c.spores, 1:lost)
     end
-    # model.current.net_rev += (model.pars.coffee_price * harvest) - model.current.costs
-    # model.current.gains += model.coffee_price * harvest * model.pars.p_density
-    model.current.prod += harvest
+    return nothing
 end
 
-function fungicide!(model::ABM)
-    model.current.costs += length(model.current.coffee_ids) * 1.0 #model.pars.fung_cost
-    model.current.fung_effect = 15
+function prune_shades!(model::SpatialRustABM, tshade::Float64)
+    if model.current.ind_shade > tshade
+        model.current.ind_shade = tshade
+    else
+        model.current.ind_shade *= 0.9
+    end
+    model.current.costs += model.mngpars.tot_prune_cost
 end
 
-# function prune!(model::ABM)
-#     # n_pruned = trunc(model.pars.prune_effort * length(model.current.shade_ids))
-#     # model.current.costs += n_pruned * model.pars.prune_cost
-#     model.current.costs += length(model.current.shade_ids) * model.pars.prune_cost
-#     # pruned = partialsort(model.current.shade_ids, 1:n_pruned, rev=true, by = x -> model[x].shade)
-#     # for pr in pruned
-#     for pr in model.current.shade_ids
-#         model[pr].shade = model.pars.target_shade
-#     end
-# end
-
-function inspect!(model::ABM)
-    n_inspected = round(Int,(model.pars.inspect_effort * length(model.current.coffee_ids)), RoundToZero)
-    # n_inspected = model.pars.n_inspected
-    cofids = sample(model.rng, model.current.coffee_ids, n_inspected, replace = false)
-
-    for c in cofids
-        # cof = model[c]
-        if model[c].hg_id != 0# && rand < model.pars.inspect_effort * (sum(model[hg_id].state[2,]) / 3)
-            #elimina las que sean > 2.5, * effort
-            rust = model[model[c].hg_id]
-            @inbounds areas = rust.state[2, 1:rust.n_lesions]
-            if any(areas .> 0.05)
-                replace!(a -> a .< 0.05 ? 0.0 : a, areas)
-                spotted = sample(model.rng, 1:rust.n_lesions, weights(areas), 5)
-                newstate = rust.state[:, Not(spotted)]
-                rust.n_lesions = size(newstate)[2]
-                rust.state = hcat(newstate, ones(4, 25 - rust.n_lesions))
+function inspect!(model::SpatialRustABM)
+    n_infected = 0
+    actv = filter(active, model.agents)
+    if model.mngpars.n_inspected < length(actv)
+        inspected = sample(model.rng, actv, model.mngpars.n_inspected, replace = false)
+    else
+        inspected = actv
+    end
+
+    rmles = model.mngpars.rm_lesions
+    for c in inspected
+        # lesion area of 0.05 means a diameter of ~0.25 cm, which is taken as minimum so grower can see it
+        nvis = sum(>(0.05), c.areas, init = 0.0)
+        # area of 0.8 means a diameter of ~1 cm
+        if nvis > 0  && (0.8 < maximum(c.areas, init = 0.0) || rand(model.rng) < nvis / 5)
+            n_infected += 1
+            spotted = unique!(sort!(sample(model.rng, 1:c.n_lesions, weights(visible.(c.areas)), rmles)))
+            deleteat!(c.ages, spotted)
+            deleteat!(c.areas, spotted)
+            deleteat!(c.spores, spotted)
+            c.n_lesions -= length(spotted)
+            if c.n_lesions == 0 && (c.deposited < 0.05)
+                c.deposited == 0.0
+                c.rusted = false
             end
-            # rust.n_lesions = round(Int, model[cof.hg_id].n_lesions * 0.1)
-            # rust.area = round(Int, model[cof.hg_id].area * 0.1)
         end
     end
+
+    model.current.costs += model.mngpars.tot_inspect_cost
+    model.current.obs_incidence = n_infected / model.mngpars.n_inspected
 end
+
+visible(a::Float64) = a > 0.05 ? a : 0.0
+
+function fungicide!(model::SpatialRustABM)
+    model.current.costs += model.mngpars.tot_fung_cost
+    model.current.fungicide = 1 # model.mngpars.fung_effect
+    model.current.fung_count += 1
+end

src/ABM/CoffeeSteps.jl

@@ -0,0 +1,65 @@
+
+function update_sunlight!(cof::Coffee, map::Matrix{Float64}, ind_shade::Float64)
+    cof.sunlight = 1.0 - @inbounds map[cof.pos...] * ind_shade
+end
+
+function veg_growth!(coffee::Coffee, pars::CoffeePars, sl::Float64)
+    photo_veg = coffee.veg * pars.photo_frac
+    PhS = pars.photo_const * (sl / (pars.k_sl + sl)) * (photo_veg / (pars.k_v + photo_veg))
+
+    coffee.veg += pars.phs_veg * PhS - pars.μ_veg * coffee.veg
+    if coffee.veg < 0.0
+        coffee.veg = 0.0
+    end
+    coffee.storage += pars.phs_sto * PhS 
+end
+
+function rep_growth!(coffee::Coffee, pars::CoffeePars, sl::Float64)
+    veg = coffee.veg
+    photo_veg = veg * pars.photo_frac
+    μ_v = pars.μ_veg * veg
+    prod = coffee.production
+
+    PhS = pars.photo_const * (sl / (pars.k_sl + sl)) *
+    (photo_veg / (pars.k_v + photo_veg))
+
+    if coffee.storage < 0.0
+        Δprod = PhS - pars.μ_prod * prod
+        if Δprod > 0.0
+            last = Δprod - μ_v
+            coffee.veg += min(last, 0.0)
+            coffee.storage += max(last, 0.0)
+        else
+            coffee.veg -= μ_v
+            coffee.production += Δprod
+        end
+    else
+        frac_v = veg / (veg + prod)
+        d_prod = (1.0 - frac_v) * PhS - pars.μ_prod * prod
+        if d_prod < 0.0
+            coffee.veg += pars.phs_veg * frac_v * PhS - μ_v
+            coffee.storage += 0.95 * d_prod
+            coffee.production += 0.05 * d_prod
+        else
+            coffee.veg += pars.phs_veg * (d_prod + frac_v * PhS) - μ_v
+            # coffee.veg += d_prod + frac_v * pars.phs_veg * PhS - μ_v
+        end
+    end
+
+    if coffee.veg < 0.0
+        coffee.veg = 0.0
+    end
+    if coffee.production < 0.0
+        coffee.production = 0.0
+    end
+end
+
+function regrow!(cof::Coffee, sl::Float64, map::Matrix{Int})
+    cof.veg = 2.0
+    cof.storage = init_storage(sl)
+    cof.exh_countdown = 0
+    @inbounds map[cof.pos...] = 1
+end
+
+
+init_storage(sl::Float64) = 75.5 * exp(-5.5 * sl) + 2.2