### 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

# TODO: check names and order of enum
# TODO: note where values come from
## IMPORTANT: do not change the order of this enum, or initlandscape() will break!
"The soil type"
@enum SoilType begin
    soiltype_nodata = 0       # 0
    soiltype_sand             # 1
    soiltype_loamy_sand       # 2
    soiltype_sandy_loam       # 3
    soiltype_loam             # 4
    soiltype_silt_loam        # 5
    soiltype_silt             # 6
    soiltype_sandy_clay_loam  # 7
    soiltype_clay_loam        # 8
    soiltype_silty_clay_loam  # 9
    soiltype_sandy_clay       # 10
    soiltype_silty_clay       # 11
    soiltype_clay             # 12
    soiltype_unknown          # 13  TODO: invented this to get a type 13
end

"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, soiltype_nodata, 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])
            (iszero(ff)) && (ff = missing)
            landscape[x,y] = Pixel(lcv, 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