text:
"Base Entity Class from whence Ships, Dropoffs and Shipyards inherit."
abstract type Entity end

"Dropoff class for housing dropoffs."
struct Dropoff <: Entity
    owner::Int
    id::Int
    position::Position
end

"Shipyard class to house shipyards."
struct Shipyard <: Entity
    owner::Int
    id::Int
    position::Position
end

"Return a move to spawn a new ship."
spawn(shipyard::Shipyard) = Commands.GENERATE

"Ship class to house ship entities."
struct Ship <: Entity
    owner::Int
    id::Int
    position::Position
    halite_amount::Int
end

"Is this ship at max halite capacity?"
is_full(ship::Ship, max_halite::Int) = ship.halite_amount >= max_halite

"Return a move to transform this ship into a dropoff."
make_dropoff(ship::Ship) = "$(Commands.CONSTRUCT) $(ship.id)"

"Return a move to move this ship in a direction without checking for collisions."
move(ship::Ship, direction::Direction) = "$(Commands.MOVE) $(ship.id) $(convert(direction))"

"Don't move this ship."
stay_still(ship::Ship) = "$(Commands.MOVE) $(ship.id) $(Commands.STAY_STILL)"

"Player object containing all items/metadata pertinent to the player."
mutable struct Player
    id::Int
    shipyard::Shipyard
    halite_amount::Int
    ships::Dict{Int, Ship}
    dropoffs::Dict{Int, Dropoff}
end

"""
Returns a singular ship mapped by the ship id.
# Arguments
- `player::Player` : The player whom owns the ship.
- `ship_id::Int` : The ship id of the ship you wish to return.
"""
get_ship(player::Player, ship_id::Int) = player.ships[ship_id]

"Returns all ship objects for given player in an array."
get_ships(player::Player) = collect(values(player.ships))

"Return a random Ship from all ships for given player."
function random_ship(player::Player)
    ships = get_ships(player)
    if length(ships) == 0
        # TODO: unsafe if player has no ships what to return?
        return @info "Tried to get a random ship but player has no ships."
    else
        return rand(ships)
    end
end

"""
Returns a singular dropoff mapped by its id.
# Arguments
- `player::Player` : The player whom owns the dropoff.
- `dropoff_id::Int` : The dropoff id to return.
"""
get_dropoff(player::Player, dropoff_id::Int) = player.dropoffs[dropoff_id]

"Returns all dropoff objects for given player in an array."
get_dropoffs(player::Player) = collect(values(player.dropoffs))

"""
Check whether the player has a ship with a given ID. Returns true if and only if the ship exists.
Useful if you track ships via IDs elsewhere and want to make sure the ship still exists.
# Arguments
- `player::Player` : The player whom owns the ship.
- `ship_id::Int` : The ID to check.
"""
has_ship(player::Player, ship_id::Int) = haskey(player.ships, ship_id)

"Creates a player object from the input given by the game engine."
function generate_players()
    player, shipyard_x, shipyard_y = parse.(Int, split(readline()))
    Player(player, Shipyard(player, -1, Position(shipyard_x, shipyard_y)), 0, Dict(), Dict())
end

"""
Updates this player object considering the input from the game engine for the current specific turn.
# Arguments
- `player::Player` : The player to update.
- `num_ships::Int` : The number of ships this player has this turn.
- `num_dropoffs::Int` : The number of dropoffs this player has this turn.
- `halite::Int` : How much halite the player has in total.
"""
function update_player!(player::Player, num_ships::Int, num_dropoffs::Int, halite::Int)
    player.ships = Dict(update_ships!(player) for _=1:num_ships)
    player.dropoffs = Dict(update_dropoffs!(player) for _=1:num_dropoffs)
    player.halite_amount = halite
    player
end

"A cell on the game map."
mutable struct MapCell
    position::Position
    halite_amount::Int
    ship::Union{Ship, Nothing}
    structure::Union{Entity, Nothing}
end

"Return whether this cell has any ships."
is_occupied(cell::MapCell) = cell.ship != nothing

"Return whether this cell has any structures."
has_structure(cell::MapCell) = cell.structure != nothing

"Return when this cell has no ships or structures."
is_empty(cell::MapCell) = !is_occupied(cell) && !has_structure(cell)

"Return what the structure type in this cell."
structure_type(cell::MapCell) = has_structure(cell) ? typeof(cell.structure) : typeof(nothing)

"""
Mark this cell as unsafe (occupied) for navigation.
Use in conjunction with GameMap.naive_navigate.
"""
function mark_unsafe!(cell::MapCell, ship::Ship)
    cell.ship = ship
    cell
end

"""
The game map.
Can be indexed by a contained entity.
Coordinates start at 0. Coordinates are normalized for you
"""
struct GameMap
    width::Int
    height::Int
    cells::Dict{Position, MapCell}
end

"""
Returns a normalized position object fitting within the bounds of the toroidal map.
i.e.: Takes a point which may or may not be within width and height bounds, and places
it within those bounds considering wraparound.
# Arguments
- `game_map::GameMap` : The game map.
- `position::Position` : A position object.
"""
normalize(game_map::GameMap, position::Position) = Position(mod(position.x, game_map.width), mod(position.y, game_map.height))

# index game_map[x, y] or game_map[Position()] without calling game_map.cells
Base.getindex(G::GameMap, i::Vararg{Int, 2}) = G.cells[normalize(G, Position(i[1], i[2]))]
Base.getindex(G::GameMap, p::Position) = G.cells[normalize(G, p)]

"""
Compute the Manhattan distance between two locations.
Accounts for wrap-around.
# Arguments
- `game_map::GameMap` : The game map.
- `pos_1::Position` : The source from where to calculate.
- `pos_2::Position` : The target to where calculate.
"""
function calculate_distance(game_map::GameMap, pos_1::Position, pos_2::Position)
    source = normalize(game_map, pos_1)
    target = normalize(game_map, pos_2)
    resulting_position = abs(source - target)
    min(resulting_position.x, game_map.width - resulting_position.x) + min(resulting_position.y, game_map.height - resulting_position.y)
end

"""
Returns where in the cardinality spectrum the target is from source. e.g.: North, East; South, West; etc.
A tuple containing the target Direction. A tuple item (or both) could be None if within same coords
NOTE: Ignores toroid
# Arguments
- `source::Position` : The source position.
- `target::Position` : The target position.
"""
function get_target_direction(source::Position, target::Position)
    y = target.y > source.y ? Direction("South") : target.y < source.y ? Direction("North") : nothing
    x = target.x > source.x ? Direction("East") : target.x < source.x ? Direction("West") : nothing
    y, x
end

"""
Returns a list of valid (closest) Directions towards your target.
Return the Direction(s) to move closer to the target point, or empty if the points are the same.
This move mechanic does not account for collisions. The multiple directions are if both directional movements
are viable.
# Arguments
- `game_map::GameMap` : The game map.
- `source::Position` : The starting position.
- `destination::Position` : The destination towards which you wish to move your object.
"""
function get_unsafe_moves(game_map::GameMap, source::Position, destination::Position)
    source = normalize(game_map, source)
    destination = normalize(game_map, destination)
    distance = abs(destination - source)
    y_cardinality, x_cardinality = get_target_direction(source, destination)
    x = distance.x < game_map.width / 2 ? x_cardinality : invert(x_cardinality)
    y = distance.y < game_map.height / 2 ? y_cardinality : invert(y_cardinality)
    [x, y]
end

"""
Returns a singular safe move towards the destination.
# Arguments
- `game_map::GameMap` : The game map.
- `ship::Ship` : The ship to move.
- `destination::Position` : Ending position.
"""
function naive_navigate(game_map::GameMap, ship::Ship, destination::Position)
    # No need to normalize destination, since get_unsafe_moves does that.
    for direction in get_unsafe_moves(game_map, ship.position, destination)
        if direction == nothing continue end
        target_pos = normalize(game_map, directional_offset(ship.position, direction))
        if !is_occupied(game_map[target_pos])
            mark_unsafe!(game_map[target_pos], ship)
            return direction
        end
    end
    Direction("Still")
end

"Creates a map object from the input given by the game engine."
function generate_game_map()
    map_width, map_height = parse.(Int, split(readline()))
    game_map = Dict{Position, MapCell}()
    for y_position in 0:map_height-1
        cells = parse.(Int, split(readline()))
        for x_position in 0:map_width-1
            pos = Position(x_position, y_position)
            i = x_position + 1
            game_map[pos] = MapCell(pos, cells[i], nothing, nothing)
        end
    end
    GameMap(map_width, map_height, game_map)
end

"Updates the map object from the input given by the game engine."
function update_game_map!(game_map::GameMap)
    for y in 0:game_map.height-1
        for x in 0:game_map.width-1
            game_map[Position(x, y)].ship = nothing
        end
    end

    for _ in 1:parse(Int, readline())
        cell_x, cell_y, cell_energy = parse.(Int, split(readline()))
        game_map[Position(cell_x, cell_y)].halite_amount = cell_energy
    end
    game_map
end

"Updates the dropoff objects Dict from the input given by the game engine for a given player."
function update_dropoffs!(player::Player)
    dropoff_id, x_position, y_position = parse.(Int, split(readline()))
    dropoff_id => Dropoff(player.id, dropoff_id, Position(x_position, y_position))
end

"Updates the ship objects Dict from the input given by the game engine for a given player."
function update_ships!(player::Player)
    ship_id, x_position, y_position, halite = parse.(Int, split(readline()))
    ship_id => Ship(player.id, ship_id, Position(x_position, y_position), halite)
end

max_stars_repo_name: johnnykwwang/Halite-III

max_issues_repo_name: johnnykwwang/Halite-III

max_forks_repo_name: johnnykwwang/Halite-III


text:
using Test, Random, FileIO, QuartzImageIO, ColorTypes
using FixedPointNumbers, TestImages, ImageAxes
using ImageMetadata
using OffsetArrays

# Saving notes:
# autumn_leaves and toucan fail as of November 2015. The "edges" of the
# leaves are visibly different after a save+load cycle. Not sure if the
# reader or the writer is to blame. Probably an alpha channel issue.
# Mri-stack and multichannel timeseries OME are both image stacks,
# but the save code only saves the first frame at the moment.

@testset "Local" begin
    imagedir = joinpath(dirname(@__FILE__), "images")
    images = readdir(imagedir)
    @testset "$image" for image in images
        img = load(joinpath(imagedir, image))
        @test isa(img, Array)
    end
end

mydir = tempdir() * "/QuartzImages"
ispath(mydir) || mkdir(mydir)

@testset "TestImages" begin
    @testset "Autumn leaves" begin
        name = "autumn_leaves"
        img = testimage(name)
        @test ndims(img) == 2
        @test eltype(img) == RGBA{N0f16}
        @test size(img) == (105, 140)
        out_name = joinpath(mydir, name * ".png")
        save(out_name, img)
        oimg = load(out_name)
        @test size(oimg) == size(img)
        @test eltype(oimg) == eltype(img)
        # Note: around 10% of pixel values are not identical. Precision?
        # @test oimg == img
    end
    @testset "Camerman" begin
        name = "cameraman"
        img = testimage(name)
        @test ndims(img) == 2
        @test eltype(img) == Gray{N0f8}
        @test size(img) == (512, 512)
        out_name = joinpath(mydir, name * ".tif")
        save(out_name, img)
        oimg = load(out_name)
        @test size(oimg) == size(img)
        @test eltype(oimg) == eltype(img)
        @test oimg == img
    end
    @testset "Earth Apollo" begin
        name = "earth_apollo17"
        img = testimage(name)
        @test ndims(img) == 2
        @test eltype(img) == RGB4{N0f8}
        @test size(img) == (3002, 3000)
        out_name = joinpath(mydir, name * ".png")
        save(out_name, img)
        oimg = load(out_name)
        @test size(oimg) == size(img)
        @test eltype(oimg) == eltype(img)
        @test oimg == img
    end
    @testset "Fabio" begin
        name = "fabio"
        img = testimage(name)
        @test ndims(img) == 2
        @test eltype(img) == RGB4{N0f8}
        @test size(img) == (256, 256)
        out_name = joinpath(mydir, name * ".png")
        save(out_name, img)
        oimg = load(out_name)
        @test size(oimg) == size(img)
        @test eltype(oimg) == eltype(img)
        @test oimg == img
    end
    @testset "House" begin
        # Note: this is a GrayA image, not supported by Apple
        # So saving as a Gray for now.
        name = "house"
        img = testimage(name)
        @test ndims(img) == 2
        @test eltype(img) == GrayA{N0f8}
        @test size(img) == (512, 512)
        out_name = joinpath(mydir, name * ".tif")
        save(out_name, img)
        oimg = load(out_name)
        @test size(oimg) == size(img)
        @test eltype(oimg) == Gray{N0f8}
        @test oimg == Gray.(img)
    end
    @testset "Jetplane" begin
        # Note: this is a GrayA image, not supported by Apple
        # So saving as a Gray for now.
        name = "jetplane"
        img = testimage(name)
        @test ndims(img) == 2
        @test eltype(img) == GrayA{N0f8}
        @test size(img) == (512, 512)
        out_name = joinpath(mydir, name * ".tif")
        save(out_name, img)
        oimg = load(out_name)
        @test size(oimg) == size(img)
        @test eltype(oimg) == Gray{N0f8}
        @test oimg == Gray.(img)
    end
    @testset "Lighthouse" begin
        name = "lighthouse"
        img = testimage(name)
        @test ndims(img) == 2
        @test eltype(img) == RGB4{N0f8}
        @test size(img) == (512, 768)
        out_name = joinpath(mydir, name * ".png")
        save(out_name, img)
        oimg = load(out_name)
        @test size(oimg) == size(img)
        @test eltype(oimg) == eltype(img)
        @test oimg == img
    end
    @testset "Mandrill" begin
        name = "mandrill"
        img = testimage(name)
        @test ndims(img) == 2
        @test eltype(img) == RGB{N0f8}
        @test size(img) == (512, 512)
        out_name = joinpath(mydir, name * ".tiff")
        save(out_name, img)
        oimg = load(out_name)
        @test size(oimg) == size(img)
        @test eltype(oimg) == eltype(img)
        @test oimg == img
    end
    @testset "Moonsurface" begin
        name = "moonsurface"
        img = testimage(name)
        @test ndims(img) == 2
        @test eltype(img) == Gray{N0f8}
        @test size(img) == (256, 256)
        out_name = joinpath(mydir, name * ".tif")
        save(out_name, img)
        oimg = load(out_name)
        @test size(oimg) == size(img)
        @test eltype(oimg) == eltype(img)
        @test oimg == img
    end
    @testset "Mountainstream" begin
        name = "mountainstream"
        img = testimage(name)
        @test ndims(img) == 2
        @test eltype(img) == RGB4{N0f8}
        @test size(img) == (512, 768)
        out_name = joinpath(mydir, name * ".png")
        save(out_name, img)
        oimg = load(out_name)
        @test size(oimg) == size(img)
        @test eltype(oimg) == eltype(img)
        @test oimg == img
    end
    @testset "MRI Stack" begin
        name = "mri-stack"
        img = testimage(name)
        @test isa(img, AxisArray)
        @test map(step, axisvalues(img)) == (1, 1, 5)
        @test ndims(img) == 3
        @test eltype(img) == Gray{N0f8}
        @test size(img) == (226, 186, 27)
        out_name = joinpath(mydir, name * ".tif")
        # This TestImage has a special case, labeling the axes.  Pop it out.
        save(out_name, img.data)
        oimg = load(out_name)
        @test size(oimg) == size(img)
        @test eltype(oimg) == eltype(img)
        @test oimg == img
    end
    @testset "M51" begin
        name = "m51"
        img = testimage(name)
        @test ndims(img) == 2
        @test eltype(img) == Gray{N0f16}
        @test size(img) == (510, 320)
        out_name = joinpath(mydir, name * ".tif")
        save(out_name, img)
        oimg = load(out_name)
        @test size(oimg) == size(img)
        @test eltype(oimg) == eltype(img)
        @test oimg == img
    end
    @testset "HeLa cells" begin
        name = "hela-cells"
        img = testimage(name)
        @test ndims(img) == 2
        @test eltype(img) == RGB{N0f16}
        @test size(img) == (512, 672)
        out_name = joinpath(mydir, name * ".tif")
        save(out_name, img)
        oimg = load(out_name)
        @test size(oimg) == size(img)
        @test eltype(oimg) == eltype(img)
        @test oimg == img
    end
    @testset "Blobs GIF" begin
        name = "blobs"
        img = testimage(name)
        @test ndims(img) == 2
        @test eltype(img) == RGB4{N0f8}
        @test size(img) == (254, 256)
        # Unclear why this needs to be a .png
        out_name = joinpath(mydir, name * ".png")
        save(out_name, img)
        oimg = load(out_name)
        @test oimg == img
    end
    @testset "Multichannel timeseries OME" begin
        # Note: the existence of OMETIFF.jl breaks this test,
        # even though we don't explicitly import it here.
        # So, we need to do some extra work to call the image
        # with our specific loader here.
        mcifile = "multi-channel-time-series.ome.tif"
        name = joinpath(TestImages.imagedir, mcifile)
        if !isfile(name)
            REPO_URL = "https://github.com/JuliaImages/TestImages.jl/blob/gh-pages/images/"
            download(REPO_URL*mcifile*"?raw=true", name)
        end
        img = QuartzImageIO.load(name)
        @test ndims(img) == 3
        @test eltype(img) == Gray{N0f8}
        @test size(img) == (167, 439, 21)
        out_name = joinpath(mydir, name * ".tif")
        save(out_name, img)
        oimg = load(out_name)
        @test size(oimg) == size(img)
        @test eltype(oimg) == eltype(img)
        @test oimg == img
    end
end

@testset "ImageMeta" begin
    # https://github.com/sisl/PGFPlots.jl/issues/5
    img = ImageMeta(rand(RGB{N0f8}, 3, 5))
    out_name = joinpath(mydir, "imagemeta.png")
    save(out_name, img)
    oimg = load(out_name)
    @test oimg == img
end

@testset "Saving" begin
    @testset "RGB" begin
        imgc = rand(RGB{Float32}, 40, 30)
        out_name = joinpath(mydir, "float32.png")
        save(out_name, imgc)
        inimg = load(out_name)
        @test size(imgc) == size(inimg)
    end

    @testset "HSV" begin
        # issue 37
        Random.seed!(42)
        img = rand(HSV{Float32}, 8, 10)
        out_name = joinpath(mydir, "hsv.png")
        save(out_name, img)
        oimg = load(out_name)
        @test HSV.(oimg)[1].h ≈ img[1].h rtol=0.02 atol=0.003
        @test HSV.(oimg)[1].s ≈ img[1].s rtol=0.02 atol=0.003
        @test HSV.(oimg)[1].v ≈ img[1].v rtol=0.02 atol=0.003
    end

    @testset "HSL" begin
        Random.seed!(42)
        img = rand(HSL{Float32}, 8, 10)
        out_name = joinpath(mydir, "hsl.png")
        save(out_name, img)
        oimg = load(out_name)
        @test HSL.(oimg)[1].h ≈ img[1].h rtol=0.04 atol=0.003
        @test HSL.(oimg)[1].s ≈ img[1].s rtol=0.02 atol=0.005
        @test HSL.(oimg)[1].l ≈ img[1].l rtol=0.02 atol=0.003
    end

    @testset "Lab" begin
        Random.seed!(42)
        img = rand(Lab{Float32}, 8, 10)
        out_name = joinpath(mydir, "Lab.png")
        save(out_name, img)
        oimg = Lab.(load(out_name))
        for I in eachindex(img)
            # Comparing oimg[I] against img[I] may fail because
            # information was lost when saving the Lab values as
            # RGB. Instead, we just verify that oimg[I] = Lab(RGB(img[I]))
            expected = Lab(RGB(img[I]))
            @test oimg[I].l ≈ expected.l rtol=0.07
            @test oimg[I].a ≈ expected.a rtol=0.07 atol=0.35
            @test oimg[I].b ≈ expected.b rtol=0.07 atol=0.35
        end
    end
end

@testset "Streams" begin
    name = "lighthouse"
    img = testimage(name)
    out_name = joinpath(mydir, name * ".png")
    @testset "saving" begin
        open(out_name, "w") do io
            QuartzImageIO.save(Stream(format"PNG", io), img)
        end
        imgcmp = load(out_name)
        @test imgcmp == img
    end
    @testset "loading" begin
        imgcmp = open(out_name) do io
            QuartzImageIO.load(Stream(format"PNG", io))
        end
        @test imgcmp == img
    end
end

@testset "OffsetArrays" begin
    img = OffsetArray(Gray{N0f8}[1 0; 0 1], 0:1, 3:4)
    fn = joinpath(mydir, "indices.png")
    save(fn, img)
    imgr = load(fn)
    @test imgr == parent(img)
end

rm(mydir, recursive=true)

@test !@isdefined ImageMagick

ext: jl


text:
############################################################
## joAbstractOperator - outer constructors
############################################################

# this file is currently not included in joAbstractOperator.jl
# common constructor exports are in respective */constructors.jl files


text:
function DividedDifferenceInterpolator()
    return DividedDifferenceInterpolator(())
end

function interpolate(obj::DividedDifferenceInterpolator, arg0::Vector{jdouble}, arg1::Vector{jdouble})
    return jcall(obj, "interpolate", PolynomialFunctionNewtonForm, (Vector{jdouble}, Vector{jdouble}), arg0, arg1)
end


text:
const CPLE_None = 0
const CPLE_AppDefined = 1
const CPLE_OutOfMemory = 2
const CPLE_FileIO = 3
const CPLE_OpenFailed = 4
const CPLE_IllegalArg = 5
const CPLE_NotSupported = 6
const CPLE_AssertionFailed = 7
const CPLE_NoWriteAccess = 8
const CPLE_UserInterrupt = 9
const CPLE_ObjectNull = 10
const CPLE_HttpResponse = 11
const CPLE_AWSBucketNotFound = 12
const CPLE_AWSObjectNotFound = 13
const CPLE_AWSAccessDenied = 14
const CPLE_AWSInvalidCredentials = 15
const CPLE_AWSSignatureDoesNotMatch = 16
@enum CPLErr::UInt32 CE_None = 0 CE_Debug = 1 CE_Warning = 2 CE_Failure = 3 CE_Fatal = 4
const CPLErrorNum = Cint
@enum CPLXMLNodeType::UInt32 CXT_Element = 0 CXT_Text = 1 CXT_Attribute = 2 CXT_Comment = 3 CXT_Literal = 4

mutable struct CPLXMLNode
    eType::CPLXMLNodeType
    pszValue::Cstring
    psNext::Ptr{Void}
    psChild::Ptr{Void}
end

const CPL_HAS_GINT64 = 1
const CPL_FRMT_GB_WITHOUT_PREFIX = "ll"
const CPL_IS_LSB = 1
const FALSE = 0
const TRUE = 1
const GInt32 = Cint
const GUInt32 = UInt32
const GInt16 = Int16
const GUInt16 = UInt16
const GByte = Cuchar
const GBool = Cint
const GIntBig = Clonglong
const GUIntBig = Culonglong
const GInt64 = GIntBig
const GUInt64 = GUIntBig
const GPtrDiff_t = GIntBig
@enum CPLVirtualMemAccessMode::UInt32 VIRTUALMEM_READONLY = 0 VIRTUALMEM_READONLY_ENFORCED = 1 VIRTUALMEM_READWRITE = 2
const VSI_STAT_EXISTS_FLAG = 0x01
const VSI_STAT_NATURE_FLAG = 0x02
const VSI_STAT_SIZE_FLAG = 0x04
const VSI_STAT_SET_ERROR_FLAG = 0x08
const VSIStatBuf = Void
const vsi_l_offset = GUIntBig
const VSILFILE = FILE
@enum VSIRangeStatus::UInt32 VSI_RANGE_STATUS_UNKNOWN = 0 VSI_RANGE_STATUS_DATA = 1 VSI_RANGE_STATUS_HOLE = 2
const RASTERIO_EXTRA_ARG_CURRENT_VERSION = 1
const GDALMD_AREA_OR_POINT = "AREA_OR_POINT"
const GDALMD_AOP_AREA = "Area"
const GDALMD_AOP_POINT = "Point"
const GDAL_DMD_LONGNAME = "DMD_LONGNAME"
const GDAL_DMD_HELPTOPIC = "DMD_HELPTOPIC"
const GDAL_DMD_MIMETYPE = "DMD_MIMETYPE"
const GDAL_DMD_EXTENSION = "DMD_EXTENSION"
const GDAL_DMD_CONNECTION_PREFIX = "DMD_CONNECTION_PREFIX"
const GDAL_DMD_EXTENSIONS = "DMD_EXTENSIONS"
const GDAL_DMD_CREATIONOPTIONLIST = "DMD_CREATIONOPTIONLIST"
const GDAL_DMD_OPENOPTIONLIST = "DMD_OPENOPTIONLIST"
const GDAL_DMD_CREATIONDATATYPES = "DMD_CREATIONDATATYPES"
const GDAL_DMD_CREATIONFIELDDATATYPES = "DMD_CREATIONFIELDDATATYPES"
const GDAL_DMD_SUBDATASETS = "DMD_SUBDATASETS"
const GDAL_DCAP_OPEN = "DCAP_OPEN"
const GDAL_DCAP_CREATE = "DCAP_CREATE"
const GDAL_DCAP_CREATECOPY = "DCAP_CREATECOPY"
const GDAL_DCAP_VIRTUALIO = "DCAP_VIRTUALIO"
const GDAL_DCAP_RASTER = "DCAP_RASTER"
const GDAL_DCAP_VECTOR = "DCAP_VECTOR"
const GDAL_DCAP_GNM = "DCAP_GNM"
const GDAL_DCAP_NOTNULL_FIELDS = "DCAP_NOTNULL_FIELDS"
const GDAL_DCAP_DEFAULT_FIELDS = "DCAP_DEFAULT_FIELDS"
const GDAL_DCAP_NOTNULL_GEOMFIELDS = "DCAP_NOTNULL_GEOMFIELDS"
const GDAL_OF_READONLY = 0x00
const GDAL_OF_UPDATE = 0x01
const GDAL_OF_ALL = 0x00
const GDAL_OF_RASTER = 0x02
const GDAL_OF_VECTOR = 0x04
const GDAL_OF_GNM = 0x08
const GDAL_OF_KIND_MASK = 0x1e
const GDAL_OF_SHARED = 0x20
const GDAL_OF_VERBOSE_ERROR = 0x40
const GDAL_OF_INTERNAL = 0x80
const GDAL_OF_DEFAULT_BLOCK_ACCESS = 0
const GDAL_OF_ARRAY_BLOCK_ACCESS = 0x0100
const GDAL_OF_HASHSET_BLOCK_ACCESS = 0x0200
const GDAL_OF_RESERVED_1 = 0x0300
const GDAL_OF_BLOCK_ACCESS_MASK = 0x0300
const GDAL_DS_LAYER_CREATIONOPTIONLIST = "DS_LAYER_CREATIONOPTIONLIST"
const GMF_ALL_VALID = 0x01
const GMF_PER_DATASET = 0x02
const GMF_ALPHA = 0x04
const GMF_NODATA = 0x08
const GDAL_DATA_COVERAGE_STATUS_UNIMPLEMENTED = 0x01
const GDAL_DATA_COVERAGE_STATUS_DATA = 0x02
const GDAL_DATA_COVERAGE_STATUS_EMPTY = 0x04
@enum GDALDataType::UInt32 GDT_Unknown = 0 GDT_Byte = 1 GDT_UInt16 = 2 GDT_Int16 = 3 GDT_UInt32 = 4 GDT_Int32 = 5 GDT_Float32 = 6 GDT_Float64 = 7 GDT_CInt16 = 8 GDT_CInt32 = 9 GDT_CFloat32 = 10 GDT_CFloat64 = 11 GDT_TypeCount = 12
@enum GDALAsyncStatusType::UInt32 GARIO_PENDING = 0 GARIO_UPDATE = 1 GARIO_ERROR = 2 GARIO_COMPLETE = 3 GARIO_TypeCount = 4
@enum GDALAccess::UInt32 GA_ReadOnly = 0 GA_Update = 1
@enum GDALRWFlag::UInt32 GF_Read = 0 GF_Write = 1
@enum GDALRIOResampleAlg::UInt32 GRIORA_NearestNeighbour = 0 GRIORA_Bilinear = 1 GRIORA_Cubic = 2 GRIORA_CubicSpline = 3 GRIORA_Lanczos = 4 GRIORA_Average = 5 GRIORA_Mode = 6 GRIORA_Gauss = 7

mutable struct GDALRasterIOExtraArg
    nVersion::Cint
    eResampleAlg::GDALRIOResampleAlg
    pfnProgress::GDALProgressFunc
    pProgressData::Ptr{Void}
    bFloatingPointWindowValidity::Cint
    dfXOff::Cdouble
    dfYOff::Cdouble
    dfXSize::Cdouble
    dfYSize::Cdouble
end

@enum GDALColorInterp::UInt32 GCI_Undefined = 0 GCI_GrayIndex = 1 GCI_PaletteIndex = 2 GCI_RedBand = 3 GCI_GreenBand = 4 GCI_BlueBand = 5 GCI_AlphaBand = 6 GCI_HueBand = 7 GCI_SaturationBand = 8 GCI_LightnessBand = 9 GCI_CyanBand = 10 GCI_MagentaBand = 11 GCI_YellowBand = 12 GCI_BlackBand = 13 GCI_YCbCr_YBand = 14 GCI_YCbCr_CbBand = 15 GCI_YCbCr_CrBand = 16
const GCI_Max = GCI_YCbCr_CrBand
@enum GDALPaletteInterp::UInt32 GPI_Gray = 0 GPI_RGB = 1 GPI_CMYK = 2 GPI_HLS = 3
const GSpacing = GIntBig

mutable struct GDAL_GCP
    pszId::Cstring
    pszInfo::Cstring
    dfGCPPixel::Cdouble
    dfGCPLine::Cdouble
    dfGCPX::Cdouble
    dfGCPY::Cdouble
    dfGCPZ::Cdouble
end

mutable struct GDALRPCInfo
    dfLINE_OFF::Cdouble
    dfSAMP_OFF::Cdouble
    dfLAT_OFF::Cdouble
    dfLONG_OFF::Cdouble
    dfHEIGHT_OFF::Cdouble
    dfLINE_SCALE::Cdouble
    dfSAMP_SCALE::Cdouble
    dfLAT_SCALE::Cdouble
    dfLONG_SCALE::Cdouble
    dfHEIGHT_SCALE::Cdouble
    adfLINE_NUM_COEFF::NTuple{20, Cdouble}
    adfLINE_DEN_COEFF::NTuple{20, Cdouble}
    adfSAMP_NUM_COEFF::NTuple{20, Cdouble}
    adfSAMP_DEN_COEFF::NTuple{20, Cdouble}
    dfMIN_LONG::Cdouble
    dfMIN_LAT::Cdouble
    dfMAX_LONG::Cdouble
    dfMAX_LAT::Cdouble
end

mutable struct GDALColorEntry
    c1::Int16
    c2::Int16
    c3::Int16
    c4::Int16
end

@enum GDALRATFieldType::UInt32 GFT_Integer = 0 GFT_Real = 1 GFT_String = 2
@enum GDALRATFieldUsage::UInt32 GFU_Generic = 0 GFU_PixelCount = 1 GFU_Name = 2 GFU_Min = 3 GFU_Max = 4 GFU_MinMax = 5 GFU_Red = 6 GFU_Green = 7 GFU_Blue = 8 GFU_Alpha = 9 GFU_RedMin = 10 GFU_GreenMin = 11 GFU_BlueMin = 12 GFU_AlphaMin = 13 GFU_RedMax = 14 GFU_GreenMax = 15 GFU_BlueMax = 16 GFU_AlphaMax = 17 GFU_MaxCount = 18
@enum GDALTileOrganization::UInt32 GTO_TIP = 0 GTO_BIT = 1 GTO_BSQ = 2
const GDAL_GTI2_SIGNATURE = "GTI2"

mutable struct GDALTransformerInfo
    abySignature::NTuple{4, GByte}
    pszClassName::Cstring
    pfnTransform::GDALTransformerFunc
    pfnCleanup::Ptr{Void}
    pfnSerialize::Ptr{Void}
    pfnCreateSimilar::Ptr{Void}
end

mutable struct OGRContourWriterInfo
    hLayer::Ptr{Void}
    adfGeoTransform::NTuple{6, Cdouble}
    nElevField::Cint
    nIDField::Cint
    nNextID::Cint
end

@enum GDALGridAlgorithm::UInt32 GGA_InverseDistanceToAPower = 1 GGA_MovingAverage = 2 GGA_NearestNeighbor = 3 GGA_MetricMinimum = 4 GGA_MetricMaximum = 5 GGA_MetricRange = 6 GGA_MetricCount = 7 GGA_MetricAverageDistance = 8 GGA_MetricAverageDistancePts = 9 GGA_Linear = 10 GGA_InverseDistanceToAPowerNearestNeighbor = 11

mutable struct GDALGridInverseDistanceToAPowerOptions
    dfPower::Cdouble
    dfSmoothing::Cdouble
    dfAnisotropyRatio::Cdouble
    dfAnisotropyAngle::Cdouble
    dfRadius1::Cdouble
    dfRadius2::Cdouble
    dfAngle::Cdouble
    nMaxPoints::GUInt32
    nMinPoints::GUInt32
    dfNoDataValue::Cdouble
end

mutable struct GDALGridInverseDistanceToAPowerNearestNeighborOptions
    dfPower::Cdouble
    dfRadius::Cdouble
    dfSmoothing::Cdouble
    nMaxPoints::GUInt32
    nMinPoints::GUInt32
    dfNoDataValue::Cdouble
end

mutable struct GDALGridMovingAverageOptions
    dfRadius1::Cdouble
    dfRadius2::Cdouble
    dfAngle::Cdouble
    nMinPoints::GUInt32
    dfNoDataValue::Cdouble
end

mutable struct GDALGridNearestNeighborOptions
    dfRadius1::Cdouble
    dfRadius2::Cdouble
    dfAngle::Cdouble
    dfNoDataValue::Cdouble
end

mutable struct GDALGridDataMetricsOptions
    dfRadius1::Cdouble
    dfRadius2::Cdouble
    dfAngle::Cdouble
    nMinPoints::GUInt32
    dfNoDataValue::Cdouble
end

mutable struct GDALGridLinearOptions
    dfRadius::Cdouble
    dfNoDataValue::Cdouble
end

mutable struct GDALTriFacet
    anVertexIdx::NTuple{3, Cint}
    anNeighborIdx::NTuple{3, Cint}
end

mutable struct GDALTriBarycentricCoefficients
    dfMul1X::Cdouble
    dfMul1Y::Cdouble
    dfMul2X::Cdouble
    dfMul2Y::Cdouble
    dfCstX::Cdouble
    dfCstY::Cdouble
end

mutable struct GDALTriangulation
    nFacets::Cint
    pasFacets::Ptr{GDALTriFacet}
    pasFacetCoefficients::Ptr{GDALTriBarycentricCoefficients}
end

const OGRERR_NONE = 0
const OGRERR_NOT_ENOUGH_DATA = 1
const OGRERR_NOT_ENOUGH_MEMORY = 2
const OGRERR_UNSUPPORTED_GEOMETRY_TYPE = 3
const OGRERR_UNSUPPORTED_OPERATION = 4
const OGRERR_CORRUPT_DATA = 5
const OGRERR_FAILURE = 6
const OGRERR_UNSUPPORTED_SRS = 7
const OGRERR_INVALID_HANDLE = 8
const OGRERR_NON_EXISTING_FEATURE = 9
const wkb25DBit = 0x80000000
const ogrZMarker = 0x21125711
const ALTER_NAME_FLAG = 0x01
const ALTER_TYPE_FLAG = 0x02
const ALTER_WIDTH_PRECISION_FLAG = 0x04
const ALTER_NULLABLE_FLAG = 0x08
const ALTER_DEFAULT_FLAG = 0x10
const ALTER_ALL_FLAG = (((ALTER_NAME_FLAG | ALTER_TYPE_FLAG) | ALTER_WIDTH_PRECISION_FLAG) | ALTER_NULLABLE_FLAG) | ALTER_DEFAULT_FLAG
const OGR_F_VAL_NULL = 0x00000001
const OGR_F_VAL_GEOM_TYPE = 0x00000002
const OGR_F_VAL_WIDTH = 0x00000004
const OGR_F_VAL_ALLOW_NULL_WHEN_DEFAULT = 0x00000008
const OGR_F_VAL_ALLOW_DIFFERENT_GEOM_DIM = 0x00000010
const OGRNullFID = -1
const OGRUnsetMarker = -21121
const OGRNullMarker = -21122
const OLCRandomRead = "RandomRead"
const OLCSequentialWrite = "SequentialWrite"
const OLCRandomWrite = "RandomWrite"
const OLCFastSpatialFilter = "FastSpatialFilter"
const OLCFastFeatureCount = "FastFeatureCount"
const OLCFastGetExtent = "FastGetExtent"
const OLCCreateField = "CreateField"
const OLCDeleteField = "DeleteField"
const OLCReorderFields = "ReorderFields"
const OLCAlterFieldDefn = "AlterFieldDefn"
const OLCTransactions = "Transactions"
const OLCDeleteFeature = "DeleteFeature"
const OLCFastSetNextByIndex = "FastSetNextByIndex"
const OLCStringsAsUTF8 = "StringsAsUTF8"
const OLCIgnoreFields = "IgnoreFields"
const OLCCreateGeomField = "CreateGeomField"
const OLCCurveGeometries = "CurveGeometries"
const OLCMeasuredGeometries = "MeasuredGeometries"
const ODsCCreateLayer = "CreateLayer"
const ODsCDeleteLayer = "DeleteLayer"
const ODsCCreateGeomFieldAfterCreateLayer = "CreateGeomFieldAfterCreateLayer"
const ODsCCurveGeometries = "CurveGeometries"
const ODsCTransactions = "Transactions"
const ODsCEmulatedTransactions = "EmulatedTransactions"
const ODsCMeasuredGeometries = "MeasuredGeometries"
const ODsCRandomLayerRead = "RandomLayerRead"
const ODsCRandomLayerWrite = "RandomLayerWrite "
const ODrCCreateDataSource = "CreateDataSource"
const ODrCDeleteDataSource = "DeleteDataSource"
const OLMD_FID64 = "OLMD_FID64"

mutable struct OGREnvelope
    MinX::Cdouble
    MaxX::Cdouble
    MinY::Cdouble
    MaxY::Cdouble
end

mutable struct OGREnvelope3D
    MinX::Cdouble
    MaxX::Cdouble
    MinY::Cdouble
    MaxY::Cdouble
    MinZ::Cdouble
    MaxZ::Cdouble
end

const OGRErr = Cint
const OGRBoolean = Cint
@enum OGRwkbGeometryType::UInt32 wkbUnknown = 0 wkbPoint = 1 wkbLineString = 2 wkbPolygon = 3 wkbMultiPoint = 4 wkbMultiLineString = 5 wkbMultiPolygon = 6 wkbGeometryCollection = 7 wkbCircularString = 8 wkbCompoundCurve = 9 wkbCurvePolygon = 10 wkbMultiCurve = 11 wkbMultiSurface = 12 wkbCurve = 13 wkbSurface = 14 wkbPolyhedralSurface = 15 wkbTIN = 16 wkbTriangle = 17 wkbNone = 100 wkbLinearRing = 101 wkbCircularStringZ = 1008 wkbCompoundCurveZ = 1009 wkbCurvePolygonZ = 1010 wkbMultiCurveZ = 1011 wkbMultiSurfaceZ = 1012 wkbCurveZ = 1013 wkbSurfaceZ = 1014 wkbPolyhedralSurfaceZ = 1015 wkbTINZ = 1016 wkbTriangleZ = 1017 wkbPointM = 2001 wkbLineStringM = 2002 wkbPolygonM = 2003 wkbMultiPointM = 2004 wkbMultiLineStringM = 2005 wkbMultiPolygonM = 2006 wkbGeometryCollectionM = 2007 wkbCircularStringM = 2008 wkbCompoundCurveM = 2009 wkbCurvePolygonM = 2010 wkbMultiCurveM = 2011 wkbMultiSurfaceM = 2012 wkbCurveM = 2013 wkbSurfaceM = 2014 wkbPolyhedralSurfaceM = 2015 wkbTINM = 2016 wkbTriangleM = 2017 wkbPointZM = 3001 wkbLineStringZM = 3002 wkbPolygonZM = 3003 wkbMultiPointZM = 3004 wkbMultiLineStringZM = 3005 wkbMultiPolygonZM = 3006 wkbGeometryCollectionZM = 3007 wkbCircularStringZM = 3008 wkbCompoundCurveZM = 3009 wkbCurvePolygonZM = 3010 wkbMultiCurveZM = 3011 wkbMultiSurfaceZM = 3012 wkbCurveZM = 3013 wkbSurfaceZM = 3014 wkbPolyhedralSurfaceZM = 3015 wkbTINZM = 3016 wkbTriangleZM = 3017 wkbPoint25D = 2147483649 wkbLineString25D = 2147483650 wkbPolygon25D = 2147483651 wkbMultiPoint25D = 2147483652 wkbMultiLineString25D = 2147483653 wkbMultiPolygon25D = 2147483654 wkbGeometryCollection25D = 2147483655
@enum OGRwkbVariant::UInt32 wkbVariantOldOgc = 0 wkbVariantIso = 1 wkbVariantPostGIS1 = 2
@enum OGRwkbByteOrder::UInt32 wkbXDR = 0 wkbNDR = 1
@enum OGRFieldType::UInt32 OFTInteger = 0 OFTIntegerList = 1 OFTReal = 2 OFTRealList = 3 OFTString = 4 OFTStringList = 5 OFTWideString = 6 OFTWideStringList = 7 OFTBinary = 8 OFTDate = 9 OFTTime = 10 OFTDateTime = 11 OFTInteger64 = 12 OFTInteger64List = 13
const OFTMaxType = OFTInteger64List
@enum OGRFieldSubType::UInt32 OFSTNone = 0 OFSTBoolean = 1 OFSTInt16 = 2 OFSTFloat32 = 3
const OFSTMaxSubType = OFSTFloat32
@enum OGRJustification::UInt32 OJUndefined = 0 OJLeft = 1 OJRight = 2
const OGRField = Void
@enum OGRSTClassId::UInt32 OGRSTCNone = 0 OGRSTCPen = 1 OGRSTCBrush = 2 OGRSTCSymbol = 3 OGRSTCLabel = 4 OGRSTCVector = 5
@enum OGRSTUnitId::UInt32 OGRSTUGround = 0 OGRSTUPixel = 1 OGRSTUPoints = 2 OGRSTUMM = 3 OGRSTUCM = 4 OGRSTUInches = 5
@enum OGRSTPenParam::UInt32 OGRSTPenColor = 0 OGRSTPenWidth = 1 OGRSTPenPattern = 2 OGRSTPenId = 3 OGRSTPenPerOffset = 4 OGRSTPenCap = 5 OGRSTPenJoin = 6 OGRSTPenPriority = 7 OGRSTPenLast = 8
@enum OGRSTBrushParam::UInt32 OGRSTBrushFColor = 0 OGRSTBrushBColor = 1 OGRSTBrushId = 2 OGRSTBrushAngle = 3 OGRSTBrushSize = 4 OGRSTBrushDx = 5 OGRSTBrushDy = 6 OGRSTBrushPriority = 7 OGRSTBrushLast = 8
@enum OGRSTSymbolParam::UInt32 OGRSTSymbolId = 0 OGRSTSymbolAngle = 1 OGRSTSymbolColor = 2 OGRSTSymbolSize = 3 OGRSTSymbolDx = 4 OGRSTSymbolDy = 5 OGRSTSymbolStep = 6 OGRSTSymbolPerp = 7 OGRSTSymbolOffset = 8 OGRSTSymbolPriority = 9 OGRSTSymbolFontName = 10 OGRSTSymbolOColor = 11 OGRSTSymbolLast = 12
@enum OGRSTLabelParam::UInt32 OGRSTLabelFontName = 0 OGRSTLabelSize = 1 OGRSTLabelTextString = 2 OGRSTLabelAngle = 3 OGRSTLabelFColor = 4 OGRSTLabelBColor = 5 OGRSTLabelPlacement = 6 OGRSTLabelAnchor = 7 OGRSTLabelDx = 8 OGRSTLabelDy = 9 OGRSTLabelPerp = 10 OGRSTLabelBold = 11 OGRSTLabelItalic = 12 OGRSTLabelUnderline = 13 OGRSTLabelPriority = 14 OGRSTLabelStrikeout = 15 OGRSTLabelStretch = 16 OGRSTLabelAdjHor = 17 OGRSTLabelAdjVert = 18 OGRSTLabelHColor = 19 OGRSTLabelOColor = 20 OGRSTLabelLast = 21

mutable struct CPLSharedFileInfo
    fp::Ptr{FILE}
    nRefCount::Cint
    bLarge::Cint
    pszFilename::Cstring
    pszAccess::Cstring
end

@enum GDALResampleAlg::UInt32 GRA_NearestNeighbour = 0 GRA_Bilinear = 1 GRA_Cubic = 2 GRA_CubicSpline = 3 GRA_Lanczos = 4 GRA_Average = 5 GRA_Mode = 6 GRA_Max = 8 GRA_Min = 9 GRA_Med = 10 GRA_Q1 = 11 GRA_Q3 = 12
@enum GWKAverageOrModeAlg::UInt32 GWKAOM_Average = 1 GWKAOM_Fmode = 2 GWKAOM_Imode = 3 GWKAOM_Max = 4 GWKAOM_Min = 5 GWKAOM_Quant = 6

mutable struct GDALWarpOptions
    papszWarpOptions::Ptr{Cstring}
    dfWarpMemoryLimit::Cdouble
    eResampleAlg::GDALResampleAlg
    eWorkingDataType::GDALDataType
    hSrcDS::GDALDatasetH
    hDstDS::GDALDatasetH
    nBandCount::Cint
    panSrcBands::Ptr{Cint}
    panDstBands::Ptr{Cint}
    nSrcAlphaBand::Cint
    nDstAlphaBand::Cint
    padfSrcNoDataReal::Ptr{Cdouble}
    padfSrcNoDataImag::Ptr{Cdouble}
    padfDstNoDataReal::Ptr{Cdouble}
    padfDstNoDataImag::Ptr{Cdouble}
    pfnProgress::GDALProgressFunc
    pProgressArg::Ptr{Void}
    pfnTransformer::GDALTransformerFunc
    pTransformerArg::Ptr{Void}
    papfnSrcPerBandValidityMaskFunc::Ptr{GDALMaskFunc}
    papSrcPerBandValidityMaskFuncArg::Ptr{Ptr{Void}}
    pfnSrcValidityMaskFunc::GDALMaskFunc
    pSrcValidityMaskFuncArg::Ptr{Void}
    pfnSrcDensityMaskFunc::GDALMaskFunc
    pSrcDensityMaskFuncArg::Ptr{Void}
    pfnDstDensityMaskFunc::GDALMaskFunc
    pDstDensityMaskFuncArg::Ptr{Void}
    pfnDstValidityMaskFunc::GDALMaskFunc
    pDstValidityMaskFuncArg::Ptr{Void}
    pfnPreWarpChunkProcessor::Ptr{Void}
    pPreWarpProcessorArg::Ptr{Void}
    pfnPostWarpChunkProcessor::Ptr{Void}
    pPostWarpProcessorArg::Ptr{Void}
    hCutline::Ptr{Void}
    dfCutlineBlendDist::Cdouble
end

const VRT_NODATA_UNSET = -1234.56
const SRS_WKT_WGS84 = "GEOGCS[\"WGS 84\",DATUM[\"WGS_1984\",SPHEROID[\"WGS 84\",6378137,298.257223563,AUTHORITY[\"EPSG\",\"7030\"]],AUTHORITY[\"EPSG\",\"6326\"]],PRIMEM[\"Greenwich\",0,AUTHORITY[\"EPSG\",\"8901\"]],UNIT[\"degree\",0.0174532925199433,AUTHORITY[\"EPSG\",\"9122\"]],AUTHORITY[\"EPSG\",\"4326\"]]"
const SRS_PT_ALBERS_CONIC_EQUAL_AREA = "Albers_Conic_Equal_Area"
const SRS_PT_AZIMUTHAL_EQUIDISTANT = "Azimuthal_Equidistant"
const SRS_PT_CASSINI_SOLDNER = "Cassini_Soldner"
const SRS_PT_CYLINDRICAL_EQUAL_AREA = "Cylindrical_Equal_Area"
const SRS_PT_BONNE = "Bonne"
const SRS_PT_ECKERT_I = "Eckert_I"
const SRS_PT_ECKERT_II = "Eckert_II"
const SRS_PT_ECKERT_III = "Eckert_III"
const SRS_PT_ECKERT_IV = "Eckert_IV"
const SRS_PT_ECKERT_V = "Eckert_V"
const SRS_PT_ECKERT_VI = "Eckert_VI"
const SRS_PT_EQUIDISTANT_CONIC = "Equidistant_Conic"
const SRS_PT_EQUIRECTANGULAR = "Equirectangular"
const SRS_PT_GALL_STEREOGRAPHIC = "Gall_Stereographic"
const SRS_PT_GAUSSSCHREIBERTMERCATOR = "Gauss_Schreiber_Transverse_Mercator"
const SRS_PT_GEOSTATIONARY_SATELLITE = "Geostationary_Satellite"
const SRS_PT_GOODE_HOMOLOSINE = "Goode_Homolosine"
const SRS_PT_IGH = "Interrupted_Goode_Homolosine"
const SRS_PT_GNOMONIC = "Gnomonic"
const SRS_PT_HOTINE_OBLIQUE_MERCATOR_AZIMUTH_CENTER = "Hotine_Oblique_Mercator_Azimuth_Center"
const SRS_PT_HOTINE_OBLIQUE_MERCATOR = "Hotine_Oblique_Mercator"
const SRS_PT_HOTINE_OBLIQUE_MERCATOR_TWO_POINT_NATURAL_ORIGIN = "Hotine_Oblique_Mercator_Two_Point_Natural_Origin"
const SRS_PT_LABORDE_OBLIQUE_MERCATOR = "Laborde_Oblique_Mercator"
const SRS_PT_LAMBERT_CONFORMAL_CONIC_1SP = "Lambert_Conformal_Conic_1SP"
const SRS_PT_LAMBERT_CONFORMAL_CONIC_2SP = "Lambert_Conformal_Conic_2SP"
const SRS_PT_LAMBERT_CONFORMAL_CONIC_2SP_BELGIUM = "Lambert_Conformal_Conic_2SP_Belgium"
const SRS_PT_LAMBERT_AZIMUTHAL_EQUAL_AREA = "Lambert_Azimuthal_Equal_Area"
const SRS_PT_MERCATOR_1SP = "Mercator_1SP"
const SRS_PT_MERCATOR_2SP = "Mercator_2SP"
const SRS_PT_MERCATOR_AUXILIARY_SPHERE = "Mercator_Auxiliary_Sphere"
const SRS_PT_MILLER_CYLINDRICAL = "Miller_Cylindrical"
const SRS_PT_MOLLWEIDE = "Mollweide"
const SRS_PT_NEW_ZEALAND_MAP_GRID = "New_Zealand_Map_Grid"
const SRS_PT_OBLIQUE_STEREOGRAPHIC = "Oblique_Stereographic"
const SRS_PT_ORTHOGRAPHIC = "Orthographic"
const SRS_PT_POLAR_STEREOGRAPHIC = "Polar_Stereographic"
const SRS_PT_POLYCONIC = "Polyconic"
const SRS_PT_ROBINSON = "Robinson"
const SRS_PT_SINUSOIDAL = "Sinusoidal"
const SRS_PT_STEREOGRAPHIC = "Stereographic"
const SRS_PT_SWISS_OBLIQUE_CYLINDRICAL = "Swiss_Oblique_Cylindrical"
const SRS_PT_TRANSVERSE_MERCATOR = "Transverse_Mercator"
const SRS_PT_TRANSVERSE_MERCATOR_SOUTH_ORIENTED = "Transverse_Mercator_South_Orientated"
const SRS_PT_TRANSVERSE_MERCATOR_MI_21 = "Transverse_Mercator_MapInfo_21"
const SRS_PT_TRANSVERSE_MERCATOR_MI_22 = "Transverse_Mercator_MapInfo_22"
const SRS_PT_TRANSVERSE_MERCATOR_MI_23 = "Transverse_Mercator_MapInfo_23"
const SRS_PT_TRANSVERSE_MERCATOR_MI_24 = "Transverse_Mercator_MapInfo_24"
const SRS_PT_TRANSVERSE_MERCATOR_MI_25 = "Transverse_Mercator_MapInfo_25"
const SRS_PT_TUNISIA_MINING_GRID = "Tunisia_Mining_Grid"
const SRS_PT_TWO_POINT_EQUIDISTANT = "Two_Point_Equidistant"
const SRS_PT_VANDERGRINTEN = "VanDerGrinten"
const SRS_PT_KROVAK = "Krovak"
const SRS_PT_IMW_POLYCONIC = "International_Map_of_the_World_Polyconic"
const SRS_PT_WAGNER_I = "Wagner_I"
const SRS_PT_WAGNER_II = "Wagner_II"
const SRS_PT_WAGNER_III = "Wagner_III"
const SRS_PT_WAGNER_IV = "Wagner_IV"
const SRS_PT_WAGNER_V = "Wagner_V"
const SRS_PT_WAGNER_VI = "Wagner_VI"
const SRS_PT_WAGNER_VII = "Wagner_VII"
const SRS_PT_QSC = "Quadrilateralized_Spherical_Cube"
const SRS_PT_AITOFF = "Aitoff"
const SRS_PT_WINKEL_I = "Winkel_I"
const SRS_PT_WINKEL_II = "Winkel_II"
const SRS_PT_WINKEL_TRIPEL = "Winkel_Tripel"
const SRS_PT_CRASTER_PARABOLIC = "Craster_Parabolic"
const SRS_PT_LOXIMUTHAL = "Loximuthal"
const SRS_PT_QUARTIC_AUTHALIC = "Quartic_Authalic"
const SRS_PT_SCH = "Spherical_Cross_Track_Height"
const SRS_PP_CENTRAL_MERIDIAN = "central_meridian"
const SRS_PP_SCALE_FACTOR = "scale_factor"
const SRS_PP_STANDARD_PARALLEL_1 = "standard_parallel_1"
const SRS_PP_STANDARD_PARALLEL_2 = "standard_parallel_2"
const SRS_PP_PSEUDO_STD_PARALLEL_1 = "pseudo_standard_parallel_1"
const SRS_PP_LONGITUDE_OF_CENTER = "longitude_of_center"
const SRS_PP_LATITUDE_OF_CENTER = "latitude_of_center"
const SRS_PP_LONGITUDE_OF_ORIGIN = "longitude_of_origin"
const SRS_PP_LATITUDE_OF_ORIGIN = "latitude_of_origin"
const SRS_PP_FALSE_EASTING = "false_easting"
const SRS_PP_FALSE_NORTHING = "false_northing"
const SRS_PP_AZIMUTH = "azimuth"
const SRS_PP_LONGITUDE_OF_POINT_1 = "longitude_of_point_1"
const SRS_PP_LATITUDE_OF_POINT_1 = "latitude_of_point_1"
const SRS_PP_LONGITUDE_OF_POINT_2 = "longitude_of_point_2"
const SRS_PP_LATITUDE_OF_POINT_2 = "latitude_of_point_2"
const SRS_PP_LONGITUDE_OF_POINT_3 = "longitude_of_point_3"
const SRS_PP_LATITUDE_OF_POINT_3 = "latitude_of_point_3"
const SRS_PP_RECTIFIED_GRID_ANGLE = "rectified_grid_angle"
const SRS_PP_LANDSAT_NUMBER = "landsat_number"
const SRS_PP_PATH_NUMBER = "path_number"
const SRS_PP_PERSPECTIVE_POINT_HEIGHT = "perspective_point_height"
const SRS_PP_SATELLITE_HEIGHT = "satellite_height"
const SRS_PP_FIPSZONE = "fipszone"
const SRS_PP_ZONE = "zone"
const SRS_PP_LATITUDE_OF_1ST_POINT = "Latitude_Of_1st_Point"
const SRS_PP_LONGITUDE_OF_1ST_POINT = "Longitude_Of_1st_Point"
const SRS_PP_LATITUDE_OF_2ND_POINT = "Latitude_Of_2nd_Point"
const SRS_PP_LONGITUDE_OF_2ND_POINT = "Longitude_Of_2nd_Point"
const SRS_PP_PEG_POINT_LATITUDE = "peg_point_latitude"
const SRS_PP_PEG_POINT_LONGITUDE = "peg_point_longitude"
const SRS_PP_PEG_POINT_HEADING = "peg_point_heading"
const SRS_PP_PEG_POINT_HEIGHT = "peg_point_height"
const SRS_UL_METER = "Meter"
const SRS_UL_FOOT = "Foot (International)"
const SRS_UL_FOOT_CONV = "0.3048"
const SRS_UL_US_FOOT = "Foot_US"
const SRS_UL_US_FOOT_CONV = "0.3048006096012192"
const SRS_UL_NAUTICAL_MILE = "Nautical Mile"
const SRS_UL_NAUTICAL_MILE_CONV = "1852.0"
const SRS_UL_LINK = "Link"
const SRS_UL_LINK_CONV = "0.20116684023368047"
const SRS_UL_CHAIN = "Chain"
const SRS_UL_CHAIN_CONV = "20.116684023368047"
const SRS_UL_ROD = "Rod"
const SRS_UL_ROD_CONV = "5.02921005842012"
const SRS_UL_LINK_Clarke = "Link_Clarke"
const SRS_UL_LINK_Clarke_CONV = "0.2011661949"
const SRS_UL_KILOMETER = "Kilometer"
const SRS_UL_KILOMETER_CONV = "1000."
const SRS_UL_DECIMETER = "Decimeter"
const SRS_UL_DECIMETER_CONV = "0.1"
const SRS_UL_CENTIMETER = "Centimeter"
const SRS_UL_CENTIMETER_CONV = "0.01"
const SRS_UL_MILLIMETER = "Millimeter"
const SRS_UL_MILLIMETER_CONV = "0.001"
const SRS_UL_INTL_NAUT_MILE = "Nautical_Mile_International"
const SRS_UL_INTL_NAUT_MILE_CONV = "1852.0"
const SRS_UL_INTL_INCH = "Inch_International"
const SRS_UL_INTL_INCH_CONV = "0.0254"
const SRS_UL_INTL_FOOT = "Foot_International"
const SRS_UL_INTL_FOOT_CONV = "0.3048"
const SRS_UL_INTL_YARD = "Yard_International"
const SRS_UL_INTL_YARD_CONV = "0.9144"
const SRS_UL_INTL_STAT_MILE = "Statute_Mile_International"
const SRS_UL_INTL_STAT_MILE_CONV = "1609.344"
const SRS_UL_INTL_FATHOM = "Fathom_International"
const SRS_UL_INTL_FATHOM_CONV = "1.8288"
const SRS_UL_INTL_CHAIN = "Chain_International"
const SRS_UL_INTL_CHAIN_CONV = "20.1168"
const SRS_UL_INTL_LINK = "Link_International"
const SRS_UL_INTL_LINK_CONV = "0.201168"
const SRS_UL_US_INCH = "Inch_US_Surveyor"
const SRS_UL_US_INCH_CONV = "0.025400050800101603"
const SRS_UL_US_YARD = "Yard_US_Surveyor"
const SRS_UL_US_YARD_CONV = "0.914401828803658"
const SRS_UL_US_CHAIN = "Chain_US_Surveyor"
const SRS_UL_US_CHAIN_CONV = "20.11684023368047"
const SRS_UL_US_STAT_MILE = "Statute_Mile_US_Surveyor"
const SRS_UL_US_STAT_MILE_CONV = "1609.347218694437"
const SRS_UL_INDIAN_YARD = "Yard_Indian"
const SRS_UL_INDIAN_YARD_CONV = "0.91439523"
const SRS_UL_INDIAN_FOOT = "Foot_Indian"
const SRS_UL_INDIAN_FOOT_CONV = "0.30479841"
const SRS_UL_INDIAN_CHAIN = "Chain_Indian"
const SRS_UL_INDIAN_CHAIN_CONV = "20.11669506"
const SRS_UA_DEGREE = "degree"
const SRS_UA_DEGREE_CONV = "0.0174532925199433"
const SRS_UA_RADIAN = "radian"
const SRS_PM_GREENWICH = "Greenwich"
const SRS_DN_NAD27 = "North_American_Datum_1927"
const SRS_DN_NAD83 = "North_American_Datum_1983"
const SRS_DN_WGS72 = "WGS_1972"
const SRS_DN_WGS84 = "WGS_1984"
const SRS_WGS84_SEMIMAJOR = 6.378137e6
const SRS_WGS84_INVFLATTENING = 298.257223563
@enum OGRAxisOrientation::UInt32 OAO_Other = 0 OAO_North = 1 OAO_South = 2 OAO_East = 3 OAO_West = 4 OAO_Up = 5 OAO_Down = 6
@enum OGRDatumType::UInt32 ODT_HD_Min = 1000 ODT_HD_Classic = 1001 ODT_HD_Geocentric = 1002 ODT_HD_Max = 1999 ODT_VD_Min = 2000 ODT_VD_Orthometric = 2001 ODT_VD_Ellipsoidal = 2002 ODT_VD_AltitudeBarometric = 2003 ODT_VD_Normal = 2004 ODT_VD_GeoidModelDerived = 2005 ODT_VD_Depth = 2006 ODT_VD_Max = 2999 ODT_LD_Min = 10000 ODT_LD_Max = 32767
const ODT_HD_Other = ODT_HD_Min
const ODT_VD_Other = ODT_VD_Min


text:
#=
This file is auto-generated. Do not edit.
=#
"""
    mutable struct OneDOneQMachine <: Machine
        R::Float64
        Xd::Float64
        Xq::Float64
        Xd_p::Float64
        Xq_p::Float64
        Td0_p::Float64
        Tq0_p::Float64
        ext::Dict{String, Any}
        states::Vector{Symbol}
        n_states::Int
        internal::InfrastructureSystemsInternal
    end

Parameters of 4-states synchronous machine: Simplified Marconato model
 The derivative of stator fluxes (ψd and ψq) is neglected and ωψd = ψd and
 ωψq = ψq is assumed (i.e. ω=1.0). This is standard when
 transmission network dynamics is neglected.

# Arguments
- `R::Float64`: Resistance after EMF in machine per unit, validation range: `(0, nothing)`
- `Xd::Float64`: Reactance after EMF in d-axis per unit, validation range: `(0, nothing)`
- `Xq::Float64`: Reactance after EMF in q-axis per unit, validation range: `(0, nothing)`
- `Xd_p::Float64`: Transient reactance after EMF in d-axis per unit, validation range: `(0, nothing)`
- `Xq_p::Float64`: Transient reactance after EMF in q-axis per unit, validation range: `(0, nothing)`
- `Td0_p::Float64`: Time constant of transient d-axis voltage, validation range: `(0, nothing)`
- `Tq0_p::Float64`: Time constant of transient q-axis voltage, validation range: `(0, nothing)`
- `ext::Dict{String, Any}`
- `states::Vector{Symbol}`: The states are:
	eq_p: q-axis transient voltage,
	ed_p: d-axis transient voltage
- `n_states::Int`: OneDOneQMachine has 2 states
- `internal::InfrastructureSystemsInternal`: power system internal reference, do not modify
"""
mutable struct OneDOneQMachine <: Machine
    "Resistance after EMF in machine per unit"
    R::Float64
    "Reactance after EMF in d-axis per unit"
    Xd::Float64
    "Reactance after EMF in q-axis per unit"
    Xq::Float64
    "Transient reactance after EMF in d-axis per unit"
    Xd_p::Float64
    "Transient reactance after EMF in q-axis per unit"
    Xq_p::Float64
    "Time constant of transient d-axis voltage"
    Td0_p::Float64
    "Time constant of transient q-axis voltage"
    Tq0_p::Float64
    ext::Dict{String, Any}
    "The states are:
	eq_p: q-axis transient voltage,
	ed_p: d-axis transient voltage"
    states::Vector{Symbol}
    "OneDOneQMachine has 2 states"
    n_states::Int
    "power system internal reference, do not modify"
    internal::InfrastructureSystemsInternal
end

function OneDOneQMachine(R, Xd, Xq, Xd_p, Xq_p, Td0_p, Tq0_p, ext=Dict{String, Any}(), )
    OneDOneQMachine(R, Xd, Xq, Xd_p, Xq_p, Td0_p, Tq0_p, ext, [:eq_p, :ed_p], 2, InfrastructureSystemsInternal(), )
end

function OneDOneQMachine(; R, Xd, Xq, Xd_p, Xq_p, Td0_p, Tq0_p, ext=Dict{String, Any}(), states=[:eq_p, :ed_p], n_states=2, internal=InfrastructureSystemsInternal(), )
    OneDOneQMachine(R, Xd, Xq, Xd_p, Xq_p, Td0_p, Tq0_p, ext, states, n_states, internal, )
end

# Constructor for demo purposes; non-functional.
function OneDOneQMachine(::Nothing)
    OneDOneQMachine(;
        R=0,
        Xd=0,
        Xq=0,
        Xd_p=0,
        Xq_p=0,
        Td0_p=0,
        Tq0_p=0,
        ext=Dict{String, Any}(),
    )
end

"""Get [`OneDOneQMachine`](@ref) `R`."""
get_R(value::OneDOneQMachine) = value.R
"""Get [`OneDOneQMachine`](@ref) `Xd`."""
get_Xd(value::OneDOneQMachine) = value.Xd
"""Get [`OneDOneQMachine`](@ref) `Xq`."""
get_Xq(value::OneDOneQMachine) = value.Xq
"""Get [`OneDOneQMachine`](@ref) `Xd_p`."""
get_Xd_p(value::OneDOneQMachine) = value.Xd_p
"""Get [`OneDOneQMachine`](@ref) `Xq_p`."""
get_Xq_p(value::OneDOneQMachine) = value.Xq_p
"""Get [`OneDOneQMachine`](@ref) `Td0_p`."""
get_Td0_p(value::OneDOneQMachine) = value.Td0_p
"""Get [`OneDOneQMachine`](@ref) `Tq0_p`."""
get_Tq0_p(value::OneDOneQMachine) = value.Tq0_p
"""Get [`OneDOneQMachine`](@ref) `ext`."""
get_ext(value::OneDOneQMachine) = value.ext
"""Get [`OneDOneQMachine`](@ref) `states`."""
get_states(value::OneDOneQMachine) = value.states
"""Get [`OneDOneQMachine`](@ref) `n_states`."""
get_n_states(value::OneDOneQMachine) = value.n_states
"""Get [`OneDOneQMachine`](@ref) `internal`."""
get_internal(value::OneDOneQMachine) = value.internal

"""Set [`OneDOneQMachine`](@ref) `R`."""
set_R!(value::OneDOneQMachine, val) = value.R = val
"""Set [`OneDOneQMachine`](@ref) `Xd`."""
set_Xd!(value::OneDOneQMachine, val) = value.Xd = val
"""Set [`OneDOneQMachine`](@ref) `Xq`."""
set_Xq!(value::OneDOneQMachine, val) = value.Xq = val
"""Set [`OneDOneQMachine`](@ref) `Xd_p`."""
set_Xd_p!(value::OneDOneQMachine, val) = value.Xd_p = val
"""Set [`OneDOneQMachine`](@ref) `Xq_p`."""
set_Xq_p!(value::OneDOneQMachine, val) = value.Xq_p = val
"""Set [`OneDOneQMachine`](@ref) `Td0_p`."""
set_Td0_p!(value::OneDOneQMachine, val) = value.Td0_p = val
"""Set [`OneDOneQMachine`](@ref) `Tq0_p`."""
set_Tq0_p!(value::OneDOneQMachine, val) = value.Tq0_p = val
"""Set [`OneDOneQMachine`](@ref) `ext`."""
set_ext!(value::OneDOneQMachine, val) = value.ext = val


text:
using
  FourierFlows,
  FourierFlows.Diffusion,
  FFTW,
  LinearAlgebra,
  Printf,
  JLD2,
  Test

using FourierFlows: parsevalsum2

using LinearAlgebra: mul!, ldiv!, norm

# the devices on which tests will run
devices = (CPU(),)
@has_cuda devices = (CPU(), GPU())
@has_cuda using CuArrays

const rtol_fft = 1e-12
const rtol_output = 1e-12
const rtol_utils = 1e-13
const rtol_timesteppers = 1e-12

const steppers = [
  "ForwardEuler",
  "RK4",
  "ETDRK4",
  "AB3",
  "FilteredForwardEuler",
  "FilteredRK4",
  "FilteredETDRK4",
  "FilteredAB3"
]

# Run tests
include("createffttestfunctions.jl")

for dev in devices

  println("testing on "*string(typeof(dev)))

  @time @testset "Grid tests" begin
    include("test_grid.jl")

    nx, Lx =  6, 2π
    ny, Ly =  8, 4π
    nz, Lz = 10, 3.0

    # Test 1D grid
    g₁ = OneDGrid(dev, nx, Lx)
    @test testnx(g₁, nx)
    @test testdx(g₁)
    @test testdk(g₁)
    @test testx(g₁)
    @test testk(g₁)
    @test testkr(g₁)
    @test testdealias(g₁)
    @test testmakefilter(dev, g₁)

    # Test 2D rectangular grid
    g₂ = TwoDGrid(dev, nx, Lx, ny, Ly)
    @test testnx(g₂, nx)
    @test testny(g₂, ny)
    @test testdx(g₂)
    @test testdy(g₂)
    @test testdk(g₂)
    @test testdl(g₂)
    @test testx(g₂)
    @test testy(g₂)
    @test testk(g₂)
    @test testkr(g₂)
    @test testl(g₂)
    @test testgridpoints(g₂)
    @test testdealias(g₂)
    @test testmakefilter(dev, g₂)

    # Test 3D parallelogram grid
    g₃ = ThreeDGrid(dev, nx, Lx, ny, Ly, nz, Lz)
    @test testnx(g₃, nx)
    @test testny(g₃, ny)
    @test testnz(g₃, nz)
    @test testdx(g₃)
    @test testdy(g₃)
    @test testdz(g₃)
    @test testdk(g₃)
    @test testdl(g₃)
    @test testdm(g₃)
    @test testx(g₃)
    @test testy(g₃)
    @test testz(g₃)
    @test testk(g₃)
    @test testkr(g₃)
    @test testl(g₃)
    @test testm(g₃)
    @test testgridpoints(g₃)
    @test testdealias(g₃)
    @test testmakefilter(dev, g₃)

    # Test typed grids
    T = Float32
    @test testtypedonedgrid(dev, nx, Lx; T=T)
    @test testtypedtwodgrid(dev, nx, Lx, ny, Ly; T=T)
    @test testtypedthreedgrid(dev, nx, Lx, ny, Ly, nz, Lz; T=T)
  end

  @time @testset "FFT tests" begin
    include("test_fft.jl")

    # Test 1D grid
    nx = 32             # number of points
    Lx = 2π             # Domain width
    g1 = OneDGrid(dev, nx, Lx)

    @test test_fft_cosmx(g1)
    @test test_rfft_cosmx(g1)
    @test test_rfft_mul_cosmx(g1)

    # Test 2D rectangular grid
    nx, ny = 32, 64     # number of points
    Lx, Ly = 2π, 3π     # Domain width
    g2 = TwoDGrid(dev, nx, Lx, ny, Ly)

    @test test_fft_cosmxcosny(g2)
    @test test_rfft_cosmxcosny(g2)
    @test test_rfft_mul_cosmxcosny(g2)
    @test test_fft_sinmxny(g2)
    @test test_rfft_sinmxny(g2)
    @test test_rfft_mul_sinmxny(g2)

    # Test 3D parallelogram grid
    nx, ny, nz = 32, 30, 16     # number of points
    Lx, Ly, Lz = 2π, 3π, 4.0    # Domain width
    g3 = ThreeDGrid(dev, nx, Lx, ny, Ly, nz, Lz)

    @test test_fft_cosmxcosny(g3)
    @test test_rfft_cosmxcosny(g3)
    @test test_rfft_mul_cosmxcosny(g3)
    @test test_fft_sinmxny(g3)
    @test test_rfft_sinmxny(g3)
    @test test_rfft_mul_sinmxny(g3)
  end

  @time @testset "IFFT tests" begin
    include("test_ifft.jl")

    # Test 1D grid
    nx = 32             # number of points
    Lx = 2π             # Domain width
    g1 = OneDGrid(dev, nx, Lx)

    @test test_ifft_cosmx(g1)
    @test test_irfft_cosmx(g1)
    @test test_irfft_mul_cosmx(g1)

    # Test 2D rectangular grid
    nx, ny = 32, 64     # number of points
    Lx, Ly = 2π, 3π     # Domain width
    g2 = TwoDGrid(dev, nx, Lx, ny, Ly)

    @test test_ifft_cosmxcosny(g2)
    @test test_irfft_cosmxcosny(g2)
    @test test_irfft_mul_cosmxcosny(g2)
    @test test_ifft_sinmxny(g2)
    @test test_irfft_sinmxny(g2)
    @test test_irfft_mul_sinmxny(g2)

    # Test 3D parallelogram grid
    nx, ny, nz = 32, 30, 16     # number of points
    Lx, Ly, Lz = 2π, 3π, 4.0    # Domain width
    g3 = ThreeDGrid(dev, nx, Lx, ny, Ly, nz, Lz)

    @test test_ifft_cosmxcosny(g3)
    @test test_irfft_cosmxcosny(g3)
    @test test_irfft_mul_cosmxcosny(g3)
    @test test_ifft_sinmxny(g3)
    @test test_irfft_sinmxny(g3)
    @test test_irfft_mul_sinmxny(g3)
  end

  @time @testset "Timestepper tests" begin
    include("test_timesteppers.jl")
    for stepper in steppers
      @test constantdiffusiontest(stepper, dev=dev)
      @test varyingdiffusiontest(stepper, dev=dev)
    end
  end

  @time @testset "Utils tests" begin
    include("test_utils.jl")

    @test test_fltype()
    @test test_cxtype()
    @test test_innereltype()
    @test test_superzeros()
    @test test_supertuplezeros()
    @test test_supersize()
    @test test_arraytype(dev)
    @test test_arraytypeTdim(dev, Float32, 2)

    # Test on a rectangular grid
    nx, ny = 64, 128   # number of points
    Lx, Ly = 2π, 3π    # Domain width
    g = TwoDGrid(dev, nx, Lx, ny, Ly)
    x, y = gridpoints(g)
    k0, l0 = 2π/Lx, 2π/Ly

    # Real and complex-valued functions
    σ = 0.5
    f1 = @. exp(-(x^2 + y^2)/(2σ^2))
    f2 = @. (cos(2k0*x + 3l0*y^2) + im*sin(2k0*x + 3l0*y^2)) * (exp(-(x^2 + y^2)/(2σ^2)) + 2im*exp(-(x^2 + y^2)/(5σ^2)))

    # Sine/Exp waves
    k1, l1 = 2*k0,  6*l0
    k2, l2 = 3*k0, -3*l0
    sinkl1 = @. sin(k1*x + l1*y)
    sinkl2 = @. sin(k2*x + l2*y)
    expkl1 = @. cos(k1*x + l1*y) + im*sin(k1*x + l1*y)
    expkl2 = @. cos(k2*x + l2*y) + im*sin(k2*x + l2*y)

    # Analytical expression for the Jacobian of sin1 and sin2 and of exp1 and exp2
    Jsinkl1sinkl2 = @. (k1*l2-k2*l1)*cos(k1*x + l1*y)*cos(k2*x + l2*y)
    Jexpkl1expkl2 = @. (k2*l1-k1*l2)*(cos((k1+k2)*x + (l1+l2)*y)+im*sin((k1+k2)*x + (l1+l2)*y))

    @test test_parsevalsum(f1, g; realvalued=true)   # Real valued f with rfft
    @test test_parsevalsum(f1, g; realvalued=false)  # Real valued f with fft
    @test test_parsevalsum(f2, g; realvalued=false)  # Complex valued f with fft
    @test test_parsevalsum2(f1, g; realvalued=true)  # Real valued f with rfft
    @test test_parsevalsum2(f1, g; realvalued=false) # Real valued f with fft
    @test test_parsevalsum2(f2, g; realvalued=false) # Complex valued f with fft

    @test test_jacobian(sinkl1, sinkl1, 0*sinkl1, g)      # Test J(a, a) = 0
    @test test_jacobian(sinkl1, sinkl2, Jsinkl1sinkl2, g) # Test J(sin1, sin2) = Jsin1sin2
    @test test_jacobian(expkl1, expkl2, Jexpkl1expkl2, g) # Test J(exp1, exp2) = Jexp1exps2

    @test test_zeros()
    @test test_domainaverage(dev, 32)
    @test test_varsexpression_fields(g)
    @test test_varsexpression_specs(g)

    # Radial spectrum tests. Note that ahρ = ∫ ah ρ dθ.
    n = 128; δ = n/10                 # Parameters
    ahkl(k, l) = exp(-(k^2+l^2)/2δ^2) # a  = exp(-ρ²/2δ²)
        ahρ(ρ) = 2π*ρ*exp(-ρ^2/2δ^2)  # aᵣ = 2π ρ exp(-ρ²/2δ²)
    @test test_radialspectrum(dev, n, ahkl, ahρ)

    ahkl(k, l) = exp(-(k^2+l^2)/2δ^2) * k^2/(k^2+l^2) # a  = exp(-ρ²/2δ²)*cos(θ)²
        ahρ(ρ) = π*ρ*exp(-ρ^2/2δ^2)                   # aᵣ = π ρ exp(-ρ²/2δ²)
    @test test_radialspectrum(dev, n, ahkl, ahρ)
  end

  @time @testset "Diagnostics tests" begin
    include("test_diagnostics.jl")

    @test test_diagnosticsteps(dev, freq=1)
    @test test_diagnosticsteps(dev, freq=2)
    @test test_diagnosticsteps(dev, nsteps=100, freq=9, ndata=20)
    # @test test_basicdiagnostics(dev)
    @test test_scalardiagnostics(dev, freq=1)
    @test test_scalardiagnostics(dev, freq=2)
  end

  @time @testset "Output tests" begin
    include("test_output.jl")

    @test test_withoutjld2()
    @test test_uniquepath()
    @test test_outputconstructor(dev)
    @test test_getindex(dev)
    @test test_saveproblem_saveoutput(dev)
    @test test_saveproblemTwoDGrid(dev)
    @test test_savediagnostic(dev)
  end

end # end loop over devices


text:
############################################################################################
# get parameter identifier
############################################################################################

identifier(sem_fit::SemFit) = identifier(sem_fit.model)
identifier(model::AbstractSemSingle) = identifier(model.imply)
identifier(model::SemEnsemble) = model.identifier

############################################################################################
# construct identifier
############################################################################################

identifier(ram_matrices::RAMMatrices) =
    Dict{Symbol, Int64}(ram_matrices.parameters .=> 1:length(ram_matrices.parameters))
function identifier(partable::ParameterTable)
    _, _, identifier = get_par_npar_identifier(partable)
    return identifier
end

############################################################################################
# get indices of a Vector of parameter labels
############################################################################################

get_identifier_indices(parameters, identifier::Dict{Symbol, Int}) =
    [identifier[par] for par in parameters]

get_identifier_indices(parameters, obj::Union{SemFit, AbstractSemSingle, SemEnsemble, SemImply}) =
    get_identifier_indices(parameters, identifier(obj))

function get_identifier_indices(parameters, obj::Union{ParameterTable, RAMMatrices})
    @warn "You are trying to find parameter indices from a ParameterTable or RAMMatrices object. \n
           If your model contains user-defined types, this may lead to wrong results. \n
           To be on the safe side, try to reference parameters by labels or query the indices from
           the constructed model (`get_identifier_indices(parameters, model)`)." maxlog=1
    return get_identifier_indices(parameters, identifier(obj))
end

############################################################################################
# documentation
############################################################################################
"""
    get_identifier_indices(parameters, model)

Returns the indices of `parameters`.

# Arguments
- `parameters::Vector{Symbol}`: parameter labels
- `model`: either a SEM or a fitted SEM

# Examples
```julia
parameter_indices = get_identifier_indices([:λ₁, λ₂], my_fitted_sem)

values = solution(my_fitted_sem)[parameter_indices]
```
"""
function get_identifier_indices end


text:
export jaccard_distance, dice_distance, union_intersection, intersection_distance

function union_intersection(a::T, b::T) where {T <: AbstractVector}
    len_a::Int = length(a)
    len_b::Int = length(b)
    ia::Int = 1
    ib::Int = 1
    intersection_size::Int = 0
    c::Int = 0
    @inbounds while ia <= len_a && ib <= len_b
        c = cmp(a[ia], b[ib])
        if c == 0
            ia += 1
            ib += 1
            intersection_size += 1
        elseif c < 0
            ia += 1
        else
            ib += 1
        end
    end

    return len_a + len_b - intersection_size, intersection_size
end

"""
sim_jaccard computes the Jaccard's coefficient between two sets
represented as sorted arrays.
Note: 1.0 - sim_jaccard computes the Jaccard's distance
"""

function jaccard_distance(a::T, b::T)::Float64 where {T <: AbstractVector}
    u, i = union_intersection(a, b)
    return 1.0 - i / u
end

function dice_distance(a::T, b::T)::Float64 where {T <: AbstractVector}
    u, i = union_intersection(a, b)
    return 1.0 - 2 * i / (length(a) + length(b))
end

mutable struct IntersectionDistance
    calls::Int
    IntersectionDistance() = new(0)
end

function intersection_distance(a::T, b::T)::Float64 where {T <: AbstractVector}
    u, i = union_intersection(a, b)

    return 1.0 - i / min(length(a), length(b))
end


text:
using Dates
using DataStructures

using .Networks

mutable struct PSCOPFContext <: AbstractContext
    network::Networks.Network
    target_timepoints::Vector{Dates.DateTime}
    management_mode::ManagementMode

    generators_initial_state::SortedDict{String,GeneratorState}

    #uncertainties
    uncertainties::Uncertainties

    #AssessmentUncertainties
    assessment_uncertainties

    horizon_timepoints::Vector{Dates.DateTime}

    market_schedule::Schedule
    tso_schedule::Schedule

    tso_actions::TSOActions

    #flows: branch, ts, s
    market_flows::SortedDict{Tuple{String, DateTime, String}, Float64}
    tso_flows::SortedDict{Tuple{String, DateTime, String}, Float64}

    out_dir::Union{String,Nothing}
end

function PSCOPFContext(network::Networks.Network, target_timepoints::Vector{Dates.DateTime},
                    management_mode::ManagementMode,
                    generators_initial_state::SortedDict{String,GeneratorState}=SortedDict{String,GeneratorState}(),
                    uncertainties::Uncertainties=Uncertainties(),
                    assessment_uncertainties=nothing,
                    out_dir=nothing
                    )
    market_schedule = Schedule(Utilitary(), Dates.DateTime(0))
    init!(market_schedule, network, target_timepoints, get_scenarios(uncertainties))
    tso_schedule = Schedule(Utilitary(), Dates.DateTime(0))
    init!(tso_schedule, network, target_timepoints, get_scenarios(uncertainties))
    return PSCOPFContext(network, target_timepoints, management_mode,
                        generators_initial_state,
                        uncertainties, assessment_uncertainties,
                        Vector{Dates.DateTime}(),
                        market_schedule,
                        tso_schedule,
                        TSOActions(),
                        SortedDict{Tuple{String, DateTime, String}, Float64}(),
                        SortedDict{Tuple{String, DateTime, String}, Float64}(),
                        out_dir)
end

function get_network(context::PSCOPFContext)
    return context.network
end

function get_target_timepoints(context::PSCOPFContext)
    return context.target_timepoints
end

function get_horizon_timepoints(context::PSCOPFContext)
    return context.horizon_timepoints
end

function set_horizon_timepoints(context::PSCOPFContext, horizon_timepoints::Vector{Dates.DateTime})
    context.horizon_timepoints = horizon_timepoints
end

function get_management_mode(context::PSCOPFContext)
    return context.management_mode
end

function get_generators_initial_state(context::PSCOPFContext)
    return context.generators_initial_state
end

function get_uncertainties(context::PSCOPFContext)::Uncertainties
    return context.uncertainties
end

function get_uncertainties(context::PSCOPFContext, ech::Dates.DateTime)::UncertaintiesAtEch
    return get_uncertainties(context)[ech]
end

function get_scenarios(context::PSCOPFContext, ech::Dates.DateTime)::Vector{String}
    uncertainties_at_ech = get_uncertainties(context, ech)
    if isnothing(uncertainties_at_ech)
        return Vector{String}()
    else
        return get_scenarios(uncertainties_at_ech)
    end
end

function get_scenarios(context::PSCOPFContext)::Vector{String}
    uncertainties = get_uncertainties(context)
    if isnothing(uncertainties)
        return Vector{String}()
    else
        return get_scenarios(uncertainties)
    end
end

function get_assessment_uncertainties(context::PSCOPFContext)
    return context.assessment_uncertainties
end

function set_current_ech!(context_p::PSCOPFContext, ech::Dates.DateTime)
    context_p.current_ech = ech
end

function get_tso_actions(context::PSCOPFContext)
    return context.tso_actions
end

function get_tso_schedule(context_p::PSCOPFContext)
    return context_p.tso_schedule
end

function get_market_schedule(context_p::PSCOPFContext)
    return context_p.market_schedule
end

function get_market_flows(context::PSCOPFContext)
    return context.market_flows
end

function get_tso_flows(context::PSCOPFContext)
    return context.tso_flows
end

function get_limitables_ids(context_p::PSCOPFContext)
    limitables = Networks.get_generators_of_type(get_network(context_p), Networks.LIMITABLE)
    limitables_ids = map(lim_gen->Networks.get_id(lim_gen), limitables)
    return limitables_ids
end

function get_initial_state(initial_states::SortedDict{String,GeneratorState}, generator::Generator)
    if Networks.get_p_min(generator) < 1e-09
        return ON
    else
        return initial_states[Networks.get_id(generator)]
    end
end

function update_market_flows!(context::PSCOPFContext)
    flows = compute_flows(context, get_market_schedule(context))
    context.market_flows = flows
end

function update_tso_flows!(context::PSCOPFContext)
    flows = compute_flows(context, get_tso_schedule(context))
    context.tso_flows = flows
end

max_stars_repo_name: klorel/pscopf

max_issues_repo_name: klorel/pscopf

max_forks_repo_name: klorel/pscopf


text:
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#
# Description
# ==============================================================================
#
#   Tests of compact type strings.
#
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #

@testset "Compact types" begin
    # Dictionary
    # --------------------------------------------------------------------------

    data = Dict(:a => Int64(1), missing => Int64(2), :c => missing)

    expected = """
┌─────────┬─────────┐
│    Keys │  Values │
│ Symbol? │  Int64? │
├─────────┼─────────┤
│       a │       1 │
│       c │ missing │
│ missing │       2 │
└─────────┴─────────┘
"""
    result = pretty_table(String, data, sortkeys = true)
    @test result == expected

    # Tables.jl API
    # --------------------------------------------------------------------------

    table = (
        a = [missing, Int64(1), Int64(2), Int64(3)],
        b = [nothing, Int64(1), missing, Int64(1)]
    )

    expected = """
┌─────────┬────────────────────┐
│       a │                  b │
│  Int64? │ U{Nothing, Int64}? │
├─────────┼────────────────────┤
│ missing │            nothing │
│       1 │                  1 │
│       2 │            missing │
│       3 │                  1 │
└─────────┴────────────────────┘
"""
    result = pretty_table(String, table)
    @test result == expected
end


text:
using Documenter, GigaFV
#https://juliadocs.github.io/Documenter.jl/stable/man/guide/index.html
# Run this by julia --project=. --color=yes docs/make.jl in the root
# of this repository.
makedocs(
    modules = [GigaFV],
    sitename="GigaFV Documentation")


text:
export solve_SCPshooting!, solve_SCP!

function solve_SCPshooting!(TOS::TrajectoryOptimizationSolution, TOP::TrajectoryOptimizationProblem, solve_method!, init_method, solver="Mosek"; max_iter=30, kwarg...)
	robot = TOP.PD.robot
	model = TOP.PD.model

	SCPP = SCPProblem(TOP)

	# Create initialization trajectory
	traj_init = init_method(TOP)

	# Run SCP with initialization trajectory
	TOS.SCPS = SCPSolution(SCPP, traj_init)
	SCPS = TOS.SCPS
	SP = ShootingProblem(TOP, SCPS)
	TOS.SS = ShootingSolution(SP, deepcopy(traj_init))
	SS = TOS.SS

	solve_method!(SCPS, SCPP, solver, 1; kwarg...)
	push!(SS.J_true, SCPS.J_true[1])

	# Until shooting method succeeds or maximum SCP iterations is reached
	ss_sol = nothing
	while (!SCPS.converged && SCPS.iterations < max_iter)
		# Attempt shooting method
		SP = ShootingProblem(TOP, SCPS)
		ss_sol = solve!(SS, SP)

		# If successful shooting runs have converged, exit
		conv_iter_spread = 2
		if SCPS.iterations > conv_iter_spread && sum(SS.convergence_measure[end-conv_iter_spread+1:end]) <= SCPS.param.convergence_threshold
			SS.converged = true
			# TODO: Check inequality constraints
			copy!(TOS.traj, SS.traj)
			TOS.total_time = SCPS.total_time + sum(SS.iter_elapsed_times)
			return
		end
		# Run another iteration of SCP
		solve_method!(SCPS, SCPP, solver, 1; kwarg...)
	end

	copy!(TOS.traj, SCPS.traj)
	TOS.total_time = SCPS.total_time + sum(SS.iter_elapsed_times)
end

function solve_SCP!(TOS::TrajectoryOptimizationSolution, TOP::TrajectoryOptimizationProblem, solve_method!, init_method, solver="Mosek"; max_iter=30, force=false, kwarg...)
	robot = TOP.PD.robot
	model = TOP.PD.model

	SCPP = SCPProblem(TOP)

	# Create initialization trajectory
	traj_init = init_method(TOP)

	# Run SCP with initialization trajectory
	SCPS = SCPSolution(SCPP, traj_init)
	TOS.traj, TOS.SCPS = SCPS.traj, SCPS
	solve_method!(SCPS, SCPP, solver, max_iter, force; kwarg...)
end

function solve_SCP!(TOS::TrajectoryOptimizationSolution, TOP::TrajectoryOptimizationProblem, solve_method!, traj_init::Trajectory, solver="Mosek"; max_iter=30, force=false, kwarg...)
	robot = TOP.PD.robot
	model = TOP.PD.model

	SCPP = SCPProblem(TOP)

	# Run SCP with initialization trajectory
	SCPS = SCPSolution(SCPP, traj_init)
	TOS.traj, TOS.SCPS = SCPS.traj, SCPS
	solve_method!(SCPS, SCPP, solver, max_iter, force; kwarg...)
end

function convergence_metric(traj::Trajectory, traj_prev::Trajectory, OAP::OptAlgorithmProblem)
  # normalized maximum relative error between iterations
  max_num, max_den = -Inf, -Inf
  for k in 1:OAP.N
    val = norm(traj.X[:,k]-traj_prev.X[:,k])
    max_num = val > max_num ? val : max_num

    val = norm(traj.X[:,k])
    max_den = val > max_den ? val : max_den
  end
  return max_num/max_den
end


text:
function glGenOne(glGenFn)
    id::Ptr{GLuint} = GLuint[0]
    glGenFn(1, id)
    glCheckError("generating a buffer, array, or texture")
    unsafe_load(id)
end
glGenBuffer() = glGenOne(glGenBuffers)
glGenVertexArray() =  glGenOne(glGenVertexArrays)
glGenTexture() =  glGenOne(glGenTextures)

function getInfoLog(obj::GLuint)
	# Return the info log for obj, whether it be a shader or a program.
	isShader = glIsShader(obj)
	getiv   = isShader == GL_TRUE ? glGetShaderiv      : glGetProgramiv
	getInfo = isShader == GL_TRUE ? glGetShaderInfoLog : glGetProgramInfoLog

	# Get the maximum possible length for the descriptive error message
	int::Ptr{GLint} = GLint[0]
	getiv(obj, GL_INFO_LOG_LENGTH, int)
	maxlength = unsafe_load(int)

	# TODO: Create a macro that turns the following into the above:
	# maxlength = @glPointer getiv(obj, GL_INFO_LOG_LENGTH, GLint)

	# Return the text of the message if there is any
	if maxlength > 0
		buffer = zeros(GLchar, maxlength)
		sizei::Ptr{GLsizei} = GLsizei[0]
		getInfo(obj, maxlength, sizei, buffer)
		length = unsafe_load(sizei)
		bytestring(pointer(buffer), length)
	else
		""
	end
end

function validateShader(shader)
	success::Ptr{GLint} = GLint[0]
	glGetShaderiv(shader, GL_COMPILE_STATUS, success)
	unsafe_load(success) == GL_TRUE
end

function glErrorMessage()
	# Return a string representing the current OpenGL error flag, or the empty string if there's no error.
	err = glGetError()
	err == GL_NO_ERROR ? "" :
	err == GL_INVALID_ENUM ? "GL_INVALID_ENUM: An unacceptable value is specified for an enumerated argument. The offending command is ignored and has no other side effect than to set the error flag." :
	err == GL_INVALID_VALUE ? "GL_INVALID_VALUE: A numeric argument is out of range. The offending command is ignored and has no other side effect than to set the error flag." :
	err == GL_INVALID_OPERATION ? "GL_INVALID_OPERATION: The specified operation is not allowed in the current state. The offending command is ignored and has no other side effect than to set the error flag." :
	err == GL_INVALID_FRAMEBUFFER_OPERATION ? "GL_INVALID_FRAMEBUFFER_OPERATION: The framebuffer object is not complete. The offending command is ignored and has no other side effect than to set the error flag." :
	err == GL_OUT_OF_MEMORY ? "GL_OUT_OF_MEMORY: There is not enough memory left to execute the command. The state of the GL is undefined, except for the state of the error flags, after this error is recorded." : "Unknown OpenGL error with error code $err."
end

function glCheckError(actionName="")
	message = glErrorMessage()
	if length(message) > 0
		if length(actionName) > 0
			error("Error ", actionName, ": ", message)
		else
			error("Error: ", message)
		end
	end
end

function createShader(source, typ)
	# Create the shader
	shader = glCreateShader(typ)::GLuint
	if shader == 0
		error("Error creating shader: ", glErrorMessage())
	end

	# Compile the shader
	glShaderSource(shader, 1, convert(Ptr{Uint8}, pointer([convert(Ptr{GLchar}, pointer(source))])), C_NULL)
	glCompileShader(shader)

	# Check for errors
	!validateShader(shader) && error("Shader creation error: ", getInfoLog(shader))
	shader
end

function createShaderProgram(f, vertexShader, fragmentShader)
	# Create, link then return a shader program for the given shaders.

	# Create the shader program
	prog = glCreateProgram()
	if prog == 0
		error("Error creating shader program: ", glErrorMessage())
	end

	# Attach the vertex shader
	glAttachShader(prog, vertexShader)
	glCheckError("attaching vertex shader")

	# Attach the fragment shader
	glAttachShader(prog, fragmentShader)
	glCheckError("attaching fragment shader")

	f(prog)

	# Finally, link the program and check for errors.
	glLinkProgram(prog)
	status::Ptr{GLint} = GLint[0]
	glGetProgramiv(prog, GL_LINK_STATUS, status)
	if unsafe_load(status) == GL_FALSE then
		glDeleteProgram(prog)
		error("Error linking shader: ", glGetInfoLog(prog))
	end

	prog
end
createShaderProgram(vertexShader, fragmentShader) = createShaderProgram(prog->0, vertexShader, fragmentShader)

function printGLInfo()
    println("GLSL version: ",    bytestring(glGetString(GL_SHADING_LANGUAGE_VERSION)))
    println("OpenGL version: ",  bytestring(glGetString(GL_VERSION)))
    println("OpenGL vendor: ",   bytestring(glGetString(GL_VENDOR)))
    println("OpenGL renderer: ", bytestring(glGetString(GL_RENDERER)))
end

max_stars_repo_name: yurivish/Julia-Playground

max_issues_repo_name: yurivish/Julia-Playground

max_forks_repo_name: yurivish/Julia-Playground


text:
# Model
include_model("hopper_impedance")
model_ft = free_time_model(model)

# reference impedance
k0 = 2.5

# Horizon
T = 51

# Time step
tf = 1.0
h = tf / (T - 1)

# Bounds
_uu = Inf * ones(model_ft.m)
_uu[model_ft.idx_u] = model_ft.uU
_uu[3] = 5.0
_uu[end] = 1.0 * h
_ul = zeros(model_ft.m)
_ul[model_ft.idx_u] = model_ft.uL
_ul[end] = 1.0 * h
_ul[3] = 0.0
ul, uu = control_bounds(model_ft, T, _ul, _uu)

# Initial and final states
z_h = 0.5
q1 = [0.0, 0.5 + z_h, 0.0, 0.5 * model.r0]

xl, xu = state_bounds(model_ft, T,
		[model_ft.qL; model_ft.qL],
		[model_ft.qU; model_ft.qU],
        x1 = [q1; Inf * ones(model.nq)])

# Objective
include_objective(["velocity", "nonlinear_stage"])
obj_tracking = quadratic_time_tracking_objective(
    [Diagonal(zeros(model_ft.n)) for t = 1:T],
    [Diagonal([1.0, 1.0, 1.0, 10.0, 10.0, zeros(model_ft.m - model_ft.nu - 1 - 2)..., 0.0]) for t = 1:T-1],
    [zeros(model_ft.n) for t = 1:T],
    [[0.0, 0.0, k0, zeros(model_ft.m - model_ft.nu - 1)..., 0.0] for t = 1:T],
    1.0)
obj_contact_penalty = PenaltyObjective(1.0e5, model_ft.m - 1)
obj_velocity = velocity_objective(
    [Diagonal(1.0 * ones(model_ft.nq)) for t = 1:T-1],
    model_ft.nq,
    h = h,
    idx_angle = collect([3]))

function get_q⁺(x)
	view(x, model.nq .+ (1:model.nq))
end

function l_stage_body_foot(x, u, t)
	k = kinematics(model, get_q⁺(x))
	return 10.0 * (k[1] - k[3])^2.0
end
l_terminal_body_foot(x) = 0.0
obj_tf = nonlinear_stage_objective(l_stage_body_foot, l_terminal_body_foot)

obj = MultiObjective([obj_tracking, obj_contact_penalty, obj_velocity])#, obj_tf])

# Constraints
include_constraints(["free_time", "contact", "loop"])
con_free_time = free_time_constraints(T)
con_contact = contact_constraints(model_ft, T)
con_loop = loop_constraints(model, (1:model_ft.n), 1, T)
con = multiple_constraints([con_free_time, con_contact, con_loop])

# Problem
prob = trajectory_optimization_problem(model_ft,
               obj,
               T,
               xl = xl,
               xu = xu,
               ul = ul,
               uu = uu,
               con = con)

# Trajectory initialization
x0 = [[q1; q1] for t = 1:T] # linear interpolation on state
u0 = [[1.0e-3 * rand(2); k0; 1.0e-3 * rand(model_ft.m - 3 - 1); h] for t = 1:T-1] # random controls

# Pack trajectories into vector
z0 = pack(x0, u0, prob)

#NOTE: may need to run examples multiple times to get good trajectories
# Solve nominal problem

if true
	# include_snopt()
	@time z̄ , info = solve(prob, copy(z0),
		nlp = :ipopt,
		tol = 1.0e-5, c_tol = 1.0e-5, mapl = 5,
		time_limit = 60)
	@show check_slack(z̄, prob)
	x̄, ū = unpack(z̄, prob)
	tf, t, h̄ = get_time(ū)

	# projection
	# Q = [Diagonal(ones(model_ft.n)) for t = 1:T]
	# R = [Diagonal(0.1 * ones(model_ft.m)) for t = 1:T-1]
	# x_proj, u_proj = lqr_projection(model_ft, x̄, ū, h̄[1], Q, R)

	# @show tf
	# @show h̄[1]
	# @save joinpath(pwd(), "examples/trajectories/hopper_vertical_gait.jld2") x̄ ū h̄ x_proj u_proj
else
	# @load joinpath(pwd(), "examples/trajectories/hopper_vertical_gait.jld2") x̄ ū h̄ x_proj u_proj
end

#
# using Plots
# plot(hcat(ū...)[1:2, :]',
#     linetype = :steppost,
#     label = "",
#     color = :red,
#     width = 2.0)
#
# plot!(hcat(u_proj...)[1:2, :]',
#     linetype = :steppost,
#     label = "", color = :black)
#
# plot(hcat(state_to_configuration(x̄)...)',
#     color = :red,
#     width = 2.0,
# 	label = "")
#
# plot!(hcat(state_to_configuration(x_proj)...)',
#     color = :black,
# 	label = "")

include(joinpath(pwd(), "models/visualize.jl"))
vis = Visualizer()
render(vis)
visualize!(vis, model_ft,
	state_to_configuration([x̄..., ([x̄[2:end] for i = 1:5]...)...]),
	Δt = h̄[1])

z_traj = [x[2] for x in x̄]
k_traj = [u[3] for u in ū]
λ_traj = [u[model_ft.idx_λ] for u in ū]
f_traj = [u[2] for u in ū]

plot(hcat(k_traj...)', linetype = :steppost, label = "impedance")
plot!(h̄[1] * hcat(λ_traj...)', linetype = :steppost, label = ["" "λ1"])
plot!(h̄[1] * hcat(f_traj...)', linetype = :steppost, label = "force")

plot(hcat(z_traj...)', linetype = :steppost, label = "z")

# @time z̄ , info = solve(prob, copy(z̄),
# 	nlp = :ipopt,
# 	tol = 1.0e-2, c_tol = 1.0e-2, mapl = 0,
# 	time_limit = 60)
#
# obj_track = quadratic_tracking_objective(
#     [Diagonal(100.0 * ones(model.n)) for t = 1:T],
#     [Diagonal(ones(model.m)) for t = 1:T-1],
#     [x̄[t] for t = 1:T],
#     [ū[t][1:end-1] for t = 1:T-1])
#
#
# prob_track = trajectory_optimization_problem(model,
#                obj_track,
#                T,
#                xl = xl,
#                xu = xu,
#                ul = [ul[t][1:end-1] for t = 1:T-1],
#                uu = [uu[t][1:end-1] for t = 1:T-1],
#                con = con_contact)
#
# z_track = pack(x̄, [ū[t][1:end-1] for t = 1:T-1], prob_track) #+ 1.0e-2 * ones(prob_track.num_var)
# # z_track[1:model.n] .+= 0.25
# # prob_track.primal_bounds[1][1:model.n] = z_track[1:model.n]
# # prob_track.primal_bounds[2][1:model.n] = z_track[1:model.n]
#
# @time z_sol , info = solve(prob_track, copy(z_track),
# 	nlp = :ipopt,
# 	tol = 1.0e-5, c_tol = 1.0e-5, mapl = 0,
# 	time_limit = 60)
#
# x_sol, u_sol = unpack(z_sol, prob_track)
# plot(hcat(state_to_configuration(x̄)...)')
# plot!(hcat(state_to_configuration(x_sol)...)')
#
# plot(hcat([ū[t][1:end-1] for t = 1:T-1]...)', color = :red, width = 2.0, label = "")
# plot!(hcat(u_sol...)', color = :black, width = 1.0)
#
# x_sol[1] - z_track[1:model.n]


text:
using Test

o2 = Dataset("/net/fluo/data2/data/ABSCO_CS_Database/v5.1_final/o2_v51.hdf")
cs_table    = o2["Gas_07_Absorption"]
wavenumber  = o2["Wavenumber"][:]
Temperature = o2["Temperature"][:]
Pressure    = o2["Pressure"][:]
Broadener   = o2["Broadener_01_VMR"][:]

i = 10
j = 64
PP = Pressure[j]/100
TT = Temperature[i,j]

hitran_data      = read_hitran(artifact("O2"), iso=1);
absorption_model = make_hitran_model(hitran_data, Voigt(), CEF=HumlicekWeidemann32SDErrorFunction(), vmr=0.0, wing_cutoff=150)
hitran_cs = absorption_cross_section(absorption_model, wavenumber, PP, TT)

a = vSmartMOM.loadAbsco("/net/fluo/data2/data/ABSCO_CS_Database/v5.1_final/o2_v51.hdf")
# Define interpolation table grid:
ν_grid = a.ν[1]:0.01:a.ν[end]
pressures = PP-10:10:PP+10.1
temperatures = TT-10:10:TT+10.1

model_interp = make_interpolation_model(a,  ν_grid, pressures, temperatures)
abscoInter_cs = absorption_cross_section(model_interp, ν_grid, PP, TT)

@test 1e20*abscoInter_cs ≈ 1e20*cs_table[:,1,i,j]


text:
using SparseArrays
using DataFrames

#
# Reading file as a strings
const Stop_Words = open("./stopwords-en.txt", "r") do file
  read(file, String) |>
    s -> split(s, "\n") |>
    a -> Dict{String, Bool}(w => true for w in a)
end

const P_REXPR1 = r"[\.,:;\"'\-?!\+\=_\*#%(){}\[\]$&~/0-9\|@><\\\^`]"

remove_punctuation(text::AbstractString) = replace(text, P_REXPR1 => " ") |>
  strip |> string

remove_punctuation(vtext::Vector{String}) = remove_punctuation.(vtext)

cleanup(text::AbstractString) = replace(text, P_REXPR1 => " ") |>
  strip |> string |> s -> lowercase.(s)

cleanup(vtext::Vector{String}) = cleanup.(vtext)

function word_count(text::String)
    hwc = Dict{String, Int}()
    for w ∈ split(text)
      w = lowercase(w)
      hwc[w] = get(hwc, w, 0) + 1
    end
    hwc
  end

## no need to worry about Vector{Union{Missing, String}}...
word_count(vtext::Vector{String}) = word_count.(vtext)

function gen_vocab(text::String; trim_len=3, trim_count=2)
  lowercase.(split(text)) |>
    a_ -> filter(w -> !haskey(Stop_Words, w), a_) |>
    vtext -> reduce(vcat, vtext; init=String[]) |>
    unique
end

function gen_token_ix(vtext::Vector{String})
  reduce(vcat, vtext) |>
     unique |>
     words -> filter(w -> length(w) > 2, words) |>
     vocab -> Dict{String, Int}(w => ix for (ix, w) in enumerate(vocab))
end

function h_encode(text::String; token_ix=token_ix)
  @assert !isnothing(token_ix)
  v = spzeros(Int8, length(token_ix)) #
  # v = zeros(Int8, length(token_ix))

  for w ∈ split(text)
    ix = get(token_ix, lowercase(w), -1)
    ix == -1 && continue
    v[ix] += one(Int8)
  end

  v
end

h_encode(vtext::Vector{String}; token_ix=token_ix) = h_encode.(vtext; token_ix)

function get_data(df::DF, features, output) where {DF <: Union{DataFrame, SubDataFrame}}
  s_features = [Symbol(f) for f ∈ features]
  df[:, :intercept] .= 1.0
  s_features = [:intercept, s_features...]
  X_matrix = convert(Matrix, select(df, s_features))
  y = df[!, output]

  (X_matrix, y)
end


text:
i1 = indicator(:stationary,1,1,Float64)
i2 = indicator(:stationary,4,2,Float32)
i3 = indicator(:stationary,10,10,Float64)

@test numproposals(i1) == 1 && numsamples(i1) == 1 && length(i1.samples) == 2 && length(i1.stationary) == 2
@test numproposals(i2) == 4 && numsamples(i2) == 2 && length(i2.samples) == 3 && length(i2.stationary) == 5
@test numproposals(i3) == 10 && numsamples(i3) == 10 && length(i3.samples) == 11 && length(i3.stationary) == 11

u2 = log([0.5f0])
u3 = log([0.8])

v1 = log([0.6])
v1c = vcat(v1,[0.0])
v2 = log([0.1f0,0.2f0,0.5f0,1.5f0])
v2c = vcat(v2+u2[1],[0.0f0])
v3 = log([1.0,1.5,2.0,0.5,0.3,0.1,2.0,0.3,1.5,0.5])
v3c = vcat(v3+u3[1],[0.0])

a1 = transitionprobability!(i1,v1)
a2 = transitionprobability!(i2,u2,v2)
a3 = transitionprobability!(i3,u3,v3)

@test_approx_eq sum(a1) 1.0
@test_approx_eq a1 exp(v1c-maximum(v1c))./sum(exp(v1c-maximum(v1c)))
@test_approx_eq sum(a2) 1.0f0
@test_approx_eq a2 exp(v2c-maximum(v2c))./sum(exp(v2c-maximum(v2c)))
@test_approx_eq sum(a3) 1.0
@test_approx_eq a3 exp(v3c-maximum(v3c))./sum(exp(v3c-maximum(v3c)))

@test (srand(567) ; sampleindicator!(i1)) == (srand(567) ; vcat([2],[rand(Distributions.Categorical(a1)) for i=1:numsamples(i1)]))
@test (srand(567) ; sampleindicator!(i2)) == (srand(567) ; vcat([5],[rand(Distributions.Categorical(a2)) for i=1:numsamples(i2)]))
@test (srand(567) ; sampleindicator!(i3)) == (srand(567) ; vcat([11],[rand(Distributions.Categorical(a3)) for i=1:numsamples(i3)]))

@test accepted(i1) == sum(i1.samples[2:end].!=i1.samples[1:end-1])
@test accepted(i2) == sum(i2.samples[2:end].!=i2.samples[1:end-1])
@test accepted(i3) == sum(i3.samples[2:end].!=i3.samples[1:end-1])

###test the show functions
println()
println("====================")
println("Test show() function")
println("====================")
show(i1)
show(i2)
show(i3)
println("====================")
println("End  show() function")
println("====================")
println()


text:
export marginalrule

@marginalrule typeof(dot)(:in1_in2) (m_out::UnivariateNormalDistributionsFamily, m_in1::PointMass{ <: AbstractVector }, m_in2::MultivariateNormalDistributionsFamily) = begin
    m, V = mean(m_out), cov(m_out)
    x = mean(m_in1)
    q_in2 = prod(ProdAnalytical(), m_in2, MvNormalWeightedMeanPrecision(x * weightedmean(m_out), x * precision(m_out) * x'))
    return (in1 = m_in1, in2 = q_in2)
end


text:
@testset "callback" begin
    a_cb = AvgStatsCallback()

    @test a_cb.order == 0
end


text:
#!jl # # EFTfitter.jl - BLUE Example
#jl # EFTfitter.jl - BLUE Example

#~When using multiple measurements of a single observable and a uniform prior for
#~the parameter representing the combined value, the combination of measurements performed with EFTfitter.jl
#~yields the same results as the Best Linear Unbiased Estimator (BLUE) method.

#~Here, we demonstrate this by using the examples of the paper
#~"How to combine correlated estimates of a single physical quantity"
#~by L. Lyons, D. Gibaut and P. Clifford (https://www.sciencedirect.com/science/article/pii/0168900288900186).
#~All numbers are taken from the example on charm particle lifetime experiments in section 5.
#~A factor of 10^13 is applied for convenience.

using EFTfitter
using BAT
using IntervalSets
using Statistics
using StatsBase
using LinearAlgebra
using Plots

#~We need one parameter for the best estimator and choose
#~a uniform distribution in the range 8 to 14 as prior:
parameters = BAT.NamedTupleDist(
    τ = 8..14,
)

#~When combining multiple measurements of the same observable,
#~only a function returning the combination parameter is needed:
estimator(params) = params.τ

#~In Eq. (17') of the reference paper the following covariance matrix is given:
covariance = [2.74 1.15 0.86 1.31;
              1.15 1.67 0.82 1.32;
              0.86 0.82 2.12 1.05;
              1.31 1.32 1.05 2.93]

#~For using this in EFTfitter.jl, we first need to convert the covariance matrix
#~into a correlation matrix and the corresponding uncertainty values:
corr, unc = EFTfitter.cov_to_cor(covariance)

measurements = (
    τ1 = Measurement(estimator,  9.5, uncertainties = (stat=unc[1],) ),
    τ2 = Measurement(estimator, 11.9, uncertainties = (stat=unc[2],) ),
    τ3 = Measurement(estimator, 11.1, uncertainties = (stat=unc[3],) ),
    τ4 = Measurement(estimator,  8.9, uncertainties = (stat=unc[4],) ),
)

correlations = (
    stat = Correlation(corr),
)

#~construct an `EFTfitterModel`:
model = EFTfitterModel(parameters, measurements, correlations)
posterior = PosteriorDensity(model);

#~sample the posterior with BAT.jl:
algorithm = MCMCSampling(mcalg =MetropolisHastings(), nsteps = 10^6, nchains = 4)
samples = bat_sample(posterior, algorithm).result

#~plot the posterior distribution for the combination parameter τ:
plot(samples, :τ, mean=true)
#md # ![blue plots](plots/plot_blue.png)

#~print numerical results of combination:
println("Mode: $(mode(samples).τ)")
println("Mean: $(mean(samples).τ) ± $(std(samples).τ)")
#md ```
#md Mode: 11.15985
#md Mean: 11.15471 ± 0.80180
#md ```

#!jl # ### Comparison with BLUE method
#jl # comparison with BLUE method:
blue = BLUE(model)
println("BLUE: $(blue.value) ± $(blue.unc)")
println("BLUE weights: $(blue.weights)")
#md ```
#md BLUE: 11.15983 ± 1.28604
#md BLUE weights: [0.145, 0.470, 0.347, 0.038]
#md ```


text:
#=

This file contains the `model_from_parameters` function, which computes all derived information
like optical thicknesses, from the input parameters. Produces a vSmartMOM_Model object.

=#

"Generate default set of parameters for Radiative Transfer calculations (from ModelParameters/)"
default_parameters() = parameters_from_yaml(joinpath(dirname(pathof(RadiativeTransfer)), "vSmartMOM", "ModelParameters", "DefaultParameters.yaml"))

"Take the parameters specified in the vSmartMOM_Parameters struct, and calculate derived attributes into a vSmartMOM_Model"
function model_from_parameters(params::vSmartMOM_Parameters)

    # Number of total bands and aerosols (for convenience)
    n_bands = length(params.spec_bands)
    n_aer = isnothing(params.scattering_params) ? 0 : length(params.scattering_params.rt_aerosols)

    # Create observation geometry
    obs_geom = ObsGeometry(params.sza, params.vza, params.vaz, params.obs_alt)

    # Create truncation type
    truncation_type = Scattering.δBGE{params.float_type}(params.l_trunc, params.Δ_angle)

    # Set quadrature points for streams
    quad_points = rt_set_streams(params.quadrature_type, params.l_trunc, obs_geom, params.polarization_type, array_type(params.architecture))

    # Get AtmosphericProfile from parameters
    vmr = isnothing(params.absorption_params) ? Dict() : params.absorption_params.vmr
    p_full, p_half, vmr_h2o, vcd_dry, vcd_h2o, new_vmr = compute_atmos_profile_fields(params.T, params.p, params.q, vmr)
    profile = AtmosphericProfile(params.T, p_full, params.q, p_half, vmr_h2o, vcd_dry, vcd_h2o, new_vmr)

    # Reduce the profile to the number of target layers (if specified)
    if params.profile_reduction_n != -1
        profile = vSmartMOM.reduce_profile(params.profile_reduction_n, profile);
    end

    # Rayleigh optical properties calculation
    greek_rayleigh = Scattering.get_greek_rayleigh(params.depol)
    τ_rayl = [zeros(params.float_type, length(params.T)) for i=1:n_bands];

    # τ_abs[iBand][iSpec,iZ]
    τ_abs     = [zeros(params.float_type, length(params.spec_bands[i]), length(profile.p_full)) for i in 1:n_bands]

    # Loop over all bands:
    for i_band=1:n_bands

        # i'th spectral band (convert from cm⁻¹ to μm)
        curr_band_λ = 1e4 ./ params.spec_bands[i_band]

        # Compute Rayleigh properties per layer for `i_band` band center
        τ_rayl[i_band]   = getRayleighLayerOptProp(profile.p_half[end], (maximum(curr_band_λ) + minimum(curr_band_λ)/2), params.depol, profile.vcd_dry);

        # If no absorption, continue to next band
        isnothing(params.absorption_params) && continue

        # Loop over all molecules in this band, obtain profile for each, and add them up
        for molec_i in 1:length(params.absorption_params.molecules[i_band])

            # Obtain hitran data for this molecule
            @timeit "Read HITRAN"  hitran_data = read_hitran(artifact(params.absorption_params.molecules[i_band][molec_i]), iso=1)

            println("Computing profile for $(params.absorption_params.molecules[i_band][molec_i]) with vmr $(profile.vmr[params.absorption_params.molecules[i_band][molec_i]]) for band #$(i_band)")

            # Calculate absorption profile
            @timeit "Absorption Coeff"  compute_absorption_profile!(τ_abs[i_band], hitran_data, params.absorption_params.broadening_function, params.absorption_params.wing_cutoff, params.absorption_params.CEF, params.architecture, profile.vmr[params.absorption_params.molecules[i_band][molec_i]], params.spec_bands[i_band], profile);
        end
    end

    # aerosol_optics[iBand][iAer]
    aerosol_optics = [Array{AerosolOptics}(undef, (n_aer)) for i=1:length(0)];

    FT2 = isnothing(params.scattering_params) ? params.float_type : typeof(params.scattering_params.rt_aerosols[1].τ_ref)

    # τ_aer[iBand][iAer,iZ]
    τ_aer = [zeros(FT2, n_aer, length(profile.p_full)) for i=1:n_bands];

    # Loop over aerosol type
    for i_aer=1:n_aer

        # Get curr_aerosol
        curr_aerosol = params.scattering_params.rt_aerosols[i_aer].aerosol

        # Create Aerosol size distribution for each aerosol species
        size_distribution = curr_aerosol.size_distribution

        # Create a univariate aerosol distribution
        mie_aerosol = Aerosol(size_distribution, curr_aerosol.nᵣ, curr_aerosol.nᵢ)
        #mie_aerosol = make_mie_aerosol(size_distribution, curr_aerosol.nᵣ, curr_aerosol.nᵢ, params.scattering_params.r_max, params.scattering_params.nquad_radius) #Suniti: why is the refractive index needed here?

        # Create the aerosol extinction cross-section at the reference wavelength:
        mie_model      = make_mie_model(params.scattering_params.decomp_type, mie_aerosol, params.scattering_params.λ_ref, params.polarization_type, truncation_type, params.scattering_params.r_max, params.scattering_params.nquad_radius)
        k_ref          = compute_ref_aerosol_extinction(mie_model, params.float_type)

        # Loop over bands
        for i_band=1:n_bands

            # i'th spectral band (convert from cm⁻¹ to μm)
            curr_band_λ = 1e4 ./ params.spec_bands[i_band]

            # Create the aerosols:
            mie_model      = make_mie_model(params.scattering_params.decomp_type, mie_aerosol, (maximum(curr_band_λ)+minimum(curr_band_λ))/2, params.polarization_type, truncation_type, params.scattering_params.r_max, params.scattering_params.nquad_radius)

            # Compute raw (not truncated) aerosol optical properties (not needed in RT eventually)
            # @show FT2
            @timeit "Mie calc"  aerosol_optics_raw = compute_aerosol_optical_properties(mie_model, FT2);

            # Compute truncated aerosol optical properties (phase function and fᵗ), consistent with Ltrunc:
            @show i_aer, i_band
            aerosol_optics[i_band][i_aer] = Scattering.truncate_phase(truncation_type, aerosol_optics_raw; reportFit=false)

            # Compute nAer aerosol optical thickness profiles
            τ_aer[i_band][i_aer,:] = params.scattering_params.rt_aerosols[i_aer].τ_ref[] * (aerosol_optics[i_band][i_aer].k/k_ref) * vSmartMOM.getAerosolLayerOptProp(1.0, params.scattering_params.rt_aerosols[i_aer].p₀, params.scattering_params.rt_aerosols[i_aer].σp, profile.p_full)

        end
    end

    # Check the floating-type output matches specified FT

    # Return the model
    return vSmartMOM_Model(params, aerosol_optics,  greek_rayleigh, quad_points, τ_abs, τ_rayl, τ_aer, obs_geom, profile)

end