### Persefone.jl - a model of agricultural landscapes and ecosystems in Europe.
###
### This file manages the landscape maps that underlie the model.
###

using Printf

## IMPORTANT: do not change the order of this enum, or initlandscape() will break!
"The land cover classes encoded in the Mundialis Sentinel data."
@enum LandCover nodata forest grass water builtup soil agriculture

"The soil type of a Pixel or FarmPlot"
@enum SoilType begin
    nosoildata = 0   # 0
    sand             # 1
    loamy_sand       # 2
    sandy_loam       # 3
    loam             # 4
    silt_loam        # 5
    silt             # 6
    sandy_clay_loam  # 7
    clay_loam        # 8
    silty_clay_loam  # 9
    sandy_clay       # 10
    silty_clay       # 11
    clay             # 12
    nosoil           # 13
end

# TODO: move this into a CSV file and read it from there
"""Bodenatlas soil type id, corresponding Persefone soil type, and numbers to Persefone SoilType enum and the
   original Bodenatlas description of the soil type"""
const df_soiltypes_bodenatlas = DataFrame(
    NamedTuple{(:bodenatlas_id, :bodenatlas_name, :persefone_soiltype), Tuple{Int, String, SoilType}}.([
        ( 1, "Abbauflächen",      nosoil),
        ( 2, "Gewässer",          nosoil),
        ( 3, "Lehmsande (ls)",    loamy_sand),
        ( 4, "Lehmschluffe (lu)", silt_loam),
        ( 5, "Moore",             nosoil),
        ( 6, "Normallehme (ll)",  loam),
        ( 7, "Reinsande (ss)",    sand),
        ( 8, "Sandlehme (sl)",    sandy_loam),
        ( 9, "Schluffsande (us)", sandy_loam),
        (10, "Schlufftone (ut)",  silty_clay),
        (11, "Siedlung",          nosoil),
        (12, "Tonlehme (tl)",     clay_loam),
        (13, "Tonschluffe (tu)",  silty_clay_loam),
        (14, "Watt",              nosoil),
    ])
)

"Map a Bodenatlas soil type integer to a Persefone `SoilType` enum"
const soiltype_bodenatlas_to_persefone =
    Dict{Int,SoilType}(row.bodenatlas_id => row.persefone_soiltype for row in eachrow(df_soiltypes_bodenatlas))


"The types of management event that can be simulated"
@enum Management tillage sowing fertiliser pesticide harvesting


"""
    Pixel

A pixel is a simple data structure to combine land use and ownership information
in a single object. The model landscape consists of a matrix of pixels.
(Note: further landscape information may be added here in future.)
"""
mutable struct Pixel
    landcover::LandCover          # land cover class at this position
    soiltype::SoilType            # soil type class at this position
    fieldid::Union{Missing,Int64} # ID of the farmplot (if any) at this position
    events::Vector{Management}    # management events that have been applied to this pixel
    animals::Vector{Int64}        # IDs of animals currently at this position
    territories::Vector{String}   # IDs of animals that claim this pixel as part of their territory
end

Pixel(landcover::LandCover, soiltype::SoilType, fieldid::Union{Missing,Int64}) =
    Pixel(landcover, soiltype, fieldid, Vector{Management}(), Vector{Int64}(), Vector{String}())
Pixel(landcover::LandCover, fieldid::Union{Missing,Int64}) =
    Pixel(landcover, nosoildata, fieldid)
Pixel(landcover::LandCover) = Pixel(landcover, missing)

showlandscape(mat::T) where T <: AbstractMatrix{Pixel} = showlandscape(stdout, mat)

function showlandscape(io::IO, mat::T) where T <: AbstractMatrix{Pixel}
    println(io, "Matrix{Pixel}:")

    println(io, "  LandCover:")
    nrow, ncol = size(mat)
    max_fieldid = maximum(skipmissing(map(x -> getfield(x, :fieldid), mat)); init=0)
    for i in axes(mat, 1)
        print(io, "    ")
        for j in axes(mat, 2)
            charid = uppercase(first(string(mat[i, j].landcover)))
            fieldid = if ismissing(mat[i, j].fieldid)
                repeat(" ", ndigits(max_fieldid))
            else
                @sprintf("%*s", ndigits(max_fieldid), mat[i, j].fieldid)
            end
            print(io, charid, fieldid, " ")
        end
        println(io)
    end
    # print legend
    legend = join(map(x -> "$(uppercase(first(string(x)))) = $x", instances(Persefone.LandCover)), ", ")
    println(io, "Legend: ", legend)

    soiltype_to_str(s::SoilType) = replace(string(s), r"^soiltype_" => "")
    soiltype_to_char(s::SoilType) = uppercase(first(soiltype_to_str(s)))
    println(io)
    println(io, "  Soil types:")
    for i in axes(mat, 1)
        print(io, "    ")
        for j in axes(mat, 2)
            charid = soiltype_to_char(mat[i, j].soiltype)
            fieldid = if ismissing(mat[i, j].fieldid)
                repeat(" ", ndigits(max_fieldid))
            else
                @sprintf("%*s", ndigits(max_fieldid), mat[i, j].fieldid)
            end
            print(io, charid, fieldid, " ")
        end
        println(io)
    end
    # print legend
    legend = join(map(x -> "$(soiltype_to_char(x)) = $(soiltype_to_str(x))", instances(Persefone.SoilType)), ", ")
    println(io, "Legend: ", legend)

    # print number of unique animals, events, territories
    println(io)
    nanimals = length(unique(reduce(vcat, map(x -> getfield(x, :animals), mat))))
    nevents = length(unique(reduce(vcat, map(x -> getfield(x, :events), mat))))
    nterritories = length(unique(reduce(vcat, map(x -> getfield(x, :territories), mat))))
    println(io, "Num. unique animals    : ", nanimals)
    println(io, "Num. unique events     : ", nevents)
    print(io, "Num. unique territories: ", nterritories)  # no newline on last line of output
    return nothing
end


"""
    FarmEvent

A data structure to define a landscape event, giving its type,
spatial extent, and duration.
"""
mutable struct FarmEvent
    management::Management
    pixels::Vector{Tuple{Int64,Int64}}
    duration::Int64
end

"""
    initlandscape(directory, landcovermap, farmfieldsmap, soiltypesmap)

Initialise the model landscape based on the map files specified in the
configuration. Returns a matrix of pixels.
"""
function initlandscape(directory::String, landcovermap::String, farmfieldsmap::String, soiltypesmap::String)
    @debug "Initialising landscape"
    landcovermap = joinpath(directory, landcovermap)
    farmfieldsmap = joinpath(directory, farmfieldsmap)
    soiltypesmap = joinpath(directory, soiltypesmap)
    !(isfile(landcovermap)) && error("Landcover map $(landcovermap) doesn't exist.")
    !(isfile(farmfieldsmap)) && error("Field map $(farmfieldsmap) doesn't exist.")
    !(isfile(soiltypesmap)) && error("Soil type map $(soiltypesmap) doesn't exist.")
    # read in TIFFs, replacing missing values with 0
    landcover = coalesce(GeoArrays.read(landcovermap), 0)
    farmfields = coalesce(GeoArrays.read(farmfieldsmap), 0)
    soiltypes = coalesce(GeoArrays.read(soiltypesmap), 0)
    if size(landcover) != size(farmfields) || size(farmfields) != size(soiltypes)
        error("Input map sizes don't match.")
    end
    width, height = size(landcover)[1:2]
    landscape = Matrix{Pixel}(undef, width, height)
    for x in 1:width
        for y in 1:height
            # convert the numeric landcover value to LandCover, then create the Pixel object
            lcv = LandCover(Int(landcover[x,y][1]/10))
            ff = Int64(farmfields[x,y][1])
            st = soiltype_bodenatlas_to_persefone[Int(soiltypes[x, y][1])]
            (iszero(ff)) && (ff = missing)
            landscape[x,y] = Pixel(lcv, st, ff)
        end
    end
    return landscape
end

"""
    updateevents!(model)

Cycle through the list of events, removing those that have expired.
"""
function updateevents!(model::SimulationModel)
    expiredevents = []
    for e in 1:length(model.events)
        event = model.events[e]
        event.duration -= 1
        if event.duration <= 0
            push!(expiredevents, e)
            for p in event.pixels
                i = findnext(x -> x == event.management, model.landscape[p...].events, 1)
                deleteat!(model.landscape[p...].events, i)
            end
        end
    end
    deleteat!(model.events, expiredevents)
end

"""
    createevent!(model, pixels, name, duration=1)

Add a farm event to the specified pixels (a vector of position tuples) for a given duration.
"""
function createevent!(model::SimulationModel, pixels::Vector{Tuple{Int64,Int64}},
                     name::Management, duration::Int64=1)
    push!(model.events, FarmEvent(name, pixels, duration))
    for p in pixels
        push!(model.landscape[p...].events, name)
    end
end

"""
    landcover(position, model)

Return the land cover class at this position (utility wrapper).
"""
function landcover(pos::Tuple{Int64,Int64}, model::SimulationModel)
    model.landscape[pos...].landcover
end

"""
    farmplot(position, model)

Return the farm plot at this position, or nothing if there is none (utility wrapper).
"""
function farmplot(pos::Tuple{Int64,Int64}, model::SimulationModel)
    id = model.landscape[pos...].fieldid
    ismissing(id) ? nothing : model.farmplots[id]
end

"""
    directionto(pos, model, habitatdescriptor)

Calculate the direction from the given location to the closest location matching the
habitat descriptor function. Returns a coordinate tuple (target - position), or nothing
if no matching habitat is found. Caution: can be computationally expensive!
"""
function directionto(pos::Tuple{Int64,Int64}, model::SimulationModel, habitatdescriptor::Function)
    #XXX split this up into a findnearest() and a directionto() function?
    # search in expanding rectangles around the origin
    (habitatdescriptor(pos, model)) && (return (0,0))
    dist = 1
    width, height = size(model.landscape)
    while dist <= width || dist <= height #XXX is there a quicker method?
        # check the upper and lower bounds of the enclosing rectangle
        y1 = pos[2] - dist
        y2 = pos[2] + dist
        minx = maximum((pos[1]-dist, 1))
        maxx = minimum((pos[1]+dist, width))
        #XXX alternately, we could check individually whether y1 and y2 are in bounds, and loop
        # over x for each separately (not sure if that is quicker, but it may be)
        for x in minx:maxx
            (y1 > 0 && habitatdescriptor((x,y1), model)) && (return (x, y1).-pos)
            (y2 <= height && habitatdescriptor((x,y2), model)) && (return (x, y2).-pos)
        end
        # check the left and right bounds of the enclosing rectangle
        x1 = pos[1] - dist
        x2 = pos[1] + dist
        miny = maximum((pos[2]-dist+1, 1))
        maxy = minimum((pos[2]+dist-1, height))
        for y in miny:maxy
            (x1 > 0 && habitatdescriptor((x1,y), model)) && (return (x1,y).-pos)
            (x2 <= width && habitatdescriptor((x2,y), model)) && (return (x2,y).-pos)
        end
        dist += 1
    end
    return nothing
end

"""
    directionto(pos, model, habitattype)

Calculate the direction from the given location to the closest habitat of the specified type.
Returns a coordinate tuple (target - position), or nothing if no matching habitat is found.
Caution: can be computationally expensive!
"""
function directionto(pos::Tuple{Int64,Int64}, model::SimulationModel, habitattype::LandCover)
    # can't use @habitat here because macros.jl is loaded later than this file
    directionto(pos, model, function(p,m) landcover(p,m) == habitattype end)
end

"""
    distanceto(pos, model, habitatdescriptor)

Calculate the distance from the given location to the closest location matching the
habitat descriptor function. Caution: can be computationally expensive!
"""
function distanceto(pos::Tuple{Int64,Int64}, model::SimulationModel, habitatdescriptor::Function)
    #XXX this is very imprecise, as diagonal distances are not calculated trigonometrically
    target = directionto(pos, model, habitatdescriptor)
    isnothing(target) && return Inf
    return maximum(abs.(target))*@param(world.mapresolution)
end

"""
    distanceto(pos, model, habitattype)

Calculate the distance from the given location to the closest habitat of the specified type.
Caution: can be computationally expensive!
"""
function distanceto(pos::Tuple{Int64,Int64}, model::SimulationModel, habitattype::LandCover)
    # can't use @habitat here because macros.jl is loaded later than this file
    distanceto(pos, model, function(p,m) landcover(p,m) == habitattype end)
end

"""
    distancetoedge(pos, model)

Calculate the distance from the given location to the closest neighbouring habitat.
Caution: can be computationally expensive!
"""
function distancetoedge(pos::Tuple{Int64,Int64}, model::SimulationModel)
    # can't use @habitat here because macros.jl is loaded later than this file
    distanceto(pos, model, function(p,m) landcover(p,m) != landcover(pos, model) end)
end

"""
    randompixel(position, model, range, habitatdescriptor)

Find a random pixel within a given `range` of the `position` that matches the
habitatdescriptor (create this using [`@habitat`](@ref)).
"""
function randompixel(pos::Tuple{Int64,Int64}, model::SimulationModel, range::Length,
                     habitatdescriptor::Function=(pos,model)->nothing)
    range = Int(floor(range / @param(world.mapresolution)))
    for x in @shuffle!(collect((pos[1]-range):(pos[1]+range)))
        for y in @shuffle!(collect((pos[2]-range):(pos[2]+range)))
            !inbounds((x,y), model) && continue
            habitatdescriptor((x,y), model) && return (x,y)
        end
    end
    nothing
end

"""
    randomdirection(model, distance)

Get a random direction coordinate tuple within the specified distance.
"""
function randomdirection(model::SimulationModel, distance::Length)
    range = Int(floor(distance / @param(world.mapresolution)))
    Tuple(@rand(-range:range, 2))
end

"""
    inbounds(pos, model)

Is the given position within the bounds of the model landscape?
"""
function inbounds(pos, model)
    dims = size(model.landscape)
    pos[1] > 0 && pos[1] <= dims[1] && pos[2] > 0 && pos[2] <= dims[2]    
end

"""
    safebounds(pos, model)

Make sure that a given position is within the bounds of the model landscape.
"""
function safebounds(pos::Tuple{Int64,Int64}, model::SimulationModel)
    width, height = size(model.landscape)
    (bounds(pos[1], max=width, min=1), bounds(pos[2], max=height, min=1))
end