Skip to content

Instantly share code, notes, and snippets.

@bkamins
Last active August 20, 2018 14:54
Show Gist options
  • Save bkamins/243348d888f6751005a67feb33a54edc to your computer and use it in GitHub Desktop.
Save bkamins/243348d888f6751005a67feb33a54edc to your computer and use it in GitHub Desktop.
JuliaLang @ SSC2018
using Random, Statistics, DataFrames
try
using PyPlot
global const HAS_PLOTS = true
catch
global const HAS_PLOTS = false
end
@enum Tree None Green Red Brown
function setup(density, sizex=251, sizey=251)
[x == 1 ? Red : rand() < density ? Green : None for x in 1:sizex, y in 1:sizey]
end
function go(grid)
if any(isequal(Red), grid)
redtrees = findall(isequal(Red), grid)
for pos in shuffle!(redtrees)
x, y = pos[1], pos[2]
for (dx, dy) in ((0, 1), (0, -1), (1, 0), (-1, 0))
nx, ny = x + dx, y + dy
if 1 ≤ nx ≤ size(grid, 1) && 1 ≤ ny ≤ size(grid, 2)
if grid[nx, ny] == Green
grid[nx, ny] = Red
end
end
end
grid[pos] = Brown
end
return false
end
return true
end
function grid2mat(grid)
mat = zeros(Int, size(grid, 1), size(grid, 2), 3)
for i in 1:size(grid, 1), j in 1:size(grid, 2)
if grid[i,j] == Green
mat[i,j,2] = 255 # it will be green
elseif grid[i,j] == Red
mat[i,j,3] = 255 # it will be blue
elseif grid[i,j] == Brown
mat[i,j,1] = 100 # it will be brown
end
end
mat
end
function runsim(density, doplot=false)
grid = setup(density)
init_green = count(==(Green), grid)
if doplot && HAS_PLOTS
img = imshow(grid2mat(grid))
title("Step 0")
end
t = 0
while true
res = go(grid)
t += 1
if doplot && HAS_PLOTS
img[:set_data](grid2mat(grid))
title("Step $t")
show()
sleep(0.0001)
end
if res
return count(==(Brown), grid[2:end, :]) / init_green * 100
end
end
end
runsim(0.65, true)
@time runsim(0.55)
function collectsim()
Random.seed!(1)
df = DataFrame(rep=Int[], density=Float64[], pct=Float64[])
for density in 0.4:0.01:0.8
println("density: ", density)
for rep in 1:100
push!(df, (rep, density, runsim(density)))
end
end
by(df, :density, x -> DataFrame(meanpct = mean(x.pct),
sdpct = sqrt(var(x.pct)/nrow(x))))
end
res = collectsim()
plt[:fill_between](res.density, res.meanpct-2res.sdpct, res.meanpct+2res.sdpct)
plot(res.density, res.meanpct, color=:red)
using Random, Statistics, Distributions, DataFrames
try
using PyPlot
global const HAS_PLOTS = true
catch
global const HAS_PLOTS = false
end
mutable struct Sheep
angle::Float64
x::Float64
y::Float64
energy::Float64
end
function move(s::Sheep)
s.angle += rand(Uniform(0,π/2)) + rand(Uniform(-π/2, 0))
s.x += cos(s.angle)
s.y += sin(s.angle)
end
function patch(s::Sheep, grid)
floor(Int, mod(s.y, size(grid, 1))) + 1,
floor(Int, mod(s.x, size(grid, 2))) + 1
end
function tick(s::Sheep, grid::Matrix{Int})
move(s)
s.energy -= 15 # parameter; in NetLogo we have 10
x, y = patch(s, grid)
if grid[x,y] ≥ 20 # parameter
s.energy += 20
grid[x,y] -= 20
end
end
function tick(sheep_set::Set{Sheep})
for s in collect(sheep_set)
if s.energy > 2000
s.energy -= 1000 # parameter
push!(sheep_set, Sheep(s.angle, s.x, s.y, 1000))
end
if s.energy ≤ 0
pop!(sheep_set, s)
end
end
end
function tick(grid::Matrix{Int})
grid .= min.(grid .+ 10, 100) # parameter
end
function setup(nsheep=700, sizex=35, sizey=35) # parameter
rand(0:100, sizex, sizey),
Set([Sheep(2π*rand(), sizey*rand(), sizex*rand(),
1000) for i in 1:nsheep])
end
function stats(grid, sheep_set)
sheepcount = zeros(Int, size(grid));
for s in sheep_set
(x, y) = patch(s, grid)
sheepcount[x, y] += 1
end
sum(grid), length(sheep_set), cor(vec(sheepcount), vec(grid))
end
function tick(grid, sheep_set)
for s in sheep_set
tick(s, grid)
end
tick(sheep_set)
tick(grid)
stats(grid, sheep_set)
end
function runsim()
grid, sheep_set = setup()
gs, ss, cs = stats(grid, sheep_set)
df = DataFrame(gs=gs, ss=ss, cs=cs)
for i in 1:2000 # parameter; change to 10000 for performance testing
push!(df, tick(grid, sheep_set))
end
df, grid, sheep_set
end
@time runsim();
@time res, grid, sheep_set = runsim();
if HAS_PLOTS
xpos = mod.([s.x for s in sheep_set], size(grid, 2))
ypos = mod.([s.y for s in sheep_set], size(grid, 1))
plt[:subplot](1, 3, 1)
plot(res.ss)
title("Sheep count")
plt[:subplot](1, 3, 2)
plot(res.cs)
title("Sheep-grass correlation")
plt[:subplot](1, 3, 3)
imshow(grid)
scatter(xpos .- 0.5, ypos .- 0.5, color=:red)
title("Landscape overview")
end
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment