text:
# This file is a part of Julia. License is MIT: http://julialang.org/license

# timing

# time() in libc.jl

# high-resolution relative time, in nanoseconds

"""
    time_ns()

Get the time in nanoseconds. The time corresponding to 0 is undefined, and wraps every 5.8 years.
"""
time_ns() = ccall(:jl_hrtime, UInt64, ())

# This type must be kept in sync with the C struct in src/gc.h
immutable GC_Num
    allocd      ::Int64 # GC internal
    deferred_alloc::Int64 # GC internal
    freed       ::Int64 # GC internal
    malloc      ::UInt64
    realloc     ::UInt64
    poolalloc   ::UInt64
    bigalloc    ::UInt64
    freecall    ::UInt64
    total_time  ::UInt64
    total_allocd::UInt64 # GC internal
    since_sweep ::UInt64 # GC internal
    collect     ::Csize_t # GC internal
    pause       ::Cint
    full_sweep  ::Cint
end

gc_num() = ccall(:jl_gc_num, GC_Num, ())

# This type is to represent differences in the counters, so fields may be negative
immutable GC_Diff
    allocd      ::Int64 # Bytes allocated
    malloc      ::Int64 # Number of GC aware malloc()
    realloc     ::Int64 # Number of GC aware realloc()
    poolalloc   ::Int64 # Number of pool allocation
    bigalloc    ::Int64 # Number of big (non-pool) allocation
    freecall    ::Int64 # Number of GC aware free()
    total_time  ::Int64 # Time spent in garbage collection
    pause       ::Int64 # Number of GC pauses
    full_sweep  ::Int64 # Number of GC full collection
end

gc_total_bytes(gc_num::GC_Num) =
    (gc_num.allocd + gc_num.deferred_alloc +
     Int64(gc_num.collect) + Int64(gc_num.total_allocd))

function GC_Diff(new::GC_Num, old::GC_Num)
    # logic from `src/gc.c:jl_gc_total_bytes`
    old_allocd = gc_total_bytes(old)
    new_allocd = gc_total_bytes(new)
    return GC_Diff(new_allocd - old_allocd,
                   Int64(new.malloc       - old.malloc),
                   Int64(new.realloc      - old.realloc),
                   Int64(new.poolalloc    - old.poolalloc),
                   Int64(new.bigalloc     - old.bigalloc),
                   Int64(new.freecall     - old.freecall),
                   Int64(new.total_time   - old.total_time),
                   new.pause              - old.pause,
                   new.full_sweep         - old.full_sweep)
end

function gc_alloc_count(diff::GC_Diff)
    diff.malloc + diff.realloc + diff.poolalloc + diff.bigalloc
end

# total time spend in garbage collection, in nanoseconds
gc_time_ns() = ccall(:jl_gc_total_hrtime, UInt64, ())

# total number of bytes allocated so far
gc_bytes() = ccall(:jl_gc_total_bytes, Int64, ())

"""
    tic()

Set a timer to be read by the next call to [`toc`](:func:`toc`) or [`toq`](:func:`toq`). The
macro call `@time expr` can also be used to time evaluation.
"""
function tic()
    t0 = time_ns()
    task_local_storage(:TIMERS, (t0, get(task_local_storage(), :TIMERS, ())))
    return t0
end

"""
    toq()

Return, but do not print, the time elapsed since the last [`tic`](:func:`tic`). The
macro calls `@timed expr` and `@elapsed expr` also return evaluation time.
"""
function toq()
    t1 = time_ns()
    timers = get(task_local_storage(), :TIMERS, ())
    if timers === ()
        error("toc() without tic()")
    end
    t0 = timers[1]::UInt64
    task_local_storage(:TIMERS, timers[2])
    (t1-t0)/1e9
end

"""
    toc()

Print and return the time elapsed since the last [`tic`](:func:`tic`). The macro call
`@time expr` can also be used to time evaluation.
"""
function toc()
    t = toq()
    println("elapsed time: ", t, " seconds")
    return t
end

# print elapsed time, return expression value
const _mem_units = ["byte", "KB", "MB", "GB", "TB", "PB"]
const _cnt_units = ["", " k", " M", " G", " T", " P"]
function prettyprint_getunits(value, numunits, factor)
    if value == 0 || value == 1
        return (value, 1)
    end
    unit = ceil(Int, log(value) / log(factor))
    unit = min(numunits, unit)
    number = value/factor^(unit-1)
    return number, unit
end

function padded_nonzero_print(value,str)
    if value != 0
        blanks = "                "[1:18-length(str)]
        println("$str:$blanks$value")
    end
end

function time_print(elapsedtime, bytes, gctime, allocs)
    @printf("%10.6f seconds", elapsedtime/1e9)
    if bytes != 0 || allocs != 0
        bytes, mb = prettyprint_getunits(bytes, length(_mem_units), Int64(1024))
        allocs, ma = prettyprint_getunits(allocs, length(_cnt_units), Int64(1000))
        if ma == 1
            @printf(" (%d%s allocation%s: ", allocs, _cnt_units[ma], allocs==1 ? "" : "s")
        else
            @printf(" (%.2f%s allocations: ", allocs, _cnt_units[ma])
        end
        if mb == 1
            @printf("%d %s%s", bytes, _mem_units[mb], bytes==1 ? "" : "s")
        else
            @printf("%.3f %s", bytes, _mem_units[mb])
        end
        if gctime > 0
            @printf(", %.2f%% gc time", 100*gctime/elapsedtime)
        end
        print(")")
    elseif gctime > 0
        @printf(", %.2f%% gc time", 100*gctime/elapsedtime)
    end
    println()
end

function timev_print(elapsedtime, diff::GC_Diff)
    allocs = gc_alloc_count(diff)
    time_print(elapsedtime, diff.allocd, diff.total_time, allocs)
    print("elapsed time (ns): $elapsedtime\n")
    padded_nonzero_print(diff.total_time,   "gc time (ns)")
    padded_nonzero_print(diff.allocd,       "bytes allocated")
    padded_nonzero_print(diff.poolalloc,    "pool allocs")
    padded_nonzero_print(diff.bigalloc,     "non-pool GC allocs")
    padded_nonzero_print(diff.malloc,       "malloc() calls")
    padded_nonzero_print(diff.realloc,      "realloc() calls")
    padded_nonzero_print(diff.freecall,     "free() calls")
    padded_nonzero_print(diff.pause,        "GC pauses")
    padded_nonzero_print(diff.full_sweep,   "full collections")
end

"""
    @time

A macro to execute an expression, printing the time it took to execute, the number of
allocations, and the total number of bytes its execution caused to be allocated, before
returning the value of the expression.

See also [`@timev`](:func:`@timev`), [`@timed`](:func:`@timed`), [`@elapsed`](:func:`@elapsed`), and
[`@allocated`](:func:`@allocated`).
"""
macro time(ex)
    quote
        local stats = gc_num()
        local elapsedtime = time_ns()
        local val = $(esc(ex))
        elapsedtime = time_ns() - elapsedtime
        local diff = GC_Diff(gc_num(), stats)
        time_print(elapsedtime, diff.allocd, diff.total_time,
                   gc_alloc_count(diff))
        val
    end
end

"""
    @timev

This is a verbose version of the `@time` macro. It first prints the same information as
`@time`, then any non-zero memory allocation counters, and then returns the value of the
expression.

See also [`@time`](:func:`@time`), [`@timed`](:func:`@timed`), [`@elapsed`](:func:`@elapsed`), and
[`@allocated`](:func:`@allocated`).
"""
macro timev(ex)
    quote
        local stats = gc_num()
        local elapsedtime = time_ns()
        local val = $(esc(ex))
        elapsedtime = time_ns() - elapsedtime
        timev_print(elapsedtime, GC_Diff(gc_num(), stats))
        val
    end
end

"""
    @elapsed

A macro to evaluate an expression, discarding the resulting value, instead returning the
number of seconds it took to execute as a floating-point number.

See also [`@time`](:func:`@time`), [`@timev`](:func:`@timev`), [`@timed`](:func:`@timed`),
and [`@allocated`](:func:`@allocated`).
"""
macro elapsed(ex)
    quote
        local t0 = time_ns()
        local val = $(esc(ex))
        (time_ns()-t0)/1e9
    end
end

# measure bytes allocated without *most* contamination from compilation
# Note: This reports a different value from the @time macros, because
# it wraps the call in a function, however, this means that things
# like:  @allocated y = foo()
# will not work correctly, because it will set y in the context of
# the local function made by the macro, not the current function
"""
    @allocated

A macro to evaluate an expression, discarding the resulting value, instead returning the
total number of bytes allocated during evaluation of the expression. Note: the expression is
evaluated inside a local function, instead of the current context, in order to eliminate the
effects of compilation, however, there still may be some allocations due to JIT compilation.
This also makes the results inconsistent with the `@time` macros, which do not try to adjust
for the effects of compilation.

See also [`@time`](:func:`@time`), [`@timev`](:func:`@timev`), [`@timed`](:func:`@timed`),
and [`@elapsed`](:func:`@elapsed`).
"""
macro allocated(ex)
    quote
        let
            local f
            function f()
                b0 = gc_bytes()
                $(esc(ex))
                gc_bytes() - b0
            end
            f()
        end
    end
end

"""
    @timed

A macro to execute an expression, and return the value of the expression, elapsed time,
total bytes allocated, garbage collection time, and an object with various memory allocation
counters.

See also [`@time`](:func:`@time`), [`@timev`](:func:`@timev`), [`@elapsed`](:func:`@elapsed`), and
[`@allocated`](:func:`@allocated`).
"""
macro timed(ex)
    quote
        local stats = gc_num()
        local elapsedtime = time_ns()
        local val = $(esc(ex))
        elapsedtime = time_ns() - elapsedtime
        local diff = GC_Diff(gc_num(), stats)
        val, elapsedtime/1e9, diff.allocd, diff.total_time/1e9, diff
    end
end

function fftw_vendor()
    if Base.libfftw_name == "libmkl_rt"
        return :mkl
    else
        return :fftw
    end
end

## printing with color ##

function with_output_color(f::Function, color::Union{Int, Symbol}, io::IO, args...)
    buf = IOBuffer()
    have_color && print(buf, get(text_colors, color, color_normal))
    try f(IOContext(buf, io), args...)
    finally
        have_color && print(buf, color_normal)
        print(io, String(take!(buf)))
    end
end

"""
    print_with_color(color::Union{Symbol, Int}, [io], strings...)

Print strings in a color specified as a symbol.

`color` may take any of the values $(Base.available_text_colors_docstring)
or an integer between 0 and 255 inclusive. Note that not all terminals support 256 colors.
"""
print_with_color(color::Union{Int, Symbol}, io::IO, msg::AbstractString...) =
    with_output_color(print, color, io, msg...)
print_with_color(color::Union{Int, Symbol}, msg::AbstractString...) =
    print_with_color(color, STDOUT, msg...)
println_with_color(color::Union{Int, Symbol}, io::IO, msg::AbstractString...) =
    with_output_color(println, color, io, msg...)
println_with_color(color::Union{Int, Symbol}, msg::AbstractString...) =
    println_with_color(color, STDOUT, msg...)

## warnings and messages ##

"""
    info(msg...; prefix="INFO: ")

Display an informational message.
Argument `msg` is a string describing the information to be displayed.
The `prefix` keyword argument can be used to override the default
prepending of `msg`.

```jldoctest
julia> info("hello world")
INFO: hello world

julia> info("hello world"; prefix="MY INFO: ")
MY INFO: hello world
```
"""
function info(io::IO, msg...; prefix="INFO: ")
    println_with_color(info_color(), io, prefix, chomp(string(msg...)))
end
info(msg...; prefix="INFO: ") = info(STDERR, msg..., prefix=prefix)

# print a warning only once

const have_warned = Set()

warn_once(io::IO, msg...) = warn(io, msg..., once=true)
warn_once(msg...) = warn(STDERR, msg..., once=true)

function warn(io::IO, msg...;
              prefix="WARNING: ", once=false, key=nothing, bt=nothing,
              filename=nothing, lineno::Int=0)
    str = chomp(string(msg...))
    if once
        if key === nothing
            key = str
        end
        (key in have_warned) && return
        push!(have_warned, key)
    end
    print_with_color(warn_color(), io, prefix, str)
    if bt !== nothing
        show_backtrace(io, bt)
    end
    if filename !== nothing
        print(io, "\nwhile loading $filename, in expression starting on line $lineno")
    end
    println(io)
    return
end

"""
    warn(msg)

Display a warning. Argument `msg` is a string describing the warning to be displayed.
"""
warn(msg...; kw...) = warn(STDERR, msg...; kw...)

warn(io::IO, err::Exception; prefix="ERROR: ", kw...) =
    warn(io, sprint(buf->showerror(buf, err)), prefix=prefix; kw...)

warn(err::Exception; prefix="ERROR: ", kw...) =
    warn(STDERR, err, prefix=prefix; kw...)

info(io::IO, err::Exception; prefix="ERROR: ", kw...) =
    info(io, sprint(buf->showerror(buf, err)), prefix=prefix; kw...)

info(err::Exception; prefix="ERROR: ", kw...) =
    info(STDERR, err, prefix=prefix; kw...)

function julia_cmd(julia=joinpath(JULIA_HOME, julia_exename()))
    opts = JLOptions()
    cpu_target = unsafe_string(opts.cpu_target)
    image_file = unsafe_string(opts.image_file)
    compile = if opts.compile_enabled == 0
                  "no"
              elseif opts.compile_enabled == 2
                  "all"
              elseif opts.compile_enabled == 3
                  "min"
              else
                  "yes"
              end
    depwarn = if opts.depwarn == 0
                  "no"
              elseif opts.depwarn == 2
                  "error"
              else
                  "yes"
              end
    `$julia -C$cpu_target -J$image_file --compile=$compile --depwarn=$depwarn`
end

julia_exename() = ccall(:jl_is_debugbuild,Cint,())==0 ? "julia" : "julia-debug"

"""
    securezero!(o)

`securezero!` fills the memory associated with an object `o` with zeros.
Unlike `fill!(o,0)` and similar code, which might be optimized away by
the compiler for objects about to be discarded, the `securezero!` function
will always be called.
"""
function securezero! end
@noinline securezero!{T<:Number}(a::AbstractArray{T}) = fill!(a, 0)
securezero!(s::String) = securezero!(s.data)
@noinline unsafe_securezero!{T}(p::Ptr{T}, len::Integer=1) =
    ccall(:memset, Ptr{T}, (Ptr{T}, Cint, Csize_t), p, 0, len*sizeof(T))
unsafe_securezero!(p::Ptr{Void}, len::Integer=1) = Ptr{Void}(unsafe_securezero!(Ptr{UInt8}(p), len))

if is_windows()
function getpass(prompt::AbstractString)
    print(prompt)
    flush(STDOUT)
    p = Array{UInt8}(128) # mimic Unix getpass in ignoring more than 128-char passwords
                          # (also avoids any potential memory copies arising from push!)
    try
        plen = 0
        while true
            c = ccall(:_getch, UInt8, ())
            if c == 0xff || c == UInt8('\n') || c == UInt8('\r')
                break # EOF or return
            elseif c == 0x00 || c == 0xe0
                ccall(:_getch, UInt8, ()) # ignore function/arrow keys
            elseif c == UInt8('\b') && plen > 0
                plen -= 1 # delete last character on backspace
            elseif !iscntrl(Char(c)) && plen < 128
                p[plen += 1] = c
            end
        end
        return unsafe_string(pointer(p), plen) # use unsafe_string rather than String(p[1:plen])
                                               # to be absolutely certain we never make an extra copy
    finally
        securezero!(p)
    end

    return ""
end
else
getpass(prompt::AbstractString) = unsafe_string(ccall(:getpass, Cstring, (Cstring,), prompt))
end

# Windows authentication prompt
if is_windows()
    immutable CREDUI_INFO
        cbSize::UInt32
        parent::Ptr{Void}
        pszMessageText::Ptr{UInt16}
        pszCaptionText::Ptr{UInt16}
        banner::Ptr{Void}
    end

    const CREDUIWIN_GENERIC                 = 0x0001
    const CREDUIWIN_IN_CRED_ONLY            = 0x0020
    const CREDUIWIN_ENUMERATE_CURRENT_USER  = 0x0200

    const CRED_PACK_GENERIC_CREDENTIALS     = 0x0004

    const ERROR_SUCCESS                     = 0x0000
    const ERROR_CANCELLED                   = 0x04c7

    function winprompt(message, caption, default_username; prompt_username = true)
        # Step 1: Create an encrypted username/password bundle that will be used to set
        #         the default username (in theory could also provide a default password)
        credbuf = Array{UInt8,1}(1024)
        credbufsize = Ref{UInt32}(sizeof(credbuf))
        succeeded = ccall((:CredPackAuthenticationBufferW, "credui.dll"), stdcall, Bool,
            (UInt32, Cwstring, Cwstring, Ptr{UInt8}, Ptr{UInt32}),
             CRED_PACK_GENERIC_CREDENTIALS, default_username, "", credbuf, credbufsize)
        @assert succeeded

        # Step 2: Create the actual dialog
        #      2.1: Set up the window
        messageArr = Base.cwstring(message)
        captionArr = Base.cwstring(caption)
        pfSave = Ref{Bool}(false)
        cred = Ref{CREDUI_INFO}(CREDUI_INFO(sizeof(CREDUI_INFO), C_NULL, pointer(messageArr), pointer(captionArr), C_NULL))
        dwflags = CREDUIWIN_GENERIC | CREDUIWIN_ENUMERATE_CURRENT_USER
        if !prompt_username
            # Disable setting anything other than default_username
            dwflags |= CREDUIWIN_IN_CRED_ONLY
        end
        authPackage = Ref{Culong}(0)
        outbuf_data = Ref{Ptr{Void}}(C_NULL)
        outbuf_size = Ref{Culong}(0)

        #      2.2: Do the actual request
        code = ccall((:CredUIPromptForWindowsCredentialsW, "credui.dll"), stdcall, UInt32, (Ptr{CREDUI_INFO}, UInt32, Ptr{Culong},
            Ptr{Void}, Culong, Ptr{Ptr{Void}}, Ptr{Culong}, Ptr{Bool}, UInt32), cred, 0, authPackage, credbuf, credbufsize[],
            outbuf_data, outbuf_size, pfSave, dwflags)

        #      2.3: If that failed for any reason other than the user canceling, error out.
        #           If the user canceled, just return a nullable
        if code == ERROR_CANCELLED
            return Nullable{Tuple{String,String}}()
        elseif code != ERROR_SUCCESS
            error(Base.Libc.FormatMessage(code))
        end

        # Step 3: Convert encrypted credentials back to plain text
        passbuf = Array{UInt16,1}(1024)
        passlen = Ref{UInt32}(length(passbuf))
        usernamebuf = Array{UInt16,1}(1024)
        usernamelen = Ref{UInt32}(length(usernamebuf))
        # Need valid buffers for domain, even though we don't care
        dummybuf = Array{UInt16,1}(1024)
        succeeded = ccall((:CredUnPackAuthenticationBufferW, "credui.dll"), Bool,
            (UInt32, Ptr{Void}, UInt32, Ptr{UInt16}, Ptr{UInt32}, Ptr{UInt16}, Ptr{UInt32}, Ptr{UInt16}, Ptr{UInt32}),
            0, outbuf_data[], outbuf_size[], usernamebuf, usernamelen, dummybuf, Ref{UInt32}(1024), passbuf, passlen)
        if !succeeded
            error(Base.Libc.FormatMessage())
        end

        # Step 4: Free the encrypted buffer
        # ccall(:SecureZeroMemory, Ptr{Void}, (Ptr{Void}, Csize_t), outbuf_data[], outbuf_size[]) - not an actual function
        unsafe_securezero!(outbuf_data[], outbuf_size[])
        ccall((:CoTaskMemFree, "ole32.dll"), Void, (Ptr{Void},), outbuf_data[])

        # Done.
        passbuf_ = passbuf[1:passlen[]-1]
        result = Nullable((String(transcode(UInt8, usernamebuf[1:usernamelen[]-1])),
            String(transcode(UInt8, passbuf_))))
        securezero!(passbuf_)
        securezero!(passbuf)

        return result
    end

end

"""
    crc32c(data, crc::UInt32=0x00000000)
Compute the CRC-32c checksum of the given `data`, which can be
an `Array{UInt8}` or a `String`.  Optionally, you can pass
a starting `crc` integer to be mixed in with the checksum.
(Technically, a little-endian checksum is computed.)
"""
function crc32c end
crc32c(a::Array{UInt8}, crc::UInt32=0x00000000) = ccall(:jl_crc32c, UInt32, (UInt32, Ptr{UInt8}, Csize_t), crc, a, sizeof(a))
crc32c(s::String, crc::UInt32=0x00000000) = crc32c(s.data, crc)

max_stars_repo_name: habemus-papadum/julia

max_issues_repo_name: habemus-papadum/julia

max_forks_repo_name: habemus-papadum/julia


text:
using DocumentFormat
using Documenter

makedocs(;
    modules=[DocumentFormat],
    authors="Julia VSCode",
    repo="https://github.com/julia-vscode/DocumentFormat.jl/blob/{commit}{path}#L{line}",
    sitename="DocumentFormat.jl",
    format=Documenter.HTML(;
        prettyurls = prettyurls = get(ENV, "CI", nothing) == "true",
        # canonical="https://www.julia-vscode.org/DocumentFormat.jl",
        # assets=String[],
    ),
    pages=[
        "Home" => "index.md",
        "Syntax Reference" => "syntax.md",
    ],
)

deploydocs(;
    repo="github.com/julia-vscode/DocumentFormat.jl",
)

ext: jl


text:
#===============================================================================
  Preconditioning

  Todo: maybe add an inverse?
===============================================================================#

struct NoPrecon end
struct HasPrecon end
function initial_preconditioner(method, x, ::NoPrecon)
    nothing
end
function initial_preconditioner(method, x, ::HasPrecon)
    method.P(x)
end
update_preconditioner(method, x, P) = update_preconditioner(method, x, P, hasprecon(method))
update_preconditioner(method, x, P, ::NoPrecon) = nothing
update_preconditioner(method, x, P::Nothing, ::HasPrecon) = method.P(x)
update_preconditioner(method, x, P, ::HasPrecon) = method.P(x, P)
apply_preconditioner(mstyle::InPlace, ::Nothing, Pg, ∇f) = copyto!(Pg, ∇f)
apply_preconditioner(mstyle::InPlace, P, Pg, ∇f) = mul!(Pg, P, ∇f)
apply_preconditioner(mstyle::OutOfPlace, ::Nothing, Pg, ∇f) = ∇f
apply_preconditioner(mstyle::OutOfPlace, P, Pg, ∇f) = P * ∇f
apply_inverse_preconditioner(mstyle::InPlace, ::Nothing, Pg, ∇f) = copyto!(Pg, ∇f)
apply_inverse_preconditioner(mstyle::InPlace, P, Pg, ∇f) = ldiv!(Pg, factorize(P), ∇f) # add special for Diagonal?
function apply_inverse_preconditioner(mstyle::InPlace, P::Diagonal, Pg, ∇f)
    Pg .= P.diag .\ ∇f
end
apply_inverse_preconditioner(mstyle::OutOfPlace, ::Nothing, Pg, ∇f) = ∇f
apply_inverse_preconditioner(mstyle::OutOfPlace, P, pg, ∇f) = P \ ∇f

max_stars_repo_name: JuliaNLSolvers/NLSolvers

max_issues_repo_name: JuliaNLSolvers/NLSolvers

max_forks_repo_name: JuliaNLSolvers/NLSolvers


text:
module Exchange

#
# Need these packages to run the code
#
using CompEcon: brent, gridmake, qnwnorm
using Interpolations
using JLD

#
# Abstract and Typle alias declarations
#
typealias VFITP Tuple{ScaledInterpolation, ScaledInterpolation}

# Abstract Model
abstract AbstractModel

Base.writemime{TM<:AbstractModel}(io::IO, ::MIME"text/plain", m::TM) =
    print(io, "The model of type $TM")

"Identifies how many exogenous variables the model has"
n_exog(m::AbstractModel) = n_exog(m.exog)

"Evaluates certainty equivalents for agent 1"
eval_lμ1(m::AbstractModel, gp::Vector, jp::Vector) =
    log(_CE(m.agent1, gp, jp, m.Π))

"Evaluates certainty equivalents for agent 1"
eval_lμ2(m::AbstractModel, gp::Vector, up::Vector) =
    log(_CE(m.agent2, gp, up, m.Π))

"Evaluate value function for agent 1 given consumption and cert equiv"
eval_lj(m::AbstractModel, c1, μ1) = _agg_lv(m.agent1, c1, μ1)

"Evaluate value function for agent 2 given consumption and cert equiv"
eval_lu(m::AbstractModel, c2, μ2) = _agg_lv(m.agent2, c2, μ2)

"This holds allocation values at each state on our grid"
type OtherItems
    a1::Vector{Float64}
    b2::Vector{Float64}
    c1::Vector{Float64}
    c2::Vector{Float64}
end

# Abstract Exog
abstract AbstractExog

"""
Simulates the model's exogenous variables. Needs an `exog` type that
describes the exogenous processes, `init` which is an initial vector of
levels for exogenous variables, and `ϵ` which are the standard normal
i.i.d. shocks.

Note: In order for this function to work, a method named `exog_step` must
be defined for the `exog::AbstractExog` type that is passed in.
"""
function simulate(exog::AbstractExog, init::Vector, ϵ::Matrix)
    capT = size(ϵ, 2)

    # allocate
    out = similar(ϵ)
    out[:, 1] = init

    N = length(init)

    # step forward
    for t in 2:capT
        next = exog_step(exog, out[:, t-1], ϵ[:, t])
        for i in 1:N
            out[i, t] = next[i]
        end
    end
    out
end

"""
Simulates exogenous variables. Calls other `simulate` method but
initializes the state of exogenous variables.
"""
simulate(exog::AbstractExog, ϵ::Matrix) = simulate(exog, init(exog), ϵ)

"""
Given a model, construct an intitial guess for the value functions on the grid.

If the keyword arugment `infile` it should point to a `.jld` file containing at
least three top level (not in HDF5 groups) variables: `lj`, `lu`, `m_ls`
specifying log J, log U and linspaces forming the grid over (lλ☆, lzh, v);
respectively.

If the file is not valid or not given, then the guess that `lj = c1` and `lu=c2`
is used.
"""
function initialize_vf(m::AbstractModel;
                       infile="solutions/$(jld_prefix(m))_1.jld")

    # Pull out useful data
    Ngrid, Nst = size(m.grid)

    # Check for the file
    try  # try to use infile. If it fails, just set l, j to c
        # Read in solution data
        # Read in solution data
        f = jldopen(infile, "r")
        lj = read(f, "lj")
        lu = read(f, "lu")
        m_ls = read(f, "m_ls")
        close(f)

        # Create interpolants
        lj_itp = scale(interpolate(lj, BSpline(Linear()), OnGrid()), m.m_ls...)
        lu_itp = scale(interpolate(lu, BSpline(Linear()), OnGrid()), m.m_ls...)

        # Allocate space to fill
        lj_grid, lu_grid = Array(Float64, Ngrid), Array(Float64, Ngrid)

        # Fill grids
        for i=1:Ngrid
            lj_grid[i] = lj_itp[m.grid_transpose[:, i]...]
            lu_grid[i] = lu_itp[m.grid_transpose[:, i]...]
        end
        lj_grid, lu_grid
    catch e
        # Otherwise guess it is the same as consumption
        warn("Failed to use infile with error:\n$e")
        lj_grid = eval_lj(m, m.other.c1, m.other.c1)
        lu_grid = eval_lu(m, m.other.c2, m.other.c2)

        lj_grid, lu_grid
    end
end

#
# Bring in our files
#
include("HelperCode/alloc_bounds_resids.jl")
include("HelperCode/agents.jl")
include("HelperCode/model_cv.jl")
include("HelperCode/model_sv.jl")
include("HelperCode/analysis.jl")

# ------------------------------------------------------------------- #
# Helper Functions
# ------------------------------------------------------------------- #
"Checks if agents have same sigma parameter"
is_same_sigma(m::AbstractModel) = abs(m.agent1.σ-m.agent1.σ) < 1e-14

"""
Finds a level of (a1, b2) for each possible pair of (lλ☆, lzh)

Takes a lower and upper bound for both lλ☆ and lzh and builds a very
fine grid over both of them. It then uses various first order conditions
to get the allocations (a1, b1, c1, a2, b2, c2) for each possible
combination of lλ☆ and lzh.

## Parameters

- `m::AbstractModel` : The type that describes our model
- `lλ_bds::Tuple{Float64, Float64}` : A lower and upper bound for lλ☆
- `lzh_bds::Tuple{Float64, Float64}` : A lower and upper bound for lzh
- `nlλ☆::Int` : The number of points to include on lλ☆ grid
- `nlzh::Int` : The number of points to include on lzh grid

## Returns

- `a1_itp::Interpolant` : An interpolating function that takes (lλ☆, lzh)
    and returns a value for a1.
- `b2_itp::Interpolant` : An interpolating function that takes (lλ☆, lzh)
    and returns a value for b2.

"""
function solve_for_a1b2(m::AbstractModel, lλ_bds, lzh_bds; nlλ☆=1500, nlzh=250)

    # Unpack parameters
    ρ1, α1, β1, ω1, σ1 = _unpack(m.agent1)
    ρ2, α2, β2, ω2, σ2 = _unpack(m.agent2)

    # Create our fine grids
    lλ☆_fg = linspace(lλ_bds[1]*1.08, lλ_bds[2]*1.08, nlλ☆)
    lzh_fg = linspace(lzh_bds[1], lzh_bds[2], nlzh)

    # Allocate space
    a1_vals = Array(Float64, nlλ☆, nlzh)
    b2_vals = Array(Float64, nlλ☆, nlzh)

    for (i_λ, lλ_val) in enumerate(lλ☆_fg)
        for (i_zh, lzh_val) in enumerate(lzh_fg)
            zh = exp(lzh_val)

            # This function evaluates the residual from combination of
            # focs
            function a1b2_resid(a1)
                # Get allocations corresponding to a1
                b2 = eval_b2(m, a1, lzh_val)
                b1 = 1/zh - b2
                a2 = zh - a1

                # Get consumption
                c1 = _agg_c(m.agent1, a1, b1)
                c2 = _agg_c(m.agent2, b2, a2)

                rhs = log((1-β1)/(1-β2)) + (ρ1-σ1)*log(c1) - (ρ2-σ2)*log(c2)
                rhs += log((1-ω1)/ω2) + (σ1-1)*log(a1) - (σ2-1)*log(a2)

                return lλ_val - rhs
            end

            # Bisect over possible values of a1 and get corresponding b2
            a1_vals[i_λ, i_zh] = brent(a1b2_resid, 1e-6, zh-1e-6)
            b2_vals[i_λ, i_zh] = eval_b2(m, a1_vals[i_λ, i_zh], lzh_val)
        end
    end

    # Build linear interpolants for a1 and b2
    a1_itp = scale(interpolate(a1_vals, BSpline(Linear()), OnGrid()), lλ☆_fg, lzh_fg)
    b2_itp = scale(interpolate(b2_vals, BSpline(Linear()), OnGrid()), lλ☆_fg, lzh_fg)

    return a1_itp, b2_itp
end

"""
Gives us the omega that corresponds with a changed sigma

Takes a value for σ (and optionally the share of imports, s) and returns a
value for omega that should keep the import share to be about the
specified level.
"""
function given_sigma_return_omega(σ::Float64; s::Float64=0.1)
    ratio = ((1-s)/s)^(1-σ)
    ω = 1 / (1 + ratio)
    return ω
end

"""
Takes consumptions and certainty equivalents and evaluates the log of
the value function for those values
"""
function eval_vfs(m::AbstractModel, c1, c2, μ1, μ2)
    lj_new = eval_lj(m, c1, μ1)
    lu_new = eval_lu(m, c2, μ2)
    return lj_new, lu_new
end

"""
Takes two vectors of certainty equivalents (one element for each state
on the grid) and evaluates the value functions using the consumption
values on the grid that are stored in the model
"""
eval_vfs(m::AbstractModel, μ1, μ2) =
    eval_vfs(m, m.other.c1, m.other.c2, μ1, μ2)

"""
Takes value functions at points on the grid and returns a tuple of
interpolants that can evaluate the value function
"""
function update_vf_itps(m::AbstractModel, lj, lu; itptype=Linear())
    Nss = map(length, m.m_ls)
    lj_rs = reshape(lj, Nss...)
    lu_rs = reshape(lu, Nss...)
    lj_itp = scale(interpolate(lj_rs, BSpline(itptype), OnGrid()), m.m_ls...)
    lu_itp = scale(interpolate(lu_rs, BSpline(itptype), OnGrid()), m.m_ls...)
    itps = (lj_itp, lu_itp)
end

# ------------------------------------------------------------------- #
# Solution methods
# ------------------------------------------------------------------- #
"""
Outputs the implied residual when we give it an lk and a lλ☆p

The lowest level of our solution method. It takes a state today (i), a
state tomorrow (j), a value for ratio of certainty equivalents (lk), and
a value for the ratio of Pareto weights tomorrow (lλ☆). It uses these to
evaluate the residual of the foc w.r.t. U_{t+1}.

## Parameters

- `m::AbstractModel` : The model
- `i::Int` : The current period t state we are solving for on grid
- `j::Int` : The current period t+1 we are evaluating residual at
- `vf_itps` : A tuple of interpolants that can evaluate value functions
- `lk::Float64` : Value for log(μ2^(ρ2-α2) / μ1^(ρ1-α1))
- `lλ☆p::Float64` : Value for Pareto weight ratio tomorrow

## Returns

- `::Float64` : The residual of foc when lk and lλ☆p take the given values
"""
function solve_given_lλ☆p(m::AbstractModel, i::Int, j::Int, vf_itps::VFITP,
                          lk::Float64, lλ☆p::Float64)

    # Unpack parameters
    ρ1, α1, β1, ω1, σ1 = _unpack(m.agent1)
    ρ2, α2, β2, ω2, σ2 = _unpack(m.agent2)

    # Get period t state
    lλ☆ = m.grid[i, 1]

    # t+1 exogenous states
    exogp = get_exog_tomorrow(m, i, j)
    lgp = exogp[end]

    # Evaluate value function at t+1
    ljp = vf_itps[1][lλ☆p, exogp[1:end-1]...]
    lup = vf_itps[2][lλ☆p, exogp[1:end-1]...]

    return up_residual(m, ljp, lup, lλ☆p, lgp, lλ☆, lk)
end

"""
Solves for lλ☆ given a specific lk and t+1 value functions

We solve for the policy rules, lλ☆p, as if we already knew what the
value for Lars' k (`lk`) was. We bisect over the `solve_given_lλ☆p`
function in order to find the correct state by state policy choice.
"""
function solve_given_lk(m::AbstractModel, i::Int, vf_itps::VFITP, lk::Float64)

    # Unpack parameters
    ρ1, α1, β1, ω1, σ1 = _unpack(m.agent1)
    ρ2, α2, β2, ω2, σ2 = _unpack(m.agent2)
    Nϵ, Nshocks = size(m.ϵ)

    # Allocate space
    ljp, lup, lλ☆p = map(n->Array(Float64, n), fill(Nϵ, 3))

    # Need to decide how far to go out for bisection. We can't go too
    # far or we move to very inaccurate regions because of extrapolation.
    # The values below seem to do well enough.
    lλ_step = max(0.15, m.m_ls[1][end] / 3)
    lλ☆_big, lλ☆_small = m.grid[i, 1] + lλ_step, m.grid[i, 1] - lλ_step

    # Get relevant state
    exogp = get_exog_tomorrow(m, i)
    lgp_i = exogp[end]
    gp_i = exp(lgp_i)

    # Iterate over every possible state. Create an object that allows
    # us to know whether bisection failed at any state
    it_warn = fill(false, Nϵ)
    for j=1:Nϵ

        vals = map(x->solve_given_lλ☆p(m, i, j, vf_itps, lk, x),
                   (lλ☆_small, lλ☆_big))

        # Check if bisection will work. If so, bisect
        if prod(vals) .<= 0
            it_warn[j] = false
            lλ☆p[j] = brent(x->solve_given_lλ☆p(m, i, j, vf_itps, lk, x),
                            lλ☆_small, lλ☆_big)

        # If not, then choose a bogus policy and change state to a warning
        else
            it_warn[j] = true
            lλ☆p[j] = m.grid[i, 1]
        end

        # Fill in values for value function
        curr_exogp_state = map(x->x[j], exogp[1:end-1])
        ljp[j] = vf_itps[1][lλ☆p[j], curr_exogp_state...]
        lup[j] = vf_itps[2][lλ☆p[j], curr_exogp_state...]
    end

    if any(it_warn)
        warn("Brent warning $(sum(it_warn)) at $i")
    end

    # Get certainty equivalents
    lμ1 = eval_lμ1(m, gp_i, exp(ljp))
    lμ2 = eval_lμ2(m, gp_i, exp(lup))

    return lk_residual(m, lμ1, lμ2, lk), lλ☆p, lμ1, lμ2
end

"""
Updates value functions using a specified value function

We want to be able to use Howard Policy Iteration so we must be able to
update value functions for a given policy. This method updates for a
specific state `i` and manages to skip any expensive bisection and
root-finding which makes it significantly faster than calling
`solve_recursiv_case`. We return the Lars k and two certainty
equivalents
"""
function solve_given_policy(m::AbstractModel, i::Int, lλ☆p::Vector{Float64},
                            vf_itps::VFITP)
    # Unpack parameters
    ρ1, α1, β1, ω1, σ1 = _unpack(m.agent1)
    ρ2, α2, β2, ω2, σ2 = _unpack(m.agent2)
    Nϵ, Nshocks = size(m.ϵ)

    # Allocate space for t+1 value functions
    ljp, lup = map(n->Array(Float64, n), fill(Nϵ, 2))

    # Get relevant states
    exogp = get_exog_tomorrow(m, i)
    lgp_i = exogp[end]
    gp_i = exp(lgp_i)

    for j=1:Nϵ
        curr_exogp_state = map(x->x[j], exogp[1:end-1])
        ljp[j] = vf_itps[1][lλ☆p[j], curr_exogp_state...]
        lup[j] = vf_itps[2][lλ☆p[j], curr_exogp_state...]
    end

    # Get certainty equivalents
    lμ1_i = eval_lμ1(m, gp_i, exp(ljp))
    lμ2_i = eval_lμ2(m, gp_i, exp(lup))

    # If we use lk=0.0 then lk_residual gives the value for lk
    lk_i = lk_residual(m, lμ1_i, lμ2_i, 0.0)

    return lk_i, lμ1_i, lμ2_i
end

"""
Finds an optimally given policy and gives necessary output to update
value functions for recursive utility

If an agent has recursive utility then the planner will find it optimal
to give him varying Pareto weights. This function solves for these
policies at a specific state (i) and for given value functions (vf_itps).
It returns Lars log k, updated policy, and certainty equivalents.
"""
function solve_recursiv_case(m::AbstractModel, i::Int, vf_itps::VFITP)

    # Unpack parameters
    ρ1, α1, β1, ω1, σ1 = _unpack(m.agent1)
    ρ2, α2, β2, ω2, σ2 = _unpack(m.agent2)
    Nϵ, Nshocks = size(m.ϵ)

    lk_lb, lk_ub = log_lars_k_bounds(m, i, vf_itps)

    # Use brent to get solution of lk
    lk = brent(x->solve_given_lk(m, i, vf_itps, x)[1], lk_lb, lk_ub)
    # Grab remaining solution items
    _, lλ☆p, lμ1, lμ2 = solve_given_lk(m, i, vf_itps, lk)

    return lk, lλ☆p, lμ1, lμ2
end

"""
Finds an optimally given policy and gives necessary output to update
value functions for additive utility

It turns out that if the agents have additive utility then we are able
to analytically solve for the optimal policy (constant Pareto weights).
This function uses these policies to call the `solve_given_policy`
function and returns Lars log k, the (constant) policies, and the
certainty equivalents
"""
function solve_additive_case(m::AbstractModel, i::Int, vf_itps::VFITP)
    # Unpack parameters
    ρ1, α1, β1, ω1, σ1 = _unpack(m.agent1)
    ρ2, α2, β2, ω2, σ2 = _unpack(m.agent2)
    Nϵ, Nshocks = size(m.ϵ)

    # Know policy function (constant relative Pareto weights)
    lλ☆p = fill(m.grid[i, 1], Nϵ)

    # Can get lk, lμ1, lμ2 from `solve_given_policy`
    lk, lμ1, lμ2 = solve_given_policy(m, i, lλ☆p, vf_itps)

    return lk, lλ☆p, lμ1, lμ2
end

"""
Solves a specified model with various options.

## Parameters

- `m::AbstractModel` : Model we want to solve
- `capT::Int` : Maximum number of periods to recurse backwards
- `tol::Float64` : Stopping tolerance for value functions
- `infile::String` : File with starting guess for value functions
- `outfile::String` : Where to save the value funtion and policies
- `print_skip::Int` : How many periods between each print (or save)
- `mapper::Function` : Function that specifies how to solve the states
    `pmap` is for parallel and `map` is for sequential

## Returns

- `lj::Vector` : Agent 1 value function at each grid point
- `lu::Vecto` : Agent 2 value function at each grid point
- `lλ☆p::Matrix` : Policy at each grid point and possible state tomorrow
- `lμ1::Vector` : Certainty equivalent for agent 1
- `lμ2::Vector` : Certainty equivalent for agent 2
- `t::Int` : Number of backward steps taken while solving
- `err::Float64` : The difference in value functions between t+1 and t.
"""
function solve_model(m::AbstractModel; capT::Int=500, tol::Float64=1e-7,
                     infile="solutions/$(jld_prefix(m))_1.jld",
                     outfile="solutions/$(jld_prefix(m))_1.jld",
                     print_skip::Int=1, mapper::Function=pmap)

    # Decide if the model has time/state additive or recursive preferences
    if (m.agent1.α ≈ m.agent1.ρ) && (m.agent2.α ≈ m.agent2.ρ)
        state_solver = solve_additive_case
    else
        state_solver = solve_recursiv_case
    end

    # make sure the directories for the solutions exist
    for fn in (infile, outfile)
        # if we have a directory prefix of fn and it isn't already a directory,
        # create one
        dir = dirname(fn)
        !isempty(dir) && !isdir(dir) && mkdir(dir)
    end

    # Get sizes
    Nss = map(length, m.m_ls)
    Ngrid, Nstates = size(m.grid)

    # Allocate Space
    lk = zeros(Float64, Ngrid)
    lμ1 = zeros(Float64, Ngrid)
    lμ2 = zeros(Float64, Ngrid)
    lλ☆p = zeros(Float64, length(m.Π), Ngrid)
    lλ☆p_old = zeros(Float64, length(m.Π), Ngrid)
    err_info = Array(Float64, 3, capT)

    # Initialize VF by setting μ = c or reading in old solution
    lj, lu = initialize_vf(m; infile=infile)

    local t, err
    t_start = time()
    itptype = Linear()
    for t in capT:-1:1

        # Update coefficients
        itptype =  t > capT+1 ? Linear() : Quadratic(Line())
        itps = update_vf_itps(m, lj, lu, itptype=itptype)

        # Update the policies once
        solns = mapper(i -> state_solver(m, i, itps), 1:Ngrid)
        for i=1:Ngrid
            lk[i], lλ☆p[:, i], lμ1[i], lμ2[i] = solns[i]
        end

        # Update value functions
        lj_new, lu_new = eval_vfs(m, exp(lμ1), exp(lμ2))

        # Compute error and copy new vfs into the old vfs (Need to change vf here! otherwise is the same as previous)
        err = max(maxabs(lj-lj_new), maxabs(lu-lu_new))
        pol_err = maxabs(lλ☆p - lλ☆p_old)
        copy!(lj, lj_new)
        copy!(lu, lu_new)
        copy!(lλ☆p_old, lλ☆p)
        err_info[:, t] = [t, err, pol_err]

        if mod(t, print_skip) == 0
            msg = @sprintf("%-3i %-2.5e  %-2.5e %-16s  %11.05f\n",
                           t, err, pol_err, "Total Time", time()-t_start)
            print(msg)

            # Get relevant info from model
            m_ls = m.m_ls
            Nss = map(length, m_ls)

            # Want to save objects needed to create an interpolant
            jldopen(outfile, "w") do file
                write(file, "lj", reshape(lj, Nss...))
                write(file, "lu", reshape(lu, Nss...))
                write(file, "llambdastar", lλ☆p)
                write(file, "err_info", err_info)
                write(file, "m_ls", m_ls)
                write(file, "m", m)
            end
        end

        # Iterate for at least hpi_maxit steps or until close enough
        if abs(pol_err) > 1e-15
            hpi_steps = abs(round(Int, min(0.0, log10(pol_err)))) * 2
        else
            hpi_steps = 25
        end

        for i=1:hpi_steps
            # Update vfs from last step
            vf_itps = update_vf_itps(m, lj, lu; itptype=itptype)

            # Use policy to update lk and CE values
            foo = mapper(i -> solve_given_policy(m, i, lλ☆p[:, i], vf_itps), 1:Ngrid)
            for i=1:Ngrid
                lk[i], lμ1[i], lμ2[i] = foo[i]
            end

            # Update value functions
            lj_new, lu_new = eval_vfs(m, exp(lμ1), exp(lμ2))

            # Check the HPI conditions
            copy!(lj, lj_new)
            copy!(lu, lu_new)

        end

        # Check for whether to break or print stuff
        if err < tol
            println(err)
            break
        end

    end

    return lj, lu, lλ☆p, lμ1, lμ2, t, err
end

function run_all_models{MT<:AbstractModel}(::Type{MT},
                                           #           α      ρ     σ
                                           params = [(-1.0,   -2.0,  0.0),
                                                     (-1.0,  -1.0,  0.0),
                                                     (-9.0,  -1.0,  0.0),
                                                     (-49.0, -1.0,  0.0),
                                                     (-9.0,  -0.01, 0.0),
                                                     (-9.0,   1/3,  0.0),
                                                     (-9.0,  -1.0,  -0.5),
                                                     (-9.0,  -1.0,  0.5)]
                                          )

    for (i, (α, ρ, σ)) in enumerate(params)
        ω = given_sigma_return_omega(σ)

        # Create model
        m = MT(α1=α, α2=α, ρ1=ρ, ρ2=ρ, σ1=σ, σ2=σ, ω1=ω, ω2=ω)

        # Create file names
        # infile = outfile  # This sets the infile to the previous solution
        outfile = "../solutions/$(jld_prefix(m))_$(i).jld"
        infile = "../solutions/$(jld_prefix(m))_$(i).jld"

        try
            println("Solving model with α=$α, ρ=$ρ, σ=$σ")
            Exchange.solve_model(m;infile=infile, outfile=outfile, mapper=map)
        catch e
            rethrow(e)
        end

        break

    end
end

end  # Module end


text:
################################################################################
#  Copyright 2020, Marta Vanin, Tom Van Acker                                  #
################################################################################
# PowerModelsDistributionStateEstimation.jl                                    #
# An extention package of PowerModelsDistribution.jl for Static Power System   #
# State Estimation.                                                            #
################################################################################
mutable struct ReducedIVRUPowerModel <: _PMD.AbstractUnbalancedIVRModel _PMD.@pmd_fields end

"only total current variables defined over the bus_arcs in PMD are considered: with no shunt admittance, these are
equivalent to the series current defined over the branches."
function variable_mc_branch_current(pm::ReducedIVRUPowerModel; nw::Int=_IM.nw_id_default, bounded::Bool=true, report::Bool=true, kwargs...)
    _PMD.variable_mc_branch_current_real(pm, nw=nw, bounded=bounded, report=report; kwargs...)
    _PMD.variable_mc_branch_current_imaginary(pm, nw=nw, bounded=bounded, report=report; kwargs...)
end

"if the formulation is not reduced, it is delegated back to PMD"
function variable_mc_branch_current(pm::_PMD.IVRUPowerModel; nw::Int=_IM.nw_id_default, bounded::Bool=true, report::Bool=true, kwargs...)
    _PMD.variable_mc_branch_current(pm, nw=nw, bounded=bounded, report=report; kwargs...)
end

"constraint_mc_generator_power is re-defined here because in PMD the same function only accepts pm::_PMD.IVRUPowerModel.
The content of the function is otherwise identical"
function constraint_mc_generator_power(pm::ReducedIVRUPowerModel, id::Int; nw::Int=_IM.nw_id_default, report::Bool=true, bounded::Bool=true)

    generator = _PMD.ref(pm, nw, :gen, id)
    bus = _PMD.ref(pm, nw,:bus, generator["gen_bus"])

    N = length(generator["connections"])
    pmin = get(generator, "pmin", fill(-Inf, N))
    pmax = get(generator, "pmax", fill( Inf, N))
    qmin = get(generator, "qmin", fill(-Inf, N))
    qmax = get(generator, "qmax", fill( Inf, N))

    if get(generator, "configuration", _PMD.WYE) == _PMD.WYE
        _PMDSE.constraint_mc_generator_power_wye(pm, nw, id, bus["index"], generator["connections"], pmin, pmax, qmin, qmax; report=report, bounded=bounded)
    else
        _PMDSE.constraint_mc_generator_power_delta(pm, nw, id, bus["index"], generator["connections"], pmin, pmax, qmin, qmax; report=report, bounded=bounded)
    end
end

"constraint_mc_generator_power_wye is re-defined here because in PMD the same function only accepts pm::_PMD.IVRUPowerModel.
The content of the function is otherwise identical"
function constraint_mc_generator_power_wye(pm::ReducedIVRUPowerModel, nw::Int, id::Int, bus_id::Int, connections::Vector{Int}, pmin::Vector{<:Real}, pmax::Vector{<:Real}, qmin::Vector{<:Real}, qmax::Vector{<:Real}; report::Bool=true, bounded::Bool=true)

    vr = _PMD.var(pm, nw, :vr, bus_id)
    vi = _PMD.var(pm, nw, :vi, bus_id)
    crg = _PMD.var(pm, nw, :crg, id)
    cig = _PMD.var(pm, nw, :cig, id)

    pg = Vector{JuMP.NonlinearExpression}([])
    qg = Vector{JuMP.NonlinearExpression}([])

    for (idx, c) in enumerate(connections)
        push!(pg, JuMP.@NLexpression(pm.model,  vr[c]*crg[c]+vi[c]*cig[c]))
        push!(qg, JuMP.@NLexpression(pm.model, -vr[c]*cig[c]+vi[c]*crg[c]))
    end

    if bounded
        for (idx,c) in enumerate(connections)
            if pmin[idx]>-Inf
                JuMP.@constraint(pm.model, pmin[idx] .<= vr[c]*crg[c]  + vi[c]*cig[c])
            end
            if pmax[idx]< Inf
                JuMP.@constraint(pm.model, pmax[idx] .>= vr[c]*crg[c]  + vi[c]*cig[c])
            end
            if qmin[idx]>-Inf
                JuMP.@constraint(pm.model, qmin[idx] .<= vi[c]*crg[c]  - vr[c]*cig[c])
            end
            if qmax[idx]< Inf
                JuMP.@constraint(pm.model, qmax[idx] .>= vi[c]*crg[c]  - vr[c]*cig[c])
            end
        end
    end

   _PMD.var(pm, nw, :pg)[id] = JuMP.Containers.DenseAxisArray(pg, connections)
   _PMD.var(pm, nw, :qg)[id] = JuMP.Containers.DenseAxisArray(qg, connections)

end
"constraint_mc_generator_power_delta is re-defined here because in PMD the same function only accepts pm::_PMD.IVRUPowerModel.
The content of the function is otherwise identical"
function constraint_mc_generator_power_delta(pm::ReducedIVRUPowerModel, nw::Int, id::Int, bus_id::Int, connections::Vector{Int}, pmin::Vector{<:Real}, pmax::Vector{<:Real}, qmin::Vector{<:Real}, qmax::Vector{<:Real}; report::Bool=true, bounded::Bool=true)
    vr = _PMD.var(pm, nw, :vr, bus_id)
    vi = _PMD.var(pm, nw, :vi, bus_id)
    crg = _PMD.var(pm, nw, :crg, id)
    cig = _PMD.var(pm, nw, :cig, id)

    nph = length(pmin)

    prev = Dict(c=>connections[(idx+nph-2)%nph+1] for (idx,c) in enumerate(connections))
    next = Dict(c=>connections[idx%nph+1] for (idx,c) in enumerate(connections))

    vrg = Dict()
    vig = Dict()
    for c in connections
        vrg[c] = JuMP.@NLexpression(pm.model, vr[c]-vr[next[c]])
        vig[c] = JuMP.@NLexpression(pm.model, vi[c]-vi[next[c]])
    end

    pg = Vector{JuMP.NonlinearExpression}([])
    qg = Vector{JuMP.NonlinearExpression}([])
    for c in connections
        push!(pg, JuMP.@NLexpression(pm.model,  vrg[c]*crg[c]+vig[c]*cig[c]))
        push!(qg, JuMP.@NLexpression(pm.model, -vrg[c]*cig[c]+vig[c]*crg[c]))
    end

    if bounded
        JuMP.@NLconstraint(pm.model, [i in 1:nph], pmin[i] <= pg[i])
        JuMP.@NLconstraint(pm.model, [i in 1:nph], pmax[i] >= pg[i])
        JuMP.@NLconstraint(pm.model, [i in 1:nph], qmin[i] <= qg[i])
        JuMP.@NLconstraint(pm.model, [i in 1:nph], qmax[i] >= qg[i])
    end

    crg_bus = Vector{JuMP.NonlinearExpression}([])
    cig_bus = Vector{JuMP.NonlinearExpression}([])
    for c in connections
        push!(crg_bus, JuMP.@NLexpression(pm.model, crg[c]-crg[prev[c]]))
        push!(cig_bus, JuMP.@NLexpression(pm.model, cig[c]-cig[prev[c]]))
    end

   _PMD.var(pm, nw, :crg_bus)[id] = JuMP.Containers.DenseAxisArray(crg_bus, connections)
   _PMD.var(pm, nw, :cig_bus)[id] = JuMP.Containers.DenseAxisArray(cig_bus, connections)
   _PMD.var(pm, nw, :pg)[id] = JuMP.Containers.DenseAxisArray(pg, connections)
   _PMD.var(pm, nw, :qg)[id] = JuMP.Containers.DenseAxisArray(qg, connections)

    if report
       _PMD.sol(pm, nw, :gen, id)[:crg_bus] = JuMP.Containers.DenseAxisArray(crg_bus, connections)
       _PMD.sol(pm, nw, :gen, id)[:cig_bus] = JuMP.Containers.DenseAxisArray(cig_bus, connections)
       _PMD.sol(pm, nw, :gen, id)[:pg] = JuMP.Containers.DenseAxisArray(pg, connections)
       _PMD.sol(pm, nw, :gen, id)[:qg] = JuMP.Containers.DenseAxisArray(qg, connections)
    end
end

"Simplified version of constraint_mc_current_balance_se, to perform state
estimation with the reduced IVR formulation: no shunts, no storage elements, no active switches"
function constraint_mc_current_balance_se(pm::_PMDSE.ReducedIVRUPowerModel, i::Int; nw::Int=_IM.nw_id_default)

    bus = _PMD.ref(pm, nw, :bus, i)
    bus_arcs = _PMD.ref(pm, nw, :bus_arcs_conns_branch, i)
    bus_arcs_trans = _PMD.ref(pm, nw, :bus_arcs_conns_transformer, i)
    bus_gens = _PMD.ref(pm, nw, :bus_conns_gen, i)
    bus_loads = _PMD.ref(pm, nw, :bus_conns_load, i)

    cr    = get(_PMD.var(pm, nw),    :cr, Dict()); _PMD._check_var_keys(cr, bus_arcs, "real current", "branch")
    ci    = get(_PMD.var(pm, nw),    :ci, Dict()); _PMD._check_var_keys(ci, bus_arcs, "imaginary current", "branch")
    crd   = get(_PMD.var(pm, nw),   :crd, Dict()); _PMD._check_var_keys(crd, bus_loads, "real current", "load")
    cid   = get(_PMD.var(pm, nw),   :cid, Dict()); _PMD._check_var_keys(cid, bus_loads, "imaginary current", "load")
    crg   = get(_PMD.var(pm, nw),   :crg, Dict()); _PMD._check_var_keys(crg, bus_gens, "real current", "generator")
    cig   = get(_PMD.var(pm, nw),   :cig, Dict()); _PMD._check_var_keys(cig, bus_gens, "imaginary current", "generator")
    crt   = get(_PMD.var(pm, nw),   :crt, Dict()); _PMD._check_var_keys(crt, bus_arcs_trans, "real current", "transformer")
    cit   = get(_PMD.var(pm, nw),   :cit, Dict()); _PMD._check_var_keys(cit, bus_arcs_trans, "imaginary current", "transformer")

    terminals = bus["terminals"]
    grounded =  bus["grounded"]
    ungrounded_terminals = [(idx,t) for (idx,t) in enumerate(terminals) if !grounded[idx]]

    for (idx, t) in ungrounded_terminals
        JuMP.@constraint(pm.model,  sum(cr[a][t] for (a, conns) in bus_arcs if t in conns)
                                    + sum(crt[a_trans][t] for (a_trans, conns) in bus_arcs_trans if t in conns)
                                    ==
                                      sum(crg[g][t]         for (g, conns) in bus_gens if t in conns)
                                    - sum(crd[d][t]         for (d, conns) in bus_loads if t in conns)
                                    )

        JuMP.@constraint(pm.model, sum(ci[a][t] for (a, conns) in bus_arcs if t in conns)
                                    + sum(cit[a_trans][t] for (a_trans, conns) in bus_arcs_trans if t in conns)
                                    ==
                                      sum(cig[g][t]         for (g, conns) in bus_gens if t in conns)
                                    - sum(cid[d][t]         for (d, conns) in bus_loads if t in conns)
                                    )
    end
end

"Simplified version of constraint_mc_current_balance, to perform power flow
calculations with the reduced IVR formulation: no shunts, no storage elements, no active switches"
function constraint_mc_current_balance(pm::_PMDSE.ReducedIVRUPowerModel, i::Int; nw::Int=_IM.nw_id_default)

    bus = _PMD.ref(pm, nw, :bus, i)
    bus_arcs = _PMD.ref(pm, nw, :bus_arcs_conns_branch, i)
    bus_arcs_trans = _PMD.ref(pm, nw, :bus_arcs_conns_transformer, i)
    bus_gens = _PMD.ref(pm, nw, :bus_conns_gen, i)
    bus_loads = _PMD.ref(pm, nw, :bus_conns_load, i)

    cr    = get(_PMD.var(pm, nw),    :cr, Dict()); _PMD._check_var_keys(cr, bus_arcs, "real current", "branch")
    ci    = get(_PMD.var(pm, nw),    :ci, Dict()); _PMD._check_var_keys(ci, bus_arcs, "imaginary current", "branch")
    crd   = get(_PMD.var(pm, nw),   :crd, Dict()); _PMD._check_var_keys(crd, bus_loads, "real current", "load")
    cid   = get(_PMD.var(pm, nw),   :cid, Dict()); _PMD._check_var_keys(cid, bus_loads, "imaginary current", "load")
    crg   = get(_PMD.var(pm, nw),   :crg, Dict()); _PMD._check_var_keys(crg, bus_gens, "real current", "generator")
    cig   = get(_PMD.var(pm, nw),   :cig, Dict()); _PMD._check_var_keys(cig, bus_gens, "imaginary current", "generator")
    crt   = get(_PMD.var(pm, nw),   :crt, Dict()); _PMD._check_var_keys(crt, bus_arcs_trans, "real current", "transformer")
    cit   = get(_PMD.var(pm, nw),   :cit, Dict()); _PMD._check_var_keys(cit, bus_arcs_trans, "imaginary current", "transformer")

    terminals = bus["terminals"]
    grounded =  bus["grounded"]
    ungrounded_terminals = [(idx,t) for (idx,t) in enumerate(terminals) if !grounded[idx]]

    for (idx, t) in ungrounded_terminals
        JuMP.@constraint(pm.model,  sum(cr[a][t] for (a, conns) in bus_arcs if t in conns)
                                    + sum(crt[a_trans][t] for (a_trans, conns) in bus_arcs_trans if t in conns)
                                    ==
                                      sum(crg[g][t]         for (g, conns) in bus_gens if t in conns)
                                    - sum(crd[d][t]         for (d, conns) in bus_loads if t in conns)
                                    )

        JuMP.@constraint(pm.model, sum(ci[a][t] for (a, conns) in bus_arcs if t in conns)
                                    + sum(cit[a_trans][t] for (a_trans, conns) in bus_arcs_trans if t in conns)
                                    ==
                                      sum(cig[g][t]         for (g, conns) in bus_gens if t in conns)
                                    - sum(cid[d][t]         for (d, conns) in bus_loads if t in conns)
                                    )
    end
end

"Constraint that defines the voltage drop along a branch for the ReducedIVRUPowerModel.
Conceptually similar to the same function in PowerModelsDistribution, but written in
function of the total current instead of the series current.
This is because in the ReducedIVRUPowerModel series current is not defined as it is identical to
the total current due to the absence of shunt admittance."
function constraint_mc_bus_voltage_drop(pm::ReducedIVRUPowerModel, i::Int; nw::Int=_IM.nw_id_default)

    branch = _PMD.ref(pm, nw, :branch, i)
    f_bus = branch["f_bus"]
    t_bus = branch["t_bus"]
    f_idx = (i, f_bus, t_bus)

    r = branch["br_r"]
    x = branch["br_x"]

    vr_fr = [_PMD.var(pm, nw, :vr, f_bus)[c] for c in branch["f_connections"]]
    vi_fr = [_PMD.var(pm, nw, :vi, f_bus)[c] for c in branch["f_connections"]]

    vr_to = [_PMD.var(pm, nw, :vr, t_bus)[c] for c in branch["t_connections"]]
    vi_to = [_PMD.var(pm, nw, :vi, t_bus)[c] for c in branch["t_connections"]]

    cr_fr = [_PMD.var(pm, nw, :cr, f_idx)[c] for c in branch["f_connections"]]
    ci_fr = [_PMD.var(pm, nw, :ci, f_idx)[c] for c in branch["f_connections"]]

    JuMP.@constraint(pm.model, vr_to .== vr_fr - r*cr_fr + x*ci_fr)
    JuMP.@constraint(pm.model, vi_to .== vi_fr - r*ci_fr - x*cr_fr)
end

"This constraint makes sure that the current entering and exiting a branch are equivalent.
It is used for state estimation and power flow calculations with the ReducedIVRUPowerModel."
function constraint_current_to_from(pm::ReducedIVRUPowerModel, i::Int; nw::Int=_IM.nw_id_default)

    branch = _PMD.ref(pm, nw, :branch, i)
    f_bus = branch["f_bus"]
    t_bus = branch["t_bus"]
    t_arc = (i, t_bus, f_bus)
    f_arc = (i, f_bus, t_bus)

    cr_fr = [_PMD.var(pm, nw, :cr, f_arc)[c] for c in branch["f_connections"]]
    ci_fr = [_PMD.var(pm, nw, :ci, f_arc)[c] for c in branch["f_connections"]]

    cr_to = [_PMD.var(pm, nw, :cr, t_arc)[c] for c in branch["t_connections"]]
    ci_to = [_PMD.var(pm, nw, :ci, t_arc)[c] for c in branch["t_connections"]]

    JuMP.@constraint(pm.model, cr_fr .== -cr_to)
    JuMP.@constraint(pm.model, ci_fr .== -ci_to)

end


text:
__precompile__()

module mbRK
export rk1, rk4

import Base.length

type ButcherTableau
	A::Array{Float64,2}
	b::Array{Float64,1}
	c::Array{Float64,1}
end

length(BT::ButcherTableau) = length(BT.b)

"""
	rk(butcherTableau)

returns a runge-kutta method using the respective tableau:

    function(f, t0, y0, h, pred[, init, callback])

e.g. `rk1`, or `rk4`. Examines the ode `f` starting from `t0`, `y0` with fixed stepsize
`h` until `pred` evalutes to `false`.
"""
function rk(BT::ButcherTableau)
	return function rk(f::Function, t₀::Float64, y₀::Vector{Float64}, h::Float64, pred::Function;
		init=Void, callback=Void)

		local A,b,c,s,t,y,K
		A,b,c = BT.A, BT.b, BT.c
		s = length(BT)
		K = cell(s)

		t,y = deepcopy(t₀), deepcopy(y₀)
		K[1] = f(t,y)

		init≠Void && init()

		while pred(t,y,K[1])
			for i in 2:s
				K[i] = f( t + h*c[i], y + h*(A[i,1:i-1]*K[1:i-1])[1] )
			end

			t += h
			y += h*sum(map(*, K, b))
			K[1] = f(t,y)

			callback≠Void && callback(t,y,K)
		end

		return Void
	end
end

BTrk1 = ButcherTableau(
	zeros(1,1),
	ones(1),
	zeros(1)
)
rk1 = rk(BTrk1)

BTrk4 = ButcherTableau(
	[
		.0	.0	.0	.0
		.5	.0	.0	.0
		.0	.5	.0	.0
		.0	.0	1.0	.0
	],
	[1/6; 1/3; 1/3; 1/6],
	[.0; .5; .5; 1.0]
)
rk4 = rk(BTrk4)

end


text:
struct Match{R}
    vertices::Vector{Int}
    half_edges::Vector{Int} # half edge id
end

function Base.show(io::IO, m::Match)
    print(io, "Match:")
    if !isempty(m.vertices)
        println(io)
        print(io, "  vertices = ", m.vertices)
    end

    if !isempty(m.half_edges)
        println(io)
        print(io, "  half_edge = ", m.half_edges)
    end
    return
end

"""
    Rule{R}

A struct for rules. `R` can be `:string_genus`, `yang_baxter_star`, `yang_baxter_triangle`,
`charge_rm_v`, `charge_rm_f`, `z_fusion`, `x_fusion`, `perm_rz`, `identity`, `genus_fusion`,
`swap_genus`.
"""
struct Rule{T} end
Rule(r::Symbol) = Rule{r}()

Base.match(r::Rule{R}, tait::Tait{P}) where {R, P} = match!(Match{R}[], r, tait)

function match!(matches, ::Rule{:string_genus}, tait::Tait)
    for (v, _) in tait.g.vs_isolated
        v in tait.genuses && push!(matches, Match{:string_genus}([v], Int[]))
    end
    return matches
end

function match!(matches, ::Rule{:yang_baxter_star}, tait::Tait)
    for v in vertices(tait)
        (is_open_vertex(tait, v) || is_genus(tait, v)) && continue
        hes = trace_vertex(tait, v)
        if length(hes) == 3
            nbs = Set(dst(tait, he) for he in hes)
            if length(nbs) == 3 && all(x->!is_open_vertex(tait, dst(tait, x)), hes)
                push!(matches, Match{:yang_baxter_star}([v], hes))
            end
        end
    end
    return matches
end

function match!(matches, ::Rule{:yang_baxter_triangle}, tait::Tait)
    for f in faces(tait)
        f == 0 && continue
        hes = trace_face(tait, f)
        length(hes) == 3 || continue
        all(x->!is_open_vertex(tait, dst(tait, x)), hes) || continue
        push!(matches, Match{:yang_baxter_triangle}([], hes))
    end
    return matches
end

function match!(matches, ::Rule{:charge_rm_v}, tait::Tait)
    for v in vertices(tait)
        is_open_vertex(tait, v) && continue
        hes = trace_vertex(tait, v)
        if length(hes) == 1
            if !is_open_half_edge(tait, hes[1])
                p = quon_param(tait, hes[1])
                if is_para_pi(p)
                    push!(matches, Match{:charge_rm_v}([v], hes))
                end
            end
            continue
        end

        i1 = findfirst(hes) do x
            if is_open_half_edge(tait, x)
                return true
            else
                return !is_para_pi(quon_param(tait, x))
            end
        end
        if i1 === nothing
            !isempty(hes) && push!(matches, Match{:charge_rm_v}([v], hes))
            continue
        end

        hes = permute!(hes, [collect(i1:length(hes)); collect(1:(i1-1))])
        hes_match = Int[]
        for i = eachindex(hes)
            is_match = false
            he = hes[i]
            if !is_open_half_edge(tait, he)
                p = quon_param(tait, he)
                if is_para_pi(p)
                    push!(hes_match, he)
                    is_match = true
                    if he == hes[end]
                        (length(hes_match) > 1) && push!(matches, Match{:charge_rm_v}([v], hes_match))
                    end
                end
            end
            if !is_match
                if length(hes_match) > 1
                    push!(matches, Match{:charge_rm_v}([v], hes_match))
                end
                !isempty(hes_match) && (hes_match = Int[])
            end
        end
    end
    return matches
end

function match!(matches, ::Rule{:charge_rm_f}, tait::Tait)
    for f in faces(tait)
        f == 0 && continue
        hes = trace_face(tait, f)

        i1 = findfirst(hes) do x
            if is_open_half_edge(tait, x)
                return true
            else
                return !is_para_pi(quon_param(tait, x))
            end
        end
        if i1 === nothing
            println(f)
            (length(hes) > 1) && push!(matches, Match{:charge_rm_f}([], hes))
            continue
        end

        hes = permute!(hes, [collect(i1:length(hes)); collect(1:(i1-1))])
        hes_match = Int[]
        for i = eachindex(hes)
            is_match = false
            he = hes[i]
            if !is_open_half_edge(tait, he)
                p = quon_param(tait, he)
                if is_orth_pi(p)
                    push!(hes_match, he)
                    is_match = true
                    if he == hes[end]
                        (length(hes_match) > 1) && push!(matches, Match{:charge_rm_f}([], hes_match))
                    end
                end
            end
            if !is_match
                if length(hes_match) > 1
                    push!(matches, Match{:charge_rm_f}([], hes_match))
                end
                !isempty(hes_match) && (hes_match = Int[])
            end
        end
    end
    return matches
end

function match!(matches, ::Rule{:z_fusion}, tait::Tait)
    for f in faces(tait)
        f == 0 && continue
        hes = trace_face(tait, f)
        length(hes) == 2 || continue
        all(he -> !is_open_half_edge(tait, he), hes) || continue
        push!(matches, Match{:z_fusion}([], hes))
    end
    return matches
end

function match!(matches, ::Rule{:x_fusion}, tait::Tait)
    for v in vertices(tait)
        is_genus(tait, v) && continue
        is_open_vertex(tait, v) && continue
        hes = trace_vertex(tait, v)
        length(hes) == 2 || continue
        all(he -> !is_open_half_edge(tait, he), hes) || continue
        push!(matches, Match{:x_fusion}([v], hes))
    end
    return matches
end

function match!(matches, ::Rule{:perm_rz}, tait::Tait)
    for v in vertices(tait)
        is_open_vertex(tait, v) && continue
        hes = trace_vertex(tait, v)
        # length(hes) >= 3 || continue
        ids = findall(hes) do he
            is_genus(tait, dst(tait, he))
        end
        # TODO: Should we check isolated genuses?
        length(ids) == 0 && continue
        for idx in ids
            he_genus = hes[idx]
            # find half_edge on the other face that has a genus
            he_out = σ_inv(tait, he_genus)
            for _ in 2:(length(hes) - 1)
                he_out = σ_inv(tait, he_out)
                he = he_out
                while !is_genus(tait, src(tait, he))
                    he = next(tait, he)
                    he == he_out && break
                end
                if he != he_out && src(tait, he) != dst(tait, he_genus)
                    push!(matches, Match{:perm_rz}([v, src(tait, he)], [he_genus, he]))
                end
            end
        end
    end
    return matches
end

function match!(matches, ::Rule{:identity}, tait::Tait)
    for (he_id, theta) in tait.quon_params
        if is_para_zero(theta) || is_orth_zero(theta)
            push!(matches, Match{:identity}([], [he_id]))
        end
    end
    return matches
end

function match!(matches, ::Rule{:genus_fusion}, tait::Tait)
    for g in genuses(tait)
        is_isolated(tait, g) && continue
        hes_g = trace_vertex(tait, g)
        for he_g in hes_g
            face(tait, he_g) == 0 && continue
            he = next(tait, he_g)
            while he != he_g
                v = src(tait, he)
                if is_genus(tait, v)
                    push!(matches, Match{:genus_fusion}([g, v], [he_g, he]))
                end
                he = next(tait, he)
            end
        end
    end
    return matches
end

function match!(matches, ::Rule{:swap_genus}, tait::Tait)
    for v0 in vertices(tait)
        is_genus(tait, v0) && continue
        hes_v0 = trace_vertex(tait, v0)
        length(hes_v0) == 4 || continue
        vs = [dst(tait, he) for he in hes_v0]

        faces_contain_4_edges = true
        for he in hes_v0
            f = face(tait, he)
            f_hes = trace_face(tait, f)
            if length(f_hes) != 4
                faces_contain_4_edges = false
                break
            end
        end
        faces_contain_4_edges || continue

        vertices_contain_4_neighbors = true
        for v in vs
            v_hes = trace_vertex(tait, v)
            if is_genus(tait, v) || length(v_hes) != 4
                vertices_contain_4_neighbors = false
                break
            end
        end
        vertices_contain_4_neighbors || continue

        # check parameters
        check_swap_parameters(tait, v0) || continue
        for he in hes_v0
            v = dst(tait, he)
            check_swap_parameters(tait, v) || continue
        end

        has_4_connected_genuses = true
        connected_genuses = Int[]
        for he in hes_v0
            he_twin = twin(tait, he)
            he_g = σ_inv(tait, σ_inv(tait, he_twin))
            v_g = dst(tait, he_g)
            if is_genus(tait, v_g)
                push!(connected_genuses, v_g)
            else
                has_4_connected_genuses = false
                break
            end
        end
        has_4_connected_genuses || continue

        # find genus has open edge
        idx = findfirst(x->!is_genus_connected_to_open_edge(tait, x), connected_genuses)
        idx === nothing && continue

        push!(matches,
            Match{:swap_genus}(
                [connected_genuses[idx:end]; connected_genuses[1:idx-1];
                    vs[idx:end]; vs[1:idx-1]; v0
                ],
                [hes_v0[idx:end]; hes_v0[1:idx-1]]
            )
        )
    end
    return matches
end

function check_swap_parameters(tait::Tait, v)
    hes = trace_vertex(tait, v)
    length(hes) == 4 || return false
    any(is_open_half_edge(tait, he) for he in hes) && return false
    (p1, p2, p3, p4) = [quon_param(tait, he) for he in hes]
    if is_para_half_pi(p1) && is_para_half_pi(p3) &&
        is_orth_half_pi(p2) && is_orth_half_pi(p4)
        return true
    elseif is_para_half_pi(p2) && is_para_half_pi(p4) &&
        is_orth_half_pi(p1) && is_orth_half_pi(p3)
        return true
    end
    return false
end

function is_genus_connected_to_open_edge(tait::Tait, genus)
    for he in trace_vertex(tait, genus)
        if is_open_vertex(tait, dst(tait, he))
            return true
        end
    end
    return false
end


text:
using PyCall, PyPlot
@pyimport pandas as pd
@pyimport seaborn as sns

# This script relies on data that is generated by running smc-abc-script.jl

emu_plot_data = pd.DataFrame(data= Dict( :param_1=>emu_out.population[3][:,1],:param_2=>emu_out.population[3][:,2],:param_3=>emu_out.population[3][:,3]))
sim_plot_data = pd.DataFrame(data= Dict( :param_1=>sim_out.population[3][:,1],:param_2=>sim_out.population[3][:,2],:param_3=>sim_out.population[3][:,3]))
PyPlot.ion()
fig, ax = PyPlot.subplots(figsize=(15,15),ncols=3, nrows=3)
     left   =  0.125  # the left side of the subplots of the figure
     right  =  0.9    # the right side of the subplots of the figure
     bottom =  0.1    # the bottom of the subplots of the figure
     top    =  0.9    # the top of the subplots of the figure
     wspace =  .5     # the amount of width reserved for blank space between subplots
     hspace =  .5    # the amount of height reserved for white space between subplots
     PyPlot.subplots_adjust(
         left    =  left,
         bottom  =  bottom,
         right   =  right,
         top     =  top,
         wspace  =  wspace,
         hspace  =  hspace
     )
     PyPlot.ioff()
     cmap = sns.cubehelix_palette(as_cmap=true, dark=0, light=1, reverse=true)
     # marginals
     sns.distplot(sim_plot_data[:param_1],bins=15, kde= true, hist=true,rug= false ,color = "cornflowerblue",ax=ax[1,1])
     sns.distplot(sim_plot_data[:param_2],bins=15, kde= true, hist=true,rug= false ,   color = "cornflowerblue", ax=ax[2,2])
     sns.distplot(sim_plot_data[:param_3],bins=15, kde= true, hist=true,rug= false ,  color = "cornflowerblue", ax=ax[3,3])
     sns.distplot(emu_plot_data[:param_1],bins=15, kde= true, hist=true,rug= false ,   color = "salmon", label = "lisi",axlabel = "Parameter 1",ax=ax[1,1])
    sns.distplot(emu_plot_data[:param_2],bins=15, kde= true, hist=true,rug= false ,  color = "salmon", label = "lisi",axlabel = "Parameter 2",ax=ax[2,2])
    sns.distplot(emu_plot_data[:param_3],bins=15, kde= true, hist=true,rug= false ,   color = "salmon", label = "lisi",axlabel = "Parameter 3",ax=ax[3,3])
    #joint
     sns.jointplot(x="param_1", y="param_2", data=emu_plot_data, kind="kde",color="lightpink", ax=ax[2,1])

     sns.jointplot(x="param_1", y="param_3", data=emu_plot_data, kind="kde", color="lightpink",ax=ax[3,1])
     sns.jointplot(x="param_2", y="param_1", data=emu_plot_data, kind="kde", color="lightpink",ax=ax[3,2])
     #HERE scatter
     sns.jointplot(y="param_1", x="param_2", data=sim_plot_data, kind="kde", color="cornflowerblue", ax=ax[1,2])
     sns.jointplot(y="param_1", x="param_3", data=sim_plot_data, kind="kde", color="cornflowerblue",ax=ax[1,3])
     sns.jointplot(y="param_2", x="param_1", data=sim_plot_data, kind="kde", color="cornflowerblue",ax=ax[2,3])
     ax[1,2][:set_xlabel]("Parameter 2")
     ax[1,2][:set_ylabel]("Parameter 1")
     #
     ax[1,3][:set_xlabel]("Parameter 3")
     ax[1,3][:set_ylabel]("Parameter 1")
     #
     ax[2,3][:set_xlabel]("Parameter 3")
     ax[2,3][:set_ylabel]("Parameter 2")
     #
     ax[2,1][:set_xlabel]("Parameter 1")
     ax[2,1][:set_ylabel]("Parameter 2")
     #
     ax[3,1][:set_xlabel]("Parameter 1")
     ax[3,1][:set_ylabel]("Parameter 3")
     #
     ax[3,2][:set_xlabel]("Parameter 2")
     ax[3,2][:set_ylabel]("Parameter 3")
PyPlot.show(fig)


text:
# setup cfg to point to test channel parameters JSON file
import YAML
cfg_folder = joinpath(dirname(@__FILE__), "..", "cfg")
cfg_file = joinpath(cfg_folder, "cfg_path.txt")

# backup existing one
restore_backup = false
if isfile(cfg_file)
	cp(cfg_file, cfg_file*".orig")
	restore_backup = true
end
if !isdir(cfg_folder)
	mkdir(cfg_folder)
end

open(cfg_file, "w") do f
	write(f, joinpath(cfg_folder, "test_measure.yml"))
end
# 	cfg = Dict{String,String}()
# 	cfg["channel_params_file"] = joinpath(dirname(@__FILE__), "ChannelParams.json")
# 	cfg["instrument_params_file"] = joinpath(dirname(@__FILE__), "Instruments.json")
# 	cfg["sequence_files_path"] = sequence_files_path
# 	JSON.print(f, cfg)
# end

using QGL
using Base.Test

try
@testset "Qubits" begin
q1 = Qubit("q1")
@test q1.label == "q1"
end

@testset "Pulse Creation" begin
q1 = Qubit("q1")
p = X90(q1)
@test p.label == "X90"
end

finally
if restore_backup
	mv(cfg_file*".orig", cfg_file, remove_destination=true)
end

end


text:
# install unregistered packages
Pkg.clone("https://github.com/tpapp/RaggedData.jl.git")
Pkg.clone("https://github.com/tpapp/ByteParsers.jl.git")
Pkg.clone("https://github.com/tpapp/DiscreteRanges.jl.git")
Pkg.clone("https://github.com/tpapp/LargeColumns.jl.git")
Pkg.clone("https://github.com/tpapp/FlexDates.jl.git")
Pkg.clone("https://github.com/tpapp/WallTimeProgress.jl.git")
# build and test
Pkg.clone(pwd())
Pkg.build("AMDB")
Pkg.test("AMDB"; coverage=true)


text:
using SeisAcoustic

# homogeneous velocity and density model
nz = 51; nx = 75;
vel = 3000 * ones(nz, nx);  # m/s
vel[26:end,:] .= 3500;      # m/s
rho = 2000 * ones(nz, nx);  # kg/m^3

# vertical and horizontal grid size
dz = 10; dx = 10;

# time step size and maximum modelling length
dt = 0.001; tmax = 0.2;  # use second as unit

# receiver location
irx = collect(1 :2 :nx);
irz = 2 * ones(Int64, length(irx));

# source location
isx1 = floor(Int64, nx/2); isz1 = 2;
isx2 = collect(5:15:nx-5); ns=length(isx2); isz2 = 2*ones(ns);

test_type = [(true , Float32, 1.0e-6 ),
             (true , Float64, 1.0e-12),
             (false, Float32, 1.0e-6 ),
             (false, Float64, 1.0e-12)
            ]

(free_surface, data_format, tol) = test_type[1]

# generate finite difference stencil
params = TdParams(rho, vel, free_surface, dz, dx, dt, tmax;
         data_format=data_format, fd_flag="taylor", order=3, npml=20, apml=900.);

# create multiple source term
source_term = get_multi_sources(isz2, isx2, params; ot=0.0*rand(ns), fdom=15.0, amp=100000);

# pre-compute source side wavefield and their boundary
dir_sourceside = joinpath(homedir(), "Desktop/source_side");
get_wavefield_bound(dir_sourceside, source_term, params);

# forward born approximation
# generate input for forward operator
m1      = randn(params.data_format, params.nz*params.nx)
path_m1 = joinpath(homedir(), "Desktop/m1.rsf");
hdr1    = RegularSampleHeader(m1);
write_RSdata(path_m1, hdr1, m1);

dir_born1 = joinpath(homedir(), "Desktop/forward_born");
born_approximation_forward(dir_born1, path_m1, irz, irx, dir_sourceside, params);

# adjoint operator
# adjoint born approximation
w   = ricker(15.0, params.dt)
hl  = floor(Int64, length(w)/2)

dir_born2 = joinpath(homedir(), "Desktop/adjoint_born");
rm(dir_born2, force=true, recursive=true)
mkdir(dir_born2)
for i = 1 : length(source_term)
    # generate data trace by trace
    rec2 = Recordings(irz, irx, params);
    for i = 1 : rec2.nr
        tmp = conv(randn(params.nt)*1000, w)
        tmp = convert(Vector{params.data_format}, tmp[hl+1:end-hl])
        rec2.p[:,i] .= tmp
    end

    file_name = join(["recordings_" "$i" ".bin"])
    path_rec  = joinpath(dir_born2, file_name)

    write_recordings(path_rec, rec2)
end


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

struct AuthApi <: SwaggerApi
    client::Swagger.Client
end

function _swaggerinternal_authRevoke(_api::AuthApi; test=nothing, _mediaType=nothing)
    _ctx = Swagger.Ctx(_api.client, "GET", Dict{String, Any}, "/auth.revoke", ["slackAuth"])
    Swagger.set_param(_ctx.query, "test", test)  # type Bool
    Swagger.set_header_accept(_ctx, ["application/json"])
    Swagger.set_header_content_type(_ctx, (_mediaType === nothing) ? ["application/x-www-form-urlencoded"] : [_mediaType])
    return _ctx
end

"""

Revokes a token.
Param: test::Bool
Return: Dict{String, Any}
"""
function authRevoke(_api::AuthApi; test=nothing, _mediaType=nothing)
    _ctx = _swaggerinternal_authRevoke(_api; test=test, _mediaType=_mediaType)
    Swagger.exec(_ctx)
end

function authRevoke(_api::AuthApi, response_stream::Channel; test=nothing, _mediaType=nothing)
    _ctx = _swaggerinternal_authRevoke(_api; test=test, _mediaType=_mediaType)
    Swagger.exec(_ctx, response_stream)
end

function _swaggerinternal_authTest(_api::AuthApi; _mediaType=nothing)
    _ctx = Swagger.Ctx(_api.client, "GET", Dict{String, Any}, "/auth.test", ["slackAuth"])
    Swagger.set_header_accept(_ctx, ["application/json"])
    Swagger.set_header_content_type(_ctx, (_mediaType === nothing) ? ["application/x-www-form-urlencoded", "application/json"] : [_mediaType])
    return _ctx
end

"""

Checks authentication & identity.
Return: Dict{String, Any}
"""
function authTest(_api::AuthApi; _mediaType=nothing)
    _ctx = _swaggerinternal_authTest(_api; _mediaType=_mediaType)
    Swagger.exec(_ctx)
end

function authTest(_api::AuthApi, response_stream::Channel; _mediaType=nothing)
    _ctx = _swaggerinternal_authTest(_api; _mediaType=_mediaType)
    Swagger.exec(_ctx, response_stream)
end

export authRevoke, authTest


text:
@testset "check softmax op's gradient" begin
    clear()
    for d in [1 2 3 (1,2) (2,3) (1,3) (1,2,3)]
        @testset "check softmax op's gradient at dim = $d" begin
            DIMS = d
            TYPE = Array{Float64};

            # [1] prepare input data and its label
            inputdims = 64;
            timeSteps = 16;
            batchsize = 32;
            x = Variable(rand(inputdims, timeSteps, batchsize), type=TYPE, keepsgrad=true);
            l = Variable(rand(inputdims, timeSteps, batchsize), type=TYPE);

            # [2] forward and backward propagation
            probs = softmax(x; dims=DIMS)
            Loss1 = crossEntropyLoss(probs, l);
            Loss2 = mseLoss(probs, l);
            COST1 = cost(0.8*Loss1 + 0.2*Loss2)
            backward();

            # [3] with a samll change of a weight
            GRAD = x.delta[1];
            DELTA = 1e-5;
            x.value[1] += DELTA;

            # [4] forward and backward propagation with a samll change of a weight
            probs = softmax(x; dims=DIMS)
            Loss1 = crossEntropyLoss(probs, l);
            Loss2 = mseLoss(probs, l);
            COST2 = cost(0.8*Loss1 + 0.2*Loss2)
            backward();

            # [5] check if the auto-grad is true or not
            dLdW = (COST2 - COST1)/DELTA;   # numerical gradient
            err  = abs((dLdW-GRAD)/(GRAD+eps(Float64)))*100;  # relative error in %
            err  = err < 1e-1 ? 0.0 : err;
            @test err<1.0
        end
    end
end

max_stars_repo_name: sonosole/Delta

max_issues_repo_name: sonosole/Delta

max_forks_repo_name: sonosole/Delta


text:
using Test
using QuantumOpticsBase

@testset "state_definitions" begin

n=30
b=FockBasis(n)
omega=40.3
T=2.3756
r=thermalstate(omega*number(b),T)
for k=1:n-1
    @test isapprox(r.data[k+1,k+1]/r.data[k,k],exp(-omega/T))
end
S=entropy_vn(r)
z=sum(exp.(-[0:n;]*omega))
s=expect(omega*number(b),r)/T+log(z)
isapprox(S,s)

alpha=rand()+im*rand()
r=coherentthermalstate(b,omega*number(b),T,alpha)
@test isapprox(expect(number(b),r),abs(alpha)^2+1/(exp(omega/T)-1), atol=1e-14)
@test isapprox(expect(destroy(b),r),alpha, atol=1e-14)
@test isapprox(entropy_vn(r),S, atol=1e-14)

rp=phase_average(r)
@test isapprox(expect(number(b),rp),abs(alpha)^2+1/(exp(omega/T)-1), atol=1e-14)
@test isapprox(expect(destroy(b),rp),0, atol=1e-14)
for k=1:n
    @test isapprox(rp.data[k,k],r.data[k,k], atol=1e-14)
end

rpas=passive_state(r)
for k=1:n-1
    @test real(rpas.data[k+1,k+1])<real(rpas.data[k,k])
end
@test isapprox(expect(number(b),rpas),1/(exp(omega/T)-1), atol=1e-14)
@test isapprox(expect(destroy(b),rpas),0, atol=1e-14)
@test isapprox(entropy_vn(rpas),S, atol=1e-14)

end # testset


text:
abstract type AbstractIntegratorRK{dType, tType} <: ODEIntegrator{dType, tType} end
abstract type AbstractIntegratorIRK{dType, tType} <: IODEIntegrator{dType, tType} end
abstract type AbstractIntegratorPRK{dType, tType} <: PODEIntegrator{dType, tType} end

IntegratorRK = Union{AbstractIntegratorRK, AbstractIntegratorIRK, AbstractIntegratorPRK}

@inline equation(integrator::IntegratorRK, i::Symbol) = integrator.params.equs[i]
@inline equations(integrator::IntegratorRK) = integrator.params.equs
@inline timestep(integrator::IntegratorRK) = integrator.params.Δt
@inline tableau(integrator::IntegratorRK)  = integrator.params.tab

@inline nstages(integrator::IntegratorRK)  = nstages(tableau(integrator))
@inline eachstage(integrator::IntegratorRK) = 1:nstages(integrator)


text:
using Base: @propagate_inbounds, @_inline_meta

"""
    Window{X}(k::Kernel, a::DenseArray, pos::CartesianIndex)

Create a stack-allocated view on `a` with interior axes `X` and cartesian indexing
relative to `pos` in the parent array.
When indexing outside `X` the returned value is determined by the extension
specified in `k`.
The total axes of the window will span the full axes of kernel `k` except for
when `k` was created with an `ExtensionNothing` in which case they will be
trimmed.
"""
function Window end

# AbstractArray interface

Base.axes(::Type{W}) where W<:Window =
    W.parameters[5] <: ExtensionArray &&
    W.parameters[5].parameters[4] != ExtensionNothing &&
    W.parameters[5].parameters[4] != ExtensionNone ?
        axes(W.parameters[4]) :
        W.parameters[3]

Base.axes(w::Window) = axes(typeof(w))
Base.size(w::Window) = length.(axes(w))

@propagate_inbounds Base.getindex(w::Window, wi::Int...) = getindex(w, CartesianIndex(wi))
@propagate_inbounds function Base.getindex(w::Window{<:Any,N}, wi::CartesianIndex{N}) where N
    # central piece to get efficient boundary handling.
    # we rely on the compiler to constant propagate this check away
    checkbounds_inner(Bool, w, wi) || return getindex_extension(w, wi)

    return parent(w)[position(w) + wi]
end

@propagate_inbounds Base.setindex!(w::Window, x, wi::Int...) = setindex!(w, x, CartesianIndex(wi))
@propagate_inbounds function Base.setindex!(w::Window{<:Any,N}, x, wi::CartesianIndex{N}) where N
    # central piece to get efficient boundary handling.
    # we rely on the compiler to constant propagate this check away
    checkbounds_inner(Bool, w, wi) || return setindex_extension!(w, x, wi)

    return parent(w)[position(w) + wi] = x
end

# Window interface

Base.eltype(::Type{<:Window{T}}) where T = T

"""
    checkbounds_inner(Bool, w::Window, i::CartesianIndex)

Return true if `i` indexes `w`'s parent in the interior and false if an
extension would be involved.
"""
@inline checkbounds_inner(::Type{Bool}, w::Window, i::CartesianIndex) =
    in(i, CartesianIndices(axes_inner(w)))

"""
    parent(w::Window)

Return reference to parent array.
"""
Base.parent(w::Window) = w.parent

"""
    axes_inner(w::Window)

Return the inner axes of the window where no extension handling is active.
"""
axes_inner(w::Window{<:Any,<:Any,X}) where X = X

"""
    position(w::Window)::CartesianIndex

Return the position of `w`, i.e. its center coordinate within the parent array.
"""
Base.position(w::Window) = w.position

"""
    Tuple(w::Window)

Create a tuple from the statically sized window `w`.

NOTE: this doesn't check bounds and thus assumes the window was properly
      created.
"""
@generated function Base.Tuple(w::Window)
    els = [ :( w[$i] ) for i in CartesianIndices(axes(w)) ]
    return :( @_inline_meta; @inbounds ($(els...),) )
end

## think twice before defining `iterate` like this, rather encourage use of
## mapreduce functions
#
#@generated function Base.iterate(w::Window, st = nothing)
#    ci = CartesianIndices(axes(w))
#    return quote
#        @_inline_meta
#        it = isnothing(st) ? iterate($ci) : iterate($ci, st)
#        isnothing(it) && return nothing
#        return @inbounds w[it[1]], it[2]
#    end
#end

# specialized mapfoldl, since Julia base has an upper limit on tuple size for
# inlining it

# since generated keyword functions may not inline, we forward to an internal
# generated non-keyword function
#   FIXME: file a julia bug
# this workaround does not seem to work when actually using non-default keyword
# arguments. explicitly specializing type arguments (F, Op) helps with the allocations
# but inlining still doesn't happen.
#   FIXME: file a julia bug

struct _DefaultValue end

@inline function Base.mapfoldl(f::F, op::Op, w::W...;
        dims = :, init = _DefaultValue()) where {F, Op, W<:Window}
    return _mapfoldl(f, op, dims, init, w...)
end

@generated function _mapfoldl(f, op, dims::Colon, init, w...)
    wax = same(axes, w...)
    weltype = same(eltype, w...)
    ws(i) = [ :( @inbounds w[$j][$i] ) for j in eachindex(w) ]
    stm_init = init === _DefaultValue ?
        :( acc = Base.mapreduce_empty(f, op, $weltype) ) :
        :( acc = init )
    stm_accs = [ :( acc = op(acc, f($(ws(i)...))) ) for i in CartesianIndices(wax) ]
    return :( @_inline_meta; $stm_init; $(stm_accs...); acc )
end

# Set of workarounds for not being able to subtype AbstractArray
#   see also `src/types.jl`
#   TODO: remove these as soon as Window <: AbstractArray

Base.ndims(w::Window) = length(axes(w))
Base.eltype(w::Window) = eltype(typeof(w))
Base.length(w::Window) = prod(size(w))
Base.CartesianIndices(w::Window) = CartesianIndices(axes(w))
Base.keys(w::Window) = CartesianIndices(w)
@inline Base.checkbounds(::Type{Bool}, w::Window, i::CartesianIndex) =
    in(i, keys(w))

# redefine inlined variations of various operations on arrays
# TODO: Julia base should define inlined dispatch function

@inline Base.mapreduce(f, op, w::Window...; dims = :, init = _DefaultValue()) =
    _mapfoldl(f, op, dims, init, w...)
@inline Base.in(x, w::Window) = _mapfoldl(==(x), |, :, _DefaultValue(), w)

@inline Base.sum(w::Window) = _mapfoldl(identity, +, :, _DefaultValue(), w)
@inline Base.sum(f, w::Window) = _mapfoldl(f, +, :, _DefaultValue(), w)
@inline Base.prod(w::Window) = _mapfoldl(identity, *, :, _DefaultValue(), w)
@inline Base.prod(f, w::Window) = _mapfoldl(f, *, :, _DefaultValue(), w)

# typemax/typemin usage here is not strictly compatible with julia base but
# simplifies the mapreduce implementation considerably
@inline Base.minimum(w::Window) = _mapfoldl(identity, min, :, typemax(eltype(w)), w)
@inline Base.minimum(f, w::Window) = _mapfoldl(f, min, :, typemax(eltype(w)), w)
@inline Base.maximum(w::Window) = _mapfoldl(identity, max, :, typemin(eltype(w)), w)
@inline Base.maximum(f, w::Window) = _mapfoldl(f, max, :, typemin(eltype(w)), w)

@inline Base.count(w::Window{Bool}) = _mapfoldl(identity, +, :, _DefaultValue(), w)
@inline Base.count(f, w::Window) = _mapfoldl(x->f(x)::Bool, +, :, _DefaultValue(), w)
@inline Base.all(w::Window{Bool}) = _mapfoldl(identity, &, :, _DefaultValue(), w)
@inline Base.all(f, w::Window) = _mapfoldl(x->f(x)::Bool, &, :, _DefaultValue(), w)
@inline Base.any(w::Window{Bool}) = _mapfoldl(identity, |, :, _DefaultValue(), w)
@inline Base.any(f, w::Window) = _mapfoldl(x->f(x)::Bool, |, :, _DefaultValue(), w)


text:
module TestLocalEval
using Base.Test
using Debug

@test_throws ErrorException @localscope

@debug function f(x)
    y = x+3
    (let y = -5, w = 4; y = -5; @localscope; end), @localscope
end

@debug @notrap function g(x)
    y = x+3
    (let y = -5, w = 4; y = -5; @localscope; end), @localscope
end

const a = 52

for func in [f, g]
    global z = 5

    si, so = func(11)
    @test so[:x] === 11 === si[:x]
    @test so[:y] === 14
    @test si[:y] === -5
    @test si[:w] === 4
    @test debug_eval(so, :(x, y, z, x*z)) === (11, 14, 5, 55)

    so[:x], so[:y], z = 1,2,3
    @test so[:x] === 1
    @test so[:y] === 2
    @test si[:y] === -5
    @test debug_eval(so, :(x, y, z, y*z)) === (1, 2, 3, 6)

    symso, symsi = Set(keys(so)), Set(keys(si))
    @test symso == Set([:x, :y])
    @test symsi == Set([:x, :y, :w])

    @test debug_eval(si, :a) === a
    # the scope only looks up local variables through indexing
    @test_throws ErrorException si[:a]
end

@debug @notrap scope = @localscope
@test debug_eval(scope, :a) === a
# the scope only looks up local variables through indexing
@test_throws ErrorException scope[:a]

end # module


text:
# returns only subtypes we can translate
# TODO (throughout) check which context should be used for this since generator is evaluated when the context of a module
function translatable_subtypes(dt::DataType)
	filter(subtypes(current_module(), dt)) do subdt
		# !(subdt <: DataGenerators.Generator) &&
		!(subdt <: Vararg)
	end
end

function transform_type_ast(ast::ASTNode)

	# extract primary datatypes and upper bound datatypes of typevars in the parsed type as well as supplemental types
	supporteddts = isempty(ast.args[:supporteddts]) ? extract_primary_datatypes(ast.args[:type]) : ast.args[:supporteddts]
	# remove Any from this list as tree is too big - instead user could specify required subtypes of Any in supplemental types
	# remove Vararg since it is handled in a special way within Tuple
	supporteddts = filter(t -> !(t in [Any; Vararg;]), supporteddts)
	# merge datatypes upwards (i.e. if A subsumes B, take just B)
	supporteddts = merge_datatypes_up(supporteddts)
	# this now gives a non-overlapping list of the datatypes we will support in the dt tree
	if isempty(supporteddts)
		error("No supportable datatypes passed in type nor supplemental types")
	end

	if !all(t -> t in DIRECTLY_SUPPORTED_TYPES, nonabstract_descendents(supporteddts; subtypefn=translatable_subtypes))
		# if some of the nonabstract descendent types of the supported types cannot be constructed directly
		# then we add Tuple and Type as supported types in order to support the signatures of constructor methods
		supporteddts = merge_datatypes_up([supporteddts; Tuple; Type])
	end

	valuenode = create_value_node()
	datatypenode = create_datatype_node()
	typenode = create_type_node()
	methodnode = create_method_node()
	dtrootnode = create_dt_root_node(supporteddts)
	ast.children = [valuenode; datatypenode; typenode; methodnode; dtrootnode;]

	add_reference(ast, :valueref, valuenode) do node
		(node.func == :method) ||
		((node.func == :cm) && !haskey(node.args, :method) && (node.args[:datatype] == Tuple))
	end

	add_reference(ast, :datatyperef, datatypenode) do node
		(node.func == :value) ||
		(node.func == :type) ||
		((node.func == :dt) && (node.args[:datatype] == Type)) ||
		((node.func == :cm) && !haskey(node.args, :method) && (node.args[:datatype] in [DataType; Union;]))
	end

	add_reference(ast, :typeref, typenode) do node
		(node.func == :datatype) ||
		(node.func == :method) ||
		((node.func == :dt) && (node.args[:datatype] == Type))
	end

	add_reference(ast, :dtref, dtrootnode) do node
		(node.func == :value) ||
		(node.func == :type)
	end

	add_reference(ast, :methodref, methodnode) do node
		(node.func == :cm) && haskey(node.args, :method)
	end

	add_choose(ast, :dt, :dtref) do node
		(node.func == :dt) && node.args[:abstract]
	end

	add_choose(ast, :cm, :cmref) do node
		(node.func == :dt) && !node.args[:abstract]
	end

  	push!(ast.refs, valuenode) # add type node as reference to root to enable reachability check

end

create_value_node() = ASTNode(:value)

create_datatype_node() = ASTNode(:datatype)

create_type_node() = ASTNode(:type)

create_method_node() = ASTNode(:method)

function create_dt_root_node(supporteddts::Vector{DataType})
	# supporteddttree is the minimal partial subtree of primary datatypes that includes the supported datatypes, and remains
	# rooted at Any
	# the handlable_subtypes function filters the subtypes to those we can handle
	supporteddttree = datatype_tree(supporteddts; subtypefn = translatable_subtypes)
	create_dt_node(Any, supporteddttree, supporteddts)
end

function create_dt_node(t::DataType, supporteddttree::Dict{DataType, Vector{DataType}}, supporteddts::Vector{DataType})
	primarydt = primary_datatype(t)
	primaryisabstract = is_abstract(primarydt)
	node = ASTNode(:dt)
	node.args[:name] = primarydt.name.name
	node.args[:datatype] = primarydt
	node.args[:abstract] = primaryisabstract
	@assert primaryisabstract == !isempty(supporteddttree[primarydt])
	if primaryisabstract
		primarysubtypes = supporteddttree[primarydt]
		node.children = map(st->create_dt_node(st, supporteddttree, supporteddts), primarysubtypes)
	else
		if !(primarydt in DIRECTLY_SUPPORTED_TYPES) # TODO could also add partially supported constructor methods even in the case of a translator constructed alternative
			append!(node.children, map(cm -> create_constructor_method_node(cm, primarydt), partially_supported_constructor_methods(primarydt, supporteddts)))
		end
		if isempty(node.children)
			push!(node.children, create_constructor_method_node(primarydt))
		end
	end
	node
end

function create_constructor_method_node(cm::SimplifiedMethod, primarydt::DataType)
	node = ASTNode(:cm)
	node.args[:method] = cm
	node.args[:datatype] = primarydt
	node
end

function create_constructor_method_node(primarydt::DataType)
	node = ASTNode(:cm)
	node.args[:datatype] = primarydt
	node
end

function add_choose(predicate::Function, node::ASTNode, childfunc::Symbol, reffunc::Symbol)
	for child in node.children
		add_choose(predicate, child, childfunc, reffunc)
	end
	if predicate(node)
		chooseablenodes = filter(child -> child.func == childfunc, node.children)
		if !isempty(chooseablenodes)
			choosenode = ASTNode(:choose)
			for chooseablenode in chooseablenodes
				refnode = ASTNode(reffunc)
				refnode.refs = [chooseablenode;]
				push!(choosenode.children, refnode)
			end
			push!(node.children, choosenode)
		end
	end
end

function add_reference(predicate::Function, node::ASTNode, func::Symbol, target::ASTNode)
	for child in node.children
		add_reference(predicate, child, func, target)
	end
	if predicate(node)
		refnode = ASTNode(func)
		refnode.refs = [target;]
		push!(node.children, refnode)
	end
end


text:
struct Policy
    a::Int
    b::Int
    letter::Char
end

struct Password
    value::AbstractString
    policy::Policy
end

"Parse a password and its policy from a string."
function parsepassword(s)
    rx = r"([0-9]+)-([0-9]+) ([a-zA-Z]): ([a-zA-Z]+)" # "1-3 a: abcdae"
    m = match(rx, s)
    pol = Policy(parse(Int, m.captures[1]), parse(Int, m.captures[2]), m.captures[3][1])
    pass = Password(m.captures[4], pol)
end

"Validate a password with its policy by letter count."
function validpasswordcount(p)
    rx = Regex(string(p.policy.letter))
    count = length(collect(eachmatch(rx, p.value)))
    count >= p.policy.a && count <= p.policy.b
end

"Validate password based on letter position."
function validpasswordpos(p)
    a = p.value[p.policy.a]
    b = p.value[p.policy.b]
    a_valid = a == p.policy.letter
    b_valid = b == p.policy.letter
    (a_valid || b_valid) && !(a_valid && b_valid)
end

parsepasswords(arr) = map(x->parsepassword(x), arr)

day02part1(arr) = arr |> parsepasswords |> x->filter(validpasswordcount, x) |> length

day02part2(arr) = arr |> parsepasswords |> x->filter(validpasswordpos, x) |> length

function day02(path)
    arr = filter(s->!isempty(s), split(read(path, String), "\n"))
    p1 = day02part1(arr)
    p2 = day02part2(arr)
    (p1, p2)
end

#=
rx = r"([0-9]+)-([0-9]+) (.): ([a-zA-Z]+)"
m = match(rx, "13-3 a: abcde")
m.captures

rx3 = r"a"
match(rx3, "aaaa").captures

rx2 = r"([^x]*x){20,23}"
match(rx2, "xxxxxxxxxxxxxxxxcxx")

p = parsepassword("2-3 a: abcdae")
validpasswordcount(p)
validpasswordpos(p)

arr = filter(s->!isempty(s), split(read("day02/input.txt", String), "\n"))
filter(x->!validpassword(x), parsepasswords(arr))

l = day02part1(arr)

day02("day02/input.txt")

t = """
1-3 a: abcde
1-3 b: cdefg
2-9 c: ccccccccc
"""
tarr = filter(s->!isempty(s), split(t, "\n"))
day02part1(tarr)

=#


text:
test = VERSION >= v"0.7-DEV" ? :(using Test) : :(using Base.Test)
eval(test)

using Logging

using ResourceBundles
import ResourceBundles: ROOT, BOTTOM, create_locid
import ResourceBundles: set_locale!, load_file

cd(abspath(pathof(ResourceBundles), "..", ".."))

# test if logger contains text in a message; finally with:clear logger.
function test_log(log::Test.TestLogger, text::AbstractString)
    try
        isempty(text) && return test_log(log)
        conmess(x::Test.LogRecord) = occursin(text, x.message)
        any(conmess.(test_log_filter(log)))
    finally
        empty!(log.logs)
    end
end

test_log(log::Test.TestLogger) = isempty(test_log_filter(log))
# ignore messages from Core module
test_log_filter(log::Test.TestLogger) = filter(lr->lr._module != Core, log.logs)

@testset "types" begin include("types.jl") end
@testset "locale" begin include("locale.jl") end
@testset "resource_bundle" begin include("resource_bundle.jl") end
@testset "string" begin include("string.jl") end
@testset "libc" begin include("libc.jl") end
# @testset "clocale" begin include("clocale.jl") end
# @testset "localetrans" begin include("localetrans.jl") end