text:
using Statistics: mean
using Oceananigans

include("utils.jl")

function impose_cooling_disk!(model::Model)
    g = model.grid
    c = model.constants

    # Parameters for generating initial surface heat flux.
    Rc = 600   # Radius of cooling disk [m].
    Ts = 20    # Surface temperature [°C].
    Q₀ = -800  # Cooling disk heat flux [W/m²].
    Q₁ = 0     # Noise added to cooling disk heat flux [W/m²].
    Ns = 5 * (c.f * Rc/g.Lz)  # Stratification or Brunt–Väisälä frequency [s⁻¹].

    αᵥ = 2.07e-4  # Volumetric coefficient of thermal expansion for water [K⁻¹].
    cₚ = 4181.3   # Isobaric mass heat capacity [J / kg·K].

    Tz = Ns^2 / (c.g * αᵥ)  # Vertical temperature gradient [K/m].

    # Center horizontal coordinates so that (x,y) = (0,0) corresponds to the center
    # of the domain (and the cooling disk).
    x₀ = g.xC .- mean(g.xC)
    y₀ = g.yC .- mean(g.yC)

    # Calculate vertical temperature profile and convert to Kelvin.
    T_ref = 273.15 .+ Ts .+ Tz .* (g.zC .- mean(Tz * g.zC))

    # Impose reference temperature profile.
    model.tracers.T.data .= repeat(reshape(T_ref, 1, 1, g.Nz), g.Nx, g.Ny, 1)

    # Add small temperature perturbations to the surface.
    model.tracers.T.data[:, :, 1] .+= 0.001*rand()

    # Generate surface heat flux field.
    Q = Q₀ .+ Q₁ * (0.5 .+ rand(g.Nx, g.Ny))

    # Set surface heat flux to zero outside of cooling disk of radius Rc.
    r₀² = @. x₀*x₀ + y₀'*y₀'
    Q[findall(r₀² .> Rc^2)] .= 0

    # Convert surface heat flux into 3D forcing term for use when calculating
    # source terms at each time step. Also convert surface heat flux [W/m²]
    # into a temperature tendency forcing [K/s].
    # @. model.forcings.FT.data[:, :, 1] = (Q / cₚ) * (g.Az / (model.eos.ρ₀ * g.V))

    # # We will impose a surface heat flux of -400 W/m² ≈ -4.5e-6 K/s in a
    # square on the surface of the domain.
    #
    # As we can only impose forcing functions right now, and imposing a cooling
    # disk takes a pretty complicated function, I'll just impose a cooling square
    # until we can figure out the best way of imposing forcings with complex
    # geometries, probably by being able to define, e.g. a forcing only at the
    # top, etc.
    @inline function cooling_disk(grid, u, v, w, T, S, i, j, k)
        if k == 1
            x = i*grid.Δx
            y = j*grid.Δy
            r² = (x - grid.Lx/2)^2 + (y - grid.Ly/2)^2
            if r² < 600^2
                return -4.5e-6
            else
                return 0
            end
        else
            return 0
        end
    end

    model.forcing = Forcing(nothing, nothing, nothing, cooling_disk, nothing)
    nothing
end

Nx, Ny, Nz = 100, 100, 50
Lx, Ly, Lz = 2000, 2000, 1000
Nt, Δt = 1080, 20   # 6 hour simulation

model = Model(N=(Nx, Ny, Nz), L=(Lx, Ly, Lz), ν=4e-2, κ=4e-2 )  # CPU Model
# model = Model(N=(Nx, Ny, Nz), L=(Lx, Ly, Lz), arch=:GPU, float_type=Float32)  # GPU Model

impose_cooling_disk!(model)

# Add a NetCDF output writer that saves NetCDF files to the current directory
# "." with a filename prefix of "deep_convection_" every 20 iterations.
nc_writer = NetCDFOutputWriter(dir=".", prefix="deep_convection_", frequency=20)
push!(model.output_writers, nc_writer)

# time_step!(model; Nt=Nt, Δt=Δt)
for i = 1:Nt
    tic = time_ns()
    time_step!(model, 1, Δt)
    println("Time: $(model.clock.time)  [$(prettytime(time_ns()-tic))]")
end

make_vertical_slice_movie(model, nc_writer, "T", Nt, Δt, 293.13, ceil(Int, Ny/2))


text:
using Images, CoordinateTransformations, ImageView

img = load("sample-images/cats-3061372_640.jpg");
tfm = LinearMap(RotMatrix(-pi/3)) # rotate by -60 degrees
img = warp(img, tfm)
imshow(img)

ext: jl


text:
export iauEcm06
"""
ICRS equatorial to ecliptic rotation matrix, IAU 2006.

This function is part of the International Astronomical Union's
SOFA (Standards of Fundamental Astronomy) software collection.

Status:  support function.

Given:
   date1,date2  double         TT as a 2-part Julian date (Note 1)

Returned:
   rm           double[3][3]   ICRS to ecliptic rotation matrix

Notes:

   1. The TT date date1+date2 is a Julian Date, apportioned in any
   convenient way between the two arguments.  For example,
   JD(TT)=2450123.7 could be expressed in any of these ways,
   among others:

         date1          date2

      2450123.7           0.0       (JD method)
      2451545.0       -1421.3       (J2000 method)
      2400000.5       50123.2       (MJD method)
      2450123.5           0.2       (date & time method)

   The JD method is the most natural and convenient to use in
   cases where the loss of several decimal digits of resolution
   is acceptable.  The J2000 method is best matched to the way
   the argument is handled internally and will deliver the
   optimum resolution.  The MJD method and the date & time methods
   are both good compromises between resolution and convenience.

   1. The matrix is in the sense

      E_ep = rm x P_ICRS,

   where P_ICRS is a vector with respect to ICRS right ascension
   and declination axes and E_ep is the same vector with respect to
   the (inertial) ecliptic and equinox of date.

   2. P_ICRS is a free vector, merely a direction, typically of unit
   magnitude, and not bound to any particular spatial origin, such
   as the Earth, Sun or SSB.  No assumptions are made about whether
   it represents starlight and embodies astrometric effects such as
   parallax or aberration.  The transformation is approximately that
   between mean J2000.0 right ascension and declination and ecliptic
   longitude and latitude, with only frame bias (always less than
   25 mas) to disturb this classical picture.

Called:
   iauObl06     mean obliquity, IAU 2006
   iauPmat06    PB matrix, IAU 2006
   iauIr        initialize r-matrix to identity
   iauRx        rotate around X-axis
   iauRxr       product of two r-matrices

This revision:  2015 December 11

SOFA release 2018-01-30

Copyright (C) 2018 IAU SOFA Board.  See notes at end.
"""

# void iauEcm06(double date1, double date2, double rm[3][3])

function iauEcm06(date1::Real, date2::Real)
   # Allocate return value
   rm = zeros(Float64, 3, 3)

   ccall((:iauEcm06, libsofa_c), Cvoid,
         (Cdouble, Cdouble,
         Ptr{Cdouble}),
         convert(Float64, date1), convert(Float64, date2),
         rm)

   return SMatrix{3,3}(rm')
end


text:
"""
    AlphaVectorPolicy(pomdp::POMDP, alphas, action_map)

Construct a policy from alpha vectors.

# Arguments
- `alphas`: an |S| x (number of alpha vecs) matrix or a vector of alpha vectors.
- `action_map`: a vector of the actions correponding to each alpha vector

    AlphaVectorPolicy{P<:POMDP, A}

Represents a policy with a set of alpha vectors.

Use `action` to get the best action for a belief, and `alphavectors` and `alphapairs` to

# Fields
- `pomdp::P` the POMDP problem
- `n_states::Int` the number of states in the POMDP
- `alphas::Vector{Vector{Float64}}` the list of alpha vectors
- `action_map::Vector{A}` a list of action corresponding to the alpha vectors
"""
struct AlphaVectorPolicy{P<:POMDP, A} <: Policy
    pomdp::P # needed for mapping states to locations in alpha vectors
    n_states::Int
    alphas::Vector{Vector{Float64}}
    action_map::Vector{A}
end

@deprecate AlphaVectorPolicy(pomdp::POMDP, alphas) AlphaVectorPolicy(pomdp, alphas, ordered_actions(pomdp))

function AlphaVectorPolicy(m::POMDP, alphas::AbstractVector, amap)
    AlphaVectorPolicy(m, length(states(m)), alphas,
                      convert(Vector{actiontype(m)}, amap))
end

# assumes alphas is |S| x (number of alpha vecs)
function AlphaVectorPolicy(p::POMDP, alphas::Matrix{Float64}, action_map)
    # turn alphas into vector of vectors
    num_actions = size(alphas, 2)
    alpha_vecs = Vector{Float64}[]
    for i = 1:num_actions
        push!(alpha_vecs, vec(alphas[:,i]))
    end

    AlphaVectorPolicy(p, length(states(p)), alpha_vecs,
                      convert(Vector{actiontype(p)}, action_map))
end

updater(p::AlphaVectorPolicy) = DiscreteUpdater(p.pomdp)

"""
Return an iterator of alpha vector-action pairs in the policy.
"""
alphapairs(p::AlphaVectorPolicy) = (p.alphas[i]=>p.action_map[i] for i in 1:length(p.alphas))

"""
Return the alpha vectors.
"""
alphavectors(p::AlphaVectorPolicy) = p.alphas

# The three functions below rely on beliefvec being implemented for the belief type
# Implementations of beliefvec are below
function value(p::AlphaVectorPolicy, b)
    bvec = beliefvec(p.pomdp, p.n_states, b)
    maximum(dot(bvec,a) for a in p.alphas)
end

function action(p::AlphaVectorPolicy, b)
    bvec = beliefvec(p.pomdp, p.n_states, b)
    num_vectors = length(p.alphas)
    best_idx = 1
    max_value = -Inf
    for i = 1:num_vectors
        temp_value = dot(bvec, p.alphas[i])
        if temp_value > max_value
            max_value = temp_value
            best_idx = i
        end
    end
    return p.action_map[best_idx]
end

function actionvalues(p::AlphaVectorPolicy, b)
    bvec = beliefvec(p.pomdp, p.n_states, b)
    num_vectors = length(p.alphas)
    max_values = -Inf*ones(length(actions(p.pomdp)))
    for i = 1:num_vectors
        temp_value = dot(bvec, p.alphas[i])
        ai = actionindex(p.pomdp, p.action_map[i])
        if temp_value > max_values[ai]
            max_values[ai] = temp_value
        end
    end
    return max_values
end

"""
    POMDPPolicies.beliefvec(m::POMDP, n_states::Int, b)

Return a vector-like representation of the belief `b` suitable for calculating the dot product with the alpha vectors.
"""
function beliefvec end

function beliefvec(m::POMDP, n, b::SparseCat)
    return sparsevec(collect(stateindex(m, s) for s in b.vals), collect(b.probs), n)
end
beliefvec(m::POMDP, n, b::DiscreteBelief) = b.b
beliefvec(m::POMDP, n, b::AbstractArray) = b

function beliefvec(m::POMDP, n_states, b)
    sup = support(b)
    bvec = zeros(n_states)
    for s in sup
        bvec[stateindex(m, s)] = pdf(b, s)
    end
    return bvec
end

function Base.push!(p::AlphaVectorPolicy, alpha::Vector{Float64}, a)
    push!(p.alphas, alpha)
    push!(p.action_map, a)
end

@deprecate beliefvec(m::POMDP, b) beliefvec(m, length(states(m)), b)

max_stars_repo_name: JuliaPOMDP/POMDPPolicies

max_issues_repo_name: JuliaPOMDP/POMDPPolicies

max_forks_repo_name: JuliaPOMDP/POMDPPolicies


text:
# timestamp.jl
#
#   Convert difference timestamps
#
# Jeremy Rüffer
# Thünen Institut
# Institut für Agrarklimaschutz
# Julia 1.5.0
# 02.04.2020
# Last Edit: 03.08.2020

using StofffluesseAG
using Test
using Dates

println("\n=======================================================")
println("============  StofffluesseAG Module Tests  ============")
println("=======================================================")

@testset "timestamp_licor Test" begin
    @test (1.481490002e9,4.9999872e7) == timestamp_licor(DateTime(2016,12,11,22,0,2,50)) || "TIMESTAMP_LICOR: timestamp_licor(DateTime(2016,12,11,22,0,2,50)) output incorrect, should be (1.481490002e9,4.9999872e7)"
    @test DateTime(2016,12,11,22,0,2,550) == timestamp_licor(1481490002,550000000) || "TIMESTAMP_LICOR: timestamp_licor(1481490002,550000000) output incorrect, should be 2016-12-11T22:00:02.55"
end
println("\n\n")

@testset "timestamp_lgr Test" begin
    @test DateTime(2017,5,15,15,13,47,994) == timestamp_lgr("15/05/2017 15:13:47.994") || "TIMESTAMP_LGR: timestamp_lgr(\"15/05/2017 15:13:47.994\") output incorrect, should be 2017-05-15T15:13:47.994"
    @test "15/05/2017 15:13:47.993" == timestamp_lgr(DateTime(2017,5,15,15,13,47,993)) || "TIMESTAMP_LGR: timestamp_lgr(DateTime(2017,5,15,15,13,47,993)) output incorrect, should be \"15/05/2017 15:13:47.993\""
end
println("\n\n")

@testset "timestamp_aerodyne Test" begin
    @test DateTime(2017,5,11,0,0,0,510) == timestamp_aerodyne(3577305600.510) || "TIMESTAMP_AERODYNE: timestamp_aerodyne(3577305600.510) output incorrect, should be 2017-05-11T00:00:00.51"
    @test 3577305600.510 == timestamp_aerodyne(DateTime(2017,5,11,0,0,0,510)) || "TIMESTAMP_AERODYNE: timestamp_aerodyne(DateTime(2017,5,11,0,0,0,510)) output incorrect, should be 3577305600.510"
end
println("\n\n")

@testset "dirlist Test" begin
    dir = joinpath(splitdir(@__FILE__)[1],"dirtest")
    targetFiles = Array{String}(undef,4)
    targetFiles[1] = joinpath(dir,"a.txt")
    targetFiles[2] = joinpath(dir,"b.txt")
    targetFiles[3] = joinpath(dir,"c.txt")
    targetFiles[4] = joinpath(dir,"subfolder","d.txt")
    targetFolders = Array{String}(undef,1)
    targetFolders[1] = joinpath(dir,"subfolder")
    @test (targetFiles[1:3], targetFolders) == dirlist(dir,recur=1) || "DIRLIST: dirlist(dir,recur=1) output incorrect"
    @test (targetFiles, targetFolders) == dirlist(dir,recur=2) || "DIRLIST: dirlist(dir,recur=2) output incorrect"
end
println("\n\n")

@testset "NaN Tests" begin
    x = [1.0,NaN,3.0,2.0]
    @test 1.0 == nanmin(x) || "NANMIN: nanmin([1.0,NaN,3.0,4.0]) output incorrect, should be 1.0"
    @test 3.0 == nanmax(x) || "NANMAX: nanmax([1.0,NaN,3.0,4.0]) output incorrect, should be 3.0"
    @test 1.0 == nanstd(x) || "NANSTD: nanstd([1.0,NaN,3.0,4.0]) output incorrect, should be 1.0"
    @test 2.0 == nanmean(x) || "NANMEAN: nanmean([1.0,NaN,3.0,4.0]) output incorrect, should be 2.0"
end
println("\n\n")

@testset "nCLD Tests" begin
    @test false == nCLD_Status(0).AUXDevice || "nCLD_Status(0).AUXDevice should be FALSE"
    @test false == nCLD_Status(0).AdditionalConvHeating || "nCLD_Status(0).AdditionalConvHeating should be FALSE"
    @test false == nCLD_Status(0).BypassPressureRegulator || "nCLD_Status(0).BypassPressureRegulator should be FALSE"
    @test false == nCLD_Status(0).CalValve || "nCLD_Status(0).CalValve should be FALSE"
    @test false == nCLD_Status(0).CalibrationActive || "nCLD_Status(0).CalibrationActive should be FALSE"
    @test false == nCLD_Status(0).HotTubingHeating || "nCLD_Status(0).HotTubingHeating should be FALSE"
    @test false == nCLD_Status(0).OzoneGenerator || "nCLD_Status(0).OzoneGenerator should be FALSE"
    @test false == nCLD_Status(0).PeltierCooling || "nCLD_Status(0).PeltierCooling should be FALSE"
    @test false == nCLD_Status(0).PowerUpPhase || "nCLD_Status(0).PowerUpPhase should be FALSE"
    @test false == nCLD_Status(0).ReactorChamberHeating || "nCLD_Status(0).ReactorChamberHeating should be FALSE"
    @test false == nCLD_Status(0).Recorder || "nCLD_Status(0).Recorder should be FALSE"
    @test false == nCLD_Status(0).RemoteMode || "nCLD_Status(0).RemoteMode should be FALSE"
    @test false == nCLD_Status(0).ScrubberHeating || "nCLD_Status(0).ScrubberHeating should be FALSE"
    @test false == nCLD_Status(0).StandbyOperation || "nCLD_Status(0).StandbyOperation should be FALSE"
    @test false == nCLD_Status(0).VacuumPump || "nCLD_Status(0).VacuumPump should be FALSE"

    @test false == nCLD_Status(42).AUXDevice || "nCLD_Status(0).AUXDevice should be FALSE"
    @test false == nCLD_Status(42).AdditionalConvHeating || "nCLD_Status(0).AdditionalConvHeating should be FALSE"
    @test false == nCLD_Status(42).BypassPressureRegulator || "nCLD_Status(0).BypassPressureRegulator should be FALSE"
    @test false == nCLD_Status(42).CalValve || "nCLD_Status(0).CalValve should be FALSE"
    @test false == nCLD_Status(42).CalibrationActive || "nCLD_Status(0).CalibrationActive should be FALSE"
    @test true == nCLD_Status(42).HotTubingHeating || "nCLD_Status(0).HotTubingHeating should be TRUE"
    @test false == nCLD_Status(42).OzoneGenerator || "nCLD_Status(0).OzoneGenerator should be FALSE"
    @test false == nCLD_Status(42).PeltierCooling || "nCLD_Status(0).PeltierCooling should be FALSE"
    @test false == nCLD_Status(42).PowerUpPhase || "nCLD_Status(0).PowerUpPhase should be FALSE"
    @test true == nCLD_Status(42).ReactorChamberHeating || "nCLD_Status(0).ReactorChamberHeating should be TRUE"
    @test false == nCLD_Status(42).Recorder || "nCLD_Status(0).Recorder should be FALSE"
    @test false == nCLD_Status(42).RemoteMode || "nCLD_Status(0).RemoteMode should be FALSE"
    @test false == nCLD_Status(42).ScrubberHeating || "nCLD_Status(0).ScrubberHeating should be FALSE"
    @test false == nCLD_Status(42).StandbyOperation || "nCLD_Status(0).StandbyOperation should be FALSE"
    @test false == nCLD_Status(42).VacuumPump || "nCLD_Status(0).VacuumPump should be FALSE"
end
println("\n\n")

@testset "findnewton Tests" begin
    data = collect(1:1000)
    @test 400 == findnewton(data,400) || "FINDNEWTON: findnewton(data,400) should return 400"
    @test 999 == findnewton(data,999) || "FINDNEWTON: findnewton(data,999) should return 999"
    @test 1 == findnewton(data,1) || "FINDNEWTON: findnewton(data,1) should return 1"
    @test 42 == findnewton(data,42) || "FINDNEWTON: findnewton(data,42) should return 42"
    @test 500 == findnewton(data,500) || "FINDNEWTON: findnewton(data,500) should return 500"
end
println("\n\n")

println("=====  StoffflüßeAG Tests Complete  =====")


text:
module AirfoilDatabase

using Printf
struct AirfoilData
    name::String
    desc::String
    x::Vector{Float64}
    y::Vector{Float64}
end

function Base.show(io::IO, airfoil::AirfoilData)
    print(io, airfoil.name, ' ', '(', airfoil.desc, ')')
end

include("database.jl")

function query_airfoil(name::AbstractString)::Vector{AirfoilData}
    que = split(lowercase(name))
    result = Vector{AirfoilData}()
    for i in 1:db_size
        idx = index[i]
        if all(q->occursin(q, idx), que)
            push!(result, database[i])
        end
    end
    return result
end

include("NACA.jl")

export AirfoilData, query_airfoil, NACA

end # module


text:
@enum JoinKind begin
    inner = 1
    left = 2
    right = 3
    outer = 4
end

"""
# Description

A logical node that represents an join operation. It has three fields:

1. `sources`: A tuple of two aliased relations that will be joined. Joins
    between more than two relations are represented as trees of binary joins.
2. `predicates`: A tuple of `Expression` objects that defines the predicates
    for the join.
3. `kind`: A enum value of type `JoinKind` indicating whether this is an inner,
    left, right or outer join.
"""
struct Join <: LogicalNode
    sources::Any
    predicates::Tuple
    kind::JoinKind
end

function Base.print(io::IO, node::Join)
    println(io, "Join LogicalNode")
    println(io, node.kind)
    println(io, type(node.source))
    for e in node.predicates
        println(io, e)
    end
end


text:
#=
Tools for normal form games.
=#

const opponents_actions_docstring = """
`opponents_actions::Union{Action,ActionProfile,Nothing}` : Profile of N-1
  opponents' actions. If N=2, then it must be a vector of reals (in which case
  it is treated as the opponent's mixed action) or a scalar of integer (in which
  case it is treated as the opponent's pure action). If N>2, then it must be a
  tuple of N-1 integers (pure actions) or N-1 vectors of reals (mixed actions).
  (For the degenerate case N=1, it must be `nothing`.)"""

# Player #

"""
    Player{N,T}

Type representing a player in an N-player normal form game.

# Fields

- `payoff_array::Array{T<:Real}` : Array representing the player's payoff
  function.
"""
struct Player{N,T<:Real}
    payoff_array::Array{T,N}

    function Player{N,T}(payoff_array::Array{T,N}) where {N,T<:Real}
        if prod(size(payoff_array)) == 0
            throw(ArgumentError("every player must have at least one action"))
        end
        return new(payoff_array)
    end
end

Player(payoff_array::Array{T,N}) where {T<:Real,N} = Player{N,T}(payoff_array)

Player{N,T}(player::Player{N,S}) where {N,T,S} =
    Player(Array{T}(player.payoff_array))

Base.convert(::Type{T}, player::Player) where {T<:Player} =
    player isa T ? player : T(player)

num_actions(p::Player) = size(p.payoff_array, 1)
num_opponents(::Player{N}) where {N} = N - 1

Base.summary(player::Player) =
    string(Base.dims2string(size(player.payoff_array)),
           " ",
           split(string(typeof(player)), ".")[end])

function Base.show(io::IO, player::Player)
    print(io, summary(player))
    println(io, ":")
    Base.print_array(io, player.payoff_array)
end

# delete_action

"""
    delete_action(player, action[, player_idx=1])

Return a new Player instance with the action(s) specified by `action` deleted
from the action set of the player specified by `player_idx`.

# Arguments

- `player::Player` : Player instance.
- `action::Union{PureAction,AbstractVector{<:PureAction}}`: The action(s) to be
  deleted.
- `player_idx::Integer` : Index of the player to delete action(s) for.

# Returns

- `::Player` : `Player` instance with the action(s) deleted as specified.

# Examples

```julia
julia> player = Player([3 0; 0 3; 1 1])
3×2 Player{2,Int64}:
 3  0
 0  3
 1  1

julia> delete_action(player, 3)
2×2 Player{2,Int64}:
 3  0
 0  3

julia> delete_action(player, 1, 2)
3×1 Player{2,Int64}:
 0
 3
 1
```
"""
function delete_action(player::Player{N,T}, action::AbstractVector{<:PureAction},
                       player_idx::Integer=1) where {N,T}
    sel = Any[Colon() for i in 1:N]
    sel[player_idx] = setdiff(axes(player.payoff_array, player_idx), action)
    payoff_array_new = player.payoff_array[sel...]::Array{T,N}
    return Player(payoff_array_new)
end

delete_action(player::Player, action::PureAction, player_idx::Integer=1) =
    delete_action(player, [action], player_idx)

# Reference
# https://stackoverflow.com/questions/53637511/how-to-delete-the-specific-row-of-n-dimensional-array-in-julia

# payoff_vector

# To resolve definition ambiguity
function payoff_vector(player::Player, opponents_actions::Tuple{})
    throw(ArgumentError("input tuple must not be empty"))
end

"""
    payoff_vector(player, opponents_actions)

Return a vector of payoff values for a Player in an N>2 player game, one for
each own action, given a tuple of the opponents' pure actions.

# Arguments

- `player::Player` : Player instance.
- `opponents_actions::PureActionProfile` : Tuple of N-1 opponents' pure actions.

# Returns

- `::Vector` : Payoff vector.
"""
function payoff_vector(player::Player, opponents_actions::PureActionProfile)
    length(opponents_actions) != num_opponents(player) &&
        throw(ArgumentError(
            "length of opponents_actions must be $(num_opponents(player))"
        ))
    payoffs = player.payoff_array
    for i in num_opponents(player):-1:1
        payoffs = _reduce_ith_opponent(payoffs, i, opponents_actions[i])
    end
    return vec(payoffs)
end

"""
    payoff_vector(player, opponents_actions)

Return a vector of payoff values for a Player in an N>2 player game, one for
each own action, given a tuple of the opponents' mixed actions.

# Arguments

- `player::Player` : Player instance.
- `opponents_actions::MixedActionProfile` : Tuple of N-1 opponents' mixed
  actions.

# Returns

- `::Vector` : Payoff vector.
"""
function payoff_vector(player::Player{N,T1},
                       opponents_actions::MixedActionProfile{T2}) where {N,T1,T2}
    length(opponents_actions) != num_opponents(player) &&
        throw(ArgumentError(
            "length of opponents_actions must be $(num_opponents(player))"
        ))
    S = promote_type(T1, T2)
    payoffs::Array{S,N} = player.payoff_array
    for i in num_opponents(player):-1:1
        payoffs = _reduce_ith_opponent(payoffs, i, opponents_actions[i])
    end
    return vec(payoffs)
end

"""
    payoff_vector(player, opponent_action)

Return a vector of payoff values for a Player in a 2-player game, one for each
own action, given the opponent's pure action.

# Arguments

- `player::Player{2}` : Player instance.
- `opponent_action::PureAction` : Opponent's pure action (integer).

# Returns

- `::Vector` : Payoff vector.
"""
function payoff_vector(player::Player{2}, opponent_action::PureAction)
    return player.payoff_array[:, opponent_action]
end

"""
    payoff_vector(player, opponent_action)

Return a vector of payoff values for a Player in a 2-player game, one for each
own action, given the opponent's mixed action.

# Arguments

- `player::Player{2}` : Player instance.
- `opponent_action::MixedAction` : Opponent's mixed action (vector of reals).

# Returns

- `::Vector` : Payoff vector.
"""
function payoff_vector(player::Player{2}, opponent_action::MixedAction)
    # player.num_opponents == 1
    return player.payoff_array * opponent_action
end

# Trivial case with player.num_opponents == 0
"""
    payoff_vector(player, opponent_action)

Return a vector of payoff values for a Player in a trivial game with 1 player,
one for each own action.

# Arguments

- `player::Player{1}` : Player instance.
- `opponent_action::Nothing`

# Returns

- `::Vector` : Payoff vector.
"""
function payoff_vector(player::Player{1}, opponent_action::Nothing)
    return player.payoff_array
end

# _reduce_ith_opponent

# Given an N-d array `payoff_array` of shape (n_{1}, ..., n_{i+1}, 1, ..., 1),
# return the N-d array of payoffs of shape (n_{1}, ..., n_{i}, 1, 1, ..., 1)
# given by fixing the i-th opponent's pure or mixed action to be `action`.
for (S, ex_mat_action) in ((PureAction, :(A[:, action])),
                           (MixedAction, :(A * action)))
    @eval function _reduce_ith_opponent(payoff_array::Array{T,N},
                                        i::Int, action::$S) where {N,T}
        shape = size(payoff_array)
        A = reshape(payoff_array, (prod(shape[1:i]), shape[i+1]))::Matrix{T}
        out = $(ex_mat_action)
        shape_new = tuple(shape[1:i]..., ones(Int, N-i)...)::NTuple{N,Int}
        return reshape(out, shape_new)
    end
end

# is_best_response

"""
    is_best_response(player, own_action, opponents_actions; tol=1e-8)

Return true if `own_action` is a best response to `opponents_actions`.

# Arguments

- `player::Player` : Player instance.
- `own_action::PureAction` : Own pure action (integer).
- $(opponents_actions_docstring)
- `tol::Real` : Tolerance to be used to determine best response actions.

# Returns

- `::Bool` : True if `own_action` is a best response to `opponents_actions`;
  false otherwise.
"""
function is_best_response(player::Player,
                          own_action::PureAction,
                          opponents_actions::Union{Action,ActionProfile,Nothing};
                          tol::Real=1e-8)
    payoffs = payoff_vector(player, opponents_actions)
    payoff_max = maximum(payoffs)
    return payoffs[own_action] >= payoff_max - tol
end

"""
    is_best_response(player, own_action, opponents_actions; tol=1e-8)

Return true if `own_action` is a best response to `opponents_actions`.

# Arguments

- `player::Player` : Player instance.
- `own_action::MixedAction` : Own mixed action (vector of reals).
- $(opponents_actions_docstring)
- `tol::Real` : Tolerance to be used to determine best response actions.

# Returns

- `::Bool` : True if `own_action` is a best response to `opponents_actions`;
  false otherwise.
"""
function is_best_response(player::Player,
                          own_action::MixedAction,
                          opponents_actions::Union{Action,ActionProfile,Nothing};
                          tol::Real=1e-8)
    payoffs = payoff_vector(player, opponents_actions)
    payoff_max = maximum(payoffs)
    return dot(own_action, payoffs) >= payoff_max - tol
end

# best_response

"""
    best_responses(player, opponents_actions; tol=1e-8)

Return all the best response actions to `opponents_actions`.

# Arguments

- `player::Player` : Player instance.
- $(opponents_actions_docstring)
- `tol::Real` : Tolerance to be used to determine best response actions.

# Returns

- `best_responses::Vector{Int}` : Vector containing all the best response
  actions.
"""
function best_responses(player::Player,
                        opponents_actions::Union{Action,ActionProfile,Nothing};
                        tol::Real=1e-8)
    payoffs = payoff_vector(player, opponents_actions)
    payoff_max = maximum(payoffs)
    best_responses = findall(x -> x >= payoff_max - tol, payoffs)
    return best_responses
end

"""
    best_response([rng=GLOBAL_RNG], player, opponents_actions;
                  tie_breaking=:smallest, tol=1e-8)

Return a best response action to `opponents_actions`.

# Arguments

- `rng::AbstractRNG=GLOBAL_RNG` : Random number generator; relevant only with
  `tie_breaking=:random`.
- `player::Player` : Player instance.
- $(opponents_actions_docstring)
- `tie_breaking::Symbol` : Control how to break a tie (see Returns for details).
- `tol::Real` : Tolerance to be used to determine best response actions.

# Returns

- `::Int` : If `tie_breaking=:smallest`, returns the best response action with
  the smallest index; if `tie_breaking=:random`, returns an action randomly
  chosen from the best response actions.
"""
function best_response(rng::AbstractRNG, player::Player,
                       opponents_actions::Union{Action,ActionProfile,Nothing};
                       tie_breaking::Symbol=:smallest,
                       tol::Real=1e-8)
    if tie_breaking == :smallest
        payoffs = payoff_vector(player, opponents_actions)
        return argmax(payoffs)
    elseif tie_breaking == :random
        brs = best_responses(player, opponents_actions; tol=tol)
        return rand(rng, brs)
    else
        throw(ArgumentError(
            "tie_breaking must be one of `:smallest` or `:random`"
        ))
    end
end

best_response(player::Player,
              opponents_actions::Union{Action,ActionProfile,Nothing};
              tie_breaking::Symbol=:smallest, tol::Real=1e-8) =
    best_response(Random.GLOBAL_RNG, player, opponents_actions,
                  tie_breaking=tie_breaking, tol=tol)

"""
    BROptions

Struct to contain options for `best_response`.

# Fields
- `tol::Real=1e-8` : Tolerance level.
- `tie_breaking::Symbol=:smallest` : `:smallest` or `:random`.
- `rng::AbstractRNG=GLOBAL_RNG` : Random number generator.
"""
@with_kw struct BROptions{T<:Real,TR<:AbstractRNG}
    tol::T = 1e-8
    tie_breaking::Symbol = :smallest
    rng::TR = Random.GLOBAL_RNG
end

"""
    best_response(player, opponents_actions, options)

Return a best response action to `opponents_actions` with options as specified
by a `BROptions` instance `options`.
"""
best_response(player::Player,
              opponents_actions::Union{Action,ActionProfile,Nothing},
              options::BROptions) =
    best_response(options.rng, player, opponents_actions;
                  tie_breaking=options.tie_breaking, tol=options.tol)

# Perturbed best response
"""
    best_response(player, opponents_actions, payoff_perturbation)

Return the perturbed best response to `opponents_actions`.

# Arguments

- `player::Player` : Player instance.
- $(opponents_actions_docstring)
- `payoff_perturbation::Vector{Float64}` : Vector of length equal to the number
  of actions of the player containing the values ("noises") to be added to the
  payoffs in determining the best response.

# Returns

- `::Int` : Best response action.
"""
function best_response(player::Player,
                       opponents_actions::Union{Action,ActionProfile,Nothing},
                       payoff_perturbation::Vector{Float64})
    length(payoff_perturbation) != num_actions(player) &&
        throw(ArgumentError(
            "length of payoff_perturbation must be $(num_actions(player))"
        ))

    payoffs = payoff_vector(player, opponents_actions) + payoff_perturbation
    return argmax(payoffs)
end

# NormalFormGame #

"""
    NormalFormGame{N,T}

Type representing an N-player normal form game.

# Fields

- `players::NTuple{N,Player{N,T<:Real}}` : Tuple of Player instances.
- `nums_actions::NTuple{N,Int}` : Tuple of the numbers of actions, one for each
  player.
"""
struct NormalFormGame{N,T<:Real}
    players::NTuple{N,Player{N,T}}
    nums_actions::NTuple{N,Int}
end

num_players(::NormalFormGame{N}) where {N} = N

function NormalFormGame(::Tuple{})  # To resolve definition ambiguity
    throw(ArgumentError("input tuple must not be empty"))
end

"""
    NormalFormGame([T=Float64], nums_actions)

Constructor of an N-player NormalFormGame, consisting of payoffs all 0.

# Arguments

- `T::Type` : Type of payoff values; defaults to `Float64` if not specified.
- `nums_actions::NTuple{N,Int}` : Numbers of actions of the N players.
"""
function NormalFormGame(T::Type, nums_actions::NTuple{N,Int}) where N
    players = Vector{Player{N,T}}(undef, N)
    for i in 1:N
        sz = ntuple(j -> nums_actions[i-1+j <= N ? i-1+j : i-1+j-N], N)
        players[i] = Player(zeros(T, sz))
    end
    return NormalFormGame(players)
end

NormalFormGame(nums_actions::NTuple{N,Int}) where {N} =
    NormalFormGame(Float64, nums_actions)

# Check that the shapes of the payoff arrays are consistent
function is_consistent(players::NTuple{N,Player{N,T}}) where {N,T}
    shape_1 = size(players[1].payoff_array)
    for i in 2:N
        shape = size(players[i].payoff_array)
        for j in 1:(N-i+1)
            shape[j] == shape_1[i-1+j] || return false
        end
        for j in (N-i+2):N
            shape[j] == shape_1[j-(N-i+1)] || return false
        end
    end
    return true
end

"""
    NormalFormGame(players)

Constructor of an N-player NormalFormGame with a tuple of N Player instances.

# Arguments

- `players::NTuple{N,Player}` : Tuple of Player instances.
"""
function NormalFormGame(players::NTuple{N,Player{N,T}}) where {N,T}
    is_consistent(players) ||
        throw(ArgumentError("shapes of payoff arrays must be consistent"))
    nums_actions = ntuple(i -> num_actions(players[i]), N)
    return NormalFormGame{N,T}(players, nums_actions)
end

"""
    NormalFormGame(players)

Constructor of an N-player NormalFormGame with a vector of N Player instances.

# Arguments

- `players::Vector{Player}` : Vector of Player instances.
"""
NormalFormGame(players::Vector{Player{N,T}}) where {N,T} =
    NormalFormGame(ntuple(i -> players[i], N))

"""
    NormalFormGame(players...)

Constructor of an N-player NormalFormGame with N Player instances.

# Arguments

- `players::Player{N,T}...` : N Player instances

# Examples

```julia
# p1, p2, and p3 are all of type `Player{3,T}` for some `T`
NormalFormGame(p1, p2, p3)
```
"""
function NormalFormGame(players::Player{N,T}...) where {N,T}
    length(players) != N && error("Need $N players")
    NormalFormGame(players)  # use constructor for Tuple of players above
end

"""
    NormalFormGame(payoffs)

Construct an N-player NormalFormGame for N>=2 with an array `payoffs` of M=N+1
dimensions, where `payoffs[a_1, a_2, ..., a_N, :]` contains a profile of N
payoff values.

# Arguments

- `payoffs::Array{T<:Real}` : Array with ndims=N+1 containing payoff profiles.
"""
function NormalFormGame(payoffs::Array{T,M}) where {T<:Real,M}
    N = M - 1
    dims = Base.front(size(payoffs))
    colons = Base.front(ntuple(j -> Colon(), M))

    size(payoffs)[end] != N && throw(ArgumentError(
        "length of the array in the last axis must be equal to
         the number of players"
    ))

    players = ntuple(
        i -> Player(permutedims(view(payoffs, colons..., i),
                                     (i:N..., 1:i-1...)::typeof(dims))
                   ),
        Val(N)
    )
    NormalFormGame(players)
end

"""
    NormalFormGame(payoffs)

Construct a symmetric 2-player NormalFormGame with a square matrix.

# Arguments

- `payoffs::Matrix{T<:Real}` : Square matrix representing each player's payoff
  matrix.
"""
function NormalFormGame(payoffs::Matrix{T}) where T<:Real
    n, m = size(payoffs)
    n != m && throw(ArgumentError(
        "symmetric two-player game must be represented by a square matrix"
    ))
    player = Player(payoffs)
    return NormalFormGame(player, player)
end

function NormalFormGame{N,T}(g::NormalFormGame{N,S}) where {N,T,S}
    players_new = ntuple(i -> Player{N,T}(g.players[i]), Val(N))
    return NormalFormGame(players_new)
end

Base.convert(::Type{T}, g::NormalFormGame) where {T<:NormalFormGame} =
    g isa T ? g : T(g)

Base.summary(g::NormalFormGame) =
    string(Base.dims2string(g.nums_actions),
           " ",
           split(string(typeof(g)), ".")[end])

# delete_action

"""
    delete_action(g, action, player_idx)

Return a new `NormalFormGame` instance with the action(s) specified by `action`
deleted from the action set of the player specified by `player_idx`.

# Arguments

- `g::NormalFormGame` : `NormalFormGame` instance.
- `action::Union{PureAction, AbstractVector{<:PureAction}}` : The action(s) to
  be deleted.
- `player_idx::Integer` : Index of the player to delete action(s) for.

# Returns

- `::NormalFormGame` : `NormalFormGame` instance with the action(s) deleted as
  specified.
"""
function delete_action(g::NormalFormGame{N},
                       action::AbstractVector{<:PureAction},
                       player_idx::Integer) where N
    players_new  = [delete_action(player, action,
                    player_idx-i+1>0 ? player_idx-i+1 : player_idx-i+1+N)
                    for (i, player) in enumerate(g.players)]
    return NormalFormGame(players_new)
end

delete_action(g::NormalFormGame, action::PureAction, player_idx::Integer) =
    delete_action(g, [action], player_idx)

# TODO: add printout of payoff arrays
function Base.show(io::IO, g::NormalFormGame)
    print(io, summary(g))
end

function Base.getindex(g::NormalFormGame{N,T},
                       index::Integer...) where {N,T}
    length(index) != N &&
        throw(DimensionMismatch("index must be of length $N"))

    payoff_profile = Array{T}(undef, N)
    for (i, player) in enumerate(g.players)
        payoff_profile[i] =
            player.payoff_array[(index[i:end]..., index[1:i-1]...)...]
    end
    return payoff_profile
end

# Trivial game with 1 player
function Base.getindex(g::NormalFormGame{1,T}, index::Integer) where T
    return g.players[1].payoff_array[index]
end

function Base.setindex!(g::NormalFormGame{N,T},
                        payoff_profile::Vector{S},
                        index::Integer...) where {N,T,S<:Real}
    length(index) != N &&
        throw(DimensionMismatch("index must be of length $N"))
    length(payoff_profile) != N &&
        throw(DimensionMismatch("assignment must be of $N-array"))

    for (i, player) in enumerate(g.players)
        player.payoff_array[(index[i:end]...,
                             index[1:i-1]...)...] = payoff_profile[i]
    end
    return payoff_profile
end

# Trivial game with 1 player
function Base.setindex!(g::NormalFormGame{1,T},
                        payoff::S,
                        index::Integer) where {T,S<:Real}
    g.players[1].payoff_array[index] = payoff
    return payoff
end

# Indexing with CartesianIndices
Base.getindex(g::NormalFormGame{N}, ci::CartesianIndex{N}) where {N} =
    g[to_indices(g, (ci,))...]
Base.setindex!(g::NormalFormGame{N}, v, ci::CartesianIndex{N}) where {N} =
    g[to_indices(g, (ci,))...] = v

# is_nash

function is_nash(g::NormalFormGame, action_profile::ActionProfile;
                 tol::Real=1e-8)
    for (i, player) in enumerate(g.players)
        own_action = action_profile[i]
        opponents_actions =
            tuple(action_profile[i+1:end]..., action_profile[1:i-1]...)
        if !(is_best_response(player, own_action, opponents_actions, tol=tol))
            return false
        end
    end
    return true
end

function is_nash(g::NormalFormGame{2}, action_profile::ActionProfile;
                 tol::Real=1e-8)
    for (i, player) in enumerate(g.players)
        own_action, opponent_action =
            action_profile[i], action_profile[3-i]
        if !(is_best_response(player, own_action, opponent_action, tol=tol))
            return false
        end
    end
    return true
end

# attach docstring for N>1 player games
@doc """
    is_nash(g, action_profile; tol=1e-8)

Return true if `action_profile` is a Nash equilibrium.

# Arguments

- `g::NormalFormGame` : Instance of N-player NormalFormGame.
- `action_profile::ActionProfile` : Tuple of N integers (pure actions) or N
  vectors of reals (mixed actions).
- `tol::Real` : Tolerance to be used to determine best response actions.

# Returns

- `::Bool`
""" is_nash

# Trivial game with 1 player
"""
    is_nash(g, action; tol=1e-8)

Return true if `action` is a Nash equilibrium of a trivial game with 1 player.

# Arguments

- `g::NormalFormGame{1}` : Instance of 1-player NormalFormGame.
- `action::Action` : Integer (pure action) or vector of reals (mixed action).
- `tol::Float64` : Tolerance to be used to determine best response actions.

# Returns

- `::Bool`
"""
is_nash(g::NormalFormGame{1}, action::Action; tol::Real=1e-8) =
    is_best_response(g.players[1], action, nothing, tol=tol)

is_nash(g::NormalFormGame{1}, action_profile::ActionProfile;
        tol::Real=1e-8) = is_nash(g, action_profile..., tol=tol)

# Utility functions

"""
    pure2mixed(num_actions, action)

Convert a pure action to the corresponding mixed action.

# Arguments

- `num_actions::Integer` : The number of the pure actions (= the length of a
  mixed action).
- `action::PureAction` : The pure action to convert to the corresponding mixed
  action.

# Returns

- `mixed_action::Vector{Float64}` : The mixed action representation of the given
  pure action.
"""
function pure2mixed(num_actions::Integer, action::PureAction)
    mixed_action = zeros(num_actions)
    mixed_action[action] = 1
    return mixed_action
end

# is_pareto_efficient & is_pareto_dominant

function pareto_inferior_to(payoff_profile1, payoff_profile2)
    all(payoff_profile2 .>= payoff_profile1) &&
    any(payoff_profile2 .> payoff_profile1)
end

function not_pareto_superior_to(payoff_profile1, payoff_profile2)
    any(payoff_profile2 .> payoff_profile1) ||
    all(payoff_profile2 .== payoff_profile1)
end

for (f, op) = ((:is_pareto_efficient, pareto_inferior_to),
               (:is_pareto_dominant, not_pareto_superior_to))
    @eval function $(f)(g::NormalFormGame,
                        action_profile::PureActionProfile)
        payoff_profile0 = g[action_profile...]
        for profile in CartesianIndices(g.nums_actions)
            if CartesianIndex(action_profile) != profile
                if ($(op)(payoff_profile0, g[profile]))
                    return false
                end
            end
        end
        return true
    end
end

@doc """
    is_pareto_efficient(g, action_profile)

Return true if `action_profile` is Pareto efficient for game `g`.

# Arguments

- `g::NormalFormGame` : Instance of N-player NormalFormGame.
- `action_profile::PureActionProfile` : Tuple of N integers (pure actions).

# Returns

- `::Bool`
""" is_pareto_efficient

@doc """
    is_pareto_dominant(g, action_profile)

Return true if `action_profile` is Pareto dominant for game `g`.

# Arguments

- `g::NormalFormGame` : Instance of N-player NormalFormGame.
- `action_profile::PureActionProfile` : Tuple of N integers (pure actions).

# Returns

- `::Bool`
""" is_pareto_dominant

# is_dominated

"""
    is_dominated(player, action; tol=1e-8,
                 lp_solver=() -> Clp.Optimizer(LogLevel=0))

Determine whether `action` is strictly dominated by some mixed action.

# Arguments

- `player::Player` : Player instance.
- `action::PureAction` : Integer representing a pure action.
- `tol::Real` : Tolerance level used in determining domination.
- `lp_solver::Union{Type{<:MathOptInterface.AbstractOptimizer},Function}` :
  Linear programming solver to be used internally. Pass a
  `MathOptInterface.AbstractOptimizer` type (such as `Clp.Optimizer`) if no
  option is needed, or a function (such as `() -> Clp.Optimizer(LogLevel=0)`) to
  supply options.

# Returns

- `::Bool` : True if `action` is strictly dominated by some mixed action; false
  otherwise.

"""
function is_dominated(
    ::Type{T}, player::Player, action::PureAction; tol::Real=1e-8,
    lp_solver::Union{Type{TO},Function}=() -> Clp.Optimizer(LogLevel=0)
) where {T<:Real,TO<:MOI.AbstractOptimizer}
    payoff_array = player.payoff_array
    m, n = size(payoff_array, 1) - 1, prod(size(payoff_array)[2:end])

    ind = trues(num_actions(player))
    ind[action] = false

    A_ub = Matrix{T}(undef, (m+1, n))  # transposed
    A_ub[1:end-1, :] .= reshape(selectdim(payoff_array, 1, action), (1, n))
    A_ub[1:end-1, :] -= reshape(selectdim(payoff_array, 1, ind), (m, n))
    A_ub[end, :] .= 1

    a_eq = ones(T, m+1)
    a_eq[end] = 0

    c = zeros(T, m+1)
    c[end] = 1

    CACHE = MOIU.UniversalFallback(MOIU.Model{Float64}())
    optimizer = MOIU.CachingOptimizer(CACHE, lp_solver())
    x = MOI.add_variables(optimizer, m+1)
    MOI.set(optimizer, MOI.ObjectiveFunction{MOI.ScalarAffineFunction{T}}(),
            MOI.ScalarAffineFunction{T}(MOI.ScalarAffineTerm{T}.(c, x), 0))
    MOI.set(optimizer, MOI.ObjectiveSense(), MOI.MAX_SENSE)
    for j in 1:n
        MOI.add_constraint(
            optimizer,
            MOI.ScalarAffineFunction{T}(
                MOI.ScalarAffineTerm{T}.(A_ub[:, j], x), 0
            ),
            MOI.LessThan{T}(0)
        )
    end
    MOI.add_constraint(
        optimizer,
        MOI.ScalarAffineFunction{T}(MOI.ScalarAffineTerm{T}.(a_eq, x), 0),
        MOI.EqualTo{T}(1)
    )
    # Nonnegativity
    for i in 1:m
        a = zeros(T, m+1)
        a[i] = -1
        MOI.add_constraint(
            optimizer,
            MOI.ScalarAffineFunction{T}(MOI.ScalarAffineTerm{T}.(a, x), 0),
            MOI.LessThan{T}(0)
        )
    end
    MOI.optimize!(optimizer)
    status = MOI.get(optimizer, MOI.TerminationStatus())

    if status == MOI.OPTIMAL
        return (MOI.get(optimizer, MOI.ObjectiveValue()) > tol)::Bool
    elseif status == MOI.INFEASIBLE
        return false
    else
        throw(ErrorException("Error: solution status $(status)"))
    end
end

function is_dominated(
    ::Type{T}, player::Player{1}, action::PureAction; tol::Real=1e-8,
    lp_solver::Union{Type{TO},Function}=() -> Clp.Optimizer(LogLevel=0)
) where {T<:Real,TO<:MOI.AbstractOptimizer}
        payoff_array = player.payoff_array
        return maximum(payoff_array) > payoff_array[action] + tol
end

is_dominated(
    player::Player, action::PureAction; tol::Real=1e-8,
    lp_solver::Union{Type{TO},Function}=() -> Clp.Optimizer(LogLevel=0)
) where {TO<:MOI.AbstractOptimizer} =
    is_dominated(Float64, player, action, tol=tol, lp_solver=lp_solver)

# dominated_actions

"""
    dominated_actions(player; tol=1e-8,
                      lp_solver=() -> Clp.Optimizer(LogLevel=0))

Return a vector of actions that are strictly dominated by some mixed actions.

# Arguments

- `player::Player` : Player instance.
- `tol::Real` : Tolerance level used in determining domination.
- `lp_solver::Union{Type{<:MathOptInterface.AbstractOptimizer},Function}` :
  Linear programming solver to be used internally. Pass a
  `MathOptInterface.AbstractOptimizer` type (such as `Clp.Optimizer`) if no
  option is needed, or a function (such as `() -> Clp.Optimizer(LogLevel=0)`) to
  supply options.

# Returns

- `out::Vector{Int}` : Vector of integers representing pure actions, each
  of which is strictly dominated by some mixed action.

"""
function dominated_actions(
    ::Type{T}, player::Player; tol::Real=1e-8,
    lp_solver::Union{Type{TO},Function}=() -> Clp.Optimizer(LogLevel=0)
) where {T<:Real,TO<:MOI.AbstractOptimizer}
    out = Int[]
    for action = 1:num_actions(player)
        if is_dominated(T, player, action, tol=tol, lp_solver=lp_solver)
            append!(out, action);
        end
    end

    return out
end

dominated_actions(
    player::Player; tol::Real=1e-8,
    lp_solver::Union{Type{TO},Function}=() -> Clp.Optimizer(LogLevel=0)
) where {TO<:MOI.AbstractOptimizer} =
    dominated_actions(Float64, player, tol=tol, lp_solver=lp_solver)


text:
using HSA
using HSA.ExtFinalization

USE_CODEGEN = false

function check(message)
	println("✓ $message")
end

rt = Runtime()
check("Initializing the hsa runtime")

gpus = HSA.all_agents(dev = HSA.DeviceTypeGpu)
if isempty(gpus)
	println("No GPU Agent found")
	exit(1)
end
check("Getting a gpu agent")

agent = gpus[1]

name = HSA.agent_info_name(agent)
check("Querying the agent name")
println("The agent name is $name.")

queue_size = HSA.agent_info_queue_max_size(agent)
check("Querying the agent maximum queue size")
println("The maximum queue size is $queue_size.")

queue = Queue(agent, queue_size, typ = HSA.QueueTypeSingle, error_callback = (s, q) -> begin
    str = HSA.status_string(s)
	println("($s) : $str")
end)
check("Creating the queue")

brig_ptr = C_NULL

if !USE_CODEGEN
	# Load the precompiled .brig file containing the
	# kernel to be executed

    script_path = dirname(@__FILE__)
    brig_path = joinpath(script_path, "vector_copy.brig")
	mod_bytes = open(readbytes, brig_path)

	brig_ptr = pointer(mod_bytes)
else # USE_CODEGEN
	# Use the HSAIL Code Generator to compile a
	# kernel function to BRIG
	import HSA.Intrinsics

	@hsa_kernel function vector_copy_kernel(a::Ptr{Int64},b::Ptr{Int64})
		idx = get_global_id(Int32(0)) + 1

		x = Base.unsafe_load(b, idx)
        Base.unsafe_store!(a, x, idx)
		return nothing
	end

	println(HSA.Compilation.src_hsail(vector_copy_kernel, Tuple{Ptr{Int64}, Ptr{Int64}}))

	brig_ptr = HSA.Compilation.brig(vector_copy_kernel, Tuple{Ptr{Int64}, Ptr{Int64}})
    assert(brig_ptr != C_NULL)
	check("Compile Kernel to BRIG")
end

program = Program()
check("Create the program")

HSA.ExtFinalization.program_add_module(program, brig_ptr)
check("Adding the brig module to the program")

isa = HSA.agent_info_isa(agent)
check("Query the agents isa")

code_object = HSA.ExtFinalization.program_finalize(program, isa, 0)
check("Finalizing the program")

finalize(program)
check("Destroying the program")

executable = Executable()
check("Create the executable")

HSA.executable_load_code_object(executable, agent, code_object)
check("Loading the code object")

HSA.executable_freeze(executable)
check("Freeze the executable")

println("Symbols in the executable:")
symbols = HSA.symbols(executable)
last_symbol_name = ""
last_symbol = ""
for (s, n) in symbols
	m = HSA.executable_symbol_info_module_name(s)
	println("$m :: $n")

	last_symbol_name = n
	last_symbol = s
end

symbol = last_symbol
try
	symbol = HSA.executable_get_symbol(executable,last_symbol_name)
catch ex
	warn("expected kernel symbol '$last_symbol_name' not found, using any symbol in the executable\n Exception: $ex")

end
check("Extract the symbol from the executable")

kernel_object = HSA.executable_symbol_info_kernel_object(symbol)
check("Extracting the symbol from the executable")
kernarg_segment_size = HSA.executable_symbol_info_kernel_kernarg_segment_size(symbol)
check("Extracting the kernarg segment size from the executable")
group_segment_size = HSA.executable_symbol_info_kernel_group_segment_size(symbol)
check("Extracting the group segment size from the executable")
private_segment_size = HSA.executable_symbol_info_kernel_private_segment_size(symbol)
check("Extracting the private segment from the executable")

signal = Signal(value = 1)
check("Creating a HSA signal")
if HSA.load(signal) != UInt64(1)
	error("incorrect signal value")
end
check("Verify initial signal value")

const N = 16

a_in = Array(Int, N)
fill!(a_in, 1)
HSA.memory_register(a_in, sizeof(a_in))
check("Registering argument memory for input parameter")

b_out = Array(Int, N)
fill!(b_out, 0)
HSA.memory_register(b_out, sizeof(b_out))
check("Registering argument memory for output parameter")

regions = HSA.regions(agent)

karg_region_idx = findfirst(r -> begin
    segment = HSA.region_info_segment(r)
    if (HSA.HSA_REGION_SEGMENT_GLOBAL != segment)
        return false
    end

    flags = HSA.region_info_global_flags(r)
    if (flags & HSA.HSA_REGION_GLOBAL_FLAG_KERNARG) != 0
        return true
	end

    return false
end, regions)

if karg_region_idx == 0
	error("No Kernarg Region found")
end
kernarg_region = regions[karg_region_idx]
check("Finding a kernarg memory region")

alloc = HSA.memory_allocate(kernarg_region, kernarg_segment_size)
kernarg_address = convert(Ptr{Ptr{Int}}, alloc.ptr)
check("Allocating kernel argument memory buffer of size $kernarg_segment_size")
unsafe_store!(kernarg_address, pointer(a_in))
unsafe_store!(kernarg_address + sizeof(pointer(a_in)), pointer(b_out))

index = HSA.load_write_index(queue) # HSA.add_write_index!(queue, UInt64(1))

dispatch_packet = KernelDispatchPacket(kernel_object,N;
    kernarg_address = kernarg_address,
    private_segment_size = private_segment_size,
    group_segment_size = group_segment_size,
    completion_signal = signal
)
dispatch_packet.header.acquire_fence_scope = HSA.HSA_FENCE_SCOPE_SYSTEM
dispatch_packet.header.release_fence_scope = HSA.HSA_FENCE_SCOPE_SYSTEM

pkg_bytes = Array(UInt8, 64)
unsafe_store!(convert(Ptr{Void}, pointer(pkg_bytes)), dispatch_packet)

println("""
    Header: $(pkg_bytes[1:2])
    Setup: $(pkg_bytes[3:4])
    Workgroup X: $(pkg_bytes[5:6])
	""")

queue[index] = dispatch_packet

HSA.store_write_index!(queue, UInt64(index + 1))
HSA.store!(queue.doorbell_signal, Int64(index))
check("Dispatching the kernel")

value = 1
while value != 0
value = wait(signal, :(<), 1, wait_state_hint = HSA.HSA_WAIT_STATE_ACTIVE)
check("Waiting for the kernel")
end

failed = 0
for (i in 1:N)
	if (b_out[i]!=a_in[i])
		failed = failed + 1

		if failed > 10
			println("...")
			break
		end

        println("Bad index: $i\n$(a_in[i]) --> $(b_out[i])")
	end
end

if(failed == 0)
	println("Passed validation.")
else
	println("VALIDATION FAILED!\n")
end

finalize(signal)
finalize(executable)
finalize(code_object)
finalize(queue)
finalize(rt)


text:
__precompile__()

module PyDSTool

using PyCall, DataStructures, DiffEqBase, RecipesBase, LinearAlgebra

const ds = PyNULL()

function __init__()
  try
    copy!(ds, pyimport_conda("PyDSTool", "pydstool", "conda-forge"))
    return
  catch err
    if err isa PyCall.PyError
      # A dirty hack to force importing PyDSTool:
      # https://github.com/JuliaDiffEq/PyDSTool.jl/issues/5
      #
      # At this point, we assume that "import PyDSTool" was failed due
      # to the bug in how it checks SciPy version number.  We
      # workaround it by monkey-patching `scipy.__version__` while
      # importing PyDSTool.  We make ./_pydstool_jl_hack.py importable
      # and execute it by importing it.
      pushfirst!(PyVector(pyimport("sys")["path"]), @__DIR__)
      pyimport("_pydstool_jl_hack")
    else
      rethrow()
    end
  end
  copy!(ds, pyimport_conda("PyDSTool", "pydstool", "conda-forge"))
end

include("constants.jl")
include("ode_construct_solve.jl")
include("bifurcation.jl")

export ds, build_args

export PYDSTOOL_CURVE_CLASSES, ALL_POINT_TYPES

export set_name,set_ics,set_pars, set_vars, set_tdata,set_fnspecs,
       set_tdomain, interpert_pts, build_ode, solve_ode,bifurcation_curve

export find_changes
end


text:
using Seismic, PyCall
@pyimport matplotlib.pyplot as pl

#Generate seismic data

nt = 256;
nx = convert(Int64,128);
h = Array{Header}(nx);
for ix = 1:nx
  h[ix] = Seismic.InitSeisHeader();
  h[ix].tracenum = ix;
  h[ix].d1 = 0.004;
  h[ix].n1 = nt;
  h[ix].imx = ix;
  h[ix].imy = 0;
  h[ix].ihx = 0;
  h[ix].ihy = 0;
end

dt=0.004;
dx=5;
f_niquest=Int(floor(1/dt))

Tmax=(nt-1)*dt

data,ext=SeisLinearEvents(nt=nt,dt=dt,nx1=nx,dx1=dx,tau=[Tmax/4., Tmax/3])#,p1=[0.0001])

#Transform to FK domain
m=fft(data,1)/sqrt(size(data,1));
m=fft(m,2)/sqrt(size(data,2))
m=fftshift(m)

m=m[Int(floor(f_niquest/2)):f_niquest,:]

# apply FK filter
data_fk= SeisFKFilter(data;dt=dt,dx=dx,va=-1500,vb=-2500,vc=2500,vd=1500)

#Transform back to TX domain.

m=fft(data_fk,1)/sqrt(size(data_fk,1));
m=fft(m,2)/sqrt(size(data_fk,2))
m=fftshift(m)

m=m[Int(floor(f_niquest/2)):f_niquest,:]

#>>>>>>>>>>>>>>>Plot<<<<<<<<<<<<<<<<

f1 = pl.figure(1)
pl.imshow(data[1:nt,:],cmap="PuOr",aspect="auto",vmin=-1,vmax=1)
pl.title("data before fk filter")
pl.xlabel("X (m)")
pl.ylabel("T (s)")
pl.show()

f2 = pl.figure(2)
pl.imshow(abs.(m),cmap="PuOr",aspect="auto",vmin=-1,vmax=1)
pl.title("fk spectrum before filter")
pl.xlabel("K (1/m)")
pl.ylabel("W (1/s)")
pl.show()

f3 = pl.figure(3)
pl.imshow(data_fk[1:nt,:],cmap="PuOr",aspect="auto",vmin=-1,vmax=1)
pl.title("data after fk filter")
pl.xlabel("X (m)")
pl.ylabel("T (s)")
pl.show()

f4 = pl.figure(4)
pl.imshow(abs.(m),cmap="PuOr",aspect="auto",vmin=-1,vmax=1)
pl.title("fk spectrum after filter")
pl.xlabel("K (1/m)")
pl.ylabel("W (1/s)")
pl.show()

pl.show()

max_stars_repo_name: fercarozzi/myseismicjulia

max_issues_repo_name: fercarozzi/myseismicjulia

max_forks_repo_name: fercarozzi/myseismicjulia


text:
using DecisionMakingPolicies
using Test

import Zygote: gradient, pullback

@testset "Stateless Softmax Tests" begin
    num_actions = 2
    T = Float32
    p = StatelessSoftmax(T, num_actions)
    @test length(p.θ) == num_actions
    @test eltype(p.θ) == T
    @test size(p.θ) == (num_actions,)

    a = rand(p())
    @test typeof(a) <: Int
    logp1 = logpdf(p(), 1)
    logp2 = logpdf(p, 1)
    @test typeof(logp1) == T
    @test typeof(logp2) == typeof(logp1)
    @test isapprox(logp1, logp2)
    g = (θ=zero(p.θ),)
    gauto1 = gradient(p->logpdf(p(),1),p)[1].θ
    gauto2 = gradient(p->logpdf(p,1),p)[1].θ

    @test eltype(gauto1) == eltype(gauto2)
    @test eltype(gauto1) == eltype(p.θ)
    @test all(isapprox.(gauto1, gauto2))

    logp3, _ = grad_logpdf!(g, p, 1)
    @test typeof(logp3) == T
    @test logp3 == logp2
    @test all(isapprox.(g.θ, gauto1))
    g2 = Array{T,1}([0.5, -0.5])
    @test all(isapprox.(g.θ, g2))

    A = zeros(Int,1)
    A2 = zeros(Int,2)
    _, logp4, _ = sample_with_trace!(g, A, p)
    @test A[1] > 0
    if A[1] == 1
        @test all(isapprox.(g.θ, g2))
    else
        @test all(isapprox.(g.θ, -g2))
    end
    sample_with_trace!(g, view(A2,1:1), p)
    @test A2[1] > 0
    @test A2[2] == 0
    a1 = A2[1]
    sample_with_trace!(g, view(A2,2:2), p)
    @test A2[1] == a1
    @test A2[2] > 0
end

@testset "Stateless Softmax Buffer Tests" begin
    num_actions = 2
    T = Float32
    p = StatelessSoftmax(T, num_actions)
    buff = SoftmaxBuffer(p)

    @test eltype(buff.action) == Int
    @test eltype(buff.p) == T
    @test eltype(buff.ψ[1]) == T
    @test size(buff.ψ[1]) == size(p.θ)
    @test size(buff.p) == size(p.θ)

    a = rand(p(buff))
    @test typeof(a) <: Int
    logp1 = logpdf(p(), 1)
    logp2 = logpdf(p(buff), 1)
    logp3 = logpdf!(buff, p, 1)
    @test typeof(logp1) == T
    @test typeof(logp2) == typeof(logp1)
    @test typeof(logp3) == typeof(logp1)
    @test isapprox(logp1, logp2)
    @test isapprox(logp1, logp3)

    gauto1 = gradient(p->logpdf(p,1),p)[1].θ
    logp3, _ = grad_logpdf!(buff, p, 1)
    @test typeof(logp3) == T
    @test logp3 == logp1
    @test all(isapprox.(buff.ψ[1], gauto1))
    g2 = Array{T,1}([0.5, -0.5])
    @test all(isapprox.(buff.ψ[1], g2))

    _, logp4, _ = sample_with_trace!(buff, p)
    @test buff.action[1] > 0
    if buff.action[1] == 1
        @test all(isapprox.(buff.ψ[1], g2))
    else
        @test all(isapprox.(buff.ψ[1], -g2))
    end
end

@testset "Linear Softmax Tests" begin
    num_actions = 2
    num_features = 3
    T = Float32
    p = LinearSoftmax(T, num_features, num_actions)

    @test eltype(p.θ) == T
    @test size(p.θ) == (num_features, num_actions)
    s = Array{T, 1}([0.0, 1.0, 2.0])
    a = rand(p(s))
    @test typeof(a) <: Int
    logp1 = logpdf(p(s), 1)
    logp2 = logpdf(p, s, 1)
    @test typeof(logp1) == T
    @test typeof(logp2) == typeof(logp1)
    @test isapprox(logp1, logp2)
    g = (θ=zero(p.θ),)

    gauto1 = gradient(p->logpdf(p(s),1),p)[1].θ
    gauto2 = gradient(p->logpdf(p,s,1),p)[1].θ

    @test eltype(gauto1) == eltype(gauto2)
    @test eltype(gauto1) == eltype(p.θ)
    @test all(isapprox.(gauto1, gauto2))

    logp3, _ = grad_logpdf!(g, p, s, 1)
    @test typeof(logp3) == T
    @test logp3 == logp2
    @test all(isapprox.(g.θ, gauto1))
    g2 = Array{T,2}([0.0 0.0; 0.5 -0.5; 1.0 -1.0])
    @test all(isapprox.(g.θ, g2))

    k = 3
    sb = reshape(collect(T, 1:num_features*k), (num_features, k))
    ab = rand(1:num_actions, (1, k))
    gb = zero(p.θ)
    logps = []
    for i in 1:k
        logp, ψ = grad_logpdf!(g, p, sb[:, i], ab[1, i])
        push!(logps, logp)
        @. gb += ψ[1]
    end
    zlogps, gbf = pullback(p->(logpdf(p,sb,ab)),p)

    A = zeros(Int,1)
    A2 = zeros(Int,2)
    _, logp4, _ = sample_with_trace!(g, A, p, s)
    @test A[1] > 0
    if A[1] == 1
        @test all(isapprox.(g.θ, g2))
    else
        @test all(isapprox.(g.θ, -g2))
    end
    sample_with_trace!(g, view(A2,1:1), p, s)
    @test A2[1] > 0
    @test A2[2] == 0
    a1 = A2[1]
    sample_with_trace!(g, view(A2,2:2), p, s)
    @test A2[1] == a1
    @test A2[2] > 0
end

@testset "Linear Softmax Buffer Tests" begin
    num_actions = 2
    num_features = 3
    T = Float32
    p = LinearSoftmax(T, num_features, num_actions)
    buff = SoftmaxBuffer(p)
    @test eltype(buff.action) == Int
    @test eltype(buff.p) == T
    @test eltype(buff.ψ[1]) == T
    @test size(buff.ψ[1]) == size(p.θ)
    @test size(buff.p) == (num_actions,)

    s = Array{T, 1}([0.0, 1.0, 2.0])
    a = rand(p(s))
    @test typeof(a) <: Int
    logp1 = logpdf(p(s), 1)
    logp2 = logpdf(p(buff, s), 1)
    logp3 = logpdf!(buff, p, s, 1)
    @test typeof(logp1) == T
    @test typeof(logp2) == typeof(logp1)
    @test typeof(logp3) == typeof(logp1)
    @test isapprox(logp1, logp2)
    @test isapprox(logp1, logp3)

    gauto1 = gradient(p->logpdf(p(s),1),p)[1].θ
    logp3, _ = grad_logpdf!(buff, p, s, 1)
    @test typeof(logp3) == T
    @test logp3 == logp2
    @test all(isapprox.(buff.ψ[1], gauto1))
    g2 = Array{T,2}([0.0 0.0; 0.5 -0.5; 1.0 -1.0])
    @test all(isapprox.(buff.ψ[1], g2))

    _, logp4, _ = sample_with_trace!(buff, p, s)
    @test buff.action[1] > 0
    if buff.action[1] == 1
        @test all(isapprox.(buff.ψ[1], g2))
    else
        @test all(isapprox.(buff.ψ[1], -g2))
    end
end


text:
using MAGMA
import LinearAlgebra.LAPACK: gels!, gesv!, getrs!, getri!
using Test
using CuArrays
using TimerOutputs

T = ComplexF64

# randomly generate a 2 by 2 matrix for testing
size_test = 1024

A = rand(T, size_test, size_test)
B = rand(T, size_test, size_test)

dA = cu(A)
dB = cu(B)

A_copy = copy(A)
B_copy = copy(B)

A_test = dA
B_test = dB

right_answer = gesv!(A, B)

magma_init()
result = magma_gesv!(A_test, B_test)
magma_finalize()

for i in 1:length(result)
    @test Array(result[i]) ≈ right_answer[i]
end


text:
using Pkg
pkg"activate ."
using NNSparseCoding; const N = NNSparseCoding;
using Test
using Plots;
theme(:dark);
using BenchmarkTools
using StatsBase,Statistics
using MultivariateStats

function scatty(M1,M2)
  scatter(M1[:], M2[:]; leg = false, ratio = 1)
  return plot!([extrema(M1)...].* 1.1, identity; linewidth = 3)
end
# function to plot one or more black/mytransfe images
function showbw(mat::Matrix{<:Real},_title="")
  heatmap(mat, color=:grays,ratio=1,cbar=false,axis=false,title=_title)
end
function showbw_many(title::String,mats...)
  mats = [mats...] # convert tuple to array
  padval=maximum(maximum,mats)
  framesize = size(mats[1],1)
  padmat = fill(padval, framesize,round(Int,0.1framesize))
  padmats = repeat([padmat,],inner=length(mats))
  matspads = permutedims(hcat(mats,padmats))[:]
  imgrow = hcat(matspads[1:end-1]...)
  showbw(imgrow , title)
end
##
#
# const dirimg = abspath(@__DIR__,"..","data")
# const sizesq=315
# all_images = N.read_natural_images(dirimg,sizesq ; verbose=true)
#
# const npatches = 1000
# const patch_size = 40
#
# patches = N.get_random_patches(npatches,all_images,patch_size)
#
#
# patchmat = N.merge_mats(patches)
#
# mean(patchmat;dims=2) |> extrema
# std(patchmat;dims=2) |> extrema
# cov(patchmat;dims=2) |> heatmap
#
#
# ##
# showbw_many("done!",sample(patches,5)...)
##

n = 80
d = n
k = 500
A,S = N.initialize_mats(n,d,k)
map!(s-> rand() < 0.95 ? 0.0 : 3s ,S,S)
X = A*S
broadcast!(-,X,X,mean(X;dims=2))

##
mytransf,Mp,Mm = N.preprocess_plus_minus(X,N.PreprocessPCA(d))

##

lambda=0.1
nsteps=10_000
stepsize=0.2
optthresh=1E-4
op = N.nnsc_optimize(Mp,d,lambda ;
  step_zero=stepsize, nstepmax=nsteps,
  opt_threshold=optthresh,verbose=true)

om = N.nnsc_optimize(Mm,d,lambda ;
  step_zero=stepsize, nstepmax=nsteps,
  opt_threshold=optthresh,verbose=true)

##

scatty(X,N.revert_plus_minus(op.AS,om.AS,mytransf))

extrema( X .- N.revert_plus_minus(op.AS,om.AS,mytransf))

###
mytransf,Mp,Mm = N.preprocess_plus_minus(X,N.PreprocessWhiten)

##

lambda=0.1
nsteps=10_000
stepsize=0.2
optthresh=1E-4
op = N.nnsc_optimize(Mp,d,lambda ;
  step_zero=stepsize, nstepmax=nsteps,
  opt_threshold=optthresh,verbose=true)

om = N.nnsc_optimize(Mm,d,lambda ;
  step_zero=stepsize, nstepmax=nsteps,
  opt_threshold=optthresh,verbose=true)

##

scatty(X,N.revert_plus_minus(op.AS,om.AS,mytransf))

extrema( X .- N.revert_plus_minus(op.AS,om.AS,mytransf))

boh = N.revert_plus_minus(op.S,zero(op.S),mytransf)


text:
# SPECTplan.jl

using SPECTrecon: SPECTplan
using BenchmarkTools: @btime
using MATLAB

function call_SPECTplan_matlab(mpath, mumap, psfs, dy)

    mat"""
    addpath($mpath)
    SPECTplan_matlab($mumap, $psfs, $dy);
    """
end

function SPECTplan_time()
    T = Float32
    nx = 64
    ny = 64
    nz = 40
    nview = 60

    mumap = rand(T, nx, ny, nz)

    nx_psf = 19
    nz_psf = 19
    dy = T(4.7952)

    psfs = rand(T, nx_psf, nz_psf, ny, nview)
    psfs = psfs .+ mapslices(reverse, psfs, dims = [1, 2])
    psfs = psfs .+ mapslices(transpose, psfs, dims = [1, 2])
    psfs = psfs ./ mapslices(sum, psfs, dims = [1, 2])

    plan = SPECTplan(mumap, psfs, dy)

    println("SPECTplanfast")
    @btime plan = SPECTplan($mumap, $psfs, $dy; mode = :fast)
    # 13.585 ms (97843 allocations: 25.42 MiB)

    println("SPECTplanmem")
    @btime plan = SPECTplan($mumap, $psfs, $dy; mode = :mem)
    # 13.578 ms (97789 allocations: 12.30 MiB)

    mpath = pwd()
    println("SPECTplan_matlab")
    println("Warning: Check if MIRT is installed")
    call_SPECTplan_matlab(mpath, mumap, psfs, dy)
    # 1.103891 seconds, real memory size: 1.30 GB
    nothing
end

# run all functions, time may vary on different machines, will alllocate ~100 MB memory.
SPECTplan_time()


text:
## Advent of Code 2021 ##
#########################

## Day 01: Sonar Sweep ##
#########################

day1_test = [199, 200, 208, 210, 200, 207, 240, 269, 260, 263]
day1_input = parse.(Int, readlines(open("input-01")))

"""Given a Vector of Int depth values, this counts how often it decreases when
going from position i to i+step_size."""
function depth_count(depth_array::Vector{Int}, step_size::Int = 1)
    counter = 0
    for i = 1:(length(depth_array)-step_size)
        if depth_array[i] < depth_array[i+step_size]
            counter = counter + 1
        end
    end
    counter
end

@assert 0 == depth_count(day1_test, 100)
@assert 7 == @show depth_count(day1_test)
@assert 1532 == @show depth_count(day1_input)

@assert 5 == @show depth_count(day1_test, 3)
@assert 1571 == @show depth_count(day1_input, 3)

## Day 02: Dive! ##
###################

day02_test = ["forward 5", "down 5", "forward 8", "up 3", "down 8", "forward 2"]
day02_input = readlines(open("input-02"))

"""Tries to parse the instruction as "direction distance" code and crashes if
that is not possible."""
function parse_swim_instruction(instruction::String)
    m = match(r"^(forward|down|up) ([0-9]+)$", instruction)
    (m[1], parse(Int, m[2]))
end

"""Given a vector of swim instructions, this method returns the final position and aim value."""
function swim(instructions::Vector{String})::Tuple{Int,Int,Int}
    aim = 0
    x = 0
    y = 0
    for instruction in instructions
        dir, t = parse_swim_instruction(instruction)
        if dir == "forward"
            x = x + t
            y = y + aim * t
        elseif dir == "down"
            aim = aim + t
        elseif dir == "up"
            aim = aim - t
        else
            throw("Direction $dir not recognized")
        end
    end
    (x, y, aim)
end

legacy_position_product(t::Tuple) = t[1] * t[3] # x * aim
position_product(t::Tuple) = t[1] * t[2] # x * y

@assert 150 == @show legacy_position_product(swim(day02_test))
@assert 2019945 == @show legacy_position_product(swim(day02_input))
@assert 900 == @show position_product(swim(day02_test))
@assert 1599311480 == @show position_product(swim(day02_input))

## Day 03: Binary Diagnostic ##
###############################

day03_test = ["00100", "11110", "10110", "10111", "10101", "01111", "00111", "11100", "10000", "11001", "00010", "01010"]
day03_input = readlines(open("input-03"))

"""Given a vector of bit strings "00110", this returns the most common character
at the given index. If they occur the same amount, returns 1."""
function most_common_digit(diagnostics::Vector{String}, index::Int)::Int
    ones_count = 0
    zeros_count = 0
    for diagnostic in diagnostics
        char = diagnostic[index]
        if char == '1'
            ones_count += 1
        elseif char == '0'
            zeros_count += 1
        else
            throw("Illegal character $(diagnostic[i]) in diagnostic code")
        end
    end
    ones_count >= zeros_count ? 1 : 0
end

@assert most_common_digit(day03_test, 1) == 1
@assert most_common_digit(day03_test, 2) == 0
@assert most_common_digit(day03_test, 3) == 1
@assert most_common_digit(day03_test, 4) == 1
@assert most_common_digit(day03_test, 5) == 0

"""Turns [1,0,0,1,0] into 18."""
function parse_binary(vec::Vector)::Int
    result = 0
    for v in vec
        result = 2 * result + v
    end
    result
end

@assert parse_binary([]) == 0
@assert parse_binary([1, 0, 0, 1, 0]) == 18

function power_consumption(diagnostics::Vector{String})
    gamma = map(i -> most_common_digit(diagnostics, i), 1:length(diagnostics[1]))
    epsilon = map(x -> 1 - x, gamma) # invert digits
    parse_binary(gamma) * parse_binary(epsilon)
end

@assert 198 == @show power_consumption(day03_test)
@assert 3959450 == @show power_consumption(day03_input) # 1565 is not the right answer

function find_rating(diagnostics::Vector{String}, invert = false, i = 1)::Int
    if length(diagnostics) > 1
        most_common = most_common_digit(diagnostics, i)
        most_common = invert ? 1 - most_common : most_common
        filtered_list = filter(d -> parse(Int, d[i]) == most_common, diagnostics)
        find_rating(filtered_list, invert, i + 1)
    elseif length(diagnostics) == 1
        parse(Int, diagnostics[1], base = 2)
    else
        throw("No unique result found")
    end
end

@assert 23 == find_rating(day03_test) # oxygen generator
@assert 10 == find_rating(day03_test, true) # co2 scrubber

life_support_rating(diagnostics) = find_rating(diagnostics) * find_rating(diagnostics, true)

@assert 230 == @show life_support_rating(day03_test)
@assert 7440311 == @show life_support_rating(day03_input)

## Day 04: Giant Squid ##
#########################

day04_test = readlines(open("test-04"))
day04_input = readlines(open("input-04"))

# A board is a Vector with 25 entries.
# A draw state is a Boolean Vector with 25 entries

function bingo(input; find_first = true)
    input = copy(input)
    call_order = parse.(Int, split(popfirst!(input), ","))

    boards = []
    while length(input) > 0
        board = popBingoBoard!(input)
        if length(board) > 0
            push!(boards, board)
        end
    end

    best_duration = find_first ? length(call_order) + 1 : 0
    score_of_best = 0
    for board in boards
        duration, drawn = bingo_duration(call_order, board)
        is_better = (find_first && duration < best_duration) || (!find_first && duration > best_duration)
        if is_better
            best_duration = duration
            score_of_best = 0
            for i = 1:25
                if drawn[i] == 0
                    score_of_best += board[i]
                end
            end
        end
    end
    score_of_best * call_order[best_duration]
end

"""Given a Bingo and a call order, this method marks one call at a time and
checks for a strait line. If there is a line, the marked pattern is returned
alongside the duration it took to get there."""
function bingo_duration(call_order::Vector{Int}, board::Vector{Int})::Tuple{Int,Vector{Int}}
    drawn = zeros(Int, 25)
    call_count = 0
    for call in call_order
        call_count += 1
        call_index = findfirst(x -> x == call, board)
        if call_index !== nothing
            drawn[call_index] = 1
            if has_line(drawn, call_index)
                break
            end
        end
    end
    call_count, drawn
end

"""Checks if there is a horizontal or vertical line through the call index."""
function has_line(drawn::Vector, call_index)
    @assert length(drawn) == 25 "Malformed drawn pattern $drawn."
    @assert 1 <= call_index <= 25 "Call index $call_index out of bounds."
    col = (call_index - 1) % 5 + 1
    row = call_index - col + 1
    if all(drawn[[col, col + 5, col + 10, col + 15, col + 20]] .== 1)
        true # Vertical
    elseif all(drawn[[row, row + 1, row + 2, row + 3, row + 4]] .== 1)
        true # Horizontal
    else
        false
    end
end

"""Pop 25 numbers and return a bingo board. Resistant against empty lines.

Returns `nothing` when there are no Bingos left.
"""
function popBingoBoard!(input)::Union{Vector{Int},Nothing}
    result = []
    while length(result) < 25
        if length(input) == 0
            if length(result) == 0
                return nothing
            else
                throw("Incomplete Bingo: $result")
            end
        end
        line = parse.(Int, split(popfirst!(input)))
        append!(result, line)
    end
    if length(result) == 25
        result
    else
        throw("Malformed Bingo: $result")
    end
end

@assert 4512 == @show bingo(day04_test)
@assert 11536 == @show bingo(day04_input)

@assert 1924 == @show bingo(day04_test; find_first = false)
@assert 1284 == @show bingo(day04_input; find_first = false)

## Day 05: Hydrothermal Venture ##
##################################

day05_test = readlines(open("test-05"))
day05_input = readlines(open("input-05"))

function vent_intersections(data; skip_diagonals = false)
    diagramm = Dict()
    for d in data
        m = match(r"([0-9]+),([0-9]+) -> ([0-9]+),([0-9]+)", d)
        if m[1] == m[3]
            x = parse(Int, m[1])
            for y in r(parse(Int, m[2]), parse(Int, m[4]))
                v = get(diagramm, (x, y), 0)
                diagramm[(x, y)] = v + 1
            end
        elseif m[2] == m[4]
            y = parse(Int, m[2])
            for x in r(parse(Int, m[1]), parse(Int, m[3]))
                v = get(diagramm, (x, y), 0)
                diagramm[(x, y)] = v + 1
            end
        elseif skip_diagonals
            continue
        else
            x1 = parse(Int, m[1])
            y1 = parse(Int, m[2])
            x2 = parse(Int, m[3])
            y2 = parse(Int, m[4])
            dx = x1 < x2 ? 1 : -1
            dy = y1 < y2 ? 1 : -1
            for i = 0:(abs(x2 - x1))
                x = x1 + i * dx
                y = y1 + i * dy
                v = get(diagramm, (x, y), 0)
                diagramm[(x, y)] = v + 1
            end
        end
    end
    intersects = 0
    for (_key, value) in diagramm
        if value > 1
            intersects += 1
        end
    end
    intersects
end

r(a, b) = a < b ? (a:b) : (b:a)

@assert 5 == @show vent_intersections(day05_test; skip_diagonals = true)
@assert 6007 == @show vent_intersections(day05_input; skip_diagonals = true)

@assert 12 == @show vent_intersections(day05_test)
@assert 19349 == @show vent_intersections(day05_input)

## Day 06: Lanternfish ##
#########################

day06_test = [3, 4, 3, 1, 2]
day06_input = parse.(Int, split(read(open("input-06"), String), ","))

function group_fish(fish_list::Vector{Int})::Dict{Int,Int}
    groups = Dict{Int,Int}()
    for fish in fish_list
        groups[fish] = get(groups, fish, 0) + 1
    end
    groups
end

@assert group_fish(day06_test) == Dict(1 => 1, 2 => 1, 3 => 2, 4 => 1)

function pass_time(fish_dict::Dict{Int,Int})::Dict{Int,Int}
    groups = Dict{Int,Int}()
    for (time, count) in fish_dict
        if time == 0
            # Special case where we spawn offspring
            groups[8] = get(groups, 8, 0) + count
            groups[6] = get(groups, 6, 0) + count
        else
            groups[time-1] = get(groups, time - 1, 0) + count
        end
    end
    groups
end

function pass_time(fish_dict::Dict{Int,Int}, repeats::Int)::Dict{Int,Int}
    @assert repeats >= 0 "Passed in a negative repeat count $repeats which is invalid."
    if repeats == 0
        fish_dict
    else
        pass_time(pass_time(fish_dict), repeats - 1)
    end
end

count_fish(fish_dict) = sum(values(fish_dict))

@assert pass_time(group_fish(day06_test)) == Dict(0 => 1, 1 => 1, 2 => 2, 3 => 1)
@assert pass_time(group_fish(day06_test), 2) == Dict(6 => 1, 8 => 1, 0 => 1, 1 => 2, 2 => 1)

@assert 5934 == @show count_fish(pass_time(group_fish(day06_test), 80))
@assert 362639 == @show count_fish(pass_time(group_fish(day06_input), 80))

@assert 26984457539 == @show count_fish(pass_time(group_fish(day06_test), 256))
@assert 1639854996917 == @show count_fish(pass_time(group_fish(day06_input), 256))

## Day 07: The Treachery of Whales ##
#####################################

day07_test = [16, 1, 2, 0, 4, 2, 7, 1, 2, 14]
day07_input = parse.(Int, split(read(open("input-07"), String), ","))

using Statistics # median

function crab_distance(crabs::Vector{Int})
    m = Int(floor(Statistics.median(crabs)))
    sum(abs.(crabs .- m))
end

@assert 37 == @show crab_distance(day07_test)
@assert 356992 == @show crab_distance(day07_input)

function crab_distance_2(crabs::Vector{Int})
    m = mean(crabs)
    fuel1 = sum(square_fuel.(abs.(crabs .- Int(floor(m)))))
    fuel2 = sum(square_fuel.(abs.(crabs .- Int(ceil(m)))))
    min(fuel1, fuel2)
end

square_fuel(d) = Int(d * (d + 1) / 2)

@assert square_fuel(0) == 0
@assert square_fuel(1) == 1
@assert square_fuel(2) == 3
@assert square_fuel(3) == 6

@assert 168 == @show crab_distance_2(day07_test)
@assert 101268110 == @show crab_distance_2(day07_input)

## Day 08: Seven Segment Search ##
##################################

day08_test = readlines(open("test-08"))
day08_input = readlines(open("input-08"))

function easy_7_digit_numbers(data)
    count = 0
    for d in data
        output = split(d, "|")[2]
        out_groups = split(output)
        for group in out_groups
            if length(group) in [2, 3, 4, 7]
                count += 1
            end
        end
    end
    count
end

@assert 26 == @show easy_7_digit_numbers(day08_test)
@assert 344 == @show easy_7_digit_numbers(day08_input)

# Part 2 seems really hard.

#  4444
# 5    1
# 5    1
#  6666
# 7    2
# 7    2
#  3333

segments_to_digit = Dict(
    Set([1, 2, 3, 4, 5, 7]) => 0, Set([1, 2]) => 1, Set([4, 1, 6, 7, 3]) => 2,
    Set([4, 1, 6, 2, 3]) => 3, Set([1, 2, 5, 6]) => 4, Set([4, 5, 6, 2, 3]) => 5,
    Set([2, 3, 4, 5, 6, 7]) => 6, Set([4, 1, 2]) => 7,
    Set([1, 2, 3, 4, 5, 6, 7]) => 8, Set([1, 2, 3, 4, 5, 6]) => 9)

# Seems like we have all 10 numbers in a random order each time. So we should be
# able to solve it the same way every time.

function solve_7_digit_puzzle(data)
    input = split(data, "|")[1]
    in_group = Set.(split(input))
    # Solve
    eight = grab!(in_group, 7)
    one = grab!(in_group, 2)
    seven = grab!(in_group, 3)
    l4 = first(setdiff(seven, one)) # Seven minus one leaves the top letter.
    six = grab_if!(in_group, x -> length(x) == 6 && length(intersect(one, x)) == 1)
    l1 = first(setdiff(one, six))
    l2 = first(intersect(one, six))
    four = grab!(in_group, 4)
    zero = grab_if!(in_group, x -> length(x) == 6 && length(intersect(four, x)) == 3)
    nine = grab_if!(in_group, x -> length(x) == 6 && length(intersect(four, x)) == 4)
    l6 = first(setdiff(eight, zero))
    l7 = first(setdiff(eight, nine))
    l5 = first(setdiff(four, Set([l1, l2, l6])))
    l3 = first(setdiff(eight, Set([l1, l2, l4, l5, l6, l7])))

    mapping = Dict(l1 => 1, l2 => 2, l3 => 3, l4 => 4, l5 => 5, l6 => 6, l7 => 7)

    total = 0

    output = split(data, "|")[2]
    out_groups = Set.(split(output))
    for group in out_groups
        g = Set([mapping[c] for c in group])
        total = 10 * total + segments_to_digit[g]
    end
    total
end

grab!(vec, l) = grab_if!(vec, x -> length(x) == l)
grab_if!(vec, test) = popat!(vec, findfirst(test, vec))

@assert 61229 == @show sum(solve_7_digit_puzzle.(day08_test))
@assert 1048410 == @show sum(solve_7_digit_puzzle.(day08_input))

## Day 09: Smoke Basin ##
#########################

day09_test = ["2199943210",
    "3987894921",
    "9856789892",
    "8767896789",
    "9899965678"]
day09_input = readlines(open("input-09"))

function height_map(input::Vector{String})::Array{Int,2}
    h = length(input)
    w = length(input[1])
    h_map = zeros(Int, (w, h))

    for y = 1:h, x = 1:w
        h_map[x, y] = parse(Int, input[y][x])
    end
    h_map
end

function get_basin_mins(h_map::Array{Int,2})::Tuple{Vector{Tuple{Int,Int}},Int}
    w, h = size(h_map)
    mins = []
    score = 0
    for x = 1:w, y = 1:h
        v = h_map[x, y]
        if x >= 2 && h_map[x-1, y] <= v
            continue
        elseif x < w && h_map[x+1, y] <= v
            continue
        elseif y >= 2 && h_map[x, y-1] <= v
            continue
        elseif y < h && h_map[x, y+1] <= v
            continue
        end
        push!(mins, (x, y))
        score += 1 + v
    end
    (mins, score)
end

@assert 15 == @show get_basin_mins(height_map(day09_test))[2]
@assert 439 == @show get_basin_mins(height_map(day09_input))[2]

# Now for part 2 I need to find connected regions that are divided by 9s.
# Sounds like I need some kind of flood-fill algorithm.

function basins(h_map::Array{Int,2})::Int
    mins = get_basin_mins(h_map)[1]
    sizes = []
    for (x, y) in mins
        b_size = flood_fill(h_map, x, y)
        push!(sizes, b_size)
    end
    prod(sort(sizes)[end-2:end])
end

function flood_fill(h_map::Array{Int,2}, x::Int, y::Int)::Int
    is_inside = zeros(Bool, size(h_map))
    flood_fill!(is_inside, h_map, x, y)
    sum(is_inside)
end

function flood_fill!(is_inside::Array{Bool,2}, h_map::Array{Int,2}, x::Int, y::Int)
    is_inside[x, y] = true
    w, h = size(h_map)
    if x >= 2 && h_map[x-1, y] != 9 && !is_inside[x-1, y]
        flood_fill!(is_inside, h_map, x - 1, y)
    end
    if x < w && h_map[x+1, y] != 9 && !is_inside[x+1, y]
        flood_fill!(is_inside, h_map, x + 1, y)
    end
    if y >= 2 && h_map[x, y-1] != 9 && !is_inside[x, y-1]
        flood_fill!(is_inside, h_map, x, y - 1)
    end
    if y < h && h_map[x, y+1] != 9 && !is_inside[x, y+1]
        flood_fill!(is_inside, h_map, x, y + 1)
    end
end

@assert 1134 == @show basins(height_map(day09_test))
@assert 900900 == @show basins(height_map(day09_input))

## Day 10: Syntax Scoring ##
############################

day10_test = readlines(open("test-10"))
day10_input = readlines(open("input-10"))

closing_brackets = Dict('(' => ')', '[' => ']', '{' => '}', '<' => '>')
bracket_score = Dict('(' => 1, '[' => 2, '{' => 3, '<' => 4, ')' => 3, ']' => 57, '}' => 1197, '>' => 25137, nothing => 0)

function corrupted(line)
    stack = []
    for char in line
        if char in "([<{"
            push!(stack, char)
            continue
        end
        if length(stack) == 0
            return char
        end
        opening = pop!(stack)
        if char != closing_brackets[opening]
            return char
        end
    end
    return nothing
end

function corruption_score(line)
    bracket_score[corrupted(line)]
end

@assert 26397 == @show sum(corruption_score.(day10_test))
@assert 315693 == @show sum(corruption_score.(day10_input))

# Part two

function closing_stack(line)
    stack = []
    for char in line
        if char in "([<{"
            push!(stack, char)
            continue
        end
        if length(stack) == 0
            throw("corrupted")
        end
        opening = pop!(stack)
        if char != closing_brackets[opening]
            throw("corrupted")
        end
    end
    stack
end

function closing_score(line)
    if corruption_score(line) > 0
        0
    else
        stack = closing_stack(line)
        v = 0
        for c in reverse(stack)
            v = 5 * v + bracket_score[c]
        end
        v
    end
end

function middle_score(lines::Vector{String})::Int
    scores = filter(x -> x > 0, closing_score.(lines))
    sort(scores)[div(1 + end, 2)]
end

@assert 288957 == @show middle_score(day10_test)
@assert 1870887234 == @show middle_score(day10_input)

## Day 11: Dumbo Octopus ##
###########################

day11_test = readlines(open("test-11"))
day11_input = readlines(open("input-11"))

function octopus_input(lines::Vector{String})::Array{Int,2}
    result = zeros(Int, 10, 10)
    for x = 1:10, y = 1:10
        result[x, y] = parse(Int, lines[y][x])
    end
    result
end

function octopus_step!(octo_grid)
    octo_grid .+= 1
    # Recursive flashing
    flashed = zeros(Bool, 10, 10)
    for x = 1:10, y = 1:10
        if flashed[x, y]
            continue
        end
        if octo_grid[x, y] > 9
            octo_flash_recurse!(octo_grid, flashed, x, y)
        end
    end
    # Reset every high octopus
    for x = 1:10, y = 1:10
        if flashed[x, y]
            octo_grid[x, y] = 0
        end
    end
    # Return the number of octopi that have flashed.
    sum(flashed)
end

function octo_flash_recurse!(octo_grid, flashed, x, y)
    flashed[x, y] = true
    triggered = octo_flash!(octo_grid, x, y)
    for octo in triggered
        if !flashed[octo[1], octo[2]]
            octo_flash_recurse!(octo_grid, flashed, octo[1], octo[2])
        end
    end
end

function octo_flash!(octo_grid, x, y)::Vector{Tuple{Int,Int}}
    triggered = []
    for dx = -1:1, dy = -1:1
        if dx == 0 && dy == 0
            continue
        elseif 1 <= x + dx <= 10 && 1 <= y + dy <= 10
            octo_grid[x+dx, y+dy] += 1
            if octo_grid[x+dx, y+dy] > 9
                push!(triggered, (x + dx, y + dy))
            end
        end
    end
    triggered
end

# a = octopus_input(day11_input)
# c = 0
# for i = 1:100
#     global c
#     c = c + @show octopus_step!(a)
# end
# @show c

# Solution = 1681

# Part 2

a = octopus_input(day11_input)
for i = 1:1000000
    c = octopus_step!(a)
    if c == 100
        @show i
        break
    end
end

# Solution = 276

## Day 12: Passage Pathing ##
#############################

day12_test_1 = readlines(open("test-12-1"))
day12_test_2 = readlines(open("test-12-2"))
day12_test_3 = readlines(open("test-12-3"))
day12_input = readlines(open("input-12"))

function cave_graph(lines::Vector{String})::Dict
    graph = Dict()
    for line in lines
        m = match(r"([a-zA-Z]+)-([a-zA-Z]+)", line)
        add_to_graph!(graph, m[1], m[2])
        add_to_graph!(graph, m[2], m[1])
    end
    graph
end

function add_to_graph!(graph, from, to)
    if !(from in keys(graph))
        graph[from] = []
    end
    push!(graph[from], to)
end

is_big_cave(cave_name) = occursin(r"^[A-Z]+$", cave_name)

function count_all_path(graph, from, with_time = false, seen = Set([from]))
    if from == "end"
        return 1
    end

    total = 0
    next_nodes = graph[from]
    for node in next_nodes
        if !(node in seen)
            if is_big_cave(node)
                total += count_all_path(graph, node, with_time, seen)
            else
                my_seen = copy(seen)
                push!(my_seen, node)
                total += count_all_path(graph, node, with_time, my_seen)
            end
        elseif with_time && node != "start" # Node is already seen but we still have time
            @assert !is_big_cave(node)
            total += count_all_path(graph, node, false, seen)
        end
    end
    total
end

@assert 10 == @show count_all_path(cave_graph(day12_test_1), "start")
@assert 19 == @show count_all_path(cave_graph(day12_test_2), "start")
@assert 226 == @show count_all_path(cave_graph(day12_test_3), "start")
@assert 3802 == @show count_all_path(cave_graph(day12_input), "start")

@assert 36 == @show count_all_path(cave_graph(day12_test_1), "start", true)
@assert 103 == @show count_all_path(cave_graph(day12_test_2), "start", true)
@assert 3509 == @show count_all_path(cave_graph(day12_test_3), "start", true)
@assert 99448 == @show count_all_path(cave_graph(day12_input), "start", true)

## Day 13: Transparent Origami ##
#################################

day13_test = readlines(open("test-13"))
day13_input = readlines(open("input-13"))

function transparent_paper_code(data)
    dots = Set()
    instructions = []
    for i = 1:length(data)
        m = match(r"^([0-9]+),([0-9]+)$", data[i])
        if m !== nothing
            push!(dots, (parse(Int, m[1]), parse(Int, m[2])))
            continue
        end
        m = match(r"^fold along ([xy])=([0-9]+)$", data[i])
        if m !== nothing
            push!(instructions, (m[1], parse(Int, m[2])))
            continue
        end
    end

    for instr in instructions
        new_dots = Set()
        if instr[1] == "x"
            for dot in dots
                if dot[1] > instr[2]
                    push!(new_dots, (2 * instr[2] - dot[1], dot[2]))
                else
                    push!(new_dots, dot)
                end
            end
        elseif instr[1] == "y"
            for dot in dots
                if dot[2] > instr[2]
                    push!(new_dots, (dot[1], 2 * instr[2] - dot[2]))
                else
                    push!(new_dots, dot)
                end
            end
        end

        dots = new_dots
    end

    # get bounding box of tuples

    matrix = zeros(Int, (maximum(t -> t[1], dots) + 1, maximum(t -> t[2], dots) + 1))
    for dot in dots
        matrix[dot[1]+1, dot[2]+1] = 1
    end
    matrix
end

# Solutions found interactively.
# transparent_paper_code(day13_test)
# @show transpose(transparent_paper_code(day13_input))

## Day 14: Extended Polymerization ##
#####################################

day14_test = readlines(open("test-14"))
day14_input = readlines(open("input-14"))

function polymer_rules(data)
    # Skip the first two lines, then read the rest into a dictionary
    rules = Dict()
    for line in data[3:end]
        rule = split(line, " -> ")
        rules[rule[1]] = rule[2][1] # String -> Char
    end
    rules
end

# Part 2 is getting harder, need some memoization I think.

function memoize_polymer(rules, max_depth)
    memory = Dict()
    # Initial knowledge: Each rule evaluates to the left letter if there there
    # are no other polymerization steps left.
    for pair in keys(rules)
        memory[(pair, 0)] = Dict(pair[1] => 1)
    end
    for depth = 1:max_depth
        for pair in keys(rules)
            # @show pair, rules
            left_pair = String([pair[1], rules[pair]])
            right_pair = String([rules[pair], pair[2]])
            values = merge(+, memory[(left_pair, depth - 1)], memory[(right_pair, depth - 1)])
            memory[(pair, depth)] = values
        end
    end
    memory
end

function polymer_score(data, depth)
    memory = memoize_polymer(polymer_rules(data), depth)
    starting_poylmer = data[1]

    total = Dict()

    left = starting_poylmer[1]
    for i = 2:length(starting_poylmer)
        right = starting_poylmer[i]
        key = (String([left, right]), depth)
        total = merge(+, total, memory[key])
        left = right
    end

    total = merge(+, total, Dict(left => 1))

    maximum(values(total)) - minimum(values(total))
end

@assert 1588 == @show polymer_score(day14_test, 10)
@assert 2740 == @show polymer_score(day14_input, 10)

@assert 2188189693529 == @show polymer_score(day14_test, 40)
@assert 2959788056211 == @show polymer_score(day14_input, 40)

## Day 15: Chiton ##
####################

# Looks like I "just" need to apply a pathfinding algorithm.

day15_test = readlines(open("test-15"))
day15_input = readlines(open("input-15"))

function square_digit_inupt(lines::Vector{String})::Array{Int,2}
    w = length(lines[1])
    h = length(lines)
    result = zeros(Int, w, h)
    for x = 1:w, y = 1:h
        result[x, y] = parse(Int, lines[y][x])
    end
    result
end

function shortest_path(risk_map)
    infinity = length(risk_map) * 9

    initial = (1, 1)
    unvisited = Set([(x, y) for x = 1:size(risk_map, 1) for y = 1:size(risk_map, 2)])
    tentative_distances = fill(infinity, size(risk_map))
    tentative_distances[initial...] = 0

    todo_map = Dict(0 => [initial])

    current = initial
    while true
        # Find a node with the lowest tentative distance that is unvisited
        while !(current in unvisited)
            # Find first entry in todo_map
            lowest_key = minimum(keys(todo_map))
            places_with_lowest_key = todo_map[lowest_key]
            current = pop!(places_with_lowest_key)
            if length(places_with_lowest_key) == 0
                delete!(todo_map, lowest_key)
            end
        end

        if current == size(risk_map)
            return tentative_distances[current...]
        end

        # Mark neiboors as visited
        for nbh in nbhd4(current, size(risk_map))
            distance_through_current = tentative_distances[current...] + risk_map[nbh...]
            if distance_through_current < tentative_distances[nbh...]
                tentative_distances[nbh...] = distance_through_current
                places_with_same_key = get!(todo_map, distance_through_current, [])
                push!(places_with_same_key, nbh)
                todo_map[distance_through_current] = places_with_same_key
            end
        end
        delete!(unvisited, current)
    end
end

function nbhd4(pos, dims)
    nbhs = []
    if pos[1] > 1
        push!(nbhs, (pos[1] - 1, pos[2]))
    end
    if pos[2] > 1
        push!(nbhs, (pos[1], pos[2] - 1))
    end
    if pos[1] < dims[1]
        push!(nbhs, (pos[1] + 1, pos[2]))
    end
    if pos[2] < dims[2]
        push!(nbhs, (pos[1], pos[2] + 1))
    end
    nbhs
end

@assert 40 == @show shortest_path(square_digit_inupt(day15_test))
@assert 707 == @show shortest_path(square_digit_inupt(day15_input))

# Increase the map size by 5 times in each direction
function bigger_risk_map(risk_map)
    result = zeros(Int, size(risk_map) .* 5)
    w = size(risk_map, 1)
    h = size(risk_map, 2)
    for x1 = 1:w, y1 = 1:h
        for x2 = 0:4, y2 = 0:4
            v = risk_map[x1, y1] + x2 + y2
            result[x1+w*x2, y1+h*y2] = v <= 9 ? v : v - 9
        end
    end
    result
end

@assert 315 == @show shortest_path(bigger_risk_map(square_digit_inupt(day15_test)))
@assert 2942 == @show shortest_path(bigger_risk_map(square_digit_inupt(day15_input)))

## Day 17: Trick Shot ##
########################

day17_test = Dict(:x1 => 20, :x2 => 30, :y1 => -10, :y2 => -5)
day17_input = Dict(:x1 => 265, :x2 => 287, :y1 => -103, :y2 => -58)

# How to stop x inside the rectange?
# x-dist(x0) = x0 * (x0+1) / 2

triangle(x) = x * (x + 1) / 2

# Test x-rest-inside for x = {6, 7}
# Input x-rest-inside for x = 23
# Firing with y = {8, 9} should make it pass through the rectangle
# ((276, 276) .- (103, 58)) ./ 23 = (7.521739130434782, 9.478260869565217)

# So I try (23, 9) how high do I get? => Would be 45, which is not the correct answer.

# Since x is now resting, I can fire y arbitrarily strong to try and end back inside.
# Should be -103 -> (23, 102).

@assert 5253 == @show triangle(102) # Day 16 part 1.

# Part 2, now with actual programming :P

# Test: x min must be at least x-rest-inside. So x >= 6 for test, x >= 23 for input.
# IMHO, the "rest inside" cases are special, because y has arbitrary time to hit.

function ends_inside(dx, dy, rect)
    x, y = (0, 0)
    while true
        # Move
        x += dx
        y += dy
        # Check
        if x >= rect[:x1] && x <= rect[:x2] && y >= rect[:y1] && y <= rect[:y2]
            return true
        elseif x > rect[:x2] || y < rect[:y1] # No chance of hitting the rectangle
            return false
        end
        # Accelerate
        if dx > 0
            dx -= 1
        elseif dx < 0
            dx += 1
        end
        dy -= 1
    end
end

# Test the ends_inside function.
include("assertEndsInside.jl")

@assert 1770 == @show sum([ends_inside(x, y, day17_input) for y = -103:102, x = 1:287])


text:
module IDSGraphs

using LazyArtifacts

# exports
export

# structure
Component, LeftRightStructure, TopBottomStructure, VerticalThirdsStructure,
HorizontalThirdsStructure, CenterNestedStructure, BottomNestedStructure, TopNestedStructure,
RightNestedStructure, BottomRightNestedStructure, BottomLeftNestedStructure,
TopRightNestedStructure, OverlapStructure,

# graph & new graph-related functions (LightGraphs not reexported)
IDSGraph, components, compounds

include("structures.jl")
include("dependencygraphs.jl")

end


text:
using ParallelAccelerator

@acc function loopIt()
  A = Array{Int8}(10000,
                  10000)
  B = Array{Int8}(10000, 10000)
  C = Array{Int8}(10000, 10000)
  D = Array{Int8}(10000, 10000)
  for i = 1:10000
    for j = 1:10000
      @fuse 2 A[i, j] = Int8(0)
    end
  end
  for i = 1:10000
    for j = 1:10000
      @fuse 1 B[i, j] = A[i, j] + Int8(1)
    end
  end
  for i = 1:10000
    for j = 1:10000
      @fuse 1 C[i, j] = B[i, j] + Int8(2)
    end
  end
  for i = 1:10000
    for j = 1:10000
      D[i, j] = C[i, j] + Int8(3)
    end
  end
  return D
end

# @acc function loopIt()
#   A = Array{Int8}(10000, 10000)
#   B = Array{Int8}(10000, 10000)
#   C = Array{Int8}(10000, 10000)
#   D = Array{Int8}(10000, 10000)
#   for i = 1:10000
#     for j = 1:10000
#       A[i, j] = Int8(0)
#     end
#   end
#   for i = 1:10000
#     for j = 1:10000
#       B[i, j] = A[i, j] + Int8(1)
#     end
#   end
#   for i = 1:10000
#     for j = 1:10000
#       C[i, j] = B[i, j] + Int8(2)
#     end
#   end
#   for i = 1:10000
#     for j = 1:10000
#       D[i, j] = C[i, j] + Int8(3)
#     end
#   end
#   return D
# end

function loopIt2()
  D = Array{Int8}(10000, 10000)
  tic()
  A = Array{Int8}(10000, 10000)
  B = Array{Int8}(10000, 10000)
  C = Array{Int8}(10000, 10000)
  for i = 1:10000
    for j = 1:10000
      A[i, j] = Int8(0)
    end
  end
  for i = 1:10000
    for j = 1:10000
      B[i, j] = A[i, j] + Int8(1)
    end
  end
  for i = 1:10000
    for j = 1:10000
      C[i, j] = B[i, j] + Int8(2)
    end
  end
  for i = 1:10000
    for j = 1:10000
      D[i, j] = C[i, j] + Int8(3)
    end
  end
  println("Julia execution time ", toq())
  return D
end

function main()
    tic()
    R_tiramisu = loopIt()
    println("Tiramisu ", toq())

    tic()
    R_regular = loopIt2()
    println("Julia ", toq())

    for i = 1:10000
      for j = 1:10000
        if R_tiramisu[i, j] != 6
          println((i, j), " ", R_tiramisu[i, j])
          break
        end
      end
    end

    println(R_tiramisu == R_regular)
    # println(D)

end

main()

max_stars_repo_name: malekbr/julia-to-tiramisu

max_issues_repo_name: malekbr/julia-to-tiramisu

max_forks_repo_name: malekbr/julia-to-tiramisu