text:
buildexpr(s::Symbol) = WSymbol(s)

function buildexpr(s::AbstractArray)
    colons = Iterators.repeated(:, ndims(s) - 1)
    nested = [s[i, colons...] for i in 1:size(s)[1]]
    WSymbol("List")(buildexpr.(nested)...)
end
buildexpr(s::Tuple) = WSymbol("List")(buildexpr.(s)...)
buildexpr(s::Dict) = WSymbol("Association")(buildexpr.(collect(s))...)
buildexpr(s::Rational) = WSymbol("Rational")(buildexpr(s.num), buildexpr(s.den))
buildexpr(s::Pair) = WSymbol("Rule")(buildexpr(s.first), buildexpr(s.second))
buildexpr(e::Any) = e

const headstype = Dict(
    "Rational" => Rational,
    "Rule" => Pair,
    "List" => Array,
    "Association" => Dict,
)

getexpr(e::Any) = e
getexpr(e::WSymbol) = Symbol(e.name)
getexpr(e::WExpr) = getexpr(Any, e)
function getexpr(::Type{Any}, e::WExpr)
    name = e.head.name
    if haskey(headstype, name)
        getexpr(headstype[name], e)
    else
        e
    end
end

getexpr(::Type{Rational}, e) = Rational(getexpr.(e.args)...)
getexpr(::Type{Pair}, e) = Pair(getexpr(e.args[1]), getexpr(e.args[2]))
getexpr(::Type{Dict}, e) = Dict(getexpr.(e.args))
function getexpr(::Type{Array}, e)
    s, ty = getsize(e)
    arr = Array{ty}(undef, prod(s))
    fillarray(arr, e, s, 0, 1)
    reshape(arr, s)
end

# Auxiliary functions for getexpr(Array, .)
getsize(a::Any) = ((), typeof(a))
function getsize(a::WExpr) where T
    dims = getsize.(a.args)
    function reducer(x, y)
        if x[1] == y[1]
            return (x[1], typejoin(x[2], y[2]))
        end
        ((zip(x[1], y[1])
         |> xs -> takewhile(p -> p[1] == p[2], xs)
         |> xs -> map(p -> p[1], xs)), Any)
    end
    subdim, ty = reduce(reducer, dims)
    ((length(a.args), subdim...), ty)
end
function fillarray(arr, e, size, acc, coef)
    for i in 1:size[1]
        if length(size) == 1
            arr[acc + coef * (i - 1) + 1] = getexpr(e.args[i])
        else
            fillarray(arr, e.args[i], size[2:end], acc + coef * (i - 1), coef * size[1])
        end
    end
end


text:
module UpdateMethods

export Heun, heun_update

struct Heun
    temppara::Real
    step::Real
end

Heun(;step, temperature) = Heun(√(2 *temperature*step), step)

function heun_update(forcefunc::Function, temppara::Real, step::Real, y::Real, t::Real)::Tuple{Real, Real}
    w = randn() * temppara
    c1 = forcefunc(y,t)
    c2 = forcefunc(y+c1*step+w,t+step)
    yn = y+step/2.0*(c1+c2)+w
    return t+step, yn
end

heun_update(f::Function, h::Heun, y::Real, t::Real)::Tuple{Real, Real} = heun_update(f, h.temppara, h.step, y, t)

end

ext: jl


text:
using Dispersal, Test

init = [1.0 4.0 7.0;
        2.0 5.0 8.0;
        3.0 6.0 9.0]

output = ArrayOutput(init; tspan=1:2)
rule = AlleeExtinction(minfounders = 8.0)
sim!(output, rule)

@test output[1] == [1.0 4.0 7.0;
                    2.0 5.0 8.0;
                    3.0 6.0 9.0]
@test output[2] == [0.0 0.0 0.0;
                    0.0 0.0 8.0;
                    0.0 0.0 9.0]


text:
using FlatBuffers
using Test
import Parameters

# test default fields
@with_kw mutable struct UltimateAnswer
    answer::Int32 = 42
    question::String
    highwaysbuilt::Int32 = 7
end

x = UltimateAnswer(;question="How many roads must a man walk down?")
@test x.answer == 42
@test FlatBuffers.default(UltimateAnswer, Int32, :answer) == 42
@test FlatBuffers.default(UltimateAnswer, Int32, :highwaysbuilt) == 7
b = FlatBuffers.Builder(UltimateAnswer)
FlatBuffers.build!(b, x)
xbytes = FlatBuffers.bytes(b)
y = FlatBuffers.read(UltimateAnswer, xbytes)

@test y.answer == x.answer
@test y.question == x.question
@test y.highwaysbuilt == x.highwaysbuilt
@test x.highwaysbuilt == 7

y = Parameters.reconstruct(x, highwaysbuilt = 0)
b = FlatBuffers.Builder(UltimateAnswer)
FlatBuffers.build!(b, y)
ybytes = FlatBuffers.bytes(b)

# check that we save bytes with default integer values
@test length(ybytes) > length(xbytes)

@test y.answer == x.answer
@test y.question == x.question
@test y.highwaysbuilt == 0

y = FlatBuffers.read(UltimateAnswer, ybytes)
@test y.highwaysbuilt == 0


text:
## --- File IO primitives

    # Open a file
    """
    ```julia
    fopen(name, mode)
    ```
    Libc `fopen` function, accessed by direct `llvmcall`.

    Returns a file pointer to a file at location specified by `name` opened for
    reading, writing, or both as specified by `mode`. Valid modes include:

      c"r": Read, from an existing file.

      c"w": Write. If the file exists, it will be overwritten.

      c"a": Append, to the end of an existing file.

    as well as `"r+"`, `c"w+"`, and `"a+"`, which enable both reading and writing.

    See also: `fclose`, `fseek`

    ## Examples
    ```julia
    julia> fp = fopen(c"testfile.txt", c"w")
    Ptr{StaticTools.FILE} @0x00007fffc92bd0b0

    julia> printf(fp, c"Here is a string")
    16

    julia> fclose(fp)
    0

    shell> cat testfile.txt
    Here is a string
    ```
    """
    fopen(name::AbstractMallocdMemory, mode::AbstractMallocdMemory) = fopen(pointer(name), pointer(mode))
    fopen(name, mode) = GC.@preserve name mode fopen(pointer(name), pointer(mode))
    function fopen(name::Ptr{UInt8}, mode::Ptr{UInt8})
        Base.llvmcall(("""
        ; External declaration of the fopen function
        declare i8* @fopen(i8*, i8*)

        define i64 @main(i64 %jlname, i64 %jlmode) #0 {
        entry:
          %name = inttoptr i64 %jlname to i8*
          %mode = inttoptr i64 %jlmode to i8*
          %fp = call i8* (i8*, i8*) @fopen(i8* %name, i8* %mode)
          %jlfp = ptrtoint i8* %fp to i64
          ret i64 %jlfp
        }

        attributes #0 = { alwaysinline nounwind ssp uwtable }
        """, "main"), Ptr{FILE}, Tuple{Ptr{UInt8}, Ptr{UInt8}}, name, mode)
    end

    # Close a file
    """
    ```julia
    fclose(fp::Ptr{FILE})
    ```
    Libc `fclose` function, accessed by direct `llvmcall`.

    Closes a file that has been previously opened by `fopen`, given a file pointer.

    See also: `fopen`, `fseek`

    ## Examples
    ```julia
    julia> fp = fopen(c"testfile.txt", c"w")
    Ptr{StaticTools.FILE} @0x00007fffc92bd0b0

    julia> printf(fp, c"Here is a string")
    16

    julia> fclose(fp)
    0

    shell> cat testfile.txt
    Here is a string
    ```
    """
    function fclose(fp::Ptr{FILE})
        Base.llvmcall(("""
        ; External declaration of the fclose function
        declare i32 @fclose(i8*)

        define i32 @main(i64 %jlfp) #0 {
        entry:
          %fp = inttoptr i64 %jlfp to i8*
          %status = call i32 (i8*) @fclose(i8* %fp)
          ret i32 %status
        }

        attributes #0 = { alwaysinline nounwind ssp uwtable }
        """, "main"), Int32, Tuple{Ptr{FILE}}, fp)
    end

    # Seek in a file
    frewind(fp::Ptr{FILE}) = fseek(fp, 0, SEEK_SET)

    """
    ```julia
    fseek(fp::Ptr{FILE}, offset::Int64, whence::Int32=SEEK_CUR)
    ```
    Libc `fseek` function, accessed by direct `llvmcall`.

    Move position within a file given a file pointer `fp` obtained from `fopen`.
    The new position will be `offset` bytes (or characters, in the event that all
    characters are non-unicode ASCII characters encoded as `UInt8`s) away from the
    position specified by `whence`.

    The position reference `whence` can take on values of either:

      SEEK_SET = Int32(0)           File start

      SEEK_CUR = Int32(1)           Current position

      SEEK_END = Int32(2)           File end

    where `SEEK_CUR` is the default value.

    See also: `fopen`, `fclose`

    ## Examples
    ```julia
    julia> fp = fopen(c"testfile.txt", c"w+")
    Ptr{StaticTools.FILE} @0x00007fffc92bd148

    julia> printf(fp, c"Here is a string!")
    17

    julia> fseek(fp, -2)
    0

    julia> Char(getc(fp))
    'g': ASCII/Unicode U+0067 (category Ll: Letter, lowercase)

    julia> fclose(fp)
    0
    ```
    """
    function fseek(fp::Ptr{FILE}, offset::Int64, whence::Int32=SEEK_CUR)
        Base.llvmcall(("""
        ; External declaration of the fseek function
        declare i32 @fseek(i8*, i64, i32)

        define i32 @main(i64 %jlfp, i64 %offset, i32 %whence) #0 {
        entry:
          %fp = inttoptr i64 %jlfp to i8*
          %status = call i32 @fseek(i8* %fp, i64 %offset, i32 %whence)
          ret i32 %status
        }

        attributes #0 = { alwaysinline nounwind ssp uwtable }
        """, "main"), Int32, Tuple{Ptr{FILE}, Int64, Int32}, fp, offset, whence)
    end
    const SEEK_SET = Int32(0)
    const SEEK_CUR = Int32(1)
    const SEEK_END = Int32(2)

## -- stdio pointers

    # Get pointer to stdout
    """
    ```julia
    stdoutp()
    ```
    Zero-argument function which returns a raw pointer to the current standard
    output filestream, `stdout`.

    ## Examples
    ```julia
    julia> stdoutp()
    Ptr{StaticTools.FILE} @0x00007fffc92b91a8

    julia> printf(stdoutp(), c"Hi there!\n")
    Hi there!
    10
    ```
    """
@static if Sys.isapple()
    function stdoutp()
        Base.llvmcall(("""
        @__stdoutp = external global i8*

        define i64 @main() #0 {
        entry:
          %ptr = load i8*, i8** @__stdoutp, align 8
          %jlfp = ptrtoint i8* %ptr to i64
          ret i64 %jlfp
        }

        attributes #0 = { alwaysinline nounwind ssp uwtable }
        """, "main"), Ptr{FILE}, Tuple{})
    end
else
    function stdoutp()
        Base.llvmcall(("""
        @stdout = external global i8*

        define i64 @main() #0 {
        entry:
          %ptr = load i8*, i8** @stdout, align 8
          %jlfp = ptrtoint i8* %ptr to i64
          ret i64 %jlfp
        }

        attributes #0 = { alwaysinline nounwind ssp uwtable }
        """, "main"), Ptr{FILE}, Tuple{})
    end
end

    # Get pointer to stderr
    """
    ```julia
    stderrp()
    ```
    Zero-argument function which returns a raw pointer to the current standard
    error filestream, `stderr`.

    ## Examples
    ```julia
    julia> stderrp()
    Ptr{StaticTools.FILE} @0x00007fffc92b9240

    julia> printf(stderrp(), c"Hi there!\n")
    Hi there!
    10
    ```
    """
@static if Sys.isapple()
    function stderrp()
        Base.llvmcall(("""
        @__stderrp = external global i8*

        define i64 @main() #0 {
        entry:
          %ptr = load i8*, i8** @__stderrp, align 8
          %jlfp = ptrtoint i8* %ptr to i64
          ret i64 %jlfp
        }

        attributes #0 = { alwaysinline nounwind ssp uwtable }
        """, "main"), Ptr{FILE}, Tuple{})
    end
else
    function stderrp()
        Base.llvmcall(("""
        @stderr = external global i8*

        define i64 @main() #0 {
        entry:
          %ptr = load i8*, i8** @stderr, align 8
          %jlfp = ptrtoint i8* %ptr to i64
          ret i64 %jlfp
        }

        attributes #0 = { alwaysinline nounwind ssp uwtable }
        """, "main"), Ptr{FILE}, Tuple{})
    end
end

    # Get pointer to stdin
    """
    ```julia
    stdinp()
    ```
    Zero-argument function which returns a raw pointer to the current standard
    input filestream, `stdin`.

    ## Examples
    ```julia
    julia> stdinp()
    Ptr{StaticTools.FILE} @0x00007fffc92b9110
    ```
    """
@static if Sys.isapple()
    function stdinp()
        Base.llvmcall(("""
        @__stdinp = external global i8*

        define i64 @main() #0 {
        entry:
          %ptr = load i8*, i8** @__stdinp, align 8
          %jlfp = ptrtoint i8* %ptr to i64
          ret i64 %jlfp
        }

        attributes #0 = { alwaysinline nounwind ssp uwtable }
        """, "main"), Ptr{FILE}, Tuple{})
    end
else
    function stdinp()
        Base.llvmcall(("""
        @stdin = external global i8*

        define i64 @main() #0 {
        entry:
          %ptr = load i8*, i8** @stdin, align 8
          %jlfp = ptrtoint i8* %ptr to i64
          ret i64 %jlfp
        }

        attributes #0 = { alwaysinline nounwind ssp uwtable }
        """, "main"), Ptr{FILE}, Tuple{})
    end
end

## --- The basics of basics: putchar/fputc

    """
    ```julia
    putchar([fp::Ptr{FILE}], c::Union{Char,UInt8})
    ```
    Libc `putchar` / `fputc` function, accessed by direct `llvmcall`.

    Prints a single character `c` (either a `Char` or a raw `UInt8`) to a file
    pointer `fp`, defaulting to the current standard output `stdout` if not
    specified.

    Returns `0` on success.

    ## Examples
    ```julia
    julia> putchar('C')
    0

    julia> putchar(0x63)
    0

    julia> putchar('\n') # Newline, flushes stdout
    Cc
    0
    ```
    """
    putchar(c::Char) = putchar(UInt8(c))
    function putchar(c::UInt8)
        Base.llvmcall(("""
        ; External declaration of the putchar function
        declare i32 @putchar(i8 nocapture) nounwind

        define i32 @main(i8 %c) #0 {
        entry:
          %status = call i32 (i8) @putchar(i8 %c)
          ret i32 0
        }

        attributes #0 = { alwaysinline nounwind ssp uwtable }
        """, "main"), Int32, Tuple{UInt8}, c)
    end
    putchar(fp::Ptr{FILE}, c::Char) = putchar(fp, UInt8(c))
    function putchar(fp::Ptr{FILE}, c::UInt8)
        Base.llvmcall(("""
        ; External declaration of the fputc function
        declare i32 @fputc(i8, i8*) nounwind

        define i32 @main(i64 %jlfp, i8 %c) #0 {
        entry:
          %fp = inttoptr i64 %jlfp to i8*
          %status = call i32 (i8, i8*) @fputc(i8 %c, i8* %fp)
          ret i32 0
        }

        attributes #0 = { alwaysinline nounwind ssp uwtable }
        """, "main"), Int32, Tuple{Ptr{FILE}, UInt8}, fp, c)
    end

    """
    ```julia
    newline([fp::Ptr{FILE}])
    ```
    Prints a single newline (`\n`) to a file pointer `fp`, defaulting  to `stdout`
    if not specified.

    Returns `0` on success.

    ## Examples
    ```julia
    julia> putchar('C')
    0

    julia> newline() # flushes stdout
    C
    0
    ```
    """
    newline() = putchar(0x0a)
    newline(fp::Ptr{FILE}) = putchar(fp, 0x0a)

## --- getchar / getc

    """
    ```julia
    getchar()
    ```
    Libc `getchar` function, accessed by direct `llvmcall`.

    Reads a single character from standard input `stdin`, returning as `UInt8`.

    ## Examples
    ```julia
    julia> getchar()
    c
    0x63
    ```
    """
    function getchar()
        Base.llvmcall(("""
        ; External declaration of the getchar function
        declare i32 @getchar()

        define i8 @main() #0 {
        entry:
          %result = call i32 @getchar()
          %c = trunc i32 %result to i8
          ret i8 %c
        }

        attributes #0 = { alwaysinline nounwind ssp uwtable }
        """, "main"), UInt8, Tuple{})
    end

    """
    ```julia
    getc(fp::Ptr{FILE})
    ```
    Libc `getc` function, accessed by direct `llvmcall`.

    Reads a single character from file pointer `fp`, returning as `Int32`
    (`-1` on EOF).

    ## Examples
    ```julia
    julia> getc(stdinp())
    c
    99
    ```
    """
    function getc(fp::Ptr{FILE})
        Base.llvmcall(("""
        ; External declaration of the fgetc function
        declare i32 @fgetc(i8*)

        define i32 @main(i64 %jlfp) #0 {
        entry:
          %fp = inttoptr i64 %jlfp to i8*
          %c = call i32 (i8*) @fgetc(i8* %fp)
          ret i32 %c
        }

        attributes #0 = { alwaysinline nounwind ssp uwtable }
        """, "main"), Int32, Tuple{Ptr{FILE}}, fp)
    end
    const EOF = Int32(-1)

## --- The old reliable: puts/fputs

    """
    ```julia
    puts([fp::Ptr{FILE}], s)
    ```
    Libc `puts`/`fputs` function, accessed by direct `llvmcall`.

    Prints a string `s` (specified either as a raw `Ptr{UInt8}` to a valid
    null-terminated string in memory or elseor else a string type such as
    `StaticString` or `MallocString` for which a valid pointer can be obtained)
    followed by a newline (`\n`) to a filestream specified by the file pointer
    `fp`, defaulting to the current standard output `stdout` if not specified.

    Returns `0` on success.

    ## Examples
    ```julia
    julia> puts(c"Hello there!")
    Hello there!
    0
    ```
    """
    puts(s::AbstractMallocdMemory) = puts(pointer(s))
    puts(s) = GC.@preserve s puts(pointer(s))
    function puts(s::Ptr{UInt8})
        Base.llvmcall(("""
        ; External declaration of the puts function
        declare i32 @puts(i8* nocapture) nounwind

        define i32 @main(i64 %jls) #0 {
        entry:
          %str = inttoptr i64 %jls to i8*
          %status = call i32 (i8*) @puts(i8* %str)
          ret i32 0
        }

        attributes #0 = { alwaysinline nounwind ssp uwtable }
        """, "main"), Int32, Tuple{Ptr{UInt8}}, s)
    end

    puts(fp::Ptr{FILE}, s::AbstractMallocdMemory) = puts(fp, pointer(s))
    puts(fp::Ptr{FILE}, s) = GC.@preserve s puts(fp, pointer(s))
    function puts(fp::Ptr{FILE}, s::Ptr{UInt8})
        Base.llvmcall(("""
        ; External declaration of the puts function
        declare i32 @fputs(i8*, i8*) nounwind

        define i32 @main(i64 %jlfp, i64 %jls) #0 {
        entry:
          %fp = inttoptr i64 %jlfp to i8*
          %str = inttoptr i64 %jls to i8*
          %status = call i32 (i8*, i8*) @fputs(i8* %str, i8* %fp)
          ret i32 0
        }

        attributes #0 = { alwaysinline nounwind ssp uwtable }
        """, "main"), Int32, Tuple{Ptr{FILE}, Ptr{UInt8}}, fp, s)
        newline(fp) # puts appends `\n`, but fputs doesn't (!)
    end
    puts(s::StringView) = (printf(s); newline())
    puts(fp::Ptr{FILE}, s::StringView) = (printf(fp, s); newline(fp))

## --- gets/fgets

    """
    ```julia
    gets!(s::MallocString, fp::Ptr{FILE}, n::Integer=length(s))
    ```
    Libc `fgets` function, accessed by direct `llvmcall`.

    Read up to `n` characters from the filestream specified by file pointer `fp`
    to the MallocString `s`. Stops when a newline is encountered, end-of-file
    is reached, or `n-1` characters have been read (whichever comes first).

    ## Examples
    ```julia
    julia> s = MallocString(undef, 100)
    m""

    julia> gets!(s, stdinp(), 3)
    Ptr{UInt8} @0x00007fb15afce550

    julia> s
    m"\n"
    ```
    """
    function gets!(s::MallocString, fp::Ptr{FILE}, n::Integer=length(s))
        Base.llvmcall(("""
        ; External declaration of the gets function
        declare i8* @fgets(i8*, i32, i8*)

        define i64 @main(i64 %jls, i64 %jlfp, i32 %n) #0 {
        entry:
          %str = inttoptr i64 %jls to i8*
          %fp = inttoptr i64 %jlfp to i8*
          %stp = call i8* (i8*, i32, i8*) @fgets(i8* %str, i32 %n, i8* %fp)
          %status = ptrtoint i8* %stp to i64
          ret i64 %status
        }

        attributes #0 = { alwaysinline nounwind ssp uwtable }
        """, "main"), Ptr{UInt8}, Tuple{Ptr{UInt8}, Ptr{FILE}, Int32}, pointer(s), fp, n % Int32)
    end
    const readline! = gets!

## --- printf/fprintf, just a string

    """
    ```julia
    printf([fp::Ptr{FILE}], [fmt], s)
    ```
    Libc `printf` function, accessed by direct `llvmcall`.

    Prints a string `s` (specified either as a raw `Ptr{UInt8}` to a valid
    null-terminated string in memory or else a string type such as `StaticString`
    or `MallocString` for which a valid pointer can be obtained) to a filestream
    specified by the file pointer `fp`, defaulting to the current standard output
    `stdout` if not specified.

    Optionally, a C-style format specifier string `fmt` may be provided as well.

    Returns the number of characters printed on success.

    ## Examples
    ```julia
    julia> printf(c"Hello there!\n")
    Hello there!
    13
    ```
    """
    printf(s::MallocString) = printf(pointer(s))
    printf(s) = GC.@preserve s printf(pointer(s))
    printf(fp::Ptr{FILE}, s::MallocString) = printf(fp, pointer(s))
    printf(fp::Ptr{FILE}, s) = GC.@preserve s printf(fp, pointer(s))
    function printf(s::Ptr{UInt8})
        Base.llvmcall(("""
        ; External declaration of the printf function
        declare i32 @printf(i8* noalias nocapture, ...)

        define i32 @main(i64 %jls) #0 {
        entry:
          %str = inttoptr i64 %jls to i8*
          %status = call i32 (i8*, ...) @printf(i8* %str)
          ret i32 %status
        }

        attributes #0 = { alwaysinline nounwind ssp uwtable }
        """, "main"), Int32, Tuple{Ptr{UInt8}}, s)
    end
    function printf(fp::Ptr{FILE}, s::Ptr{UInt8})
        Base.llvmcall(("""
        ; External declaration of the fprintf function
        declare i32 @fprintf(i8*, i8*)

        define i32 @main(i64 %jlfp, i64 %jls) #0 {
        entry:
          %fp = inttoptr i64 %jlfp to i8*
          %str = inttoptr i64 %jls to i8*
          %status = call i32 (i8*, i8*) @fprintf(i8* %fp, i8* %str)
          ret i32 %status
        }

        attributes #0 = { alwaysinline nounwind ssp uwtable }
        """, "main"), Int32, Tuple{Ptr{FILE}, Ptr{UInt8}}, fp, s)
    end
    function printf(s::StringView)
        for i ∈ eachindex(s)
            putchar(s[i])
        end
        return length(s) % Int32
    end
    function printf(fp::Ptr{FILE}, s::StringView)
        for i ∈ eachindex(s)
            putchar(fp, s[i])
        end
        return length(s) % Int32
    end

## --- printf/fprintf, with a format string, just like in C

    printf(fmt::MallocString, s::MallocString) = printf(pointer(fmt), pointer(s))
    printf(fmt, s) = GC.@preserve fmt s printf(pointer(fmt), pointer(s))
    printf(fp::Ptr{FILE}, fmt::MallocString, s::MallocString) = printf(fp::Ptr{FILE}, pointer(fmt), pointer(s))
    printf(fp::Ptr{FILE}, fmt, s) = GC.@preserve fmt s printf(fp::Ptr{FILE}, pointer(fmt), pointer(s))
    function printf(fmt::Ptr{UInt8}, s::Ptr{UInt8})
        Base.llvmcall(("""
        ; External declaration of the printf function
        declare i32 @printf(i8* noalias nocapture, ...)

        define i32 @main(i64 %jlf, i64 %jls) #0 {
        entry:
          %fmt = inttoptr i64 %jlf to i8*
          %str = inttoptr i64 %jls to i8*
          %status = call i32 (i8*, ...) @printf(i8* %fmt, i8* %str)
          ret i32 0
        }

        attributes #0 = { alwaysinline nounwind ssp uwtable }
        """, "main"), Int32, Tuple{Ptr{UInt8}, Ptr{UInt8}}, fmt, s)
    end
    function printf(fp::Ptr{FILE}, fmt::Ptr{UInt8}, s::Ptr{UInt8})
        Base.llvmcall(("""
        ; External declaration of the fprintf function
        declare i32 @fprintf(i8*, ...)

        define i32 @main(i64 %jlfp, i64 %jlf, i64 %jls) #0 {
        entry:
          %fp = inttoptr i64 %jlfp to i8*
          %fmt = inttoptr i64 %jlf to i8*
          %str = inttoptr i64 %jls to i8*
          %status = call i32 (i8*, ...) @fprintf(i8* %fp, i8* %fmt, i8* %str)
          ret i32 %status
        }

        attributes #0 = { alwaysinline nounwind ssp uwtable }
        """, "main"), Int32, Tuple{Ptr{FILE}, Ptr{UInt8}, Ptr{UInt8}}, fp, fmt, s)
    end

    printf(fmt::StaticString, n::Union{Number, Ptr}) = GC.@preserve fmt printf(pointer(fmt), n)
    printf(fmt::MallocString, n::Union{Number, Ptr}) = printf(pointer(fmt), n)
    printf(fp::Ptr{FILE}, fmt::StaticString, n::Union{Number, Ptr}) = GC.@preserve fmt printf(fp::Ptr{FILE}, pointer(fmt), n)
    printf(fp::Ptr{FILE}, fmt::MallocString, n::Union{Number, Ptr}) = printf(fp::Ptr{FILE}, pointer(fmt), n)

    # Floating point numbers
    function printf(fmt::Ptr{UInt8}, n::Float64)
        Base.llvmcall(("""
        ; External declaration of the printf function
        declare i32 @printf(i8* noalias nocapture, ...)

        define i32 @main(i64 %jlf, double %d) #0 {
        entry:
          %fmt = inttoptr i64 %jlf to i8*
          %status = call i32 (i8*, ...) @printf(i8* %fmt, double %d)
          ret i32 0
        }

        attributes #0 = { alwaysinline nounwind ssp uwtable }
        """, "main"), Int32, Tuple{Ptr{UInt8}, Float64}, fmt, n)
    end
    function printf(fp::Ptr{FILE}, fmt::Ptr{UInt8}, n::Float64)
        Base.llvmcall(("""
        ; External declaration of the printf function
        declare i32 @fprintf(i8*, ...)

        define i32 @main(i64 %jlfp, i64 %jlf, double %n) #0 {
        entry:
          %fp = inttoptr i64 %jlfp to i8*
          %fmt = inttoptr i64 %jlf to i8*
          %status = call i32 (i8*, ...) @fprintf(i8* %fp, i8* %fmt, double %n)
          ret i32 %status
        }

        attributes #0 = { alwaysinline nounwind ssp uwtable }
        """, "main"), Int32, Tuple{Ptr{FILE}, Ptr{UInt8}, Float64}, fp, fmt, n)
    end

    # Just convert all other Floats to double
    printf(fmt::Ptr{UInt8}, n::AbstractFloat) = printf(fmt::Ptr{UInt8}, Float64(n))
    printf(fp::Ptr{FILE}, fmt::Ptr{UInt8}, n::AbstractFloat) = printf(fp::Ptr{FILE}, fmt::Ptr{UInt8}, Float64(n))

    # Integers
    function printf(fmt::Ptr{UInt8}, n::T) where T <: Union{Int64, UInt64, Ptr}
        Base.llvmcall(("""
        ; External declaration of the printf function
        declare i32 @printf(i8* noalias nocapture, ...)

        define i32 @main(i64 %jlf, i64 %n) #0 {
        entry:
          %fmt = inttoptr i64 %jlf to i8*
          %status = call i32 (i8*, ...) @printf(i8* %fmt, i64 %n)
          ret i32 0
        }

        attributes #0 = { alwaysinline nounwind ssp uwtable }
        """, "main"), Int32, Tuple{Ptr{UInt8}, T}, fmt, n)
    end
    function printf(fp::Ptr{FILE}, fmt::Ptr{UInt8}, n::T) where T <: Union{Int64, UInt64, Ptr}
        Base.llvmcall(("""
        ; External declaration of the printf function
        declare i32 @fprintf(i8*, ...)

        define i32 @main(i64 %jlfp, i64 %jlf, i64 %n) #0 {
        entry:
          %fp = inttoptr i64 %jlfp to i8*
          %fmt = inttoptr i64 %jlf to i8*
          %status = call i32 (i8*, ...) @fprintf(i8* %fp, i8* %fmt, i64 %n)
          ret i32 %status
        }

        attributes #0 = { alwaysinline nounwind ssp uwtable }
        """, "main"), Int32, Tuple{Ptr{FILE}, Ptr{UInt8}, T}, fp, fmt, n)
    end

    function printf(fmt::Ptr{UInt8}, n::T) where T <: Union{Int32, UInt32}
        Base.llvmcall(("""
        ; External declaration of the printf function
        declare i32 @printf(i8* noalias nocapture, ...)

        define i32 @main(i64 %jlf, i32 %n) #0 {
        entry:
          %fmt = inttoptr i64 %jlf to i8*
          %status = call i32 (i8*, ...) @printf(i8* %fmt, i32 %n)
          ret i32 0
        }

        attributes #0 = { alwaysinline nounwind ssp uwtable }
        """, "main"), Int32, Tuple{Ptr{UInt8}, T}, fmt, n)
    end
    function printf(fp::Ptr{FILE}, fmt::Ptr{UInt8}, n::T) where T <: Union{Int32, UInt32}
        Base.llvmcall(("""
        ; External declaration of the printf function
        declare i32 @fprintf(i8*, ...)

        define i32 @main(i64 %jlfp, i64 %jlf, i32 %n) #0 {
        entry:
          %fp = inttoptr i64 %jlfp to i8*
          %fmt = inttoptr i64 %jlf to i8*
          %status = call i32 (i8*, ...) @fprintf(i8* %fp, i8* %fmt, i32 %n)
          ret i32 %status
        }

        attributes #0 = { alwaysinline nounwind ssp uwtable }
        """, "main"), Int32, Tuple{Ptr{FILE}, Ptr{UInt8}, T}, fp, fmt, n)
    end

    function printf(fmt::Ptr{UInt8}, n::T) where T <: Union{Int16, UInt16}
        Base.llvmcall(("""
        ; External declaration of the printf function
        declare i32 @printf(i8* noalias nocapture, ...)

        define i32 @main(i64 %jlf, i16 %n) #0 {
        entry:
          %fmt = inttoptr i64 %jlf to i8*
          %status = call i32 (i8*, ...) @printf(i8* %fmt, i16 %n)
          ret i32 0
        }

        attributes #0 = { alwaysinline nounwind ssp uwtable }
        """, "main"), Int32, Tuple{Ptr{UInt8}, T}, fmt, n)
    end
    function printf(fp::Ptr{FILE}, fmt::Ptr{UInt8}, n::T) where T <: Union{Int16, UInt16}
        Base.llvmcall(("""
        ; External declaration of the printf function
        declare i32 @fprintf(i8*, ...)

        define i32 @main(i64 %jlfp, i64 %jlf, i16 %n) #0 {
        entry:
          %fp = inttoptr i64 %jlfp to i8*
          %fmt = inttoptr i64 %jlf to i8*
          %status = call i32 (i8*, ...) @fprintf(i8* %fp, i8* %fmt, i16 %n)
          ret i32 %status
        }

        attributes #0 = { alwaysinline nounwind ssp uwtable }
        """, "main"), Int32, Tuple{Ptr{FILE}, Ptr{UInt8}, T}, fp, fmt, n)
    end

    function printf(fmt::Ptr{UInt8}, n::T) where T <: Union{Int8, UInt8}
        Base.llvmcall(("""
        ; External declaration of the printf function
        declare i32 @printf(i8* noalias nocapture, ...)

        define i32 @main(i64 %jlf, i8 %n) #0 {
        entry:
          %fmt = inttoptr i64 %jlf to i8*
          %status = call i32 (i8*, ...) @printf(i8* %fmt, i8 %n)
          ret i32 0
        }

        attributes #0 = { alwaysinline nounwind ssp uwtable }
        """, "main"), Int32, Tuple{Ptr{UInt8}, T}, fmt, n)
    end
    function printf(fp::Ptr{FILE}, fmt::Ptr{UInt8}, n::T) where T <: Union{Int8, UInt8}
        Base.llvmcall(("""
        ; External declaration of the printf function
        declare i32 @fprintf(i8*, ...)

        define i32 @main(i64 %jlfp, i64 %jlf, i8 %n) #0 {
        entry:
          %fp = inttoptr i64 %jlfp to i8*
          %fmt = inttoptr i64 %jlf to i8*
          %status = call i32 (i8*, ...) @fprintf(i8* %fp, i8* %fmt, i8 %n)
          ret i32 %status
        }

        attributes #0 = { alwaysinline nounwind ssp uwtable }
        """, "main"), Int32, Tuple{Ptr{FILE}, Ptr{UInt8}, T}, fp, fmt, n)
    end

## ---


text:
"""
    AbstractAlgorithm

Abstract supertype of all optimization algorithms.
"""
abstract type AbstractAlgorithm end

"""
    UnivariateAlgorithm

Abstract supertype of univariate optimization algorithms.
"""
abstract type UnivariateAlgorithm <: AbstractAlgorithm end

"""
    MultivariateAlgorithm

Abstract supertype of multivariate optimization algorithms.
"""
abstract type MultivariateAlgorithm <: AbstractAlgorithm end

"""
    Solution{Tx, Tf}

Solution type for optimization models.
"""
struct Solution{Tx, Tf, M}
    converged::Bool
    iter::Int
    minimizer::Tx
    minimum::Tf
    metadata::M
end

Solution(converged, iter, minimizer, minimum) = Solution(converged, iter, minimizer, minimum, Dict())

function Base.show(io::IO, sol::Solution)
    println(io, "Solution Summary:")
    println(io, "  • Converged: ", sol.converged)
    println(io, "  • Total Iterations: ", sol.iter)
    println(io, "  • Minimizer: ", sol.minimizer)
    print(io, "  • Minimum: ", sol.minimum)
end


text:
using SerialPorts

sensor_1 = SerialPort("/dev/ttyUSB0", 9600)
sensor_2 = SerialPort("/dev/ttyUSB1", 9600)

i = 0
while i < 50
	#println("Sensor 1:", readavailable(sensor_1))
	#println("Sensor 2:", readavailable(sensor_2))
	x = readavailable(sensor_1)
	if x != ""
		println("Sensor 1 ", x)
		i += 1
		x = ""
	end
	y = readavailable(sensor_2)
	if y != ""
		println("Sensor 2 ", y)
		i += 1
		y = ""
	end
end


text:
using SafeTestsets

@safetestset "diagrams" begin
    include("diagrams.jl")
end
@safetestset "aqua" begin
    include("aqua.jl")
end


text:
export pid, pidplots, rlocus, leadlink, laglink, leadlinkat, leadlinkcurve, stabregionPID, loopshapingPI

"""
Calculates and returns a PID controller on transfer function form.
`time` indicates whether or not the parameters are given as gains (default) or as time constants
`series` indicates  whether or not the series form or parallel form (default) is desired

`C = pid(; kp=0, ki=0; kd=0, time=false, series=false)`
"""
function pid(; kp=0., ki=0., kd=0., time=false, series=false)
    s = tf("s")
    if series
        return time ? kp*(one(kp) + one(kp)/(ki*s) + kd*s) : kp*(one(kp) + ki/s + kd*s)
    else
        return time ? kp + one(kp)/(ki*s) + kd*s : kp + ki/s + kd*s
    end
end

"""
Plots interesting figures related to closing the loop around process `P` with a PID controller
Send in a bunch of PID-parameters in any of the vectors kp, ki, kd. The vectors must be the same length.

`time` indicates whether or not the parameters are given as gains (default) or as time constants

`series` indicates  whether or not the series form or parallel form (default) is desired

Available plots are `:gof` for Gang of four, `:nyquist`, `:controller` for a bode plot of the controller TF and `:pz` for pole-zero maps

One can also supply a frequency vector ω to be used in Bode and Nyquist plots

`pidplots(P, args...; kps=0, kis=0, kds=0, time=false, series=false, ω=0)`

See also `loopshapingPI`, `stabregionPID`
"""
function pidplots(P::LTISystem, args...; kps=0, kis=0, kds=0, time=false, series=false, ω=0, grid = false)

    if grid
        kp = [i for i in kps, j in kis, k in kds][:]
        ki = [j for i in kps, j in kis, k in kds][:]
        kd = [k for i in kps, j in kis, k in kds][:]
        kps, kis, kds = kp, ki, kd
    else
        n = max(length(kps), length(kis), length(kds))
        kps = kps == 0 ? zeros(n) : kps
        kis = kis == 0 ? zeros(n) : kis
        kds = kds == 0 ? zeros(n) : kds
    end
    ω   = ω   == 0 ? exp10.(range(-3, stop=3, length=500)) : ω

    getColorSys(i)   = convert(Colors.RGB,Colors.HSV(360*((i-1)/(length(kps)))^1.5,0.9,0.8))

    gof_        = in(:gof        ,args)
    nyquist_    = in(:nyquist    ,args)
    controller_ = in(:controller ,args)
    pz_         = in(:pz         ,args)
    nichols_    = in(:nichols    ,args)

    Cs = LTISystem[]
    PCs = LTISystem[]
    Ts  = LTISystem[]
    labels = Array{String,2}(undef, 1,length(kps))
    colors =  Array{Colors.RGB{Float64},2}(undef, 1, length(kps))
    for (i,kp) = enumerate(kps)
        ki = kis[i]
        kd = kds[i]
        label = "\$k_p\$ = $(round(kp, digits=3)),      \$k_i\$ = $(round(ki, digits=3)),      \$k_d\$ = $(round(kd, digits=3))"

        C = pid(kp=kp,ki=ki,kd=kd,time=time,series=series)
        S,D,N,T = gangoffour(P,C)
        push!(Cs, C)
        push!(PCs, P*C)
        push!(Ts, T)
        labels[i] = label
    end

    if nyquist_
        nq = nyquistplot(PCs, ω, lab=labels, title="Nyquist curves")
    end
    if gof_
        bd = gangoffourplot(P, Cs, ω, lab=labels)
    end
    if pz_
        pzmap(Ts, title="Pole-zero map")
    end
    if controller_
        cplot = bodeplot(Cs, ω, lab=labels, title="Controller bode plot")
    end

end

@userplot Rlocusplot
@deprecate rlocus(args...;kwargs...) rlocusplot(args...;kwargs...)

# function getpoles(G,K)
#     poles = map(k -> pole(feedback(G,tf(k))), K)
#     cat(2,poles...)'
# end

function getpoles(G, K) # If OrdinaryDiffEq is installed, we override getpoles with an adaptive method
    P          = numpoly(G)[1]
    Q          = denpoly(G)[1]
    f          = (y,_,k) -> ComplexF64.(Polynomials.roots(k[1]*P+Q))
    prob       = OrdinaryDiffEq.ODEProblem(f,f(0.,0.,0.),(0.,K[end]))
    integrator = OrdinaryDiffEq.init(prob,OrdinaryDiffEq.Tsit5(),reltol=1e-8,abstol=1e-8)
    poleout    = Vector{Vector{ComplexF64}}()
    for i in integrator
       push!(poleout,integrator.k[1])
    end
    poleout = hcat(poleout...)'
end

"""
    rlocusplot(P::LTISystem, K)
Computes and plots the root locus of the SISO LTISystem P with
a negative feedback loop and feedback gains `K`, if `K` is not provided,
range(1e-6,stop=500,length=10000) is used.
If `OrdinaryDiffEq.jl` is installed and loaded by the user (`using OrdinaryDiffEq`), `rlocusplot` will use an adaptive step-size algorithm to
select values of `K`. A scalar `Kmax` can then be given as second argument.
"""
rlocus
@recipe function rlocus(p::Rlocusplot; K=Float64[])
    P = p.args[1]
    K = isempty(K) ? range(1e-6,stop=500,length=10000) : K
    Z = tzero(P)

    poles = getpoles(P,K)
    redata = real.(poles)
    imdata = imag.(poles)
    ylim = (max(-50,minimum(imdata)), min(50,maximum(imdata)))
    xlim = (max(-50,minimum(redata)), min(50,maximum(redata)))
    title --> "Root locus"
    xguide --> "Re(roots)"
    yguide --> "Im(roots)"
    @series begin
        legend --> false
        ylims  --> ylim
        xlims  --> xlim
        label := ""
        redata, imdata
    end
    @series begin
        seriestype := :scatter
        markershape --> :circle
        markersize --> 10
        label --> "Zeros"
        real.(Z), imag.(Z)
    end
    @series begin
        seriestype := :scatter
        markershape --> :xcross
        markersize --> 10
        label --> "Open-loop poles"
        redata[1,:], imdata[1,:]
    end
end

"""
`laglink(a, M; h=0)`

Returns a phase retarding link, the rule of thumb `a = 0.1ωc` guarantees less than 6 degrees phase margin loss. The bode curve will go from `M`, bend down at `a/M` and level out at 1 for frequencies > `a`
"""
function laglink(a, M; h=0)
    numerator = [1/a, 1]
    denominator = [M/a, 1]
    gain = M
    G = tf(gain*numerator,denominator)
    return  h <= 0 ? G : c2d(G,h)
end

"""
`leadlink(b, N, K; h=0)`

Returns a phase advancing link, the top of the phase curve is located at `ω = b√(N)` where the link amplification is `K√(N)` The bode curve will go from `K`, bend up at `b` and level out at `KN` for frequencies > `bN`

The phase advance at `ω = b√(N)` can be plotted as a function of `N` with `leadlinkcurve()`

Values of `N < 1` will give a phase retarding link.

See also `leadlinkat` `laglink`
"""
function leadlink(b, N, K; h=0)
    numerator = [1/b, 1]
    denominator = [1/(b*N), 1]
    gain = K
    G = tf(gain*numerator,denominator)
    return  h <= 0 ? G : c2d(G,h)

end

"""
`leadlinkat(ω, N, K; h=0)`

Returns a phase advancing link, the top of the phase curve is located at `ω` where the link amplification is `K√(N)` The bode curve will go from `K`, bend up at `ω/√(N)` and level out at `KN` for frequencies > `ω√(N)`

The phase advance at `ω` can be plotted as a function of `N` with `leadlinkcurve()`

Values of `N < 1` will give a phase retarding link.

See also `leadlink` `laglink`
"""
function leadlinkat(ω, N, K; h=0)
    b = ω / sqrt(N)
    return leadlink(b,N,K,h=h)
end

"""
Plot the phase advance as a function of `N` for a lead link (phase advance link)

If an input argument `s` is given, the curve is plotted from `s` to 10, else from 1 to 10.

See also `Leadlink, leadlinkat`
"""
function leadlinkcurve(start=1)
    N = range(start, stop=10, length=50)
    dph = 180/pi*map(Ni->atan(sqrt(Ni))-atan(1/sqrt(Ni)), N)
    Plots.plot(N,dph, xlabel="N", ylabel="Phase advance [deg]")
end

"""
`fig, kp, ki = stabregionPID(P, [ω]; kd=0, doplot = true)`

Segments of the curve generated by this program
is the boundary of the stability region for a
process with transfer function P(s)
The PID controller is assumed to be on the form kp +ki/s +kd s

The curve is found by analyzing
P(s)*C(s) = -1 ⟹
|PC| = |P| |C| = 1
arg(P) + arg(C) = -π

If `P` is a function (e.g. s -> exp(-sqrt(s)) ), the stability of feedback loops using PI-controllers can be analyzed for processes with models with arbitrary analytic functions

See also `stabregionPID`, `loopshapingPI`, `pidplots`
"""
function stabregionPID(P, ω = _default_freq_vector(P,:bode); kd=0, doplot = true)
    Pv  = freqresp(P,ω)[:,1,1]
    r   = abs.(Pv)
    phi = angle.(Pv)
    kp  = -cos.(phi)./r
    ki  = kd.*ω.^2 .- ω.*sin.(phi)./r
    Plots.plot(kp,ki,linewidth = 1.5, xlabel="\$k_p\$", ylabel="\$k_i\$", title="Stability region of P,     \$k_d\$ = $(round(kd, digits=4))"), kp, ki
end

function stabregionPID(P::Function, ω = exp10.(range(-3, stop=1, length=50)); kd=0, doplot = true)
    Pv      = P.(im*ω)
    r       = abs.(Pv)
    phi     = angle.(Pv)
    kp      = -cos.(phi)./r
    ki      = kd.*ω.^2 .- ω.*sin.(phi)./r
    Plots.plot(kp,ki,linewidth = 1.5, xlabel="\$k_p\$", ylabel="\$k_i\$", title="Stability region of P,     \$k_d\$ = $(round(kd, digits=4))"), kp, ki
end

"""
`kp,ki,C = loopshapingPI(P,ωp; ϕl,rl, phasemargin)`

Selects the parameters of a PI-controller such that the Nyquist curve of `P` at the frequency `ωp` is moved to `rl exp(i ϕl)`

If `phasemargin` is supplied, `ϕl` is selected such that the curve is moved to an angle of `phasemargin - 180` degrees

If no `rl` is given, the magnitude of the curve at `ωp` is kept the same and only the phase is affected, the same goes for `ϕl` if no phasemargin is given.

See also `pidplots`, `stabregionPID`
"""
function loopshapingPI(P,ωp; ϕl=0,rl=0, phasemargin = 0, doplot = false)
    Pw = P(im*ωp)[1]
    ϕp = angle(Pw)
    rp = abs.(Pw)

    if phasemargin > 0
        ϕl = deg2rad(-180+phasemargin)
    else
        ϕl = ϕl == 0 ? ϕp : ϕl
    end
    rl = rl == 0 ? rp : rl

    kp = rl/rp*cos(ϕp-ϕl)
    ki = rl*ωp/rp*sin(ϕp-ϕl)

    C = pid(kp=kp, ki=ki)

    if doplot
        gangoffourplot(P,[tf(1),C])
        nyquistplot([P, P*C])
    end
    return kp,ki,C
end


text:
module ConditionalDistributionForest

using Random
using DecisionTree

export compute_conditional_quantiles!, compute_conditional_quantiles, compute_conditional_distribution!, compute_conditional_distribution, build_distribution_forest

include("tree_tools.jl")

function build_distribution_forest(
    labels              :: Vector{T},
    features            :: Matrix{S},
    n_subfeatures       = -1,
    n_trees             = 10,
    partial_sampling    = 0.7,
    max_depth           = -1,
    min_samples_leaf    = 1,
    min_samples_split   = 2,
    min_purity_increase = 0.0;
    rng                 = Random.GLOBAL_RNG,
    return_bootstrap    = false) where {S, T}

    if return_bootstrap
        return build_forest_bootstraps(labels,
                                       features,
                                       n_subfeatures,
                                       n_trees,
                                       partial_sampling,
                                       max_depth,
                                       min_samples_leaf,
                                       min_samples_split,
                                       min_purity_increase,
                                       rng = rng)
    else
        return build_forest(labels,
                            features,
                            n_subfeatures,
                            n_trees,
                            partial_sampling,
                            max_depth,
                            min_samples_leaf,
                            min_samples_split,
                            min_purity_increase,
                            rng = rng)
    end
end

function compute_conditional_quantiles!(quant::AbstractMatrix{T},
                                        forest::DecisionTree.Ensemble{S, T},
                                        Y::Vector{T},
                                        X::Matrix{S},
                                        Xvec::Matrix{S},
                                        quantiles::Vector{U},
                                        bootstraps = nothing) where {S, T, U <:AbstractFloat}
    n_prime = size(Xvec, 1)
    n_quantiles = length(quantiles)

    weights = compute_forest_weights(forest, X, Xvec, bootstraps)

    perm_Y = sortperm(Y)
    for iquantile in 1:n_quantiles
        for i in 1:n_prime
            compt = 0
            sum_weights = 0.
            while(sum_weights <= quantiles[iquantile])
                compt += 1
                sum_weights += weights[perm_Y[compt], i]
            end

            quant[i, iquantile] = Y[perm_Y[compt]]
        end
    end

    return nothing
end

function compute_conditional_quantiles!(quant::AbstractMatrix{T},
                                        Y::Vector{T},
                                        X::Matrix{S},
                                        Xvec::Matrix{S},
                                        quantiles::Vector{U},
                                        n_subfeatures       = -1,
                                        n_trees             = 10,
                                        partial_sampling    = 0.7,
                                        max_depth           = -1,
                                        min_samples_leaf    = 1,
                                        min_samples_split   = 2,
                                        min_purity_increase = 0.0;
                                        rng                 = Random.GLOBAL_RNG) where {S, T, U <:AbstractFloat}
    forest, bootstraps = build_forest_bootstraps(Y, X, n_subfeatures, n_trees, partial_sampling, max_depth, min_samples_leaf, min_samples_split, min_purity_increase, rng = rng)

    compute_conditional_quantiles!(quant, forest, Y, X, Xvec, quantiles, bootstraps)

    return nothing
end

function compute_conditional_quantiles(forest::DecisionTree.Ensemble{S, T},
                                       Y::Vector{T},
                                       X::Matrix{S},
                                       Xvec::Matrix{S},
                                       quantiles::Vector{U},
                                       bootstraps = nothing) where {S, T, U <:AbstractFloat}
    quant = Matrix{T}(undef, size(Xvec, 1), length(quantiles))
    compute_conditional_quantiles!(quant, forest, Y, X, Xvec, quantiles, bootstraps)
    return quant
end

function compute_conditional_quantiles(Y::Vector{T},
                                       X::Matrix{S},
                                       Xvec::Matrix{S},
                                       quantiles::Vector{U},
                                       n_subfeatures       = -1,
                                       n_trees             = 10,
                                       partial_sampling    = 0.7,
                                       max_depth           = -1,
                                       min_samples_leaf    = 1,
                                       min_samples_split   = 2,
                                       min_purity_increase = 0.0;
                                       rng                 = Random.GLOBAL_RNG) where {S, T, U <:AbstractFloat}
    forest, bootstraps = build_forest_bootstraps(Y, X, n_subfeatures, n_trees, partial_sampling, max_depth, min_samples_leaf, min_samples_split, min_purity_increase, rng = rng)

    return compute_conditional_quantiles(forest, Y, X, Xvec, quantiles, bootstraps)
end

function compute_conditional_distribution!(fdr::AbstractMatrix{U},
                                           forest::DecisionTree.Ensemble{S, T},
                                           Y::Vector{T},
                                           X::Matrix{S},
                                           Xvec::Matrix{S},
                                           Ydiscr::Vector{T},
                                           bootstraps = nothing) where {S, T, U <:AbstractFloat}
    weights = compute_forest_weights(forest, X, Xvec, bootstraps)

    for iy in 1:length(Ydiscr)
        @views fdr[iy, :] .= vec(sum(weights[Y .<= Ydiscr[iy], :], dims=1))
    end

    return nothing
end

function compute_conditional_distribution!(fdr::AbstractMatrix{U},
                                           Y::Vector{T},
                                           X::Matrix{S},
                                           Xvec::Matrix{S},
                                           Ydiscr::Vector{T},
                                           n_subfeatures       = -1,
                                           n_trees             = 10,
                                           partial_sampling    = 0.7,
                                           max_depth           = -1,
                                           min_samples_leaf    = 1,
                                           min_samples_split   = 2,
                                           min_purity_increase = 0.0;
                                           rng                 = Random.GLOBAL_RNG) where {S, T, U <:AbstractFloat}
    forest, bootstraps = build_forest_bootstraps(Y, X, n_subfeatures, n_trees, partial_sampling, max_depth, min_samples_leaf, min_samples_split, min_purity_increase, rng = rng)

    compute_conditional_distribution!(fdr, forest, Y, X, Xvec, Ydiscr, bootstraps)

    return nothing
end

function compute_conditional_distribution(forest::DecisionTree.Ensemble{S, T},
                                          Y::Vector{T},
                                          X::Matrix{S},
                                          Xvec::Matrix{S},
                                          Ydiscr::Vector{T},
                                          bootstraps = nothing) where {S, T}
    fdr = Array{Float64, 2}(undef, length(Ydiscr), size(Xvec, 1))
    compute_conditional_distribution!(fdr, forest, Y, X, Xvec, Ydiscr, bootstraps)

    return fdr
end

function compute_conditional_distribution(Y::Vector{T},
                                          X::Matrix{S},
                                          Xvec::Matrix{S},
                                          Ydiscr::Vector{T},
                                          n_subfeatures       = -1,
                                          n_trees             = 10,
                                          partial_sampling    = 0.7,
                                          max_depth           = -1,
                                          min_samples_leaf    = 1,
                                          min_samples_split   = 2,
                                          min_purity_increase = 0.0;
                                          rng                 = Random.GLOBAL_RNG) where {S, T}
    forest, bootstraps = build_forest_bootstraps(Y, X, n_subfeatures, n_trees, partial_sampling, max_depth, min_samples_leaf, min_samples_split, min_purity_increase, rng = rng)

    fdr = Array{Float64, 2}(undef, length(Ydiscr), size(Xvec, 1))

    compute_conditional_distribution!(fdr, forest, Y, X, Xvec, Ydiscr, bootstraps)

    return fdr
end

function build_forest_bootstraps(
    labels              :: Vector{T},
    features            :: Matrix{S},
    n_subfeatures       = -1,
    n_trees             = 10,
    partial_sampling    = 0.7,
    max_depth           = -1,
    min_samples_leaf    = 1,
    min_samples_split   = 2,
    min_purity_increase = 0.0;
    rng                 = Random.GLOBAL_RNG) where {S, T}

    if n_trees < 1
        throw("the number of trees must be >= 1")
    end
    if !(0.0 < partial_sampling <= 1.0)
        throw("partial_sampling must be in the range (0,1]")
    end

    if n_subfeatures == -1
        n_features = size(features, 2)
        n_subfeatures = round(Int, sqrt(n_features))
    end

    t_samples = length(labels)
    n_samples = floor(Int, partial_sampling * t_samples)

    rngs = DecisionTree.mk_rng(rng)::Random.AbstractRNG
    forest = Vector{DecisionTree.LeafOrNode{S, T}}(undef, n_trees)

    bootstraps = Array{Int64, 2}(undef, n_samples, n_trees)
    Threads.@threads for i in 1:n_trees
        bootstraps[:, i]  .= rand(rngs, 1:t_samples, n_samples)
        forest[i] = build_tree(
            labels[bootstraps[:, i]],
            features[bootstraps[:, i], :],
            n_subfeatures,
            max_depth,
            min_samples_leaf,
            min_samples_split,
            min_purity_increase,
            rng = rngs)
    end

    return Ensemble{S, T}(forest), bootstraps
end

function compute_forest_weights(model, X, Xvec, bootstraps)
    n_samples = size(X, 1)
    n_prime = size(Xvec, 1)
    n_trees = length(model.trees)

    weights = zeros(n_samples, n_prime)
    for l in 1:n_trees
	    leafids = get_leaf_of_feature(model.trees[l], Xvec)
        @views repartition = get_repartition(model.trees[l], X, bootstraps[:, l])

	    for i in 1:n_prime
            val = 1 / length(repartition[leafids[i]])
	        for j in repartition[leafids[i]]
	            weights[j, i] += val
	        end
	    end
    end

    weights ./= n_trees
    return weights
end

function compute_forest_weights(model, X, Xvec, ::Nothing)
    n_samples = size(X, 1)
    n_prime = size(Xvec, 1)
    n_trees = length(model.trees)

    weights = zeros(n_samples, n_prime)
    for l in 1:n_trees
	    leafids = get_leaf_of_feature(model.trees[l], Xvec)
        repartition = get_repartition(model.trees[l], X, 1:n_samples)

	    for i in 1:n_prime
            val = 1 / length(repartition[leafids[i]])
	        for j in repartition[leafids[i]]
	            weights[j, i] += val
	        end
	    end
    end

    weights ./= n_trees
    return weights
end

end # module


text:
using Test
using TopOpt
using TopOpt.TopOptProblems: getE
using Base.Iterators
# if get(ENV, "CI", nothing) != "true"
#     import Makie
#     @eval using TopOpt.TrussTopOptProblems.TrussVisualization: visualize
# end

ins_dir = joinpath(@__DIR__, "instances", "ground_meshes");

@testset "Tim problem to solve" for (problem_dim, lc_ind) in product(["2d", "3d"], [0, 1])
    # problem_dim = "2d"
    # file_name = "tim_2d.json"
    # lc_ind = 1
    file_name = "tim_$(problem_dim).json"
    problem_file = joinpath(ins_dir, file_name)

    node_points, elements, mats, crosssecs, fixities, load_cases = parse_truss_json(problem_file);
    loads = load_cases[string(lc_ind)]

    problem = TrussProblem(Val{:Linear}, node_points, elements, loads, fixities, mats, crosssecs);

    ndim, nnodes, ncells = length(node_points[1]), length(node_points), length(elements)
    @test getE(problem) == [m.E for m in mats]

    V = 0.3 # volume fraction
    xmin = 0.001 # minimum density
    rmin = 4.0; # density filter radius

    penalty = TopOpt.PowerPenalty(1.0) # 1
    solver = FEASolver(Displacement, Direct, problem, xmin = xmin,
        penalty = penalty);
    ## call solver to trigger assemble!
    solver()

    # ##############################
    # #! buckling
    # import TrussTopOpt.TrussTopOptProblems: default_quad_order
    # einfo = ElementFEAInfo(problem, TrussTopOpt.TrussTopOptProblems.default_quad_order(problem), Val{:Static})
    # ginfo = GlobalFEAInfo(problem)

    # using TrussTopOpt.TrussTopOptProblems: buckling, get_Kσs
    # buckling(problem, solver.globalinfo, solver.elementinfo)

    # ##############################

    # TODO TopOpt.LogBarrier
    # TODO linear_elasticity, du/dx
    # * Compliance
    comp = TopOpt.Compliance(problem, solver)
    obj = comp
    volfrac = TopOpt.Volume(problem, solver)
    constr = x -> volfrac(x) - V

    mma_options = options = Nonconvex.MMAOptions(
        maxiter = 3000, tol = Nonconvex.Tolerance(kkt = 0.001),
    )
    convcriteria = Nonconvex.KKTCriteria()
    x0 = fill(V, length(solver.vars))
    optimizer = Optimizer(
        obj, constr, x0, Nonconvex.MMA87(),
        options = mma_options, convcriteria = convcriteria,
    )
    simp = SIMP(optimizer, solver, penalty.p)

    # ? 1.0 might induce an infeasible solution, which gives the optimizer a hard time to escape
    # from infeasible regions and return a result
    result = simp(x0);

    println("="^10)
    println("tim-$(problem_dim) - LC $(lc_ind) - #elements $(ncells), #dof: $(ncells*ndim): opt iter $(simp.optimizer.workspace.iter)")
    println("$(result.convstate)")

    solver()

    # if get(ENV, "CI", nothing) != "true"
    #     fig = visualize(
    #         problem, solver.u; topology = result.topology, vector_arrowsize = 0.1,
    #         vector_linewidth=0.8, default_exagg_scale=ndim == 3 ? 1.0 : 0.01,
    #         exagg_range = ndim == 3 ? 10.0 : 0.1,
    #     )
    #     Makie.display(fig)
    # end
end # end testset


text:
function construct_device!(
    optimization_container::OptimizationContainer,
    sys::PSY.System,
    model::DeviceModel{R, D},
    ::Type{S},
) where {
    R <: PSY.RenewableGen,
    D <: AbstractRenewableDispatchFormulation,
    S <: PM.AbstractPowerModel,
}
    devices = get_available_components(R, sys)

    if !validate_available_devices(R, devices)
        return
    end

    # Variables
    add_variables!(optimization_container, ActivePowerVariable, devices, D())
    add_variables!(optimization_container, ReactivePowerVariable, devices, D())

    # Constraints
    add_constraints!(
        optimization_container,
        RangeConstraint,
        ActivePowerVariable,
        devices,
        model,
        S,
        get_feedforward(model),
    )
    add_constraints!(
        optimization_container,
        RangeConstraint,
        ReactivePowerVariable,
        devices,
        model,
        S,
        get_feedforward(model),
    )
    feedforward!(optimization_container, devices, model, get_feedforward(model))

    # Cost Function
    cost_function!(optimization_container, devices, model, S)

    return
end

function construct_device!(
    optimization_container::OptimizationContainer,
    sys::PSY.System,
    model::DeviceModel{R, D},
    ::Type{S},
) where {
    R <: PSY.RenewableGen,
    D <: AbstractRenewableDispatchFormulation,
    S <: PM.AbstractActivePowerModel,
}
    devices = get_available_components(R, sys)

    if !validate_available_devices(R, devices)
        return
    end

    # Variables
    add_variables!(optimization_container, ActivePowerVariable, devices, D())

    # Constraints
    add_constraints!(
        optimization_container,
        RangeConstraint,
        ActivePowerVariable,
        devices,
        model,
        S,
        get_feedforward(model),
    )
    feedforward!(optimization_container, devices, model, get_feedforward(model))

    # Cost Function
    cost_function!(optimization_container, devices, model, S)

    return
end

function construct_device!(
    optimization_container::OptimizationContainer,
    sys::PSY.System,
    ::DeviceModel{R, FixedOutput},
    ::Type{S},
) where {R <: PSY.RenewableGen, S <: PM.AbstractPowerModel}
    devices = get_available_components(R, sys)

    if !validate_available_devices(R, devices)
        return
    end

    nodal_expression!(optimization_container, devices, S)

    return
end


text:
using Particles
using Test
using Random
using LinearAlgebra

@testset "Particles" begin
    include("component.jl")
    include("statepriors.jl")
    include("particle.jl")
    include("fearnhead.jl")
    include("nix2.jl")
    include("niw.jl")
    include("labeled.jl")
    include("clustering.jl")
end


text:
abstract type SplittingCriterion  end
struct DifferenceSplit      <: SplittingCriterion end

abstract type SortBy end
struct SortByMean 	<: SortBy end

mutable struct ModelSettings
	model_type::String # 1
	minWeight::Float64 # 2
	randomw::Float64 # 3
	crit::SplittingCriterion # fn version of 4
	maxSplittingPoints::Int # 5
	iterations::Int # 6
	learningRate::Float64 # 7 moderationfactor
	nScores::Int # 8
	adaptiveLearningRate::Float64 # 9
	prem_buffer::Int # 11
	startAtMean::Bool # 12
	writeTree::Bool # 13
	number_of_num_features::Int # 14
	spawnsmaller::Bool # 26
    boolRankOptimization::Bool # 35
	boolRandomizeOnlySplitAtTopNode::Bool # 37
	subsampling_prop::Float64 # 38 #NOTE Friedman suggests to sample without replacement see Friedman 2002 Stochastic Gradient Boosting,” Computational Statistics and Data Analysis 38(4):367-378
	subsampling_features_prop::Float64 # 39
	version::String # 40
	preppedJLDFileExists::Bool # 41
	catSortByThreshold::Int # 42
	catSortBy::SortBy # 43
	scorebandsstartingpoints::Array{Int,1} # 44
	showTimeUsedByEachIteration::Bool # 45
	smoothEstimates::String # 46
	deriveFitPerScoreFromObservedRatios::Bool # deriveFitPerScoreFromObservedRatios::Bool (if this is true then the fit per score will be based on the observed ratio per score (instead of the fitted ratio per score)) . Note that this does not influence the structure of the tree, but only the fitted ratio per score (and scoreband)
    roptForcedPremIncr::Bool # 47
    writeSasCode::Bool # 49
	writeIterationMatrix::Bool # 50
	writeResult::Bool # 51
	writeStatistics::Bool # 52
	boolCreateZipFile::Bool # 53
	writeCsharpCode::Bool
	writeVbaCode::Bool
	nDepthToStartParallelization::Int
	baggingWeightTreesError::String
	cBB_niterBoosting::Int
	cBB_niterBagging::Int
	fixedinds::Array{Int,1}
	boolNumeratorStats::Bool
	statsByVariables::Array{Int,1}
	statsRandomByVariable::Int
	saveJLDFile::Bool # if this is set to false no *.jld2 file is saved #todo we should rename this to ensure it is clear that this refers to the prepped data only!
	saveResultAsJLDFile::Bool # this refers to the model output/result
	print_details::Bool # whether julia should print details to the log or not (note to developers: warnings/errors and certain important infos should always be printed!)
	seed::Int # optional entropy
	graphvizexecutable::String # e.g. C:\Program Files (x86)\Graphviz2.38\ if this
	showProgressBar_time::Bool # bar form iterators package (or some other package)
	prroduceEstAndLeafMatrices::Bool
	write_dot_graph::Bool
	calculateGini::Bool # whether or not to calculate the gini (which requires sorting of the data (which takes time!))
	calculatePoissonError::Bool
    performanceMeasure::String
	fitForStatsAndCharts::String
	ignoreZeroDenominatorValues::Bool

	# the following are treated specially
	df_name_vector::Array{String,1}
    number_of_char_features::Int # this is calculated by Julia from number_of_num_features and ncolsdfIndata (which is defined by the data)

	chosen_apply_tree_fn::String # defined below
	moderationvector::Array{Float64,1} # 50 onwards

    function ModelSettings()
	# initialize default settings
	model_type = "build_tree" # 1
	minWeight = -0.1 # 2
	randomw = 0.0 # 3
	crit = DifferenceSplit()	 # 4
	maxSplittingPoints = 250 # 5 #arbitrary choice. Preliminary testing has shown no performance difference for different values (say 10 or 200) for this setting
	iterations = 2 # 6
	learningRate = 0.1 # 7
	nScores = 1000 # 8
	adaptiveLearningRate = 1.0 # 9
	prem_buffer = 0 # 11
	startAtMean = true # 12
	writeTree = true # 13
	number_of_num_features = -1 # 14
	spawnsmaller = true # 26
	boolRankOptimization = false # 35
	boolRandomizeOnlySplitAtTopNode = true # 37
	subsampling_prop = 1.0 # 38
	subsampling_features_prop = 1.0 # 39
	version = "not_initialized" # 40
	preppedJLDFileExists = false # 41
	catSortByThreshold = 7 # 42
	catSortBy = SortByMean() # 43
	scorebandsstartingpoints = [parse(Int, x) for x in split("1,100,200,300,400,500,600,700,800,900", ',')] # scorebandsstartingpoints=parse(Int,split("1,100,200,300,400,500,600,700,800,900",',')) #44
	showTimeUsedByEachIteration = false # true #45
	smoothEstimates = "1" # 46 this is a string for now (as there could be different smoothing methods specified by this string)
	deriveFitPerScoreFromObservedRatios = true
	roptForcedPremIncr = false # 47
        writeSasCode = false # 49
	writeIterationMatrix = false# 50
	writeResult = false # true #51
	writeStatistics = true # 52
	boolCreateZipFile = false # true #53
	writeCsharpCode = false # true
	writeVbaCode = false
	nDepthToStartParallelization = -1
	baggingWeightTreesError = "mae" # ::String
	# combined bagging&boosting model
	cBB_niterBoosting = 0
	cBB_niterBagging = 0
	fixedinds = Array{Int}(undef, 0)
	boolNumeratorStats = false
	statsByVariables = Int[]
	statsRandomByVariable = 5
	saveJLDFile = true
	saveResultAsJLDFile = false
	print_details = true
	seed = 90210
	graphvizexecutable = "" # C:\\Program Files (x86)\\Graphviz2.38\\bin\\dot.exe""
	showProgressBar_time = true
	prroduceEstAndLeafMatrices = false
	write_dot_graph = false
	calculateGini = false
	calculatePoissonError = false
        performanceMeasure = "Lift Val"
	fitForStatsAndCharts = "rawRelativities" # rawRelativities,unsmoothedPerScore,smoothedPerScore
	ignoreZeroDenominatorValues = false

	# the following are treated specially
	df_name_vector = Array{String}(undef, 0) #
        number_of_char_features = -1 #
	chosen_apply_tree_fn = "apply_tree_by_leaf" # #apply_tree_by_row does not seem to work for (certain?) boosting models
	moderationvector = [0.1] #

	return new(model_type,minWeight,randomw,crit,maxSplittingPoints,iterations,learningRate,nScores,adaptiveLearningRate,prem_buffer,startAtMean,writeTree,number_of_num_features,spawnsmaller,boolRankOptimization,boolRandomizeOnlySplitAtTopNode,subsampling_prop,subsampling_features_prop,version,preppedJLDFileExists,catSortByThreshold,catSortBy,scorebandsstartingpoints,showTimeUsedByEachIteration,smoothEstimates,deriveFitPerScoreFromObservedRatios,roptForcedPremIncr,writeSasCode,writeIterationMatrix,writeResult,writeStatistics,boolCreateZipFile,writeCsharpCode,writeVbaCode,nDepthToStartParallelization,baggingWeightTreesError,cBB_niterBoosting,cBB_niterBagging,fixedinds,boolNumeratorStats,statsByVariables,statsRandomByVariable,saveJLDFile,saveResultAsJLDFile,print_details,seed,graphvizexecutable,showProgressBar_time,prroduceEstAndLeafMatrices,write_dot_graph,calculateGini,calculatePoissonError,performanceMeasure,fitForStatsAndCharts,ignoreZeroDenominatorValues
	 ,df_name_vector,number_of_char_features,chosen_apply_tree_fn,moderationvector)
    end  # ModelSettings()

end # type definition

function writeAllFieldsToArray(s)
	nn = fieldnames(typeof(s))
	res = Array{AbstractString}(undef, length(nn), 2)
	i = 0
	for x in nn
		i += 1
        @show x
        tmp = AbstractString[string(x) string(getfield(s, x))][:]
        @show tmp
        @show size(tmp)
        @show typeof(tmp)
		res[i,:] = AbstractString[string(x) string(getfield(s, x))][:]
	end
	return res
end

mi = ModelSettings()
foo = writeAllFieldsToArray(mi)


text:
# canonical data version: 1.3.0

using Test

include("anagram.jl")

@testset "no matches" begin
    @test detect_anagrams("diaper", ["hello", "world", "zombies", "pants"]) == []
end

@testset "detects simple anagram" begin
    @test detect_anagrams("ant", ["tan", "stand", "at"]) == ["tan"]
end

@testset "does not detect false positives" begin
    @test detect_anagrams("galea", ["eagle"]) == []
end

@testset "detects two anagrams" begin
    @test detect_anagrams("master", ["stream", "pigeon", "maters"]) == ["stream", "maters"]
end

@testset "does not detect anagram subsets" begin
    @test detect_anagrams("good", ["dog", "goody"]) == []
end

@testset "detects anagram" begin
    @test detect_anagrams("listen", ["enlists", "google", "inlets", "banana"]) == ["inlets"]
end

@testset "detects three anagrams" begin
    @test detect_anagrams("allergy", ["gallery", "ballerina", "regally", "clergy", "largely", "leading"]) == ["gallery", "regally", "largely"]
end

@testset "does not detect identical words" begin
    @test detect_anagrams("corn", ["corn", "dark", "Corn", "rank", "CORN", "cron", "park"]) == ["cron"]
end

@testset "does not detect non-anagrams with identical checksum" begin
    @test detect_anagrams("mass", ["last"]) == []
end

@testset "detects anagrams case-insensitively" begin
    @test detect_anagrams("Orchestra", ["cashregister", "Carthorse", "radishes"]) == ["Carthorse"]
end

@testset "detects anagrams using case-insensitive subject" begin
    @test detect_anagrams("Orchestra", ["cashregister", "carthorse", "radishes"]) == ["carthorse"]
end

@testset "detects anagrams using case-insensitive possible matches" begin
    @test detect_anagrams("orchestra", ["cashregister", "Carthorse", "radishes"]) == ["Carthorse"]
end

@testset "does not detect a word as its own anagram" begin
    @test detect_anagrams("banana", ["Banana"]) == []
end

@testset "does not detect a anagram if the original word is repeated" begin
    @test detect_anagrams("go", ["go Go GO"]) == []
end

@testset "anagrams must use all letters exactly once" begin
    @test detect_anagrams("tapper", ["patter"]) == []
end

@testset "capital word is not own anagram" begin
    @test detect_anagrams("BANANA", ["Banana"]) == []
end


text:
include("lib_for_robot.jl")
function task1(r)
    for side in (HorizonSide(i) for i=0:3)
        putmarks(r,side)
        move_marks(r,reverse(side))
    end
    putmarker!(r)
end


text:
# stuff that is universal to all power models

export
    GenericPowerModel,
    setdata, setsolver, solve,
    run_generic_model, build_generic_model, solve_generic_model,
    ids, ref, var, ext

""
abstract type AbstractPowerFormulation end

""
abstract type AbstractConicPowerFormulation <: AbstractPowerFormulation end

"""
```
type GenericPowerModel{T<:AbstractPowerFormulation}
    model::JuMP.Model
    data::Dict{String,Any}
    setting::Dict{String,Any}
    solution::Dict{String,Any}
    var::Dict{Symbol,Any} # model variable lookup
    ref::Dict{Symbol,Any} # reference data
    ext::Dict{Symbol,Any} # user extentions
end
```
where

* `data` is the original data, usually from reading in a `.json` or `.m` (patpower) file,
* `setting` usually looks something like `Dict("output" => Dict("branch_flows" => true))`, and
* `ref` is a place to store commonly used pre-computed data from of the data dictionary,
    primarily for converting data-types, filtering out deactivated components, and storing
    system-wide values that need to be computed globally. See `build_ref(data)` for further details.

Methods on `GenericPowerModel` for defining variables and adding constraints should

* work with the `ref` dict, rather than the original `data` dict,
* add them to `model::JuMP.Model`, and
* follow the conventions for variable and constraint names.
"""
mutable struct GenericPowerModel{T<:AbstractPowerFormulation}
    model::Model

    data::Dict{String,Any}
    setting::Dict{String,Any}
    solution::Dict{String,Any}

    ref::Dict{Symbol,Any} # data reference data
    var::Dict{Symbol,Any} # JuMP variables
    con::Dict{Symbol,Any} # JuMP constraint references
    cnw::Int # current network index value

    # Extension dictionary
    # Extensions should define a type to hold information particular to
    # their functionality, and store an instance of the type in this
    # dictionary keyed on an extension-specific symbol
    ext::Dict{Symbol,Any}
end

# default generic constructor
function GenericPowerModel(data::Dict{String,Any}, T::DataType; ext = Dict{String,Any}(), setting = Dict{String,Any}(), solver = JuMP.UnsetSolver())

    # TODO is may be a good place to check component connectivity validity
    # i.e. https://github.com/lanl-ansi/PowerModels.jl/issues/131

    ref = build_ref(data) # refrence data

    var = Dict{Symbol,Any}(:nw => Dict{Int,Any}())
    con = Dict{Symbol,Any}(:nw => Dict{Int,Any}())
    for nw_id in keys(ref[:nw])
        var[:nw][nw_id] = Dict{Symbol,Any}()
        con[:nw][nw_id] = Dict{Symbol,Any}()
    end

    cnw = minimum([k for k in keys(ref[:nw])])

    pm = GenericPowerModel{T}(
        Model(solver = solver), # model
        data, # data
        setting, # setting
        Dict{String,Any}(), # solution
        ref,
        var, # vars
        con,
        cnw,
        ext # ext
    )

    return pm
end

### Helper functions for ignoring multinetwork support
ids(pm::GenericPowerModel, key::Symbol) = ids(pm, pm.cnw, key)
ids(pm::GenericPowerModel, n::Int, key::Symbol) = keys(pm.ref[:nw][n][key])

ref(pm::GenericPowerModel, key::Symbol) = ref(pm, pm.cnw, key)
ref(pm::GenericPowerModel, key::Symbol, idx) = ref(pm, pm.cnw, key, idx)
ref(pm::GenericPowerModel, n::Int, key::Symbol) = pm.ref[:nw][n][key]
ref(pm::GenericPowerModel, n::Int, key::Symbol, idx) = pm.ref[:nw][n][key][idx]

Base.var(pm::GenericPowerModel, key::Symbol) = var(pm, pm.cnw, key)
Base.var(pm::GenericPowerModel, key::Symbol, idx) = var(pm, pm.cnw, key, idx)
Base.var(pm::GenericPowerModel, n::Int, key::Symbol) = pm.var[:nw][n][key]
Base.var(pm::GenericPowerModel, n::Int, key::Symbol, idx) = pm.var[:nw][n][key][idx]

con(pm::GenericPowerModel, key::Symbol) = con(pm, pm.cnw, key)
con(pm::GenericPowerModel, key::Symbol, idx) = con(pm, pm.cnw, key, idx)
con(pm::GenericPowerModel, n::Int, key::Symbol) = pm.con[:nw][n][key]
con(pm::GenericPowerModel, n::Int, key::Symbol, idx) = pm.con[:nw][n][key][idx]

ext(pm::GenericPowerModel, key::Symbol) = ext(pm, pm.cnw, key)
ext(pm::GenericPowerModel, key::Symbol, idx) = ext(pm, pm.cnw, key, idx)
ext(pm::GenericPowerModel, n::Int, key::Symbol) = pm.ext[:nw][n][key]
ext(pm::GenericPowerModel, n::Int, key::Symbol, idx) = pm.ext[:nw][n][key][idx]

# TODO Ask Miles, why do we need to put JuMP. here?  using at top level should bring it in
function JuMP.setsolver(pm::GenericPowerModel, solver::MathProgBase.AbstractMathProgSolver)
    setsolver(pm.model, solver)
end

function JuMP.solve(pm::GenericPowerModel)
    status, solve_time, solve_bytes_alloc, sec_in_gc = @timed solve(pm.model)

    try
        solve_time = getsolvetime(pm.model)
    catch
        warn(LOGGER, "there was an issue with getsolvetime() on the solver, falling back on @timed.  This is not a rigorous timing value.");
    end

    return status, solve_time
end

""
function run_generic_model(file::String, model_constructor, solver, post_method; kwargs...)
    data = PowerModels.parse_file(file)
    return run_generic_model(data, model_constructor, solver, post_method; kwargs...)
end

""
function run_generic_model(data::Dict{String,Any}, model_constructor, solver, post_method; solution_builder = get_solution, kwargs...)
    pm = build_generic_model(data, model_constructor, post_method; kwargs...)

    solution = solve_generic_model(pm, solver; solution_builder = solution_builder)

    return solution
end

""
function build_generic_model(file::String,  model_constructor, post_method; kwargs...)
    data = PowerModels.parse_file(file)
    return build_generic_model(data, model_constructor, post_method; kwargs...)
end

""
function build_generic_model(data::Dict{String,Any}, model_constructor, post_method; multinetwork=false, kwargs...)
    # NOTE, this model constructor will build the ref dict using the latest info from the data
    pm = model_constructor(data; kwargs...)

    if !multinetwork && data["multinetwork"]
        warn(LOGGER, "building a single network model with multinetwork data, only network ($(pm.cnw)) will be used.")
    end

    post_method(pm)

    return pm
end

""
function solve_generic_model(pm::GenericPowerModel, solver; solution_builder = get_solution)
    setsolver(pm.model, solver)

    status, solve_time = solve(pm)

    return build_solution(pm, status, solve_time; solution_builder = solution_builder)
end

"""
Returns a dict that stores commonly used pre-computed data from of the data dictionary,
primarily for converting data-types, filtering out deactivated components, and storing
system-wide values that need to be computed globally.

Some of the common keys include:

* `:off_angmin` and `:off_angmax` (see `calc_theta_delta_bounds(data)`),
* `:bus` -- the set `{(i, bus) in ref[:bus] : bus["bus_type"] != 4}`,
* `:gen` -- the set `{(i, gen) in ref[:gen] : gen["gen_status"] == 1 && gen["gen_bus"] in keys(ref[:bus])}`,
* `:branch` -- the set of branches that are active in the network (based on the component status values),
* `:arcs_from` -- the set `[(i,b["f_bus"],b["t_bus"]) for (i,b) in ref[:branch]]`,
* `:arcs_to` -- the set `[(i,b["t_bus"],b["f_bus"]) for (i,b) in ref[:branch]]`,
* `:arcs` -- the set of arcs from both `arcs_from` and `arcs_to`,
* `:bus_arcs` -- the mapping `Dict(i => [(l,i,j) for (l,i,j) in ref[:arcs]])`,
* `:buspairs` -- (see `buspair_parameters(ref[:arcs_from], ref[:branch], ref[:bus])`),
* `:bus_gens` -- the mapping `Dict(i => [gen["gen_bus"] for (i,gen) in ref[:gen]])`.
* `:arcs_from_dc` -- the set `[(i,b["f_bus"],b["t_bus"]) for (i,b) in ref[:dcline]]`,
* `:arcs_to_dc` -- the set `[(i,b["t_bus"],b["f_bus"]) for (i,b) in ref[:dcline]]`,
* `:arcs_dc` -- the set of arcs from both `arcs_from_dc` and `arcs_to_dc`,
* `:bus_arcs_dc` -- the mapping `Dict(i => [(l,i,j) for (l,i,j) in ref[:arcs_dc]])`, and
* `:buspairs_dc` -- (see `buspair_parameters(ref[:arcs_from_dc], ref[:dcline], ref[:bus])`),

If `:ne_branch` exists, then the following keys are also available with similar semantics:

* `:ne_branch`, `:ne_arcs_from`, `:ne_arcs_to`, `:ne_arcs`, `:ne_bus_arcs`, `:ne_buspairs`.
"""
function build_ref(data::Dict{String,Any})
    refs = Dict{Symbol,Any}()
    nws = refs[:nw] = Dict{Int,Any}()

    if data["multinetwork"]
        nws_data = data["nw"]
    else
        nws_data = Dict{String,Any}("0" => data)
    end

    for (n,nw_data) in nws_data
        nw_id = parse(Int, n)
        ref = nws[nw_id] = Dict{Symbol,Any}()

        for (key, item) in nw_data
            if isa(item, Dict)
                item_lookup = Dict([(parse(Int, k), v) for (k,v) in item])
                ref[Symbol(key)] = item_lookup
            else
                ref[Symbol(key)] = item
            end
        end

        off_angmin, off_angmax = calc_theta_delta_bounds(nw_data)
        ref[:off_angmin] = off_angmin
        ref[:off_angmax] = off_angmax

        # filter turned off stuff
        ref[:bus] = filter((i, bus) -> bus["bus_type"] != 4, ref[:bus])
        ref[:gen] = filter((i, gen) -> gen["gen_status"] == 1 && gen["gen_bus"] in keys(ref[:bus]), ref[:gen])
        ref[:branch] = filter((i, branch) -> branch["br_status"] == 1 && branch["f_bus"] in keys(ref[:bus]) && branch["t_bus"] in keys(ref[:bus]), ref[:branch])
        ref[:dcline] = filter((i, dcline) -> dcline["br_status"] == 1 && dcline["f_bus"] in keys(ref[:bus]) && dcline["t_bus"] in keys(ref[:bus]), ref[:dcline])

        ref[:arcs_from] = [(i,branch["f_bus"],branch["t_bus"]) for (i,branch) in ref[:branch]]
        ref[:arcs_to]   = [(i,branch["t_bus"],branch["f_bus"]) for (i,branch) in ref[:branch]]
        ref[:arcs] = [ref[:arcs_from]; ref[:arcs_to]]

        ref[:arcs_from_dc] = [(i,dcline["f_bus"],dcline["t_bus"]) for (i,dcline) in ref[:dcline]]
        ref[:arcs_to_dc]   = [(i,dcline["t_bus"],dcline["f_bus"]) for (i,dcline) in ref[:dcline]]
        ref[:arcs_dc]      = [ref[:arcs_from_dc]; ref[:arcs_to_dc]]

        # maps dc line from and to parameters to arcs
        arcs_dc_param = ref[:arcs_dc_param] = Dict()
        for (l,i,j) in ref[:arcs_from_dc]
            arcs_dc_param[(l,i,j)] = Dict{String,Any}(
                "pmin" => ref[:dcline][l]["pminf"],
                "pmax" => ref[:dcline][l]["pmaxf"],
                "pref" => ref[:dcline][l]["pf"],
                "qmin" => ref[:dcline][l]["qminf"],
                "qmax" => ref[:dcline][l]["qmaxf"],
                "qref" => ref[:dcline][l]["qf"]
            )
            arcs_dc_param[(l,j,i)] = Dict{String,Any}(
                "pmin" => ref[:dcline][l]["pmint"],
                "pmax" => ref[:dcline][l]["pmaxt"],
                "pref" => ref[:dcline][l]["pt"],
                "qmin" => ref[:dcline][l]["qmint"],
                "qmax" => ref[:dcline][l]["qmaxt"],
                "qref" => ref[:dcline][l]["qt"]
            )
        end

        bus_gens = Dict([(i, []) for (i,bus) in ref[:bus]])
        for (i,gen) in ref[:gen]
            push!(bus_gens[gen["gen_bus"]], i)
        end
        ref[:bus_gens] = bus_gens

        bus_arcs = Dict([(i, []) for (i,bus) in ref[:bus]])
        for (l,i,j) in ref[:arcs]
            push!(bus_arcs[i], (l,i,j))
        end
        ref[:bus_arcs] = bus_arcs

        bus_arcs_dc = Dict([(i, []) for (i,bus) in ref[:bus]])
        for (l,i,j) in ref[:arcs_dc]
            push!(bus_arcs_dc[i], (l,i,j))
        end
        ref[:bus_arcs_dc] = bus_arcs_dc

        # a set of buses to support multiple connected components
        ref_buses = Dict()
        for (k,v) in ref[:bus]
            if v["bus_type"] == 3
                ref_buses[k] = v
            end
        end

        if length(ref_buses) == 0
            big_gen = biggest_generator(ref[:gen])
            gen_bus = big_gen["gen_bus"]
            ref_buses[gen_bus] = ref[:bus][gen_bus]
            warn(LOGGER, "no reference bus found, setting bus $(gen_bus) as reference based on generator $(big_gen["index"])")
        end

        if length(ref_buses) > 1
            warn(LOGGER, "multiple reference buses found, $(keys(ref_buses)), this can cause infeasibility if they are in the same connected component")
        end

        ref[:ref_buses] = ref_buses

        ref[:buspairs] = buspair_parameters(ref[:arcs_from], ref[:branch], ref[:bus])

        if haskey(ref, :ne_branch)
            ref[:ne_branch] = filter((i, branch) -> branch["br_status"] == 1 && branch["f_bus"] in keys(ref[:bus]) && branch["t_bus"] in keys(ref[:bus]), ref[:ne_branch])

            ref[:ne_arcs_from] = [(i,branch["f_bus"],branch["t_bus"]) for (i,branch) in ref[:ne_branch]]
            ref[:ne_arcs_to]   = [(i,branch["t_bus"],branch["f_bus"]) for (i,branch) in ref[:ne_branch]]
            ref[:ne_arcs] = [ref[:ne_arcs_from]; ref[:ne_arcs_to]]

            ne_bus_arcs = Dict([(i, []) for (i,bus) in ref[:bus]])
            for (l,i,j) in ref[:ne_arcs]
                push!(ne_bus_arcs[i], (l,i,j))
            end
            ref[:ne_bus_arcs] = ne_bus_arcs

            ref[:ne_buspairs] = buspair_parameters(ref[:ne_arcs_from], ref[:ne_branch], ref[:bus])
        end

    end

    return refs
end

"find the largest active generator in the network"
function biggest_generator(gens)
    biggest_gen = nothing
    biggest_value = -Inf
    for (k,gen) in gens
        if gen["pmax"] > biggest_value
            biggest_gen = gen
            biggest_value = gen["pmax"]
        end
    end
    assert(biggest_gen != nothing)
    return biggest_gen
end

"compute bus pair level structures"
function buspair_parameters(arcs_from, branches, buses)
    buspair_indexes = collect(Set([(i,j) for (l,i,j) in arcs_from]))

    bp_angmin = Dict([(bp, -Inf) for bp in buspair_indexes])
    bp_angmax = Dict([(bp, Inf) for bp in buspair_indexes])
    bp_branch = Dict([(bp, Inf) for bp in buspair_indexes])

    for (l,branch) in branches
        i = branch["f_bus"]
        j = branch["t_bus"]

        bp_angmin[(i,j)] = max(bp_angmin[(i,j)], branch["angmin"])
        bp_angmax[(i,j)] = min(bp_angmax[(i,j)], branch["angmax"])
        bp_branch[(i,j)] = min(bp_branch[(i,j)], l)
    end

    buspairs = Dict([((i,j), Dict(
        "branch"=>bp_branch[(i,j)],
        "angmin"=>bp_angmin[(i,j)],
        "angmax"=>bp_angmax[(i,j)],
        "rate_a"=>branches[bp_branch[(i,j)]]["rate_a"],
        "tap"=>branches[bp_branch[(i,j)]]["tap"],
        "vm_fr_min"=>buses[i]["vmin"],
        "vm_fr_max"=>buses[i]["vmax"],
        "vm_to_min"=>buses[j]["vmin"],
        "vm_to_max"=>buses[j]["vmax"]
        )) for (i,j) in buspair_indexes])

    return buspairs
end


text:
function checkrepeat(OptCpnt)
    Opt = Dict(
    "V"=>["Vx","Vz"],
    "Tii"=>["Txx","Tzz"])

    for Cpnt in ["V","Tii"]
        if Cpnt in OptCpnt
            filter!(x->!(x in Opt[Cpnt]),OptCpnt)
        end
    end
    return OptCpnt
end

#=== elastic ===#
#=== write Recordings into bin file ===#
function wDB(
    fid::IOStream,
    wf::elwf2d,
    BDnloc::Array{Int64,2},
    nwf::nelwf2d,
    OptCpnt::Array{String})

    RecLocLib = Dict(
    "P"=>(BDnloc[2,:]-1).*nwf.BDntpp[1] + BDnloc[1,:],
    "Vx"=>(BDnloc[2,:]-1).*nwf.BDnvx[1] + BDnloc[1,:],
    "Vz"=>(BDnloc[2,:]-1).*nwf.BDnvz[1] + BDnloc[1,:],
    "Txx"=>(BDnloc[2,:]-1).*nwf.BDntpp[1] + BDnloc[1,:],
    "Tzz"=>(BDnloc[2,:]-1).*nwf.BDntpp[1] + BDnloc[1,:],
    "Txz"=>(BDnloc[2,:]-1).*nwf.BDntxz[1] + BDnloc[1,:])
    #"Wy"=>(BDnloc[2,:]-1).*nwf.N_BDWy[1] + BDnloc[1,:])

    RecLib = Dict(
    "P"=>vec((wf.txx+wf.tzz)/2)[RecLocLib["P"]],
    "Vx"=>vec(wf.vx)[RecLocLib["Vx"]],
    "Vz"=>vec(wf.vz)[RecLocLib["Vz"]],
    "Txx"=>vec(wf.txx)[RecLocLib["Txx"]],
    "Txz"=>vec(wf.txz)[RecLocLib["Txz"]],
    "Tzz"=>vec(wf.tzz)[RecLocLib["Tzz"]],
    "V"=>[vec(wf.vx)[RecLocLib["Vx"]];
    vec(wf.vz)[RecLocLib["Vz"]]],
    "Tii"=>[vec(wf.txx)[RecLocLib["Txx"]];
    vec(wf.tzz)[RecLocLib["Tzz"]]])
    #"Wy"=>wf.VecBDWy)

    for OptC in OptCpnt
        write(fid,RecLib[OptC])
    end
end

#=== write Wavefields into bin file ===#
function wDB(
    fid::IOStream,
    wf::elwf2d,
    OptCpnt::Array{String},
    it::Int64,
    Slices::Union{UnitRange{Int64},StepRange{Int64,Int64},Array{Int64,1},Int64})

    WFLib = Dict(
    "P"=>vec((wf.txx+wf.tzz)/2),
    "Vx"=>vec(wf.vx),
    "Vz"=>vec(wf.vz),
    "Txx"=>vec(wf.txx),
    "Tzz"=>vec(wf.tzz),
    "Txz"=>vec(wf.txz),
    "V"=>[vec(wf.vx);vec(wf.vz)],
    "Tii"=>[vec(wf.txx);vec(wf.tzz)])
    #"Wy"=>wf.VecBDWy)
    if it in Slices
        for OptC in OptCpnt
            write(fid,WFLib[OptC])
        end
    end
end

#=== write both Recordings and Wavefields into bin file ===#
function wDB(
    fidRec::IOStream,
    fidWF::IOStream,
    wf::elwf2d,
    BDnloc::Array{Int64,2},
    nwf::nelwf2d,
    OptCpnt::Array{String},
    it::Int64,
    Slices::Union{UnitRange{Int64},StepRange{Int64,Int64},Array{Int64,1},Int64})

    RecLocLib = Dict(
    "P"=>(BDnloc[2,:]-1).*nwf.BDntpp[1] + BDnloc[1,:],
    "Vx"=>(BDnloc[2,:]-1).*nwf.BDnvx[1] + BDnloc[1,:],
    "Vz"=>(BDnloc[2,:]-1).*nwf.BDnvz[1] + BDnloc[1,:],
    "Txx"=>(BDnloc[2,:]-1).*nwf.BDntpp[1] + BDnloc[1,:],
    "Tzz"=>(BDnloc[2,:]-1).*nwf.BDntpp[1] + BDnloc[1,:],
    "Txz"=>(BDnloc[2,:]-1).*nwf.BDntxz[1] + BDnloc[1,:])
    #"Wy"=>(BDnloc[2,:]-1).*nwf.N_BDWy[1] + BDnloc[1,:])

    RecLib = Dict(
    "P"=>vec((wf.txx+wf.tzz)/2)[RecLocLib["P"]],
    "Vx"=>vec(wf.vx)[RecLocLib["Vx"]],
    "Vz"=>vec(wf.vz)[RecLocLib["Vz"]],
    "Txx"=>vec(wf.txx)[RecLocLib["Txx"]],
    "Txz"=>vec(wf.txz)[RecLocLib["Txz"]],
    "Tzz"=>vec(wf.tzz)[RecLocLib["Tzz"]],
    "V"=>[vec(wf.vx)[RecLocLib["Vx"]];
    vec(wf.vz)[RecLocLib["Vz"]]],
    "Tii"=>[vec(wf.txx)[RecLocLib["Txx"]];
    vec(wf.tzz)[RecLocLib["Tzz"]]])
    #"Wy"=>wf.VecBDWy)

    WFLib = Dict(
    "P"=>vec((wf.txx+wf.tzz)/2),
    "Vx"=>vec(wf.vx),
    "Vz"=>vec(wf.vz),
    "Txx"=>vec(wf.txx),
    "Tzz"=>vec(wf.tzz),
    "Txz"=>vec(wf.txz),
    "V"=>[vec(vx);vec(wf.vz)],
    "Tii"=>[vec(txx);vec(tzz)])
    #"Wy"=>wf.VecBDWy)

    for OptC in OptCpnt
        write(fidRec,RecLib[OptC])
        if it in Slices
            write(fid,WFLib[OptC])
        end
    end

end

#=== acoustic ===#
#=== write Recordings into bin file ===#
function wDB(
    fid::IOStream,
    wf::acwf2d,
    BDnloc::Array{Int64,2},
    nwf::nacwf2d,
    OptCpnt::Array{String})

    RecLocLib = Dict(
    "P"=>(BDnloc[2,:]-1).*nwf.BDntpp[1] + BDnloc[1,:],
    "Vx"=>(BDnloc[2,:]-1).*nwf.BDnvx[1] + BDnloc[1,:],
    "Vz"=>(BDnloc[2,:]-1).*nwf.BDnvz[1] + BDnloc[1,:],
    "Txx"=>(BDnloc[2,:]-1).*nwf.BDntpp[1] + BDnloc[1,:],
    "Tzz"=>(BDnloc[2,:]-1).*nwf.BDntpp[1] + BDnloc[1,:])
    #"Wy"=>(BDnloc[2,:]-1).*nwf.N_BDWy[1] + BDnloc[1,:])

    RecLib = Dict(
    "P"=>vec((wf.txx+wf.tzz)/2)[RecLocLib["P"]],
    "Vx"=>vec(wf.vx)[RecLocLib["Vx"]],
    "Vz"=>vec(wf.vz)[RecLocLib["Vz"]],
    "Txx"=>vec(wf.txx)[RecLocLib["Txx"]],
    "Tzz"=>vec(wf.tzz)[RecLocLib["Tzz"]],
    "V"=>[vec(wf.vx)[RecLocLib["Vx"]];
    vec(wf.vz)[RecLocLib["Vz"]]],
    "Tii"=>[vec(wf.txx)[RecLocLib["Txx"]];
    vec(wf.tzz)[RecLocLib["Tzz"]]])
    #"Wy"=>wf.VecBDWy)

    for OptC in OptCpnt
        write(fid,RecLib[OptC])
    end
end

#=== write Wavefields into bin file ===#
function wDB(
    fid::IOStream,
    wf::acwf2d,
    OptCpnt::Array{String},
    it::Int64,
    Slices::Union{UnitRange{Int64},StepRange{Int64,Int64},Array{Int64,1},Int64})

    WFLib = Dict(
    "P"=>vec((wf.txx+wf.tzz)/2),
    "Vx"=>vec(wf.vx),
    "Vz"=>vec(wf.vz),
    "Txx"=>vec(wf.txx),
    "Tzz"=>vec(wf.tzz),
    "V"=>[vec(wf.vx);vec(wf.vz)],
    "Tii"=>[vec(wf.txx);vec(wf.tzz)])
    #"Wy"=>wf.VecBDWy)
    if it in Slices
        for OptC in OptCpnt
            write(fid,WFLib[OptC])
        end
    end
end

#=== write both Recordings and Wavefields into bin file ===#
function wDB(
    fidRec::IOStream,
    fidWF::IOStream,
    wf::acwf2d,
    BDnloc::Array{Int64,2},
    nwf::nacwf2d,
    OptCpnt::Array{String},
    it::Int64,
    Slices::Union{UnitRange{Int64},StepRange{Int64,Int64},Array{Int64,1},Int64})

    RecLocLib = Dict(
    "P"=>(BDnloc[2,:]-1).*nwf.BDntpp[1] + BDnloc[1,:],
    "Vx"=>(BDnloc[2,:]-1).*nwf.BDnvx[1] + BDnloc[1,:],
    "Vz"=>(BDnloc[2,:]-1).*nwf.BDnvz[1] + BDnloc[1,:],
    "Txx"=>(BDnloc[2,:]-1).*nwf.BDntpp[1] + BDnloc[1,:],
    "Tzz"=>(BDnloc[2,:]-1).*nwf.BDntpp[1] + BDnloc[1,:])
    #"Wy"=>(BDnloc[2,:]-1).*nwf.N_BDWy[1] + BDnloc[1,:])

    RecLib = Dict(
    "P"=>vec((wf.txx+wf.tzz)/2)[RecLocLib["P"]],
    "Vx"=>vec(wf.vx)[RecLocLib["Vx"]],
    "Vz"=>vec(wf.vz)[RecLocLib["Vz"]],
    "Txx"=>vec(wf.txx)[RecLocLib["Txx"]],
    "Tzz"=>vec(wf.tzz)[RecLocLib["Tzz"]],
    "V"=>[vec(wf.vx)[RecLocLib["Vx"]];
    vec(wf.vz)[RecLocLib["Vz"]]],
    "Tii"=>[vec(wf.txx)[RecLocLib["Txx"]];
    vec(wf.tzz)[RecLocLib["Tzz"]]])
    #"Wy"=>wf.VecBDWy)

    WFLib = Dict(
    "P"=>vec((wf.txx+wf.tzz)/2),
    "Vx"=>vec(wf.vx),
    "Vz"=>vec(wf.vz),
    "Txx"=>vec(wf.txx),
    "Tzz"=>vec(wf.tzz),
    "V"=>[vec(vx);vec(wf.vz)],
    "Tii"=>[vec(txx);vec(tzz)])
    #"Wy"=>wf.VecBDWy)

    for OptC in OptCpnt
        write(fidRec,RecLib[OptC])
        if it in Slices
            write(fid,WFLib[OptC])
        end
    end

end


text:
# Adaptation of BayesLogit C code by Nicholas Polson, James Scott, Jesse Windle, 2012-2019

#------------------------------------------------------------------------------
# Constants
const _TRUNC = 0.64
const tol = 1e-8
const grid_size = 81
const ygrid = [
    0.0625,
    0.06698584,
    0.07179365,
    0.07694653,
    0.08246924,
    0.08838835,
    0.09473229,
    0.1015315,
    0.1088188,
    0.1166291,
    0.125,
    0.1339717,
    0.1435873,
    0.1538931,
    0.1649385,
    0.1767767,
    0.1894646,
    0.2030631,
    0.2176376,
    0.2332582,
    0.25,
    0.2679434,
    0.2871746,
    0.3077861,
    0.329877,
    0.3535534,
    0.3789291,
    0.4061262,
    0.4352753,
    0.4665165,
    0.5,
    0.5358867,
    0.5743492,
    0.6155722,
    0.659754,
    0.7071068,
    0.7578583,
    0.8122524,
    0.8705506,
    0.933033,
    1,
    1.071773,
    1.148698,
    1.231144,
    1.319508,
    1.414214,
    1.515717,
    1.624505,
    1.741101,
    1.866066,
    2,
    2.143547,
    2.297397,
    2.462289,
    2.639016,
    2.828427,
    3.031433,
    3.24901,
    3.482202,
    3.732132,
    4,
    4.287094,
    4.594793,
    4.924578,
    5.278032,
    5.656854,
    6.062866,
    6.498019,
    6.964405,
    7.464264,
    8,
    8.574188,
    9.189587,
    9.849155,
    10.55606,
    11.31371,
    12.12573,
    12.99604,
    13.92881,
    14.92853,
    16,
]

const vgrid = [
    -256,
    -222.8609,
    -194.0117,
    -168.897,
    -147.0334,
    -128,
    -111.4305,
    -97.00586,
    -84.4485,
    -73.51668,
    -63.99997,
    -55.71516,
    -48.50276,
    -42.22387,
    -36.75755,
    -31.99844,
    -27.85472,
    -24.24634,
    -21.10349,
    -18.36524,
    -15.97843,
    -13.89663,
    -12.07937,
    -10.49137,
    -9.101928,
    -7.884369,
    -6.815582,
    -5.875571,
    -5.047078,
    -4.315237,
    -3.667256,
    -3.092143,
    -2.580459,
    -2.124095,
    -1.716085,
    -1.350442,
    -1.022007,
    -0.7263359,
    -0.4595871,
    -0.2184366,
    0,
    0.1982309,
    0.3784427,
    0.5425468,
    0.6922181,
    0.828928,
    0.953973,
    1.068498,
    1.173516,
    1.269928,
    1.358533,
    1.440046,
    1.515105,
    1.584282,
    1.64809,
    1.706991,
    1.761401,
    1.811697,
    1.858218,
    1.901274,
    1.941143,
    1.978081,
    2.012318,
    2.044068,
    2.073521,
    2.100856,
    2.126234,
    2.149802,
    2.171696,
    2.192042,
    2.210954,
    2.228537,
    2.244889,
    2.260099,
    2.274249,
    2.287418,
    2.299673,
    2.311082,
    2.321703,
    2.331593,
    2.340804,
]

function y_eval(v::Float64)
    tol = 1e-6
    r = sqrt(abs(v))
    if (v > tol)
        return tan(r) / r
    elseif (v < -tol)
        return tanh(r) / r
    else
        #y = 1 + (1 / 3) * r^2 + (2 / 15) * r^4 + (17 / 315) * r^6
        return 1 + (1 / 3) * v + (2 / 15) * (v^2) + (17 / 315) * (v^3)
    end
end

function ydy_eval(v::Float64) #! put yp from second to first
    yp2 = y_eval(v)
    if (abs(v) >= tol)
        dyp2 = 0.5 * (yp2^2 + (1 - yp2) / v)
    else
        dyp2 = 0.5 * (yp2^2 - (1 / 3) - ((2 / 15) * v))
    end
    return yp2, dyp2
end

function f_eval(v::Float64, params::Float64)
    return y_eval(v) - params
end

function fdf_eval(v::Float64, params::Float64)
    fp2, dfp2 = ydy_eval(v)
    fp2 = fp2 - params
    return fp2, dfp2
end

function df_eval(v::Float64)
    f, df = ydy_eval(v)
    return df
end

function v_eval(y::Float64, tol::Float64, max_iter::Int64)
    ylower = ygrid[1]
    yupper = ygrid[grid_size]

    if (y < ylower)
        return -(1 / (y^2))
    elseif (y > yupper)
        v = atan(0.5 * y * __PI)
        return v^2
    elseif (y == 1)
        return 0.0
    end

    id = (_log_c(y) / _log_c(2.0) + 4.0) / 0.1

    idlow = Int64(floor(id) + 1)
    idhigh = Int64(min(idlow + 1, grid_size))
    vl = vgrid[idlow]  # lower bound
    vh = vgrid[idhigh] # upper bound

    iter = 0
    diff = tol + 1.0
    vnew = copy(vl)
    vold = copy(vl)

    while (diff > tol && iter <= max_iter)
        iter += 1
        vold = copy(vnew)
        f0, f1 = fdf_eval(vold, y)
        vnew = vold - f0 / f1
        vnew = vnew > vh ? vh : vnew
        vnew = vnew < vl ? vl : vnew
        diff = abs(vnew - vold)
    end

    if (iter > max_iter)
        println("v_eval: reached max_iter: ", iter)
    end

    return vnew
end


text:
export QBasedPolicy

using MacroTools: @forward
using Flux
using Setfield

"""
    QBasedPolicy(;learner::Q, explorer::S)

Use a Q-`learner` to generate estimations of action values.
Then an `explorer` is applied on the estimations to select an action.
"""
Base.@kwdef mutable struct QBasedPolicy{Q<:AbstractLearner,E<:AbstractExplorer} <:
                           AbstractPolicy
    learner::Q
    explorer::E
end

Flux.functor(x::QBasedPolicy) = (learner = x.learner,), y -> @set x.learner = y.learner

(π::QBasedPolicy)(env) = π(env, ActionStyle(env))
(π::QBasedPolicy)(env, ::MinimalActionSet) = get_actions(env)[π.explorer(π.learner(env))]
(π::QBasedPolicy)(env, ::FullActionSet) =
    get_actions(env)[π.explorer(π.learner(env), get_legal_actions_mask(env))]

RLBase.get_prob(p::QBasedPolicy, env) = get_prob(p, env, ActionStyle(env))
RLBase.get_prob(p::QBasedPolicy, env, ::MinimalActionSet) =
    get_prob(p.explorer, p.learner(env))
RLBase.get_prob(p::QBasedPolicy, env, ::FullActionSet) =
    get_prob(p.explorer, p.learner(env), get_legal_actions_mask(env))

@forward QBasedPolicy.learner RLBase.get_priority, RLBase.update!

function Flux.testmode!(p::QBasedPolicy, mode = true)
    testmode!(p.learner, mode)
    testmode!(p.explorer, mode)
end


text:
module FastConv

# package code goes here
using Base.Cartesian

include("./utils.jl");

end # module


text:
using DataFrames, StatsBase, StatsPlots, CSV, CovidNYT, Dates
gr(size=(700, 1000))

ProjDir = @__DIR__

function ma(x::Vector{Float64}, wind::Int)
    len = length(x)
    y = Vector{Float64}(undef, len)
    for i in 1:len
        if i < wind
            y[i] = 0.0
        else
            y[i] = mean(x[i - wind + 1:i])
        end
    end
    return y
end

df_states = CSV.read(rel_path_covidnyt("..", "data", "us-states.csv"), DataFrame)
df_states[!, :month] = Dates.month.(df_states[:, :date])
df_states[!, :day] = Dates.day.(df_states[:, :date])
#first(df_states, 5) |> display

states = ["New York", "Florida", "California", "Colorado"]
counties = ["New York City", "San Diego", "District of Columbia", "Denver"]

wind = 7

p = Vector{Plots.Plot{Plots.GRBackend}}(undef, 2*length(states))
df = Vector{DataFrame}(undef, length(states))

plot_indx = 0
for (indx, state) in enumerate(states)
  global plot_indx += 1
  df[indx] = df_states[df_states.state .== state, :]

  p[plot_indx] = plot(leg=:topleft, title=state, ylab="log10")
  plot!((df[indx][:, :date]), log10.(df[indx][:, :cases]), lab="cases")
  plot!(df[indx][:, :date], log10.(df[indx][:, :deaths]), lab="deaths")

  local c = df[indx][:, :cases]
  c[2:end] = c[2:end] - c[1:end-1]
  df[indx][!, :new_cases] = c
  local d = df[indx][:, :deaths]
  d[2:end] = d[2:end] - d[1:end-1]
  df[indx][!, :new_deaths] = d

  df[indx][!, :ma] = ma(Float64.(df[indx][:, :new_cases]), wind)

  plot_indx += 1
  p[plot_indx] = plot(leg=:topleft, title=state, ylab="new cases", )
  plot!(df[indx][:, :date], df[indx][:, :new_cases], lab="New cases", color=:lightblue)
  plot!(df[indx][:, :date], df[indx][:, :new_deaths], lab="New deaths")
  yminmax = (0.0, maximum(df[indx][:, :ma]))
  plot!(df[indx][:, :date], df[indx][:, :ma], lab="$(wind) day ma", color=:darkblue, ylims=yminmax)

end

plot(p..., layout=(4, 2))
savefig("$(ProjDir)/graphs/s_$(string(today())).png")
#savefig("$(ProjDir)/graphs/s_$(string(today())).pdf")

df_counties = CSV.read(rel_path_covidnyt("..", "data", "us-counties.csv"), DataFrame)
df_counties[!, :month] = Dates.month.(df_counties[:, :date])
df_counties[!, :day] = Dates.day.(df_counties[:, :date])
#first(df_counties, 5) |> display
max_values = [7000, 5000, 800, 800]

plot_indx = 0
for (indx, county) in enumerate(counties)
  global plot_indx += 1
  df[indx] = df_counties[df_counties.county .== county, :]

  p[plot_indx] = plot(leg=:topleft, title=county, ylab="log10")
  plot!((df[indx][:, :date]), log10.(df[indx][:, :cases]), lab="cases")
  plot!(df[indx][:, :date], log10.(df[indx][:, :deaths]), lab="deaths")

  local c = df[indx][:, :cases]
  c[2:end] = c[2:end] - c[1:end-1]
  df[indx][!, :new_cases] = c
  local d = df[indx][:, :deaths]
  d[2:end] = d[2:end] - d[1:end-1]
  df[indx][!, :new_deaths] = d

  df[indx][!, :ma] = ma(Float64.(df[indx][:, :new_cases]), wind)

  plot_indx += 1
  yminmax = (0.0, max_values[indx])
  #yminmax = (0.0, maximum(df[indx][:, :new_cases]))
  p[plot_indx] = plot(leg=:topleft, title=county, ylims=yminmax, ylab="new cases")
  plot!(df[indx][:, :date], df[indx][:, :new_cases], lab="New cases", color=:lightblue)
  plot!(df[indx][:, :date], df[indx][:, :new_deaths], lab="New deaths")
  plot!(df[indx][:, :date], df[indx][:, :ma], lab="$(wind) day ma", color=:darkblue)

end

plot(p..., layout=(4, 2))
savefig("$(ProjDir)/graphs/c_$(string(today())).png")
#savefig("$(ProjDir)/graphs/c_$(string(today())).pdf")