text:
import Base.+, Base.-, Base.*, Base.<, Base.>

export Square, SquareFile, SquareRank, SquareDelta

export distance,
    file, filefromchar, isok, rank, rankfromchar, squarefromstring, tochar, tostring

export FILE_NONE, FILE_A, FILE_B, FILE_C, FILE_D, FILE_E, FILE_F, FILE_G, FILE_H
export RANK_NONE, RANK_1, RANK_2, RANK_3, RANK_4, RANK_5, RANK_6, RANK_7, RANK_8
export FILE_MIN, FILE_MAX, RANK_MIN, RANK_MAX

export SQ_A1, SQ_A2, SQ_A3, SQ_A4, SQ_A5, SQ_A6, SQ_A7, SQ_A8
export SQ_B1, SQ_B2, SQ_B3, SQ_B4, SQ_B5, SQ_B6, SQ_B7, SQ_B8
export SQ_C1, SQ_C2, SQ_C3, SQ_C4, SQ_C5, SQ_C6, SQ_C7, SQ_C8
export SQ_D1, SQ_D2, SQ_D3, SQ_D4, SQ_D5, SQ_D6, SQ_D7, SQ_D8
export SQ_E1, SQ_E2, SQ_E3, SQ_E4, SQ_E5, SQ_E6, SQ_E7, SQ_E8
export SQ_F1, SQ_F2, SQ_F3, SQ_F4, SQ_F5, SQ_F6, SQ_F7, SQ_F8
export SQ_G1, SQ_G2, SQ_G3, SQ_G4, SQ_G5, SQ_G6, SQ_G7, SQ_G8
export SQ_H1, SQ_H2, SQ_H3, SQ_H4, SQ_H5, SQ_H6, SQ_H7, SQ_H8
export SQ_NONE

export DELTA_N, DELTA_S, DELTA_E, DELTA_W
export DELTA_NW, DELTA_NE, DELTA_SW, DELTA_SE

"""
    SquareFile

Type representing the file of a square on a chess board.

Usually, a `SquareFile` is obtained either by calling the function `file()` on a
`Square` or through one of the constants `FILE_A`, `FILE_B`, ..., `FILE_H`.
"""
struct SquareFile
    val::Int
end

"""
    SquareRank

Type representing the rank of a square on a chess board.

Usually, a `SquareRank` is obtained either by calling the function `rank()` on a
`Square` or through one of the constants `RANK_1`, `RANK_2`, ..., `RANK_8`.
"""
struct SquareRank
    val::Int
end

function Base.show(io::IO, f::SquareFile)
    if isok(f)
        print(io, "FILE_$(uppercase(tochar(f)))")
    elseif f == FILE_NONE
        print(io, "FILE_NONE")
    else
        print(io, "SquareFile($(f.val))")
    end
end

function Base.show(io::IO, r::SquareRank)
    if isok(r)
        print(io, "RANK_$(tochar(r))")
    elseif r == RANK_NONE
        print(io, "RANK_NONE")
    else
        print(io, "SquareRank($(r.val))")
    end
end

(<)(f1::SquareFile, f2::SquareFile) = f1.val < f2.val
(<)(r1::SquareRank, r2::SquareRank) = r1.val < r2.val

const FILE_NONE = SquareFile(0)
const FILE_A = SquareFile(1)
const FILE_B = SquareFile(2)
const FILE_C = SquareFile(3)
const FILE_D = SquareFile(4)
const FILE_E = SquareFile(5)
const FILE_F = SquareFile(6)
const FILE_G = SquareFile(7)
const FILE_H = SquareFile(8)
const FILE_MIN = FILE_A
const FILE_MAX = FILE_H

const RANK_NONE = SquareRank(0)
const RANK_1 = SquareRank(8)
const RANK_2 = SquareRank(7)
const RANK_3 = SquareRank(6)
const RANK_4 = SquareRank(5)
const RANK_5 = SquareRank(4)
const RANK_6 = SquareRank(3)
const RANK_7 = SquareRank(2)
const RANK_8 = SquareRank(1)
const RANK_MIN = RANK_8
const RANK_MAX = RANK_1

"""
    isok(f::SquareFile)

Tests whether a `SquareFile` contains a valid value.

# Examples

```julia-repl
julia> isok(FILE_A)
true

julia> isok(SquareFile(0))
false
```
"""
function isok(f::SquareFile)
    f ≥ FILE_MIN && f ≤ FILE_MAX
end

"""
    isok(r::SquareRank)

Tests whether a `SquareRank` contains a valid value.

# Examples

```julia-repl
julia> isok(RANK_6)
true

julia> isok(SquareRank(42))
false
```
"""
function isok(r::SquareRank)
    r ≥ RANK_MIN && r ≤ RANK_MAX
end

"""
    filefromchar(c::Char)

Tries to convert a character to a file.

The return value is a `Union{SquareFile, Nothing}`. The `nothing` is returned
in case the character does not represent a valid file.

# Examples

```julia-repl
julia> filefromchar('c')
FILE_C

julia> filefromchar('2') == nothing
true
```
"""
function filefromchar(c::Char)::Union{SquareFile,Nothing}
    result = SquareFile(Int(c) - Int('a') + 1)
    isok(result) ? result : nothing
end

"""
    rankfromchar(c::Char)

Tries to convert a character to a rank.

The return value is a `Union{SquareRank, Nothing}`. The `nothing` is returned
in case the character does not represent a valid rank.

# Examples

```julia-repl
julia> rankfromchar('2')
RANK_2

julia> rankfromchar('x') == nothing
true
```
"""
function rankfromchar(c::Char)::Union{SquareRank,Nothing}
    result = SquareRank(8 - Int(c) + Int('1'))
    isok(result) ? result : nothing
end

"""
    tochar(f::SquareFile)

Converts a `SquareFile` to a character.

# Examples

```julia-repl
julia> tochar(FILE_E)
'e': ASCII/Unicode U+0065 (category Ll: Letter, lowercase)
```
"""
function tochar(f::SquareFile)
    Char(f.val - 1 + Int('a'))
end

"""
    tochar(r::SquareRank)

Converts `SquareRank` to a character.

# Examples

```julia-repl
julia> tochar(RANK_3)
'3': ASCII/Unicode U+0033 (category Nd: Number, decimal digit)
```
"""
function tochar(r::SquareRank)
    Char(8 - r.val + Int('1'))
end

"""
    Square

Type representing a square on a chess board.

A `Square` can be constructed either with an `Int` (with the convention a8=1,
a7=2, ..., a1=8, b8=9, b7=10, ..., h1=64) or with a `SquareFile` and a
`SquareRank`. There are also constants `SQ_A1`, ..., `SQ_H8` for all 64
squares on the board.

# Examples

```julia-repl
julia> Square(FILE_G, RANK_6)
SQ_G6

julia> Square(8)
SQ_A1
```
"""
struct Square
    val::Int
end

function Base.show(io::IO, s::Square)
    if isok(s)
        print(io, "SQ_$(uppercase(tostring(s)))")
    elseif s == SQ_NONE
        print(io, "SQ_NONE")
    else
        print(io, "Square($(s.val))")
    end
end

"""
    Square(f::SquareFile, r::SquareRank)

Construct a square with the given file and rank.

# Examples

```julia-repl
julia> Square(FILE_D, RANK_5)
SQ_D5
```
"""
Square(f::SquareFile, r::SquareRank) = Square(r.val + 8 * (f.val - 1))

"""
    file(s::Square)

Compute the file of the square `s`.

# Examples

```julia-repl
julia> file(SQ_C4)
FILE_C
```
"""
function file(s::Square)
    SquareFile(fld1(s.val, 8))
end

"""
    rank(s::Square)

Compute the rank of the square `s`.

# Examples

```julia-repl
julia> rank(SQ_C4)
RANK_4
```
"""
function rank(s::Square)
    SquareRank(mod1(s.val, 8))
end

"""
    isok(s::Square)

Tests whether a `Square` has a valid value.

# Examples

```julia-repl
julia> isok(SQ_G7)
true

julia> isok(Square(42))
true

julia> isok(Square(100))
false

julia> isok(Square(0))
false
```
"""
function isok(s::Square)
    isok(file(s)) && isok(rank(s))
end

"""
    tostring(s::Square)

Converts a square to a string in standard algebraic notation. If the square
has an invalid value, the returned string is `"??"`.

# Examples

```julia-repl
julia> tostring(SQ_E4)
"e4"

julia> tostring(Square(100))
"??"
```
"""
function tostring(s::Square)
    if isok(s)
        result = IOBuffer()
        write(result, tochar(file(s)), tochar(rank(s)))
        String(take!(result))
    else
        "??"
    end
end

"""
    squarefromstring(s::String)

Tries to convert a string to a `Square`.

The return value is of type `Union{Square, Nothing}`.

If the input string is too short, or if the two first characters do not
represent a square in standard algebraic notation, returns `nothing`. If the
first two characters do represent avalid square, that square is returned, even
if there are additional characters.

# Examples

```julia-repl
julia> squarefromstring("d6")
SQ_D6

julia> squarefromstring("xy") == nothing
true

julia> squarefromstring("") == nothing
true

julia> squarefromstring("g1f3")
SQ_G1
```
"""
function squarefromstring(s::String)::Union{Square,Nothing}
    if length(s) < 2
        nothing
    else
        f = filefromchar(s[1])
        r = rankfromchar(s[2])
        if isa(f, SquareFile) && isa(r, SquareRank)
            Square(f, r)
        else
            nothing
        end
    end
end

"""
    distance(f1::SquareFile, f2::SquareFile)

The horizontal distance between two files.
"""
function distance(f1::SquareFile, f2::SquareFile)::Int
    abs(f2.val - f1.val)
end

"""
    distance(r1::SquareRank, r2::SquareRank)

The vertical distance between two ranks.
"""
function distance(r1::SquareRank, r2::SquareRank)::Int
    abs(r2.val - r1.val)
end

"""
    distance(s1::Square, s2::Square)

The distance between two squares, counted by number of king moves.
"""
function distance(s1::Square, s2::Square)::Int
    max(distance(file(s1), file(s2)), distance(rank(s1), rank(s2)))
end

"""
    SquareDelta

A type representing the delta or vector between two squares.

A `SquareDelta` value is usually obtained either through one of the constants
`DELTA_N`, `DELTA_S`, `DELTA_E`, `DELTA_W`, `DELTA_NW`, `DELTA_NE`, `DELTA_SW`,
`DELTA_SE`, or by subtracting two square values.

It is possible to add or subtract two `SquareDelta`s, to multiply a
`SquareDelta` by an integer scalar, or to add or subtract a `SquareDelta` to a
`Square`.

# Examples

```julia-repl
julia> DELTA_N + DELTA_W == DELTA_NW
true

julia> SQ_D3 - SQ_C3 == DELTA_E
true

julia> SQ_G8 - 3 * DELTA_N
SQ_G5
```
"""
struct SquareDelta
    val::Int
end

(-)(s1::Square, s2::Square) = SquareDelta(s1.val - s2.val)
(+)(s::Square, d::SquareDelta) = Square(s.val + d.val)
(-)(s::Square, d::SquareDelta) = Square(s.val - d.val)
(+)(d1::SquareDelta, d2::SquareDelta) = SquareDelta(d1.val + d2.val)
(-)(d1::SquareDelta, d2::SquareDelta) = SquareDelta(d1.val - d2.val)
(*)(i::Integer, d::SquareDelta) = SquareDelta(i * d.val)

const SQ_A1 = Square(FILE_A, RANK_1)
const SQ_A2 = Square(FILE_A, RANK_2)
const SQ_A3 = Square(FILE_A, RANK_3)
const SQ_A4 = Square(FILE_A, RANK_4)
const SQ_A5 = Square(FILE_A, RANK_5)
const SQ_A6 = Square(FILE_A, RANK_6)
const SQ_A7 = Square(FILE_A, RANK_7)
const SQ_A8 = Square(FILE_A, RANK_8)

const SQ_B1 = Square(FILE_B, RANK_1)
const SQ_B2 = Square(FILE_B, RANK_2)
const SQ_B3 = Square(FILE_B, RANK_3)
const SQ_B4 = Square(FILE_B, RANK_4)
const SQ_B5 = Square(FILE_B, RANK_5)
const SQ_B6 = Square(FILE_B, RANK_6)
const SQ_B7 = Square(FILE_B, RANK_7)
const SQ_B8 = Square(FILE_B, RANK_8)

const SQ_C1 = Square(FILE_C, RANK_1)
const SQ_C2 = Square(FILE_C, RANK_2)
const SQ_C3 = Square(FILE_C, RANK_3)
const SQ_C4 = Square(FILE_C, RANK_4)
const SQ_C5 = Square(FILE_C, RANK_5)
const SQ_C6 = Square(FILE_C, RANK_6)
const SQ_C7 = Square(FILE_C, RANK_7)
const SQ_C8 = Square(FILE_C, RANK_8)

const SQ_D1 = Square(FILE_D, RANK_1)
const SQ_D2 = Square(FILE_D, RANK_2)
const SQ_D3 = Square(FILE_D, RANK_3)
const SQ_D4 = Square(FILE_D, RANK_4)
const SQ_D5 = Square(FILE_D, RANK_5)
const SQ_D6 = Square(FILE_D, RANK_6)
const SQ_D7 = Square(FILE_D, RANK_7)
const SQ_D8 = Square(FILE_D, RANK_8)

const SQ_E1 = Square(FILE_E, RANK_1)
const SQ_E2 = Square(FILE_E, RANK_2)
const SQ_E3 = Square(FILE_E, RANK_3)
const SQ_E4 = Square(FILE_E, RANK_4)
const SQ_E5 = Square(FILE_E, RANK_5)
const SQ_E6 = Square(FILE_E, RANK_6)
const SQ_E7 = Square(FILE_E, RANK_7)
const SQ_E8 = Square(FILE_E, RANK_8)

const SQ_F1 = Square(FILE_F, RANK_1)
const SQ_F2 = Square(FILE_F, RANK_2)
const SQ_F3 = Square(FILE_F, RANK_3)
const SQ_F4 = Square(FILE_F, RANK_4)
const SQ_F5 = Square(FILE_F, RANK_5)
const SQ_F6 = Square(FILE_F, RANK_6)
const SQ_F7 = Square(FILE_F, RANK_7)
const SQ_F8 = Square(FILE_F, RANK_8)

const SQ_G1 = Square(FILE_G, RANK_1)
const SQ_G2 = Square(FILE_G, RANK_2)
const SQ_G3 = Square(FILE_G, RANK_3)
const SQ_G4 = Square(FILE_G, RANK_4)
const SQ_G5 = Square(FILE_G, RANK_5)
const SQ_G6 = Square(FILE_G, RANK_6)
const SQ_G7 = Square(FILE_G, RANK_7)
const SQ_G8 = Square(FILE_G, RANK_8)

const SQ_H1 = Square(FILE_H, RANK_1)
const SQ_H2 = Square(FILE_H, RANK_2)
const SQ_H3 = Square(FILE_H, RANK_3)
const SQ_H4 = Square(FILE_H, RANK_4)
const SQ_H5 = Square(FILE_H, RANK_5)
const SQ_H6 = Square(FILE_H, RANK_6)
const SQ_H7 = Square(FILE_H, RANK_7)
const SQ_H8 = Square(FILE_H, RANK_8)

const SQ_NONE = Square(65)

const DELTA_N = SQ_A2 - SQ_A1
const DELTA_S = -1 * DELTA_N
const DELTA_E = SQ_B1 - SQ_A1
const DELTA_W = -1 * DELTA_E
const DELTA_NW = DELTA_N + DELTA_W
const DELTA_NE = DELTA_N + DELTA_E
const DELTA_SW = DELTA_S + DELTA_W
const DELTA_SE = DELTA_S + DELTA_E


text:
module MaxHelpingHandHeatWaveNoColorGrade

using ..Ahorn, Maple

@mapdef Effect "MaxHelpingHand/HeatWaveNoColorGrade" HeatWaveNoColorGrade(only::String="*", exclude::String="", controlColorGradeWhenActive::Bool=false)

placements = HeatWaveNoColorGrade

function Ahorn.canFgBg(effect::HeatWaveNoColorGrade)
    return true, true
end

end

ext: jl

max_stars_repo_name: Kitty-Cats/MaxHelpingHand

max_issues_repo_name: Kitty-Cats/MaxHelpingHand

max_forks_repo_name: Kitty-Cats/MaxHelpingHand


text:
"""
    struct LookupResult <: Any
Hold a pointer to a `clang::LookupResult` object.
"""
struct LookupResult
    ptr::CXLookupResult
end

Base.unsafe_convert(::Type{CXLookupResult}, x::LookupResult) = x.ptr
Base.cconvert(::Type{CXLookupResult}, x::LookupResult) = x


text:
using SafeTestsets
@safetestset "Stub Tests:" begin include("stub_tests.jl") end


text:
using Zygote, Test

function tasks1(x)
  ch = Channel(Inf)
  put!(ch, x^2) + take!(ch)
end

@test gradient(tasks1, 5) == (20,)

function tasks2(x)
  ch = Channel(0)
  t = @async put!(ch, x^2)
  y = take!(ch)
  wait(t)
  return y
end

@test gradient(tasks2, 5) == (10,)

function tasks3(x)
  ch = Channel(0)
  @sync begin
    @async put!(ch, x^2)
    take!(ch)
  end
end

@test gradient(tasks3, 5) == (10,)

tasks4(x) = fetch(@async x^2)

@test gradient(tasks4, 5) == (10,)

VERSION > v"1.3-" && include("threads.jl")

@test Zygote.pullback(Array, [1f0])[1] == [1f0]

@testset "#300" begin
  t = (rand(2, 2), rand(2, 2))
  ps = Params(t)
  gs = gradient(()->sum(t[1]), ps)
  @test gs[t[1]] == ones(2, 2)
end

@testset "#594" begin
  struct A x::Float64 end
  f(a,v) = a.x + v
  g(A,V) = sum(f.(A,V))
  X = A.(randn(2))
  Y = randn(2,2)
  ∇ = gradient(g,X,Y)
  @test ∇[1] == [(x = 2.0,); (x = 2.0,)]
  @test ∇[2] == [1 1; 1 1]
end


text:
using Xpress

M = 3000
N = 2000
A = rand(M,N)
b = rand(M)
f = rand(N)
lb = zeros(N)

for i in 1:4000
    if mod(i,10) == 0
        @show i
    end
    model = Xpress.xpress_model(
        name="lp_02",
        sense=:maximize,
        f = f,
        A = A,
        b = b,
        lb = lb,
        finalize_env = true)
end


text:
using Revise

includet("src/Spring.jl")
IArray = InternalArray

# 3. Validating the solution --------------------------------------------------------------

const step = 1e-3
const nit = Int(40/step)  # round to avoid small rounding errors

# We define a problem with no friction
prob = SpringProblem(p0 = 0.8, v0 = 0, m = 1, k = 1.5)

x = eulermethod(prob, step, nit)

timeline = [prob.t0 + i*step for i in 1:length(x)]
plot(timeline, IArray(x, 1), label = "standard")

# 5. Finding the frequency of the unforced movement ---------------------------------------

# Empirically, We can see in the plot that the distance between two
# equal values (!= 0) is aproximately 5.
# Let's get the actual distance from our solution.

p0 = 1       # Reference point
P0 = -1      # Next equal point
for i = 50:length(x) # We avoid the first numbers because they are very close
    if norm(x[p0][1]-x[i][1]) < 1e-3
        P0 = i
        break
    end
end

if P0 != -1
    period = timeline[P0]-timeline[p0] # Is similar to what we observed
    freq = 1/period
end

# 6. Using this frequency to force it -----------------------------------------------------

cond_prob = SpringProblem(p0 = 0.8, v0 = 0, m = 1, k = 1.5, A = 1, ω = 2*π*freq)

y = eulermethod(cond_prob, step, nit)
plot!(timeline, IArray(y, 1), label = "conditioned")

savefig("media/spring_frequency_unforced.png")

max_stars_repo_name: useredsa/numeric-differential-equations

max_issues_repo_name: useredsa/numeric-differential-equations

max_forks_repo_name: useredsa/numeric-differential-equations


text:
#   Generating #define constants
const crypto_auth_hmacsha512256_BYTES       = 32
const crypto_auth_hmacsha512256_KEYBYTES    = 32
const crypto_auth_hmacsha512256_ref         = crypto_auth_hmacsha512256
const crypto_auth_hmacsha512256_ref_verify  = crypto_auth_hmacsha512256_verify
const crypto_auth_BYTES                     = crypto_auth_hmacsha512256_BYTES
const crypto_auth_KEYBYTES                  = crypto_auth_hmacsha512256_KEYBYTES
const crypto_auth_PRIMITIVE                 = "hmacsha512256"
const crypto_auth_hmacsha256_BYTES          = 32
const crypto_auth_hmacsha256_KEYBYTES       = 32
const crypto_auth_hmacsha256_ref            = crypto_auth_hmacsha256
const crypto_auth_hmacsha256_ref_verify     = crypto_auth_hmacsha256_verify
const crypto_box_curve25519xsalsa20poly1305_PUBLICKEYBYTES = 32
const crypto_box_curve25519xsalsa20poly1305_SECRETKEYBYTES = 32
const crypto_box_curve25519xsalsa20poly1305_BEFORENMBYTES = 32
const crypto_box_curve25519xsalsa20poly1305_NONCEBYTES = 24
const crypto_box_curve25519xsalsa20poly1305_ZEROBYTES = 32
const crypto_box_curve25519xsalsa20poly1305_BOXZEROBYTES = 16
const crypto_box_curve25519xsalsa20poly1305_MACBYTES = (crypto_box_curve25519xsalsa20poly1305_ZEROBYTES - crypto_box_curve25519xsalsa20poly1305_BOXZEROBYTES)
const crypto_box_curve25519xsalsa20poly1305_ref = crypto_box_curve25519xsalsa20poly1305
const crypto_box_curve25519xsalsa20poly1305_ref_open = crypto_box_curve25519xsalsa20poly1305_open
const crypto_box_curve25519xsalsa20poly1305_ref_keypair = crypto_box_curve25519xsalsa20poly1305_keypair
const crypto_box_curve25519xsalsa20poly1305_ref_beforenm = crypto_box_curve25519xsalsa20poly1305_beforenm
const crypto_box_curve25519xsalsa20poly1305_ref_afternm = crypto_box_curve25519xsalsa20poly1305_afternm
const crypto_box_curve25519xsalsa20poly1305_ref_open_afternm = crypto_box_curve25519xsalsa20poly1305_open_afternm
const crypto_box_PUBLICKEYBYTES             = crypto_box_curve25519xsalsa20poly1305_PUBLICKEYBYTES
const crypto_box_SECRETKEYBYTES             = crypto_box_curve25519xsalsa20poly1305_SECRETKEYBYTES
const crypto_box_BEFORENMBYTES              = crypto_box_curve25519xsalsa20poly1305_BEFORENMBYTES
const crypto_box_NONCEBYTES                 = crypto_box_curve25519xsalsa20poly1305_NONCEBYTES
const crypto_box_ZEROBYTES                  = crypto_box_curve25519xsalsa20poly1305_ZEROBYTES
const crypto_box_BOXZEROBYTES               = crypto_box_curve25519xsalsa20poly1305_BOXZEROBYTES
const crypto_box_MACBYTES                   = crypto_box_curve25519xsalsa20poly1305_MACBYTES
const crypto_box_PRIMITIVE                  = "curve25519xsalsa20poly1305"
const crypto_core_hsalsa20_OUTPUTBYTES      = 32
const crypto_core_hsalsa20_INPUTBYTES       = 16
const crypto_core_hsalsa20_KEYBYTES         = 32
const crypto_core_hsalsa20_CONSTBYTES       = 16
const crypto_core_hsalsa20_ref2             = crypto_core_hsalsa20
const crypto_core_salsa2012_OUTPUTBYTES     = 64
const crypto_core_salsa2012_INPUTBYTES      = 16
const crypto_core_salsa2012_KEYBYTES        = 32
const crypto_core_salsa2012_CONSTBYTES      = 16
const crypto_core_salsa2012_ref             = crypto_core_salsa2012
const crypto_core_salsa208_OUTPUTBYTES      = 64
const crypto_core_salsa208_INPUTBYTES       = 16
const crypto_core_salsa208_KEYBYTES         = 32
const crypto_core_salsa208_CONSTBYTES       = 16
const crypto_core_salsa208_ref              = crypto_core_salsa208
const crypto_core_salsa20_OUTPUTBYTES       = 64
const crypto_core_salsa20_INPUTBYTES        = 16
const crypto_core_salsa20_KEYBYTES          = 32
const crypto_core_salsa20_CONSTBYTES        = 16
const crypto_core_salsa20_ref               = crypto_core_salsa20
const crypto_generichash_blake2b_BYTES_MIN  = 1
const crypto_generichash_blake2b_BYTES_MAX  = 64
const crypto_generichash_blake2b_KEYBYTES_MIN = 0
const crypto_generichash_blake2b_KEYBYTES_MAX = 64
const crypto_generichash_blake2b_BLOCKBYTES = 128
const crypto_generichash_blake2b_ref        = crypto_generichash_blake2b
const crypto_generichash_blake2b_ref_init   = crypto_generichash_blake2b_init
const crypto_generichash_blake2b_ref_update = crypto_generichash_blake2b_update
const crypto_generichash_blake2b_ref_final  = crypto_generichash_blake2b_final
const crypto_generichash_BYTES_MIN          = crypto_generichash_blake2b_BYTES_MIN
const crypto_generichash_BYTES_MAX          = crypto_generichash_blake2b_BYTES_MAX
const crypto_generichash_KEYBYTES_MIN       = crypto_generichash_blake2b_KEYBYTES_MIN
const crypto_generichash_KEYBYTES_MAX       = crypto_generichash_blake2b_KEYBYTES_MAX
const crypto_generichash_BLOCKBYTES         = crypto_generichash_blake2b_BLOCKBYTES
const crypto_generichash_PRIMITIVE          = "blake2b"
const crypto_hashblocks_sha256_STATEBYTES   = 32
const crypto_hashblocks_sha256_BLOCKBYTES   = 64
const crypto_hashblocks_sha256_ref          = crypto_hashblocks_sha256
const crypto_hashblocks_sha512_STATEBYTES   = 64
const crypto_hashblocks_sha512_BLOCKBYTES   = 128
const crypto_hashblocks_sha512_ref          = crypto_hashblocks_sha512
const crypto_hash_sha512_BYTES              = 64
const crypto_hash_sha512_BLOCKBYTES         = 128
const crypto_hash_sha512_ref                = crypto_hash_sha512
const crypto_hash_BYTES                     = crypto_hash_sha512_BYTES
const crypto_hash_BLOCKBYTES                = crypto_hash_sha512_BLOCKBYTES
const crypto_hash_PRIMITIVE                 = "sha512"
const crypto_hash_sha256_BYTES              = 32
const crypto_hash_sha256_BLOCKBYTES         = 64
const crypto_hash_sha256_ref                = crypto_hash_sha256
const crypto_onetimeauth_poly1305_BYTES     = 16
const crypto_onetimeauth_poly1305_KEYBYTES  = 32
const crypto_onetimeauth_BYTES              = crypto_onetimeauth_poly1305_BYTES
const crypto_onetimeauth_KEYBYTES           = crypto_onetimeauth_poly1305_KEYBYTES
const crypto_onetimeauth_PRIMITIVE          = "poly1305"
const crypto_scalarmult_curve25519_BYTES    = 32
const crypto_scalarmult_curve25519_SCALARBYTES = 32
const crypto_scalarmult_curve25519_ref      = crypto_scalarmult_curve25519
const crypto_scalarmult_curve25519_ref_base = crypto_scalarmult_curve25519_base
const crypto_scalarmult_BYTES               = crypto_scalarmult_curve25519_BYTES
const crypto_scalarmult_SCALARBYTES         = crypto_scalarmult_curve25519_SCALARBYTES
const crypto_scalarmult_PRIMITIVE           = "curve25519"
const crypto_secretbox_xsalsa20poly1305_KEYBYTES = 32
const crypto_secretbox_xsalsa20poly1305_NONCEBYTES = 24
const crypto_secretbox_xsalsa20poly1305_ZEROBYTES = 32
const crypto_secretbox_xsalsa20poly1305_BOXZEROBYTES = 16
const crypto_secretbox_xsalsa20poly1305_ref = crypto_secretbox_xsalsa20poly1305
const crypto_secretbox_xsalsa20poly1305_ref_open = crypto_secretbox_xsalsa20poly1305_open
const crypto_secretbox_KEYBYTES             = crypto_secretbox_xsalsa20poly1305_KEYBYTES
const crypto_secretbox_NONCEBYTES           = crypto_secretbox_xsalsa20poly1305_NONCEBYTES
const crypto_secretbox_ZEROBYTES            = crypto_secretbox_xsalsa20poly1305_ZEROBYTES
const crypto_secretbox_BOXZEROBYTES         = crypto_secretbox_xsalsa20poly1305_BOXZEROBYTES
const crypto_secretbox_PRIMITIVE            = "xsalsa20poly1305"
const crypto_shorthash_siphash24_BYTES      = 8
const crypto_shorthash_siphash24_KEYBYTES   = 16
const crypto_shorthash_siphash24_ref        = crypto_shorthash_siphash24
const crypto_shorthash_BYTES                = crypto_shorthash_siphash24_BYTES
const crypto_shorthash_KEYBYTES             = crypto_shorthash_siphash24_KEYBYTES
const crypto_shorthash_PRIMITIVE            = "siphash24"
const crypto_sign_ed25519_SECRETKEYBYTES    = 64
const crypto_sign_ed25519_PUBLICKEYBYTES    = 32
const crypto_sign_ed25519_BYTES             = 64
const crypto_sign_ed25519_ref10             = crypto_sign_ed25519
const crypto_sign_ed25519_ref10_open        = crypto_sign_ed25519_open
const crypto_sign_ed25519_ref10_keypair     = crypto_sign_ed25519_keypair
const crypto_sign_ed25519_ref10_seed_keypair = crypto_sign_ed25519_seed_keypair
const crypto_sign_edwards25519sha512batch_SECRETKEYBYTES = 64
const crypto_sign_edwards25519sha512batch_PUBLICKEYBYTES = 32
const crypto_sign_edwards25519sha512batch_BYTES = 64
const crypto_sign_edwards25519sha512batch_ref = crypto_sign_edwards25519sha512batch
const crypto_sign_edwards25519sha512batch_ref_open = crypto_sign_edwards25519sha512batch_open
const crypto_sign_edwards25519sha512batch_ref_keypair = crypto_sign_edwards25519sha512batch_keypair
const crypto_sign_BYTES                     = crypto_sign_ed25519_BYTES
const crypto_sign_PUBLICKEYBYTES            = crypto_sign_ed25519_PUBLICKEYBYTES
const crypto_sign_SECRETKEYBYTES            = crypto_sign_ed25519_SECRETKEYBYTES
const crypto_sign_PRIMITIVE                 = "ed25519"
const crypto_stream_aes128ctr_KEYBYTES      = 16
const crypto_stream_aes128ctr_NONCEBYTES    = 16
const crypto_stream_aes128ctr_BEFORENMBYTES = 1408
const crypto_stream_aes128ctr_portable      = crypto_stream_aes128ctr
const crypto_stream_aes128ctr_portable_xor  = crypto_stream_aes128ctr_xor
const crypto_stream_aes128ctr_portable_beforenm = crypto_stream_aes128ctr_beforenm
const crypto_stream_aes128ctr_portable_afternm = crypto_stream_aes128ctr_afternm
const crypto_stream_aes128ctr_portable_xor_afternm = crypto_stream_aes128ctr_xor_afternm
const crypto_stream_xsalsa20_KEYBYTES       = 32
const crypto_stream_xsalsa20_NONCEBYTES     = 24
const crypto_stream_xsalsa20_ref            = crypto_stream_xsalsa20
const crypto_stream_xsalsa20_ref_xor        = crypto_stream_xsalsa20_xor
const crypto_stream_xsalsa20_ref_beforenm   = crypto_stream_xsalsa20_beforenm
const crypto_stream_xsalsa20_ref_afternm    = crypto_stream_xsalsa20_afternm
const crypto_stream_xsalsa20_ref_xor_afternm = crypto_stream_xsalsa20_xor_afternm
const crypto_stream_KEYBYTES                = crypto_stream_xsalsa20_KEYBYTES
const crypto_stream_NONCEBYTES              = crypto_stream_xsalsa20_NONCEBYTES
const crypto_stream_PRIMITIVE               = "xsalsa20"
const crypto_stream_salsa2012_KEYBYTES      = 32
const crypto_stream_salsa2012_NONCEBYTES    = 8
const crypto_stream_salsa2012_ref           = crypto_stream_salsa2012
const crypto_stream_salsa2012_ref_xor       = crypto_stream_salsa2012_xor
const crypto_stream_salsa2012_ref_beforenm  = crypto_stream_salsa2012_beforenm
const crypto_stream_salsa2012_ref_afternm   = crypto_stream_salsa2012_afternm
const crypto_stream_salsa2012_ref_xor_afternm = crypto_stream_salsa2012_xor_afternm
const crypto_stream_salsa208_KEYBYTES       = 32
const crypto_stream_salsa208_NONCEBYTES     = 8
const crypto_stream_salsa208_ref            = crypto_stream_salsa208
const crypto_stream_salsa208_ref_xor        = crypto_stream_salsa208_xor
const crypto_stream_salsa208_ref_beforenm   = crypto_stream_salsa208_beforenm
const crypto_stream_salsa208_ref_afternm    = crypto_stream_salsa208_afternm
const crypto_stream_salsa208_ref_xor_afternm = crypto_stream_salsa208_xor_afternm
const crypto_stream_salsa20_KEYBYTES        = 32
const crypto_stream_salsa20_NONCEBYTES      = 8
const crypto_stream_salsa20_ref             = crypto_stream_salsa20
const crypto_stream_salsa20_ref_xor         = crypto_stream_salsa20_xor
const crypto_stream_salsa20_ref_beforenm    = crypto_stream_salsa20_beforenm
const crypto_stream_salsa20_ref_afternm     = crypto_stream_salsa20_afternm
const crypto_stream_salsa20_ref_xor_afternm = crypto_stream_salsa20_xor_afternm
const crypto_verify_16_BYTES                = 16
const crypto_verify_16_ref                  = crypto_verify_16
const crypto_verify_32_BYTES                = 32
const crypto_verify_32_ref                  = crypto_verify_32


text:
using Pkg; Pkg.activate("..")

using OceanTurb, Printf

for stepper in (:ForwardEuler, :BackwardEuler)
    model = Diffusion.Model(N=100, stepper=stepper)

    iterate!(model, 0.1)

    # Measure memory allocation
    @printf "\nTesting Diffusion with %s...\n" stepper
    for nt = (10, 100, 1000, 10000)
        @printf "nt: % 6d" nt
        @time iterate!(model, 0.1, nt)
    end
end

for stepper in (:ForwardEuler, :BackwardEuler)
    model = KPP.Model(N=100, stepper=stepper)

    iterate!(model, 0.1)

    # Measure memory allocation
    @printf "\nTesting KPP with %s...\n" stepper
    for nt = (10, 100, 1000, 10000)
        @printf "nt: % 6d" nt
        @time iterate!(model, 0.1, nt)
    end
end


text:
#
# User facing struct that unifies a ModelDef and a ModelInstance and delegates
# function calls to one or the other as appropriate.
#
using MacroTools

"""
    modeldef(m)

Return the `ModelDef` contained by Model `m`.
"""
modeldef(m::Model) = m.md

"""
    modelinstance(m::Model)

Return the `ModelInstance` contained by Model `m`.
"""
modelinstance(m::Model) = m.mi

"""
    modelinstance_def(m::Model)

Return the `ModelDef` of the `ModelInstance` contained by Model `m`.
"""
modelinstance_def(m::Model) = modeldef(modelinstance(m))

"""
    is_built(m::Model)

Return true if Model `m` is built, otherwise return false.
"""
is_built(m::Model) = !(dirty(m.md) || modelinstance(m) === nothing)

"""
    is_built(mm::MarginalModel)

Return true if `MarginalModel` `mm` is built, meaning both its `base` and `modified`
`Model`s are built, otherwise return false.
"""
is_built(mm::MarginalModel) = (is_built(mm.base) && is_built(mm.modified))

@delegate compinstance(m::Model, name::Symbol) => mi
@delegate has_comp(m::Model, name::Symbol) => md

@delegate number_type(m::Model) => md

@delegate internal_param_conns(m::Model) => md
@delegate external_param_conns(m::Model) => md

@delegate model_params(m::Model) => md
@delegate model_param(m::Model, comp_name::Symbol, param_name::Symbol; missing_ok = false) => md
@delegate model_param(m::Model, name::Symbol; missing_ok=false) => md

@delegate connected_params(m::Model) => md
@delegate unconnected_params(m::Model) => md
@delegate nothing_params(m::Model) => md

@delegate add_connector_comps!(m::Model) => md

"""
    connect_param!(m::Model, dst_comp_name::Symbol, dst_par_name::Symbol,
                    src_comp_name::Symbol, src_var_name::Symbol,
                    backup::Union{Nothing, Array}=nothing; ignoreunits::Bool=false,
                    backup_offset::Union{Int, Nothing}=nothing)

Bind the parameter `dst_par_name` of one component `dst_comp_name` of model `m`
to a variable `src_var_name` in another component `src_comp_name` of the same model
using `backup` to provide default values and the `ignoreunits` flag to indicate the need
to check match units between the two.  The `backup_offset` argument, which is only valid
when `backup` data has been set, indicates that the backup data should be used for
a specified number of timesteps after the source component begins. ie. the value would be
`1` if the destination componentm parameter should only use the source component
data for the second timestep and beyond.
"""
@delegate connect_param!(m::Model, dst_comp_name::Symbol, dst_par_name::Symbol,
                        src_comp_name::Symbol, src_var_name::Symbol,
                        backup::Union{Nothing, Array}=nothing; ignoreunits::Bool=false,
                        backup_offset::Union{Nothing, Int} = nothing) => md

"""
    connect_param!(m::Model, comp_name::Symbol, param_name::Symbol, model_param_name::Symbol;
                   check_attributes::Bool=true, ignoreunits::Bool=false))

Connect a parameter `param_name` in the component `comp_name` of composite `obj` to
the model parameter `model_param_name`.
"""
@delegate connect_param!(m::Model, comp_name::Symbol, param_name::Symbol, model_param_name::Symbol;
                        check_attributes::Bool=true, ignoreunits::Bool = false) => md

"""
    connect_param!(m::Model, dst::Pair{Symbol, Symbol}, src::Pair{Symbol, Symbol}, backup::Array; ignoreunits::Bool=false)

Bind the parameter `dst[2]` of one component `dst[1]` of model `m`
to a variable `src[2]` in another component `src[1]` of the same model
using `backup` to provide default values and the `ignoreunits` flag to indicate the need
to check match units between the two.  The `backup_offset` argument, which is only valid
when `backup` data has been set, indicates that the backup data should be used for
a specified number of timesteps after the source component begins. ie. the value would be
`1` if the destination componentm parameter should only use the source component
data for the second timestep and beyond.

"""
function connect_param!(m::Model, dst::Pair{Symbol, Symbol}, src::Pair{Symbol, Symbol},
                           backup::Union{Nothing, Array}=nothing;
                           ignoreunits::Bool=false, backup_offset::Union{Int, Nothing} = nothing)
    connect_param!(m.md, dst[1], dst[2], src[1], src[2], backup; ignoreunits=ignoreunits, backup_offset=backup_offset)
end

"""
    disconnect_param!(m::Model, comp_name::Symbol, param_name::Symbol)

Remove any parameter connections for a given parameter `param_name` in a given component
`comp_def` in model `m`.
"""
@delegate disconnect_param!(m::Model, comp_name::Symbol, param_name::Symbol) => md

"""
    add_model_param!(m::Model, name::Symbol, value::ModelParameter)

Add an model parameter with name `name` and Model Parameter `value` to Model `m`.
"""
@delegate add_model_param!(m::Model, name::Symbol, value::ModelParameter) => md

"""
    add_model_param!(m: Model, name::Symbol, value::Number;
                    param_dims::Union{Nothing,Array{Symbol}} = nothing,
                    is_shared::Bool = false)

Create and add a model parameter with name `name` and Model Parameter `value`
to Model `m`. The Model Parameter will be created with value `value`, dimensions
`param_dims` which can be left to be created automatically from the Model Def, and
an is_shared attribute `is_shared` which defaults to false.
"""
@delegate add_model_param!(m::Model, name::Symbol,
                            value::Union{Number, AbstractArray, AbstractRange, Tuple};
                            param_dims::Union{Nothing,Array{Symbol}} = nothing) => md
"""
    add_model_param!(m::Model, name::Symbol,
                    value::Union{AbstractArray, AbstractRange, Tuple};
                    param_dims::Union{Nothing,Array{Symbol}} = nothing,
                    is_shared::Bool = false)

Create and add a model parameter with name `name` and Model Parameter `value`
to Model `m`. The Model Parameter will be created with value `value`, dimensions
`param_dims` which can be left to be created automatically from the Model Def, and
an is_shared attribute `is_shared` which defaults to false.
"""
@delegate add_model_param!(m::Model, name::Symbol,
                            value::Union{AbstractArray, AbstractRange, Tuple};
                            param_dims::Union{Nothing,Array{Symbol}} = nothing,
                            is_shared::Bool = false) => md
"""
    add_internal_param_conn(m::Model, conn::InternalParameterConnection)

Add internal parameter connection `conn` to model `m`.
"""
@delegate add_internal_param_conn!(m::Model, conn::InternalParameterConnection) => md

"""
    set_leftover_params!(m::Model, parameters::Dict)

Set all of the parameters in `Model` `m` that don't have a value and are not connected
to some other component to a value from a dictionary `parameters`. This method assumes
the dictionary keys are strings (or convertible into Strings ie. Symbols) that
match the names of unset parameters in the model, and all resulting new model
parameters will be shared parameters.

Note that this function `set_leftover_params! has been deprecated, and uses should
be transitioned to using `update_leftover_params!` with keys specific to component-parameter
pairs i.e. (comp_name, param_name) => value in the dictionary.
"""
@delegate set_leftover_params!(m::Model, parameters) => md

"""
    update_leftover_params!(m::Model, parameters::Dict)

Update all of the parameters in `Model` `m` that don't have a value and are not connected
to some other component to a value from a dictionary `parameters`. This method assumes
the dictionary keys are Tuples of Symbols (or convertible to Symbols ie. Strings)
of (comp_name, param_name) that match the component-parameter pair of
unset parameters in the model.  All resulting connected model parameters will be
unshared model parameters.
"""
@delegate update_leftover_params!(m::Model, parameters) => md

"""
    update_param!(m::Model, name::Symbol, value; update_timesteps = nothing)

Update the `value` of an model parameter in model `m`, referenced by
`name`. The update_timesteps keyword argument is deprecated, we keep it here
just to provide warnings.
"""
@delegate update_param!(m::Model, name::Symbol, value; update_timesteps = nothing) => md

"""
    update_param!(m::Model, comp_name::Symbol, param_name::Symbol, value)

Update the `value` of the unshared model parameter in Model `m`'s Model Def connected
to component `comp_name`'s parameter `param_name`.
"""
@delegate update_param!(m::Model, comp_name::Symbol, param_name::Symbol, value) => md

"""
    update_params!(m::Model, parameters::Dict; update_timesteps = nothing)

For each (k, v) in the provided `parameters` dictionary, `update_param!`
is called to update the model parameter identified by k to value v.

For updating unshared parameters, each key k must be a Tuple matching the name of a
component in `obj` and the name of an parameter in that component.

For updating shared parameters, each key k must be a symbol or convert to a symbol
matching the name of a shared model parameter that already exists in the model.
"""
@delegate update_params!(m::Model, parameters::Dict; update_timesteps = nothing) => md

"""
    add_comp!(
        m::Model, comp_id::ComponentId, comp_name::Symbol=comp_id.comp_name;
        first::NothingInt=nothing,
        last::NothingInt=nothing,
        before::NothingSymbol=nothing,
        after::NothingSymbol=nothing,
        rename::NothingPairList=nothing
    )

Add the component indicated by `comp_id` to the model indicated by `m`. The component is added
at the end of the list unless one of the keywords `before` or `after` is specified. Note
that a copy of `comp_id` is made in the composite and assigned the give name. The optional
argument `rename` can be a list of pairs indicating `original_name => imported_name`. The optional
arguments `first` and `last` indicate the times bounding the run period for the given component,
which must be within the bounds of the model and if explicitly set are fixed.  These default
to flexibly changing with the model's `:time` dimension.
"""
function add_comp!(m::Model, comp_id::ComponentId, comp_name::Symbol=comp_id.comp_name; kwargs...)
    comp_def = add_comp!(m.md, comp_id, comp_name; kwargs...)
    return ComponentReference(m.md, comp_name)
end

"""
    add_comp!(
        m::Model, comp_def::AbstractComponentDef, comp_name::Symbol=comp_id.comp_name;
        first::NothingInt=nothing,
        last::NothingInt=nothing,
        before::NothingSymbol=nothing,
        after::NothingSymbol=nothing,
        rename::NothingPairList=nothing
    )

Add the component `comp_def` to the model indicated by `m`. The component is added at
the end of the list unless one of the keywords, `first`, `last`, `before`, `after`. Note
that a copy of `comp_id` is made in the composite and assigned the give name. The optional
argument `rename` can be a list of pairs indicating `original_name => imported_name`. The optional
arguments `first` and `last` indicate the times bounding the run period for the given component,
which must be within the bounds of the model and if explicitly set are fixed.  These default
to flexibly changing with the model's `:time` dimension.
"""
function add_comp!(m::Model, comp_def::AbstractComponentDef, comp_name::Symbol=comp_def.comp_id.comp_name; kwargs...)
    return add_comp!(m, comp_def.comp_id, comp_name; kwargs...)
end

"""
    replace!(
        m::Model,
        old_new::Pair{Symbol, ComponentDef},
        before::NothingSymbol=nothing,
        after::NothingSymbol=nothing,
        reconnect::Bool=true
    )

For the pair `comp_name => comp_def` in `old_new`, replace the component with name `comp_name` in
the model `m` with the new component specified by `comp_def`. The new component is added
in the same position as the old component, unless one of the keywords `before` or `after` is
specified for a different position. The optional boolean argument `reconnect` with default value
`true` indicates whether the existing parameter connections should be maintained in the new
component. Returns a ComponentReference for the added component.
"""
function Base.replace!(m::Model, old_new::Pair{Symbol, ComponentDef}; kwargs...)
    comp_name, comp_def = old_new
    _replace!(m.md, comp_name => comp_def.comp_id; kwargs...)
    return ComponentReference(m.md, comp_name)
end

@delegate ComponentReference(m::Model, name::Symbol) => md

"""
    components(m::Model)

Return an iterator on the components in a model's model instance.
"""
@delegate components(m::Model) => mi

@delegate compdefs(m::Model) => md

@delegate compdef(m::Model, comp_name::Symbol) => md

@delegate Base.length(m::Model) => md

@delegate first_and_step(m::Model) => md

@delegate time_labels(m::Model) => md

"""
    getspan(m::Model, comp_name::Symbol)

Return the number of timesteps a given component in a model will run for.
"""
@delegate getspan(m::Model, comp_name::Symbol) => md

"""
    datumdef(comp_def::ComponentDef, item::Symbol)

Return a DatumDef for `item` in the given component `comp_def`.
"""
function datumdef(comp_def::AbstractComponentDef, item::Symbol)
    if has_variable(comp_def, item)
        return variable(comp_def, item)

    elseif has_parameter(comp_def, item)
        return parameter(comp_def, item)
    else
        error("Cannot access data item; :$item is not a variable or a parameter in component $(comp_def.comp_id).")
    end
end

"""
    datumdef(m::Model, comp_name::Symbol, item::Symbol)

Return a DatumDef for `item` in the given component `comp_name` of model `m`.
"""
datumdef(m::Model, comp_name::Symbol, item::Symbol) = datumdef(compdef(m.md, comp_name), item)

"""
    dim_names(m::Model, comp_name::Symbol, datum_name::Symbol)

Return the dimension names for the variable or parameter `datum_name`
in the given component `comp_name` in model `m`.
"""
function dim_names(m::Model, comp_name::Symbol, datum_name::Symbol)
    # the line below would work if the comp_name is in the top level of components in m's component structure
    # return dim_names(compdef(m, comp_name), datum_name)

    paths = _get_all_paths(m)
    comp_path = paths[comp_name]
    comp_def = find_comp(m, comp_path)
    return dim_names(comp_def, datum_name)
end

dim_names(mm::MarginalModel, comp_name::Symbol, datum_name::Symbol) = dim_names(mm.base, comp_name, datum_name)

@delegate dimension(m::Model, dim_name::Symbol) => md

# Allow access of the form my_model[:grosseconomy, :tfp]
@delegate Base.getindex(m::Model, comp_name::Symbol, datum_name::Symbol) => mi

"""
    dim_count(m::Model, dim_name::Symbol)

Return the size of index `dim_name` in model `m`.
"""
@delegate dim_count(m::Model, dim_name::Symbol) => md
@delegate dim_counts(m::Model, dims::Vector{Symbol}) => md
@delegate dim_count_dict(m::Model) => md

"""
    dim_keys(m::Model, dim_name::Symbol)

Return keys for dimension `dim-name` in model `m`.
"""
@delegate dim_keys(m::Model, dim_name::Symbol) => md
"""
     dim_keys(mi::ModelInstance, dim_name::Symbol)

 Return keys for dimension `dim-name` in model instance `mi`.
 """
 @delegate dim_keys(mi::ModelInstance, dim_name::Symbol) => md
"""
    dim_key_dict(m::Model)

Return a dict of dimension keys for all dimensions in model `m`.
"""
@delegate dim_key_dict(m::Model) => md
"""
    dim_values(m::Model, name::Symbol)

Return values for dimension `name` in Model `m`.
"""
@delegate dim_values(m::Model, name::Symbol) => md
"""
    dim_value_dict(m::Model)

Return a dictionary of the values of all dimensions in Model `m`.
"""
@delegate dim_value_dict(m::Model) => md
"""
    set_dimension!(m::Model, name::Symbol, keys::Union{Vector, Tuple, AbstractRange})

Set the values of `m` dimension `name` to integers 1 through `count`, if `keys`` is
an integer; or to the values in the vector or range if `keys`` is either of those types.
"""
@delegate set_dimension!(m::Model, name::Symbol, keys::Union{Int, Vector, Tuple, AbstractRange}) => md

@delegate check_parameter_dimensions(m::Model, value::AbstractArray, dims::Vector, name::Symbol) => md

@delegate parameter_names(m::Model, comp_name::Symbol) => md

@delegate parameter_dimensions(m::Model, comp_name::Symbol, param_name::Symbol) => md

@delegate parameter_unit(m::Model, comp_name::Symbol, param_name::Symbol) => md

parameter(m::Model, comp_def::ComponentDef, param_name::Symbol) = parameter(comp_def, param_name)

parameter(m::Model, comp_name::Symbol, param_name::Symbol) = parameter(m, compdef(m, comp_name), param_name)

"""
    parameters(m::Model, comp_name::Symbol)

Return a list of the parameter definitions for `comp_name` in model `m`.
"""
parameters(m::Model, comp_name::Symbol) = parameters(compdef(m, comp_name))

variable(m::Model, comp_name::Symbol, var_name::Symbol) = variable(compdef(m, comp_name), var_name)

@delegate variable_unit(m::Model, comp_path::ComponentPath, var_name::Symbol) => md

@delegate variable_dimensions(m::Model, comp_path::ComponentPath, var_name::Symbol) => md
@delegate variable_dimensions(m::Model, comp::Symbol, var_name::Symbol) => md
@delegate variable_dimensions(m::Model, comp_path::NTuple{N, Symbol} where N, var_name::Symbol) => md

"""
    variables(m::Model, comp_name::Symbol)

Return an iterator on the variable definitions for `comp_name` in model `m`.
"""
variables(m::Model, comp_name::Symbol) = variables(compdef(m, comp_name))

@delegate variable_names(m::Model, comp_name::Symbol) => md

"""
    add_shared_param!(m::Model, name::Symbol, value::Any; dims::Array{Symbol}=Symbol[], datatype::DataType=Nothing)

User-facing API function to add a shared parameter to Model `m` with name
`name` and value `value`, and an array of dimension names `dims` which dfaults to
an empty vector.  The `is_shared` attribute of the added Model Parameter will be `true`.

The `value` can by a scalar, an array, or a NamedAray. Optional keyword argument 'dims' is a list
of the dimension names of the provided data, and will be used to check that they match the
model's index labels.  This must be included if the `value` is not a scalar, and defaults
to an empty vector. Optional keyword argument `datatype` allows user to specify a datatype
to use for the shared model parameter.
"""
@delegate add_shared_param!(m::Model, name::Symbol, value::Any; dims::Array{Symbol}=Symbol[], data_type::DataType=Nothing) => md

"""
    add_model_array_param!(m::Model, name::Symbol, value::Union{AbstractArray, TimestepArray}, dims)

Add a one or two dimensional (optionally, time-indexed) array parameter `name`
with value `value` to the model `m`.
"""
@delegate add_model_array_param!(m::Model, name::Symbol, value::Union{AbstractArray, TimestepArray}, dims) => md

"""
    add_model_scalar_param!(m::Model, name::Symbol, value::Any)

Add a scalar type parameter `name` with value `value` to the model `m`.
"""
@delegate add_model_scalar_param!(m::Model, name::Symbol, value::Any) => md

"""
    delete!(m::Model, component::Symbol; deep::Bool=false)

Delete a `component` by name from a model `m`'s ModelDef, and nullify the ModelInstance.
If `deep=true` then any model model parameters connected only to
this component will also be deleted.
"""
@delegate Base.delete!(m::Model, comp_name::Symbol; deep::Bool=false) => md

"""
    delete_param!(m::Model, model_param_name::Symbol)

Delete `model_param_name` from a model `m`'s ModelDef's list of model parameters, and
also remove all external parameters connections that were connected to `model_param_name`.
"""
@delegate delete_param!(m::Model, model_param_name::Symbol) => md

"""
    set_param!(m::Model, comp_name::Symbol, param_name::Symbol, value; dims=nothing)

Set the parameter of a component `comp_name` in a model `m` to a given `value`.
The `value` can by a scalar, an array, or a NamedAray. Optional keyword argument 'dims'
is a list of the dimension names of the provided data, and will be used to check
that they match the model's index labels.
"""
@delegate set_param!(m::Model, comp_name::Symbol, param_name::Symbol, value; dims=nothing) => md

"""
    set_param!(m::Model, comp_name::Symbol, param_name::Symbol, model_param_name::Symbol, value; dims=nothing)

Set the parameter `param_name` of a component `comp_name` in a model `m` to a given `value`,
storing the value in the model's parameter list by the provided name `model_param_name`.
The `value` can by a scalar, an array, or a NamedAray. Optional keyword argument 'dims'
is a list of the dimension names of the provided data, and will be used to check
that they match the model's index labels.
"""
@delegate set_param!(m::Model, comp_name::Symbol, param_name::Symbol, model_param_name::Symbol, value; dims=nothing) => md

"""
    set_param!(m::Model, param_name::Symbol, value; dims=nothing)

Set the value of a parameter in all components of the model that have a parameter of
the specified name.
"""
@delegate set_param!(m::Model, param_name::Symbol, value; dims=nothing, ignoreunits::Bool=false) => md

@delegate import_params!(m::Model) => md

"""
    Base.run(m::Model; ntimesteps::Int=typemax(Int), rebuild::Bool=false,
            dim_keys::Union{Nothing, Dict{Symbol, Vector{T} where T <: DimensionKeyTypes}}=nothing)

Run model `m` once.
"""
function Base.run(m::Model; ntimesteps::Int=typemax(Int), rebuild::Bool=false,
                  dim_keys::Union{Nothing, Dict{Symbol, Vector{T} where T <: DimensionKeyTypes}}=nothing)
    if length(m) == 0
        error("Cannot run a model with no components.")
    end

    if (rebuild || ! is_built(m))
        build!(m)
    end

    # println("Running model...")
    mi = modelinstance(m)
    run(mi, ntimesteps, dim_keys)
    nothing
end

##
## DEPRECATIONS - Should move from warning --> error --> removal
##

# -- throw errors --

"""
    replace_comp!(
        m::Model, comp_def::ComponentDef, comp_name::Symbol=comp_id.comp_name;
        before::NothingSymbol=nothing,
        after::NothingSymbol=nothing,
        reconnect::Bool=true
    )

Deprecated function for replacing the component with name `comp_name` in model `m` with the
new component specified by `comp_def`. Use the following syntax instead:

`replace!(m, comp_name => comp_def; kwargs...)`

See docstring for `replace!` for further description of available functionality.
"""
function replace_comp!(m::Model, comp_def::ComponentDef, comp_name::Symbol=comp_def.comp_id.comp_name; kwargs...)
    error("Function `replace_comp!(m, comp_def, comp_name; kwargs...)` has been deprecated. Use `replace!(m, comp_name => comp_def; kwargs...)` instead.")
end

"""
    replace_comp!(
        m::Model, comp_id::ComponentId, comp_name::Symbol=comp_id.comp_name;
        before::NothingSymbol=nothing,
        after::NothingSymbol=nothing,
        reconnect::Bool=true
    )

Deprecated function for replacing the component with name `comp_name` in model `m` with the
new component specified by `comp_id`. Use the following syntax instead:

`replace!(m, comp_name => Mimi.compdef(comp_id); kwargs...)`

See docstring for `replace!` for further description of available functionality.
"""
function replace_comp!(m::Model, comp_id::ComponentId, comp_name::Symbol=comp_id.comp_name; kwargs...)
    error("Function `replace_comp!(m, comp_id, comp_name; kwargs...)` has been deprecated. Use `replace!(m, comp_name => Mimi.compdef(comp_id); kwargs...)` instead.")
end

# -- throw warnings --

@delegate set_external_param!(m::Model, name::Symbol,
                              value::Union{Number, AbstractArray, AbstractRange, Tuple};
                              param_dims::Union{Nothing,Array{Symbol}} = nothing) => md

@delegate set_external_param!(m::Model, name::Symbol, value::ModelParameter) => md
@delegate set_external_array_param!(m::Model, name::Symbol, value::Union{AbstractArray, TimestepArray}, dims) => md
@delegate set_external_scalar_param!(m::Model, name::Symbol, value::Any) => md
@delegate external_params(m::Model) => md
@delegate external_param(m::Model, name::Symbol; missing_ok=false) => md


text:
#*******************************************************/
#* Copyright(c) 2018 by Artelys                        */
#* This source code is subject to the terms of the     */
#* MIT Expat License (see LICENSE.md)                  */
#*******************************************************/

#++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
# This example demonstrates how to use Knitro to solve the following
# mixed-integer nonlinear optimization problem(MINLP).  This model
# is test problem 1(Synthesis of processing system) in
# M. Duran & I.E. Grossmann,  "An outer approximation algorithm for
# a class of mixed integer nonlinear programs", Mathematical
# Programming 36, pp. 307-339, 1986.  The problem also appears as
# problem synthes1 in the MacMINLP test set.
#
# min   5 x3 + 6 x4 + 8 x5 + 10 x0 - 7 x2 -18 log(x1 + 1)
#      - 19.2 log(x0 - x1 + 1) + 10
# s.t.  0.8 log(x1 + 1) + 0.96 log(x0 - x1 + 1) - 0.8 x2 >= 0
#       log(x1 + 1) + 1.2 log(x0 - x1 + 1) - x2 - 2 x5 >= -2
#       x1 - x0 <= 0
#       x1 - 2 x3 <= 0
#       x0 - x1 - 2 x4 <= 0
#       x3 + x4 <= 1
#       0 <= x0 <= 2
#       0 <= x1 <= 2
#       0 <= x2 <= 1
#       x0, x1, x2 continuous
#       x3, x4, x5 binary
#
# The solution is(1.30098, 0, 1, 0, 1, 0).
#++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

using KNITRO
using Test

function example_minlp1(; verbose=true)
    #*------------------------------------------------------------------*
    #*     FUNCTION callbackEvalFC                                      *
    #*------------------------------------------------------------------*
    # The signature of this function matches KNITRO.KN_eval_callback in knitro.h.
    # Only "obj" and "c" are set in the KNITRO.KN_eval_result structure.
    function callbackEvalFC(kc, cb, evalRequest, evalResult, userParams)
        x = evalRequest.x

        # Evaluate nonlinear objective structure
        dTmp1 = x[1] - x[2] + 1.0
        dTmp2 = x[2] + 1.0
        evalResult.obj[1] = -18.0 * log(dTmp2) - 19.2 * log(dTmp1)

        # Evaluate nonlinear constraint structure
        evalResult.c[1] = 0.8 * log(dTmp2) + 0.96 * log(dTmp1)
        evalResult.c[2] = log(dTmp2) + 1.2 * log(dTmp1)

        return 0
    end

    #*------------------------------------------------------------------*
    #*     FUNCTION callbackEvalGA                                      *
    #*------------------------------------------------------------------*
    # The signature of this function matches KNITRO.KN_eval_callback in knitro.h.
    # Only "objGrad" and "jac" are set in the KNITRO.KN_eval_result structure.
    function callbackEvalGA(kc, cb, evalRequest, evalResult, userParams)
        x = evalRequest.x

        # Evaluate gradient of nonlinear objective structure
        dTmp1 = x[1] - x[2] + 1.0
        dTmp2 = x[2] + 1.0
        evalResult.objGrad[1] = -(19.2 / dTmp1)
        evalResult.objGrad[2] = (-18.0 / dTmp2) + (19.2 / dTmp1)

        # Gradient of nonlinear structure in constraint 0.
        evalResult.jac[1] = 0.96 / dTmp1                      # wrt x0
        evalResult.jac[2] = (-0.96 / dTmp1) + (0.8 / dTmp2)   # wrt x1
        # Gradient of nonlinear structure in constraint 1.
        evalResult.jac[3] = 1.2 / dTmp1                       # wrt x0
        evalResult.jac[4] = (-1.2 / dTmp1) + (1.0 / dTmp2)    # wrt x1

        return 0
    end

    #*------------------------------------------------------------------*
    #*     FUNCTION callbackEvalH                                       *
    #*------------------------------------------------------------------*
    # The signature of this function matches KNITRO.KN_eval_callback in knitro.h.
    # Only "hess" or "hessVec" are set in the KNITRO.KN_eval_result structure.
    function callbackEvalH(kc, cb, evalRequest, evalResult, userParams)
        x = evalRequest.x
        lambda_ = evalRequest.lambda
        # Scale objective component of hessian by sigma
        sigma = evalRequest.sigma

        # Evaluate the non-zero components in the Hessian of the Lagrangian.
        # Note: Since the Hessian is symmetric, we only provide the
        #       nonzero elements in the upper triangle(plus diagonal).
        dTmp1 = x[1] - x[2] + 1.0
        dTmp2 = x[2] + 1.0
        evalResult.hess[1] = sigma * (19.2 /(dTmp1 * dTmp1)) +
                                lambda_[1] * (-0.96 / (dTmp1 * dTmp1)) + lambda_[2] * (-1.2 / (dTmp1 * dTmp1))
        evalResult.hess[2] = sigma * (-19.2 /(dTmp1 * dTmp1)) +
                                lambda_[1] * (0.96 / (dTmp1 * dTmp1)) + lambda_[2] * (1.2 / (dTmp1 * dTmp1))
        evalResult.hess[3] = sigma * ((19.2 / (dTmp1 * dTmp1)) +(18.0 / (dTmp2 * dTmp2))) +
                                lambda_[1] * ((-0.96 / (dTmp1 * dTmp1)) -(0.8 / (dTmp2 * dTmp2))) +
                                lambda_[2] * ((-1.2 / (dTmp1 * dTmp1)) -(1.0 / (dTmp2 * dTmp2)))

        return 0
    end

    #*------------------------------------------------------------------*
    #*     FUNCTION callbackProcessNode                                 *
    #*------------------------------------------------------------------*
    # The signature of this function matches KNITRO.KN_user_callback in knitro.h.
    # Argument "kcSub" is the context pointer for the last node
    # subproblem solved inside Knitro.  The application level context
    # pointer is passed in through "userParams".
    function callbackProcessNode(kc, x, lambda_, userParams)
        # The Knitro context pointer was passed in through "userParams".

        # Print info about the status of the MIP solution.
        numNodes = KNITRO.KN_get_mip_number_nodes(kc)
        relaxBound = KNITRO.KN_get_mip_relaxation_bnd(kc)
        # Note: To retrieve solution information about the node subproblem
        # we need to pass in "kcSub" here.
        nodeObj = KNITRO.KN_get_obj_value(kc)
        mip_io = KNITRO.KN_get_mip_incumbent_obj(kc)
        mip_ag = KNITRO.KN_get_mip_abs_gap(kc)
        mip_rg = KNITRO.KN_get_mip_rel_gap(kc)

        if verbose
            println("callbackProcessNode:")
            println("    Node number    = ", numNodes)
            println("    Node objective = ", nodeObj)
            println("    Current relaxation bound = ", relaxBound)
            println("    Current incumbent bound  = ", mip_io)
            println("    Absolute integrality gap = ", mip_ag)
            println("    Relative integrality gap = ", mip_rg)
        end

        # User defined termination example.
        # Uncomment below to force termination after 3 nodes.
        #if(numNodes == 3)
        #    return KNITRO.KN_RC_USER_TERMINATION

        return 0
    end

    #*------------------------------------------------------------------*
    #*     main                                                         *
    #*------------------------------------------------------------------*

    # Create a new Knitro solver instance.
    kc = KNITRO.KN_new()

    # Illustrate how to override default options.
    KNITRO.KN_set_param(kc, "mip_method", KNITRO.KN_MIP_METHOD_BB)
    KNITRO.KN_set_param(kc, "algorithm", KNITRO.KN_ALG_ACT_CG)
    KNITRO.KN_set_param(kc, "outmode", KNITRO.KN_OUTMODE_SCREEN)
    KNITRO.KN_set_param(kc, KNITRO.KN_PARAM_OUTLEV, KNITRO.KN_OUTLEV_ALL)
    KNITRO.KN_set_param(kc, KNITRO.KN_PARAM_MIP_OUTINTERVAL, 1)
    KNITRO.KN_set_param(kc, KNITRO.KN_PARAM_MIP_MAXNODES, 10000)

    # Initialize Knitro with the problem definition.

    # Add the variables and set their bounds and types.
    # Note: any unset lower bounds are assumed to be
    # unbounded below and any unset upper bounds are
    # assumed to be unbounded above.
    n = 6
    KNITRO.KN_add_vars(kc, n)
    KNITRO.KN_set_var_lobnds(kc, zeros(Float64, n))
    KNITRO.KN_set_var_upbnds(kc, [2., 2., 1., 1., 1., 1.])
    KNITRO.KN_set_var_types(kc, [KNITRO.KN_VARTYPE_CONTINUOUS,
                                KNITRO.KN_VARTYPE_CONTINUOUS,
                                KNITRO.KN_VARTYPE_CONTINUOUS,
                                KNITRO.KN_VARTYPE_BINARY,
                                KNITRO.KN_VARTYPE_BINARY,
                                KNITRO.KN_VARTYPE_BINARY])

    # Note that variables x2..x5 only appear linearly in the
    # problem.  We mark them as linear variables, which may
    # help Knitro do more extensive presolving resulting in
    # faster solves.
    for i in 2:(n-1)
        KNITRO.KN_set_var_property(kc, i, KNITRO.KN_VAR_LINEAR)
    end

    # Add the constraints and set their bounds
    KNITRO.KN_add_cons(kc, 6)
    KNITRO.KN_set_con_lobnds(kc,  [0, -2,
                                -KNITRO.KN_INFINITY,
                                -KNITRO.KN_INFINITY,
                                -KNITRO.KN_INFINITY,
                                -KNITRO.KN_INFINITY])
    KNITRO.KN_set_con_upbnds(kc, [KNITRO.KN_INFINITY, KNITRO.KN_INFINITY, 0, 0, 0, 1])

    # Add the linear structure in the objective function.
    objGradIndexVars = Int32[3, 4, 5, 0, 2]
    objGradCoefs = [5.0, 6.0, 8.0, 10.0, -7.0]
    KNITRO.KN_add_obj_linear_struct(kc, objGradIndexVars, objGradCoefs)

    # Add the constant in the objective function.
    KNITRO.KN_add_obj_constant(kc, 10.0)

    # Load the linear structure for all constraints at once.
    jacIndexCons = Int32[0, 1, 1, 2, 2, 3, 3, 4, 4, 4, 5, 5]
    jacIndexVars = Int32[2, 2, 5, 1, 0, 1, 3, 0, 1, 4, 3, 4]
    jacCoefs = [-0.8, -1.0, -2.0, 1.0, -1.0, 1.0, -2.0, 1.0, -1.0, -2.0, 1.0, 1.0]
    KNITRO.KN_add_con_linear_struct(kc, jacIndexCons, jacIndexVars, jacCoefs)

    # Add a callback function "callbackEvalFC" to evaluate the nonlinear
    # structure in the objective and first two constraints.  Note that
    # the linear terms in the objective and first two constraints were
    # added above in "KNITRO.KN_add_obj_linear_struct()" and
    # "KNITRO.KN_add_con_linear_struct()" and will not be specified in the
    # callback.
    cIndices = Int32[0, 1] # Constraint indices for callback
    cb = KNITRO.KN_add_eval_callback(kc, true, cIndices, callbackEvalFC)

    # Also add a callback function "callbackEvalGA" to evaluate the
    # gradients of all nonlinear terms specified in the callback.  If
    # not provided, Knitro will approximate the gradients using finite-
    # differencing.  However, we recommend providing callbacks to
    # evaluate the exact gradients whenever possible as this can
    # drastically improve the performance of Knitro.
    # Objective gradient non-zero structure for callback
    objGradIndexVarsCB = Int32[0, 1]
    # Constraint Jacobian non-zero structure for callback
    jacIndexConsCB = Int32[0, 0, 1, 1]
    jacIndexVarsCB = Int32[0, 1, 0, 1]
    KNITRO.KN_set_cb_grad(kc, cb, callbackEvalGA,
                        nV=length(objGradIndexVarsCB),
                        objGradIndexVars=objGradIndexVarsCB,
                        jacIndexCons=jacIndexConsCB,
                        jacIndexVars=jacIndexVarsCB)

    hessIndexVars1CB = Int32[0, 0, 1]
    hessIndexVars2CB = Int32[0, 1, 1]
    # Add a callback function "callbackEvalH" to evaluate the Hessian
    #(i.e. second derivative matrix) of the objective.  If not specified,
    # Knitro will approximate the Hessian. However, providing a callback
    # for the exact Hessian(as well as the non-zero sparsity structure)
    # can greatly improve Knitro performance and is recommended if possible.
    # Since the Hessian is symmetric, only the upper triangle is provided.
    KNITRO.KN_set_cb_hess(kc, cb, length(hessIndexVars1CB), callbackEvalH,
                        hessIndexVars1=hessIndexVars1CB,
                        hessIndexVars2=hessIndexVars2CB)

    # Specify that the user is able to provide evaluations
    #  of the Hessian matrix without the objective component.
    #  turned off by default but should be enabled if possible.
    KNITRO.KN_set_param(kc, KNITRO.KN_PARAM_HESSIAN_NO_F, KNITRO.KN_HESSIAN_NO_F_ALLOW)

    # Set minimize or maximize(if not set, assumed minimize)
    KNITRO.KN_set_obj_goal(kc, KNITRO.KN_OBJGOAL_MINIMIZE)

    # Set a callback function that performs some user-defined task
    # after completion of each node in the branch-and-bound tree.
    KNITRO.KN_set_mip_node_callback(kc, callbackProcessNode)

    kn_outlev = verbose ? KNITRO.KN_OUTLEV_ALL : KNITRO.KN_OUTLEV_NONE
    KNITRO.KN_set_param(kc, KNITRO.KN_PARAM_OUTLEV, kn_outlev)

    # Solve the problem.
    #
    # Return status codes are defined in "knitro.h" and described
    # in the Knitro manual.

    nStatus = KNITRO.KN_solve(kc)
    # An example of obtaining solution information.
    nSTatus, objSol, x, lambda_ = KNITRO.KN_get_solution(kc)

    if verbose
        println("Optimal objective value  = ", objSol)
        println("Optimal x")
        for i in 1:n
            println("  x[$i] = ", x[i])
        end
    end

    # Delete the Knitro solver instance.
    KNITRO.KN_free(kc)

    @testset "Example minlp1" begin
        @test nStatus == 0
        @test objSol ≈ 6.0097589
        @test x ≈ [1.30097589, 0., 1., 0., 1., 0.] atol=1e-5
    end
end

example_minlp1(; verbose=isdefined(Main, :KN_VERBOSE) ? KN_VERBOSE : true)


text:
using DrWatson
@quickactivate "ExercisesRepo"
include(srcdir("style.jl"))

# %% 1d climate
using DynamicalSystems, PyPlot, Roots, ForwardDiff

# %% Bifurcation versus ε
fig = figure(figsize = (figx/2, figx/3))
ax = gca()
dTdt(T, ε=0.65) = 0.5 + (0.4/π)*atan((T-263)/2) - 10*ε*(0.002T)^4
es = 0.3:0.005:0.9
for e in es
    f = (T) -> dTdt(T, e)
    roots = Roots.find_zeros(f, 50, 500)
    for (i, r) in enumerate(roots)
        c = ForwardDiff.derivative(f, r) < 0 ? "k" : "w"
        plot(e, r, marker = "o", mec = "k", mew = 1,
        markersize = 6, mfc = c )
    end
end
plot([],[], ls = "None", marker = "o", mec = "k", mfc = "k", label = "stable")
plot([],[], ls = "None", marker = "o", mec = "k", mfc = "w", label = "unstable")
legend(markerscale = 3)
xlabel("\$\\epsilon\$")
ylabel("\$T^*\$")
tight_layout()

e1, T1 = 0.48755, 256.38
e2, T2 = 0.78661, 268.075
for (e, T) in ((e1, T1), (e2, T2))
    # plot(e, T; marker="o", ms = 10, color = "k", mew = 1, mec = "C2",
    #            fillstyle = "top", markerfacecoloralt= "w")
    plot(e, T; marker="o", ms = 10, color = "C2")
end

fig.subplots_adjust(left = 0.18, bottom = 0.2, top = 0.97)
# fsave(joinpath(figdir, "bif_example"), fig)

# %% Period doubling
using DynamicalSystems, PyPlot

fig = figure(figsize = (figx, 1.2figx))
axlo = fig.add_subplot(3, 1, 1)
axro = fig.add_subplot(3, 1, 3)

zoomaxs = []
for i in 1:3
	push!(zoomaxs, fig.add_subplot(3, 3, 3+i))
end

axlo.clear()
ds = Systems.logistic(0.5)
i = 1
pvalues = range(2.5, 4; length = 2001)
n = 2000
Ttr = 2000
p_index = 1
output = orbitdiagram(ds, i, p_index, pvalues; n = n, Ttr = Ttr)

L = length(pvalues)
x = Vector{Float64}(undef, n*L)
y = copy(x)
for j in 1:L
    x[(1 + (j-1)*n):j*n] .= pvalues[j]
    y[(1 + (j-1)*n):j*n] .= output[j]
end
axlo.plot(x, y, ls = "None", ms = 0.5, color = "black", marker = "o", alpha = 0.01)
axlo.set_xlim(pvalues[1], pvalues[end]); axlo.set_ylim(0,1)
axlo.set_ylabel("\$x_\\infty\$", labelpad = 10);
axlo.text(2.6, 0.1, "logistic map", size = 36)
r∞ = 3.56995
axlo.axvline(r∞; color = "C3", alpha = 0.7, ls = "dashed")

# Zoomin boxes
zbox1 = ((3.629, 0.295), (3.640, 0.32))
zbox2 = ((3.6338, 0.304), (3.6342, 0.3065))
zbox3 = ((3.63407, 0.30485), (3.63408, 0.3058))
allboxes = (zbox1, zbox2, zbox3)
for (k, zbox) in enumerate(allboxes)
	pval = range(zbox[1][1], zbox[2][1]; length = 1001)
	nsm = n÷3
	output = orbitdiagram(ds, i, p_index, pval;
		n=nsm, Ttr, ulims = (zbox[1][2], zbox[2][2])
	)
	L = length(pval)
	x = Vector{Float64}(undef, nsm*L)
	y = copy(x)
	for j in 1:L
	    x[(1 + (j-1)*nsm):j*nsm] .= pval[j]
	    y[(1 + (j-1)*nsm):j*nsm] .= output[j]
	end
	zoomaxs[k].clear()
	zoomaxs[k].plot(x, y, ls = "None", ms = 0.5, color = "black", marker = "o", alpha = 0.01)
	zoomaxs[k].set_xlim(zbox[1][1], zbox[2][1])
	zoomaxs[k].set_ylim(zbox[1][2], zbox[2][2])
	zoomaxs[k].axis("off")
end

zbox0 = ((pvalues[1], 0.0), (pvalues[end], 1.0))
allzoomaxs = [axlo, zoomaxs...]
for i in 1:3
	origin = allzoomaxs[i]
	zoomin = allzoomaxs[i+1]
    zbox = allboxes[i]
	c = i != 3 ? "C$(i)" : "C$(i+2)"
    axis_zoomin!(zoomin, origin, zbox, zbox, c, 1.0,
	connect_lines = i > 1, lw = 4.0)
end

# arrow indicator
axlo.arrow(
	(zbox1[1][1] + zbox1[2][1])/2, 0.05, 0, zbox1[2][2] - 0.1;
	color = "C1", width = 0.005, length_includes_head = true, head_length = 0.1
)

# Roessler
ro = Systems.roessler()
pvalues = range(1.5, stop = 5.18, length = 2001)
i = 1
plane = (2, 0.0)
tf = 4000.0
p_index = 3
output = produce_orbitdiagram(ro, plane, i, p_index, pvalues;
                              tfinal = tf, Ttr = 2000.0)

for (j, p) in enumerate(pvalues)
    axro.plot(fill(p, length(output[j])), output[j], lw = 0,
    marker = "o", ms = 0.2, color = "black", alpha = 0.1)
end
axro.set_xlim(pvalues[1], pvalues[end]);
axro.set_ylim(minimum(minimum(o) for o in output), maximum(maximum(o) for o in output))
axro.set_ylabel("\$x_\\infty\$ (section)", labelpad = 25)
axro.set_xlabel("parameter (\$r\$ for logistic, \$c\$ for Rossler)")
axro.set_yticks(3:3:9)
axro.text(1.7, 7, "Rössler system", size = 36)
fig.tight_layout()
fig.subplots_adjust(bottom = 0.1, top = 0.97, left = 0.09, hspace = 0.3, wspace = 0.1, right = 0.98)
# fsave(plotsdir("orbit_diagrams"), fig)

# %% # Testing how orbit diagram of quasiperiodic looks like
using DynamicalSystems, PyPlot

hh = Systems.henonheiles([0.0, 0.1, 0.5, 0.0])
psos = poincaresos(hh, (1, 0.0), 5000.0)

figure()
for i in 0:0.001:1
    plot(fill(i, length(psos)), psos[:, 2], lw = 0,
    marker = "o", ms = 0.2, color = "black", alpha = 0.1)
end

sm = Systems.standardmap([π/2, 0.1])
n = 1000
Ttr = 2000
pvalues = range(0, 1.2; length = 2001)
output = orbitdiagram(sm, 2, 1, pvalues; n = n, Ttr = Ttr)
figure()
ax = gca()
L = length(pvalues)
x = Vector{Float64}(undef, n*L)
y = copy(x)
for j in 1:L
    x[(1 + (j-1)*n):j*n] .= pvalues[j]
    y[(1 + (j-1)*n):j*n] .= output[j]
end
ax.plot(x, y, ls = "None", ms = 0.5, color = "black", marker = "o", alpha = 0.02)

# %% Roessler trajectory -> PSOS -> Lorenz map
using InteractiveChaos
using DynamicalSystems, Makie
using AbstractPlotting.MakieLayout

N = 10000.0

ds = Systems.roessler()
tr = trajectory(ds, N; Ttr = 100.0)

scene, layout = layoutscene(resolution = (2100, 700), )
display(scene)

# Plot trajectory

trplot = layout[1,1] = LScene(scene, scenekw = (camera = cam3d!, raw = false))
lines!(trplot, columns(tr)...; color = COLORSCHEME[1], linewidth = 2.0)

# Adjust ticks and sizes of 3D plot
xticks!(trplot.scene; xtickrange = [-10, 0, 10], xticklabels = string.([-10, 0, 10]))
yticks!(trplot.scene; ytickrange = [-10, 0, 10], yticklabels = string.([-10, 0, 10]))
zticks!(trplot.scene; ztickrange = [5, 15, 25], zticklabels = string.([5, 15, 25]))
trplot.scene[Axis][:ticks][:textsize] = (15, 15, 15)
trplot.scene[Axis][:names][:textsize] = (20,20,20)

# Plot plane and section
o = Point3f0(-10, 0, 0)
w = Point3f0(25, 0, 25)
p = FRect3D(o, w)
a = RGBAf0(0,0,0,0)
c = to_color(COLORSCHEME[2])
img = AbstractPlotting.ImagePattern([c a; a c]); # This throws an error if it shows, you can ignore that
mesh!(trplot, p; color = img)

psos = poincaresos(ds, (2, 0.0), N; Ttr = 100.0)
scatter!(trplot, columns(psos)...; color = COLORSCHEME[3], markersize = 500)

# Plot section separately
psplot = layout[1, 2] = LAxis(scene)
scatter!(psplot, psos[:, 1], psos[:, 3]; color = COLORSCHEME[3])
psplot.xlabel = "xₙ"
psplot.ylabel = "zₙ"

LS = 50
TS = 40

psplot.xticklabelsize = TS
psplot.yticklabelsize = TS
psplot.xlabelsize = LS
psplot.ylabelsize = LS

# Plot lorenz map
loplot = layout[1, 3] = LAxis(scene)
scatter!(loplot, psos[1:end-1, 3], psos[2:end, 3]; color = COLORSCHEME[3])
loplot.xlabel = "zₙ"
loplot.ylabel = "zₙ₊₁"
loplot.xticklabelsize = TS
loplot.yticklabelsize = TS
loplot.xlabelsize = LS
loplot.ylabelsize = LS
display(scene)

# Makie.save(plotsdir("roessler_map.png"), scene)

# %% Bifurcation kit code
# Better check https://rveltz.github.io/BifurcationKit.jl/dev/iterator/#
using BifurcationKit, SparseArrays, LinearAlgebra
const BK = BifurcationKit
using ForwardDiff

F(T, ε=0.65) = @. 0.5 + (0.4/π)*atan((T-263)/2) - 10*ε*(0.002T)^4

function deriv(x, p)
	f = (x) -> F(x, p)
	return ForwardDiff.derivative.(Ref(f), x)
end
Jac_m = (x, p) -> diagm(0 => deriv(x, p))

# parameters for the continuation
opts = ContinuationPar(
	dsmax = 0.1, dsmin = 1e-3,
	ds = 0.01, maxSteps = 2000, pMin = 0.0, pMax = 0.9, saveSolEveryStep = 0,
	newtonOptions = BK.NewtonPar(tol = 1e-8, verbose = false),
	detectBifurcation = true
)

# we define an iterator to hold the continuation routine
p0 = 0.3
u0 = [350.0]
iter = BK.PALCIterable(F, Jac_m, u0, p0, (BK.@lens _), opts; verbosity = 2)

resp = Float64[]
resx = Float64[]
lambdas = Float64[]

# this is the PALC algorithm
for state in iter
	# we save the current solution on the branch
	push!(resx, getx(state)[1])
	push!(resp, getp(state))
	BK.computeEigenvalues!(iter, state)
	eig = state.eigvals
	push!(lambdas, real(eig[1]))
end

fig = figure(figsize = (figx/2, figx/3))
ax = gca()
i1 = findfirst(l -> l ≥ 0, lambdas)
i2 = findlast(l -> l ≥ 0, lambdas)
ax.plot(resp[1:i1], resx[1:i1], color = "C3", label = "stable")
ax.plot(resp[i1:i2], resx[i1:i2], color = "C3", ls = "--", label = "unstable")
ax.plot(resp[i2+1:end], resx[i2+1:end], color = "C3")
ax.plot([resp[i1], resp[i2]], [resx[i1], resx[i2]], ls = "None",
		marker="o", ms = 10, color = "k", mew = 1, mec = "k",
        fillstyle = "top", markerfacecoloralt= "w"
)

# PyPlot.plot(resp, resx)
ax.legend(markerscale = 3)
ax.set_xlabel("\$\\epsilon\$")
ax.set_ylabel("\$T^*\$")
fig.tight_layout()
fig.subplots_adjust(left = 0.18, bottom = 0.2, top = 0.97)
# fsave(joinpath(figdir, "bif_example"), fig)

# %% Intermittency in Roessler
using DynamicalSystems, PyPlot

axro = (figure(); gca())
ro = Systems.roessler()
pvalues = range(4., stop = 8.18, length = 2001)
i = 1
plane = (2, 0.0)
tf = 4000.0
p_index = 3
output = produce_orbitdiagram(ro, plane, i, p_index, pvalues;
                              tfinal = tf, Ttr = 2000.0)

for (j, p) in enumerate(pvalues)
    axro.plot(fill(p, length(output[j])), output[j], lw = 0,
    marker = "o", ms = 0.2, color = "black", alpha = 0.1)
end
axro.set_xlim(pvalues[1], pvalues[end]);

# %%
ro = Systems.roessler(ones(3))
cc = 5.187

cs = range(5.187; step = -0.002, length = 3)

fig, axs = subplots(length(cs), 1; sharex = true)
T = 3000
dt = 0.05
for i in 1:length(cs)
	set_parameter!(ro, 3, cs[i])
	tr = trajectory(ro, T; Ttr = 200.0, dt)
	axs[i].plot(0:dt:T, tr[:, 1], lw = 0.5)
	axs[i].set_yticks([])
	axs[i].text(1.02, 0.5, "\$c=$(cs[i])\$"; bbox = bbox,
	transform = axs[i].transAxes, va = :center)
end
axs[1].set_xlim(0, T)
fig.tight_layout()
fig.subplots_adjust(hspace = 0.1, bottom = 0.1)
# fsave(plotsdir("intermittency"), fig)

max_stars_repo_name: moritz-g/ExercisesRepo

max_issues_repo_name: moritz-g/ExercisesRepo

max_forks_repo_name: moritz-g/ExercisesRepo


text:
# Syntax patterns
function prefix(io, call, argument1)
    print(io, call)
    print(io, ' ')
    print(io, argument1)
end

function parenthesize_prefix(io, call, argument1)
    print(io, call)
    print(io, ' ')
    print(io, '(')
    print(io, argument1)
    print(io, ')')
end

function postfix(io, call, argument1)
    print(io, argument1)
    print(io, ' ')
    print(io, call)
end

function tight_infix(io, call, argument1, argument2)
    print(io, argument1)
    print(io, call)
    print(io, argument2)
end

function infix(io, call, argument1, argument2, arguments...)
    print(io, argument1)
    print(io, ' ')
    print(io, call)
    print(io, ' ')
    infix(io, call, argument2, arguments...)
end

function infix(io, call, argument1)
    print(io, argument1)
end

function head_call_tail(io, call, argument1, arguments...)
    print(io, argument1)
    print(io, ' ')
    print(io, call)
    print(io, ' ')
    join(io, arguments, ", ")
end

function call_tail_head(io, call, argument1, arguments...)
    print(io, call)
    print(io, ' ')
    join(io, arguments, ", ")
    print(io, ' ')
    print(io, argument1)
end

function when_then_else(io, condition, action, arguments...)
    print(io, "WHEN ")
    print(io, condition)
    print(io, " THEN ")
    print(io, action)
    print(io, ' ')
    when_then_else(io, arguments...)
end

function when_then_else(io, action)
    print(io, "ELSE ")
    print(io, action)
end

function case(io, call, arguments...)
    print(io, call)
    print(io, ' ')
    when_then_else(io, arguments...)
    print(io, ' ')
    print(io, "END")
end

# Use the call to determine the syntax pattern
function show(io::IO, sql_expression::SQLExpression)
    call = sql_expression.call
    arguments = sql_expression.arguments

    if call === :.
        tight_infix(io, call, arguments...)
    elseif call === :(=)
        infix(io, call, arguments...)
    elseif call === :<
        infix(io, call, arguments...)
    elseif call === Symbol("<>")
        infix(io, call, arguments...)
    elseif call === Symbol("||")
        infix(io, call, arguments...)
    elseif call === :*
        infix(io, call, arguments...)
    elseif call === :/
        infix(io, call, arguments...)
    elseif call === :+
        infix(io, call, arguments...)
    elseif call === :-
        infix(io, call, arguments...)
    elseif call === :%
        infix(io, call, arguments...)
    elseif call === :||
        infix(io, call, arguments...)
    elseif call === :AND
        infix(io, call, arguments...)
    elseif call === :AS
        infix(io, call, arguments...)
    elseif call === :CASE
        case(io, call, arguments...)
    elseif call === :DESC
        postfix(io, call, arguments...)
    elseif call === :FROM
        parenthesize_prefix(io, call, arguments...)
    elseif call === Symbol("GROUP BY")
        head_call_tail(io, call, arguments...)
    elseif call === :IN
        infix(io, call, arguments...)
    elseif call === Symbol("INNER JOIN")
        infix(io, call, arguments...)
    elseif call === Symbol("IS NOT DISTINCT FROM")
        infix(io, call, arguments...)
    elseif call === Symbol("IS NULL")
        postfix(io, call, arguments...)
    elseif call === :LIKE
        infix(io, call, arguments...)
    elseif call === :LIMIT
        infix(io, call, arguments...)
    elseif call === :OFFSET
        infix(io, call, arguments...)
    elseif call === :ON
        infix(io, call, arguments...)
    elseif call === :OR
        infix(io, call, arguments...)
    elseif call === Symbol("ORDER BY")
        head_call_tail(io, call, arguments...)
    elseif call === :NOT
        prefix(io, call, arguments...)
    elseif call === :SELECT
        call_tail_head(io, call, arguments...)
    elseif call === Symbol("SELECT DISTINCT")
        call_tail_head(io, call, arguments...)
    elseif call === :WHERE
        infix(io, call, arguments...)
    else
        print(io, call)
        print(io, '(')
        join(io, arguments, ", ")
        print(io, ')')
    end
end

has_select(something) = false
has_select(expression::SQLExpression) =
    expression.call in (:SELECT, Symbol("SELECT DISTINCT"))

seek_select(something) = false
seek_select(expression::SQLExpression) =
    has_select(expression) ||
    any(seek_select(clause) for clause in expression.arguments)

# Default to * if no columns are specifcied
function finalize(sql_expression::SQLExpression)
    if !seek_select(sql_expression)
        SQLExpression(:SELECT, sql_expression, :*)
    else
        sql_expression
    end
end

function finalize(something)
    something
end

"""
    get_sql(it)

Use `get_sql` if you would like to see the SQL code generated by an SQLite
query.

```jldoctest example
julia> using QuerySQLite

julia> database = Database(joinpath(pathof(QuerySQLite) |> dirname |> dirname, "test", "Chinook_Sqlite.sqlite"));

julia> get_sql(database.Track)
SELECT * FROM (Track)
```
"""
function get_sql(it)
    finalize(translate(it.code))
end
export get_sql


text:
import Optim

function run_optimizer(
    optimizer::Optim.AbstractOptimizer,
    wrk,
    fg!,
    info_hook,
    check_convergence!
)

    tol_options = Optim.Options(
        # just so we can instantiate `optimizer_state` before `callback`
        x_tol=get(wrk.kwargs, :x_tol, 0.0),
        f_tol=get(wrk.kwargs, :f_tol, 0.0),
        g_tol=get(wrk.kwargs, :g_tol, 1e-8),
    )
    initial_x = wrk.pulsevals
    method = optimizer
    objective = Optim.promote_objtype(method, initial_x, :finite, true, Optim.only_fg!(fg!))
    optimizer_state = Optim.initial_state(method, tol_options, objective, initial_x)
    # Instantiation of `optimizer_state` calls `fg!` and sets the value of the
    # functional and gradient for the  `initial_x` in objective.F and
    # objective.DF, respectively. The `optimizer_state` is set correspondingly:
    @assert optimizer_state.x ==
            optimizer_state.x_previous ==
            objective.x_f ==
            objective.x_df
    @assert optimizer_state.g_previous == objective.DF
    # ... but `f_x_previous` does not match the initial `x_previous`:
    @assert isnan(optimizer_state.f_x_previous)

    # update the result object and check convergence
    function callback(optimization_state::Optim.OptimizationState)
        @assert optimization_state.value == objective.F
        #if optimization_state.iteration > 0
        #    @assert norm(
        #       optimizer_state.x .-
        #       (optimizer_state.x_previous .+ optimizer_state.alpha .* optimizer_state.s)
        #    ) < 1e-14
        #end
        wrk.gradient .= optimizer_state.g_previous
        wrk.searchdirection .= optimizer_state.s  # TODO: may depend on type of optimizer_state
        iter = wrk.result.iter_start + optimization_state.iteration
        update_result!(wrk, optimization_state, optimizer_state, iter)
        #update_hook!(...) # TODO
        info_tuple = info_hook(wrk, optimization_state, optimizer_state, wrk.result.iter)
        wrk.fg_count .= 0
        (info_tuple !== nothing) && push!(wrk.result.records, info_tuple)
        check_convergence!(wrk.result)
        return wrk.result.converged
    end

    options = Optim.Options(
        callback=callback,
        iterations=wrk.result.iter_stop - wrk.result.iter_start, # TODO
        x_tol=get(wrk.kwargs, :x_tol, 0.0),
        f_tol=get(wrk.kwargs, :f_tol, 0.0),
        g_tol=get(wrk.kwargs, :g_tol, 1e-8),
        show_trace=get(wrk.kwargs, :show_trace, false),
        extended_trace=get(wrk.kwargs, :extended_trace, false),
        store_trace=get(wrk.kwargs, :store_trace, false),
        show_every=get(wrk.kwargs, :show_every, 1),
        allow_f_increases=get(wrk.kwargs, :allow_f_increases, false),
    )

    res = Optim.optimize(objective, initial_x, method, options, optimizer_state)
    return res

end

function update_result!(
    wrk::GrapeWrk,
    optimization_state::Optim.OptimizationState,
    optimizer_state::Optim.AbstractOptimizerState,
    i::Int64
)
    # TODO: make this depend only on wrk. Should not be backend-dependent
    res = wrk.result
    for (k, Ψ) in enumerate(wrk.fw_states)
        copyto!(res.states[k], Ψ)
    end
    res.J_T_prev = res.J_T
    res.J_T = optimization_state.value
    (i > 0) && (res.iter = i)
    if i >= res.iter_stop
        res.converged = true
        res.message = "Reached maximum number of iterations"
        # Note: other convergence checks are done in user-supplied
        # check_convergence routine
    end
    prev_time = res.end_local_time
    res.end_local_time = now()
    res.secs = Dates.toms(res.end_local_time - prev_time) / 1000.0
end

function finalize_result!(wrk::GrapeWrk, optim_res::Optim.MultivariateOptimizationResults)
    L = length(wrk.controls)
    res = wrk.result
    if !optim_res.ls_success
        @error "optimization failed (linesearch)"
        res.message = "Failed linesearch"
    end
    if optim_res.stopped_by.f_increased
        @error "loss of monotonic convergence (try allow_f_increases=true)"
        res.message = "Loss of monotonic convergence"
    end
    if !res.converged
        @warn "Optimization failed to converge"
    end
    res.end_local_time = now()
    ϵ_opt = reshape(Optim.minimizer(optim_res), L, :)
    for l = 1:L
        res.optimized_controls[l] = discretize(ϵ_opt[l, :], res.tlist)
    end
    res.optim_res = optim_res
end

"""Print optimization progress as a table.

This functions serves as the default `info_hook` for an optimization with
GRAPE.
"""
function print_table(
    wrk,
    optimization_state::Optim.OptimizationState,
    optimizer_state::Optim.AbstractOptimizerState,
    iteration,
    args...
)
    # TODO: make this depend only on wrk. Should not be backend-dependent
    J_T = wrk.result.J_T
    ΔJ_T = J_T - wrk.result.J_T_prev
    secs = wrk.result.secs

    iter_stop = "$(get(wrk.kwargs, :iter_stop, 5000))"
    widths = [max(length("$iter_stop"), 6), 11, 11, 11, 8, 8]

    if iteration == 0
        header = ["iter.", "J_T", "|∇J_T|", "ΔJ_T", "FG(F)", "secs"]
        for (header, w) in zip(header, widths)
            print(lpad(header, w))
        end
        print("\n")
    end

    strs = (
        "$iteration",
        @sprintf("%.2e", J_T),
        @sprintf("%.2e", optimization_state.g_norm),
        (iteration > 0) ? @sprintf("%.2e", ΔJ_T) : "n/a",
        @sprintf("%d(%d)", wrk.fg_count[1], wrk.fg_count[2]),
        @sprintf("%.1f", secs),
    )
    for (str, w) in zip(strs, widths)
        print(lpad(str, w))
    end
    print("\n")
    flush(stdout)
end


text:
@testset "10.regular-expression-matching.jl" begin
    @test is_match("aa", "a") == false
    @test is_match("aa", "a*") == true
    @test is_match("ab", ".*") == true
    @test is_match("aab", "c*a*b") == true
    @test is_match("mississippi", "mis*is*p*.") == false
    @test is_match("aab", "c.") == false
    @test is_match("", "c") == false
end


text:
@testset "contrasts" begin

    d = DataFrame(y = rand(6),
                  x = [:b, :a, :c, :a, :a, :b])

    mf = ModelFrame(@formula(y ~ x), d)

    ## testing equality of ContrastsMatrix
    should_equal = [ContrastsMatrix(DummyCoding(), [:a, :b, :c]),
                    ContrastsMatrix(DummyCoding(base=:a), [:a, :b, :c]),
                    ContrastsMatrix(DummyCoding(base=:a, levels=[:a, :b, :c]), [:a, :b, :c]),
                    ContrastsMatrix(DummyCoding(levels=[:a, :b, :c]), [:a, :b, :c])]

    for c in should_equal
        @test mf.schema[term(:x)].contrasts == c
        @test hash(mf.schema[term(:x)].contrasts) == hash(c)
    end

    should_not_equal = [ContrastsMatrix(EffectsCoding(), [:a, :b, :c]),
                        ContrastsMatrix(DummyCoding(), [:b, :c]),
                        ContrastsMatrix(DummyCoding(), [:b, :c, :a]),
                        ContrastsMatrix(DummyCoding(), [:b, :c, :a]),
                        ContrastsMatrix(DummyCoding(), [:b, :c, :a]),
                        ContrastsMatrix(DummyCoding(base=:b), [:a, :b, :c]),
                        ContrastsMatrix(DummyCoding(base=:a, levels=[:b, :a, :c]), [:a, :b, :c])]

    for c in should_not_equal
        @test mf.schema[term(:x)].contrasts != c
    end

    # Dummy coded contrasts by default:
    @test ModelMatrix(mf).m == [1  1  0
                                1  0  0
                                1  0  1
                                1  0  0
                                1  0  0
                                1  1  0]
    @test coefnames(mf) == ["(Intercept)"; "x: b"; "x: c"]

    mmm = ModelMatrix(mf).m
    setcontrasts!(mf, x = DummyCoding())
    @test ModelMatrix(mf).m == mmm

    setcontrasts!(mf, x = EffectsCoding())
    @test ModelMatrix(mf).m == [1  1  0
                                1 -1 -1
                                1  0  1
                                1 -1 -1
                                1 -1 -1
                                1  1  0]
    @test coefnames(mf) == ["(Intercept)"; "x: b"; "x: c"]

    # change base level of contrast
    setcontrasts!(mf, x = EffectsCoding(base = :b))
    @test ModelMatrix(mf).m == [1 -1 -1
                                1  1  0
                                1  0  1
                                1  1  0
                                1  1  0
                                1 -1 -1]
    @test coefnames(mf) == ["(Intercept)"; "x: a"; "x: c"]

    # change levels of contrast
    setcontrasts!(mf, x = EffectsCoding(levels = [:c, :b, :a]))
    @test ModelMatrix(mf).m == [1  1  0
                                1  0  1
                                1 -1 -1
                                1  0  1
                                1  0  1
                                1  1  0]
    @test coefnames(mf) == ["(Intercept)"; "x: b"; "x: a"]

    # change levels and base level of contrast
    setcontrasts!(mf, x = EffectsCoding(levels = [:c, :b, :a], base = :a))
    @test ModelMatrix(mf).m == [1  0  1
                                1 -1 -1
                                1  1  0
                                1 -1 -1
                                1 -1 -1
                                1  0  1]
    @test coefnames(mf) == ["(Intercept)"; "x: c"; "x: b"]

    # respect order of levels

    data = DataFrame(x = levels!(categorical(['A', 'B', 'C', 'C', 'D']), ['C', 'B', 'A', 'D']))
    f = apply_schema(@formula(x ~ 1), schema(data))
    @test modelcols(f.lhs, data) == [0 1 0; 1 0 0; 0 0 0; 0 0 0; 0 0 1]

    # Helmert coded contrasts
    setcontrasts!(mf, x = HelmertCoding())
    @test ModelMatrix(mf).m == [1  1 -1
                                1 -1 -1
                                1  0  2
                                1 -1 -1
                                1 -1 -1
                                1  1 -1]
    @test coefnames(mf) == ["(Intercept)"; "x: b"; "x: c"]

    # Mismatching types of data and contrasts levels throws an error:
    @test_throws ArgumentError setcontrasts!(mf, x = EffectsCoding(levels = ["a", "b", "c"]))

    # Missing data is handled gracefully, dropping columns when a level is lost
    dm = deepcopy(d)
    allowmissing!(dm,:x)
    dm[3, :x] = missing
    mf_missing = ModelFrame(@formula(y ~ x), dm,
                            contrasts = Dict(:x => EffectsCoding()))
    @test ModelMatrix(mf_missing).m == [1  1
                                        1 -1
                                        1 -1
                                        1 -1
                                        1  1]
    @test coefnames(mf_missing) == ["(Intercept)"; "x: b"]

    # Sequential difference coding
    setcontrasts!(mf, x = SeqDiffCoding())
    seqdiff_contr = pinv([-1 1 0
                          0 -1 1]);
    @test ModelMatrix(mf).m ≈ hcat(ones(6), seqdiff_contr[[2, 1, 3, 1, 1, 2], :])

    # Things that are bad to do:
    # Applying contrasts that only have a subset of data levels:
    @test_throws ArgumentError setcontrasts!(mf, x = EffectsCoding(levels = [:a, :b]))
    # Applying contrasts that expect levels not found in data:
    @test_throws ArgumentError setcontrasts!(mf, x = EffectsCoding(levels = [:a, :b, :c, :d]))
    # Asking for base level that's not found in data
    @test_throws ArgumentError setcontrasts!(mf, x = EffectsCoding(base = :e))

    # Manually specified contrasts
    contrasts = [0  1
                 -1 -.5
                 1  -.5]
    setcontrasts!(mf, x = StatsModels.ContrastsCoding(contrasts))
    @test ModelMatrix(mf).m == [1 -1 -.5
                                1  0  1
                                1  1 -.5
                                1  0  1
                                1  0  1
                                1 -1 -.5]

    contrasts2 = [1 0
                  1 1
                  0 1]
    setcontrasts!(mf, x = StatsModels.ContrastsCoding(contrasts2))
    @test ModelMatrix(mf).m == [1  1  1
                                1  1  0
                                1  0  1
                                1  1  0
                                1  1  0
                                1  1  1]

    hypotheses2 = pinv(contrasts2)
    # need labels for hypothesis coding
    # TODO for a future release, make this an error
    # @test_throws ArgumentError HypothesisCoding(hypotheses2)

    hyp_labels = ["2a+b-c", "-a+b+2c"]
    setcontrasts!(mf, x = HypothesisCoding(hypotheses2, labels=hyp_labels))
    @test ModelMatrix(mf).m ≈ [1  1  1
                               1  1  0
                               1  0  1
                               1  1  0
                               1  1  0
                               1  1  1]

    # different results for non-orthogonal hypotheses/contrasts:
    hypotheses3 = [1 1 0
                   0 1 1]
    hyp_labels3 = ["a+b", "b+c"]
    hc3 = HypothesisCoding(hypotheses3, labels=hyp_labels3)
    setcontrasts!(mf, x = hc3)
    @test !(ModelMatrix(mf).m ≈ [1  1  1
                                 1  1  0
                                 1  0  1
                                 1  1  0
                                 1  1  0
                                 1  1  1])

    # accepts <:AbstractMatrix
    hypotheses4 = hcat([1, 1, 0], [0, 1, 1])'
    hc4 = HypothesisCoding(hypotheses4, labels=hyp_labels3)
    @test hc4.contrasts ≈ hc3.contrasts

    # specify labels via Vector{Pair}
    hc5 = HypothesisCoding(["a_and_b" => [1, 1, 0], "b_and_c" => [0, 1, 1]])
    @test hc5.contrasts[:, 1] ≈ hc3.contrasts[:,1]
    @test hc5.contrasts[:, 2] ≈ hc3.contrasts[:,2]
    @test hc5.labels == ["a_and_b", "b_and_c"]

    # throw argument error if number of levels mismatches
    @test_throws ArgumentError setcontrasts!(mf, x = StatsModels.ContrastsCoding(contrasts[1:2, :]))
    @test_throws ArgumentError setcontrasts!(mf, x = StatsModels.ContrastsCoding(hcat(contrasts, contrasts)))

    # contrasts types must be instantiated (should throw ArgumentError, currently
    # MethodError on apply_schema)
    @test_broken setcontrasts!(mf, x = DummyCoding)

    @testset "hypothesis coding" begin

        # to get the scaling right, divide by three (because three levels)
        effects_hyp = [-1 2 -1
                       -1 -1 2] ./ 3

        @test modelmatrix(setcontrasts!(mf,
                                        x = HypothesisCoding(effects_hyp,
                                                             labels=levels(d.x)[2:end]))) ≈
            modelmatrix(setcontrasts!(mf, x = EffectsCoding()))

        d2 = DataFrame(y = rand(100),
                       x = repeat([:a, :b, :c, :d], inner=25))

        sdiff_hyp = HypothesisCoding([-1 1 0 0
                                      0 -1 1 0
                                      0 0 -1 1],
                                     labels = ["b-a", "c-b", "d-c"])

        effects_hyp = HypothesisCoding([-1 3 -1 -1
                                        -1 -1 3 -1
                                        -1 -1 -1 3] ./ 4,
                                       labels = levels(d2.x)[2:end])

        f = apply_schema(@formula(y ~ 1 + x), schema(d2))

        f_sdiff = apply_schema(f, schema(d2, Dict(:x => sdiff_hyp)))
        f_effects = apply_schema(f, schema(d2, Dict(:x => effects_hyp)))

        y_means = by(d2, :x, :y => mean).y_mean

        y, X_sdiff = modelcols(f_sdiff, d2)
        @test X_sdiff \ y ≈ [mean(y_means); diff(y_means)]

        y, X_effects = modelcols(f_effects, d2)
        @test X_effects \ y ≈ [mean(y_means); y_means[2:end] .- mean(y_means)]

        @test X_effects ≈ modelcols(apply_schema(f.rhs,
                                                 schema(d2, Dict(:x=>EffectsCoding()))),
                                    d2)
    end

    @testset "hypothesis_matrix" begin
        using StatsModels: contrasts_matrix, hypothesis_matrix, needs_intercept

        cmat = contrasts_matrix(DummyCoding(), 1, 4)
        @test needs_intercept(cmat) == true
        cmat1 = hypothesis_matrix(cmat)
        @test cmat1 ≈
            [ 1 0 0 0
             -1 1 0 0
             -1 0 1 0
             -1 0 0 1]
        @test eltype(cmat1) <: Integer
        @test eltype(hypothesis_matrix(cmat1, tolerance=0)) == eltype(cmat)
        @test eltype(hypothesis_matrix(cmat1, tolerance=0.0)) == eltype(cmat)

        # incorrect interpretation without considering intercept:
        @test hypothesis_matrix(cmat, intercept=false) ≈
            [0 1 0 0
             0 0 1 0
             0 0 0 1]

        cmat2 = contrasts_matrix(HelmertCoding(), 1, 4)
        @test needs_intercept(cmat2) == false
        hmat2 = hypothesis_matrix(cmat2)
        @test hmat2 ≈
            [-1/2   1/2   0    0
             -1/6  -1/6   1/3  0
             -1/12 -1/12 -1/12 1/4]

        @test eltype(hmat2) <: Rational

        @test hypothesis_matrix(cmat2, intercept=true) ≈
            vcat([1/4 1/4 1/4 1/4], hmat2)

        cmat3 = [-1. -1
                  1   0
                  0   1]
        @test needs_intercept(cmat3) == false
        cmat3p = copy(cmat3)
        cmat3p[1] += 1e-3
        @test needs_intercept(cmat3p) == true
    end

    @testset "levels and baselevel" begin
        using DataAPI: levels
        using StatsModels: baselevel, FullDummyCoding, ContrastsCoding

        levs = [:a, :b, :c, :d]
        base = [:a]
        for C in [DummyCoding, EffectsCoding, SeqDiffCoding, HelmertCoding]
            c = C()
            @test levels(c) == nothing
            @test baselevel(c) == nothing

            c = C(levels=levs)
            @test levels(c) == levs
            @test baselevel(c) == nothing

            c = C(base=base)
            @test levels(c) == nothing
            @test baselevel(c) == base

            c = C(levels=levs, base=base)
            @test levels(c) == levs
            @test baselevel(c) == base
        end

        c = FullDummyCoding()
        @test baselevel(c) == nothing
        @test levels(c) == nothing

        @test_throws MethodError FullDummyCoding(levels=levs)
        @test_throws MethodError FullDummyCoding(base=base)

        c = HypothesisCoding(rand(3,4))
        @test baselevel(c) == levels(c) == nothing
        c = HypothesisCoding(rand(3,4), levels=levs)
        @test baselevel(c) == nothing
        @test levels(c) == levs
        # no notion of base level for HypothesisCoding
        @test_throws MethodError HypothesisCoding(rand(3,4), base=base)

        c = ContrastsCoding(rand(4,3))
        @test baselevel(c) == levels(c) == nothing
        c = ContrastsCoding(rand(4,3), levels=levs)
        @test baselevel(c) == nothing
        @test levels(c) == levs
        # no notion of base level for ContrastsCoding
        @test_throws MethodError ContrastsCoding(rand(4,3), base=base)

    end
end


text:
s = open("input.txt") do file
    readlines(file)
end

function getCount(password, letter::Char)
    count = 0;
    for i = 1:length(password)
        if password[i] == letter
            count += 1;
        end
    end

    return count;
end

function getGoodCount(s)
    goodCount = 0;
    for line in s
        m = match(r"(\d+)-(\d+) (.): (.*)", line);

        minCount = parse(Int32, m.captures[1]);
        maxCount = parse(Int32, m.captures[2]);
        letter = Char(m.captures[3][1]);

        password = m.captures[4];

        count = getCount(password, letter);

        if minCount <= count <= maxCount
            goodCount += 1;
        end
    end

    return goodCount;
end

println(getGoodCount(s));