JuliaLang @ SSC2018

forestfire.jl

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)

sheep.jl

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