text:
using SuiteSparseMatrixCollection
using MatrixMarket
using SuiteSparseGraphBLAS
SuiteSparseGraphBLAS.gbset(SuiteSparseGraphBLAS.FORMAT, SuiteSparseGraphBLAS.BYROW)
using BenchmarkTools
using SparseArrays
using LinearAlgebra
include("tc.jl")
include("pr.jl")
graphs = [
    #"karate",
    #"com-Youtube",
    #"as-Skitter",
    #"com-LiveJournal",
    #"com-Orkut",
    "com-Friendster",
]

ssmc = ssmc_db()
matrices = filter(row -> row.name ∈ graphs, ssmc)
BenchmarkTools.DEFAULT_PARAMETERS.gcsample = true
for name ∈ graphs
    path = fetch_ssmc(matrices[matrices.name .== name, :])[1]
    G = GBMatrix(convert(SparseMatrixCSC{Float64}, MatrixMarket.mmread(joinpath(path, "$name.mtx"))))
    SuiteSparseGraphBLAS.gbset(G, SuiteSparseGraphBLAS.FORMAT, SuiteSparseGraphBLAS.BYROW)
    GC.gc()
    G[:,:, mask=G, desc=SuiteSparseGraphBLAS.S] = 1
    diag(G)
    println("$name | $(size(G)) | $(nnz(G)) edges")
    d = reduce(+, G; dims=2)
    # for centrality in [PR, TC1, TC3]
    for centrality in [TC3]
        println("Benchmarking $(string(centrality)) on $(name)")
        i = 0.0
        centrality(G, d) #warmup
        for run ∈ 1:3
            i += @elapsed centrality(G, d)
        end
        println("$(string(centrality)) on $(name) over 3 runs took an average of: $(i / 3)s")
    end
end


text:
type FidMtx
     nz :: Int64
     nx :: Int64
     ext:: Int64
     iflag :: Int64
     dz :: Float64
     dx :: Float64
     dt :: Float64
     MvzBvz :: SparseMatrixCSC{Float64,Int64}
     MvzBp  :: SparseMatrixCSC{Float64,Int64}
     MvxBvx :: SparseMatrixCSC{Float64,Int64}
     MvxBp  :: SparseMatrixCSC{Float64,Int64}
     MpzBpz :: SparseMatrixCSC{Float64,Int64}
     MpzBvz :: SparseMatrixCSC{Float64,Int64}
     MpxBpx :: SparseMatrixCSC{Float64,Int64}
     MpxBvx :: SparseMatrixCSC{Float64,Int64}
end

function DispStable!(vmax::Float64, vmin::Float64, f0::Float64, dz::Float64, dx::Float64, dt::Float64)
    h = vmin / 5 / f0
    tt = 6*dx / (7*sqrt(3)*vmax)
    if dx >= h || dz >= h || dt >= tt
       println("maximum acceptable grid size: $h")
       println("maximum acceptable time step: $tt")
       error("unstable or frequency dispersion")
    end
    return nothing
end

function InitFidMtx(nz::Int64, nx::Int64, ext::Int64, iflag::Int64, dz::Float64, dx::Float64, dt::Float64, vmax::Float64, vmin::Float64, f0::Float64, v::Array{Float64,2})
    DispStable!(vmax, vmin, f0, dz, dx, dt)
    Cpml = InitCpml(nz, nx, ext, vmax, dz, dx, iflag)
    # Cpml = InitCpml(nz, nx, ext, vmax, rho, dz, dx, iflag)
    (MvzBvz, MvzBp)  = Mvz(Cpml.vz,    dz, dt, ext, iflag)
    (MvxBvx, MvxBp)  = Mvx(Cpml.vx,    dx, dt, ext, iflag)
    (MpzBpz, MpzBvz) = Mpz(Cpml.pz, v, dz, dt, ext, iflag)
    (MpxBpx, MpxBvx) = Mpx(Cpml.px, v, dx, dt, ext, iflag)
    fidMtx = FidMtx(nz, nx, ext, iflag, dz, dx, dt, MvzBvz, MvzBp, MvxBvx, MvxBp, MpzBpz, MpzBvz, MpxBpx, MpxBvx)
    return fidMtx
end

function InitFidMtx_back(nz::Int64, nx::Int64, ext::Int64, iflag::Int64, dz::Float64, dx::Float64, dt::Float64, vmax::Float64, vmin::Float64, f0::Float64, v::Array{Float64,2})
    DispStable!(vmax, vmin, f0, dz, dx, dt)
    Cpml = InitCpml(nz, nx, ext, vmax, dz, dx, iflag)
    # Cpml = InitCpml(nz, nx, ext, vmax, rho, dz, dx, iflag)
    (MvzBvz, MvzBp)  = Mvz_back(Cpml.vz,    dz, dt, ext, iflag)
    (MvxBvx, MvxBp)  = Mvx_back(Cpml.vx,    dx, dt, ext, iflag)
    (MpzBpz, MpzBvz) = Mpz_back(Cpml.pz, v, dz, dt, ext, iflag)
    (MpxBpx, MpxBvx) = Mpx_back(Cpml.px, v, dx, dt, ext, iflag)
    fidMtx = FidMtx(nz, nx, ext, iflag, dz, dx, dt, MvzBvz, MvzBp, MvxBvx, MvxBp, MpzBpz, MpzBvz, MpxBpx, MpxBvx)
    return fidMtx
end

# function InitFidMtx(nz::Int64, nx::Int64, ext::Int64, iflag::Int64, dz::Float64, dx::Float64, dt::Float64, vmax::Float64, vmin::Float64, f0::Float64, rho::Array{Float64,2}, lambda::Array{Float64,2})
#     DispStable!(vmax, vmin, f0, dz, dx, dt)
#     Cpml = InitCpml(nz, nx, ext, vmax, rho, iflag)
#     # Cpml = InitCpml(nz, nx, ext, vmax, rho, dz, dx, iflag)
#     (MvzBvz, MvzBp)  = Mvz(Cpml.vz, rho,    dz, dt, ext, iflag)
#     (MvxBvx, MvxBp)  = Mvx(Cpml.vx, rho,    dx, dt, ext, iflag)
#     (MpzBpz, MpzBvz) = Mpz(Cpml.pz, lambda, dz, dt, ext, iflag)
#     (MpxBpx, MpxBvx) = Mpx(Cpml.px, lambda, dx, dt, ext, iflag)
#     fidMtx = FidMtx(nz, nx, ext, iflag, dz, dx, dt, MvzBvz, MvzBp, MvxBvx, MvxBp, MpzBpz, MpzBvz, MpxBpx, MpxBvx)
#     return fidMtx
# end

function vp2lambda(vp::Array{Float64,2}, rho::Array{Float64,2})
    (m , n ) = size(vp)
    (m1, n1) = size(rho)
    if m1 != m || n1 != n
       error("size not match")
    end
    lambda = Array(typeof(vp[1]), m, n)
    for j = 1 : n
        for i = 1 : m
            lambda[i,j] = rho[i,j] * vp[i,j]^2
        end
    end
    return lambda
end

function lambda2vp(lambda::Array{Float64, 2}, rho::Array{Float64, 2})
    (m, n) = size(lambda)
    (m1, n1) = size(rho)
    if m1 != m || n1 != n
       error("size not match")
    end
    vp = Array(typeof(lambda[1]), m, n)
    for j = 1 : n
        for i = 1 : m
            vp[i,j] = sqrt(lambda[i,j] / rho[i,j])
        end
    end
    return vp
end

ext: jl

max_stars_repo_name: WneleiGao/AcousticWaveLab

max_issues_repo_name: WneleiGao/AcousticWaveLab

max_forks_repo_name: WneleiGao/AcousticWaveLab


text:
## Matlab author: Seongjin Yoon, Univ. of Michigan, Jan 2015
## Translated by: Harshit Nanda

export ir_radon_zwart_powell

"""
    output = ir_radon_zwart_powell(theta, rr)
Analytic 2D Radon transform of Zwart-Powell box spline. To use this
you call the theta matrix which represents the Zwart-Powell element
by the Box-Spline directions and the distance between ray and the point.

in
- `theta`         ray angle in radian
- `rr`            distance between the point and the ray (normalized by the pixel size)

out
- `output`        radon transform of the Zwart-Powell box spline element
"""
function ir_radon_zwart_powell(theta, rr)
    dim = size(theta)
    theta = vec(theta)
    zeta = zeros(length(theta), 4)

    # setting each dimension
    zeta[:,1] = cos.(theta)
    zeta[:,2] = sin.(theta)
    zeta[:,3] .= zeta[:,1] .+ zeta[:,2]
    zeta[:,4] .= zeta[:,1] .- zeta[:,2]

    # using matlab floating point value of eps(1.0)
    cond = abs.(zeta) .>= 1.1920929e-07

    N = sum(cond, dims = 2)
    N = vec(N)
    output = BoxSp4(rr[:], zeta, cond, N) ./ factorial.(N.-1)
    output = reshape(output, dim)
    return output
end

## BoxSp0
BoxSp0(y, N::Int) = y > 0 ? y^(N-1) : zero(y)

## BoxSp1
function BoxSp1(y, zeta, cond, N)
    good = cond[:,1]
    output = (BoxSp0.(y .+ (0.5 .* zeta[:,1]), N) .- BoxSp0.(y .- (0.5 .* zeta[:,1]), N)) ./ zeta[:,1]
    output[.~good] .= BoxSp0.(y[.~good], N[.~good])
    return output
end

## BoxSp2
function BoxSp2(y, zeta, cond, N)
    good = cond[:,2]

    # Average of second column
    output = (BoxSp1(y .+ (0.5.*zeta[:,2]), zeta, cond, N) .- BoxSp1(y .- (0.5 .* zeta[:,2]), zeta, cond, N)) ./ zeta[:,2]
    output[.~good] .= BoxSp1(y[.~good], zeta[.~good,:], cond[.~good,:], N[.~good])
    return output
end

## BoxSp3
function BoxSp3(y, zeta, cond, N)
    good = cond[:,3]

    # Average of third column
    output = (BoxSp2(y .+ 0.5 .* zeta[:,3], zeta, cond, N) .- BoxSp2(y .- 0.5 .* zeta[:,3], zeta, cond, N)) ./ zeta[:,3]
    output[.~good] .= BoxSp2(y[.~good], zeta[.~good,:], cond[.~good,:], N[.~good])
    return output
end

## BoxSp4
function BoxSp4(y, zeta, cond, N)
    good = cond[:,4]

    # Average of fourth column
    output = (BoxSp3(y .+ 0.5 .* zeta[:,4], zeta, cond, N) .- BoxSp3(y .- 0.5 .* zeta[:,4], zeta, cond, N)) ./ zeta[:,4]
    output[.~good] .= BoxSp3(y[.~good], zeta[.~good,:], cond[.~good,:], N[.~good])
    return output
end


text:
using FactCheck, LibBSON, Compat, DataStructures

facts("BSONOID") do
    oid = BSONOID()
    @fact length(string(oid)) --> 24
    @fact (oid == oid) --> true
    @fact typeof(hash(oid)) --> UInt
    oid2 = BSONOID()
    @fact (oid == oid2) --> false
    @fact string(BSONOID("540b628bca2e9b0d4e7dfa61")) --> "540b628bca2e9b0d4e7dfa61"
end

facts("BSONObject") do
    bsonObject = BSONObject(@compat Dict{Any,Any}(
        "null"=>nothing,
        "bool"=>true,
        "int"=>42,
        "double"=>3.141,
        "string"=>"Hello, Jérôme",
        "anotherNull"=>:null,
        "minkey"=>:minkey,
        "maxkey"=>:maxkey,
        "regularSymbol"=>:symbol,
        "subdict"=> (@compat Dict{Any,Any}(
            "key"=>"value"
            )),
        "array"=> Any["hello", (@compat Dict{Any,Any}("foo" => Any[56, false]))]
        ))
    @fact length(bsonObject) --> 11
    @fact string(bsonObject) --> "{ \"string\" : \"Hello, Jérôme\", \"anotherNull\" : null, \"null\" : null, \"regularSymbol\" : \"symbol\", \"bool\" : true, \"int\" : 42, \"minkey\" : { \"\$minKey\" : 1 }, \"maxkey\" : { \"\$maxKey\" : 1 }, \"array\" : [ \"hello\", { \"foo\" : [ 56, false ] } ], \"subdict\" : { \"key\" : \"value\" }, \"double\" : 3.141 }"
    @fact dict(bsonObject) --> @compat Dict{Any,Any}("string"=>"Hello, Jérôme","anotherNull"=>nothing,"null"=>nothing,"regularSymbol"=>"symbol","bool"=>true,"int"=>42,"minkey"=>:minkey,"maxkey"=>:maxkey,"array"=>Any["hello",(@compat Dict{Any,Any}("foo"=>Any[56,false]))],"subdict"=>(@compat Dict{Any,Any}("key"=>"value")),"double"=>3.141)
    append(bsonObject, "int64", -57)
    @fact bsonObject["int64"] --> -57
    append(bsonObject, "int32", 0x12345678)
    @fact bsonObject["int32"] --> 305419896

    context("BSONObject Copy Constructor") do
        initialDict = Dict{Any,Any}("someText" => "hello")
        bsonObject = BSONObject(@compat initialDict)
        bsonObjectCopyBad  = bsonObject
        bsonObjectCopyGood = BSONObject(bsonObject)
        append(bsonObject, "someValue", -123.456)
        @fact dict(bsonObject)          --> not(initialDict)
        @fact dict(bsonObjectCopyBad)   --> not(initialDict)
        @fact dict(bsonObjectCopyGood)  --> initialDict
    end

    context("BSONObject with OID") do
        oid = BSONOID()
        bsonObject = BSONObject(@compat Dict{Any,Any}("oid" => oid))
        @fact (length(string(bsonObject)) > 0) --> true
        @fact (bsonObject["oid"] == oid) --> true
    end

    context("BSONObject with DateTime") do
        ts = now()
        d = Date()
        obj = BSONObject(@compat Dict("time" => ts, "date" => d))
        @fact obj["time"] --> ts
        @fact obj["date"] --> d
    end

    context("BSONObject from JSON") do
        @fact_throws BSONObject("invalid JSON")
        @fact string(BSONObject("{\"pi\": 3.141}")) --> "{ \"pi\" : 3.141 }"
    end

    context("BSONObject containing BSONObject") do
        subBSONObject = BSONObject(@compat Dict{Any,Any}("key"=>"value"))
        bsonObject = BSONObject(@compat Dict{Any,Any}("sub"=>subBSONObject))
        @fact string(bsonObject) --> "{ \"sub\" : { \"key\" : \"value\" } }"
        @fact string(bsonObject["sub"]) --> "{ \"key\" : \"value\" }"
    end
end

facts("BSONArray") do
    bsonArray = BSONArray([
        nothing,
        true,
        42,
        3.141,
        "Hello, Jérôme",
        :null,
        :minkey,
        :maxkey,
        :symbol
        ])
    @fact length(bsonArray) --> 9
    @fact string(bsonArray) --> "[ null, true, 42, 3.141, \"Hello, Jérôme\", null, { \"\$minKey\" : 1 }, { \"\$maxKey\" : 1 }, \"symbol\" ]"
    @fact vector(bsonArray) --> Any[nothing,true,42,3.141,"Hello, Jérôme",nothing,:minkey,:maxkey,"symbol"]
    append(bsonArray, BSONArray([false]))
    append(bsonArray, -67)
    append(bsonArray, ["hello", 6.7])
    @fact vector(bsonArray) --> Any[nothing,true,42,3.141,"Hello, Jérôme",nothing,:minkey,:maxkey,"symbol",Any[false],-67,["hello",6.7]]

    bsonArray = BSONArray()
    append(bsonArray, BSONOID())
    @fact length(bsonArray) --> 1

    context("Array of DateTimes") do
      ts = now()
      d = Date()
      arr = BSONArray([d,ts])
      @fact length(arr) --> 2
      @fact arr[1] --> d
      @fact arr[2] --> ts
      for i = 1:length(arr)
        @fact typeof(arr[i]) <: DateTime --> true
      end
    end
end

facts("BSONObject: get!") do
    bsonObject = BSONObject()
    @fact_throws bsonObject["foo"]
    @fact get!(bsonObject, "foo", true) --> true
    @fact get!(bsonObject, "foo", true) --> true
end

facts("BSONError") do
    err = BSONError()
    @fact typeof(string(err)) <: AbstractString --> true
end

facts("Issue 18") do
    @fact string(BSONObject(OrderedDict("a" => OrderedDict("b" => "c")))) --> "{ \"a\" : { \"b\" : \"c\" } }"
end


text:
function pull_request_build(::NewPackage,
                            pr::GitHub.PullRequest,
                            current_pr_head_commit_sha::String,
                            registry::GitHub.Repo;
                            auth::GitHub.Authorization,
                            authorized_authors::Vector{String},
                            registry_head::String,
                            registry_master::String,
                            suggest_onepointzero::Bool,
                            whoami::String)
    # first check if the PR is open, and the author is authorized - if not, then quit
    # then, delete ALL reviews by me
    # then check rules 1-5. if fail, post comment.
    # then check rules 6-7. if fail, post comment.
    # if everything passed, add an approving review by me
    # 1. Normal capitalization - name should match r"^[A-Z]\w*[a-z][0-9]?$" - i.e. starts with a capital letter, ASCII alphanumerics only, ends in lowercase
    # 2. Not too short - at least five letters - you can register names shorter than this, but doing so requires someone to approve
    # 3. Standard initial version number - one of 0.0.1, 0.1.0, 1.0.0
    # 4. Repo URL ends with /$name.jl.git where name is the package name
    # 5. Compat for all dependencies - all [deps] should also have [compat] entries (and Julia itself) - [compat] entries should have upper bounds
    # 6. Version can be installed - given the proposed changes to the registry, can we resolve and install the new version of the package?
    # 7. Version can be loaded - once it's been installed (and built?), can we load the code?
    pkg, version = parse_pull_request_title(NewPackage(), pr)
    @info("This is a new package pull request", pkg, version)
    pr_author_login = author_login(pr)
    if is_open(pr)
        if pr_author_login in authorized_authors
            my_retry(() -> delete_all_of_my_reviews!(registry, pr; auth = auth, whoami = whoami))
            my_retry(() -> GitHub.create_status(registry, current_pr_head_commit_sha; auth=auth, params=Dict("state" => "pending", "context" => "automerge/decision", "description" => "New package. Pending.")))
            g0, m0 = pr_only_changes_allowed_files(NewPackage(), registry, pr, pkg; auth = auth)
            newp_g1, newp_m1 = meets_normal_capitalization(pkg)
            newp_g2, newp_m2 = meets_name_length(pkg)
            newp_g3, newp_m3 = meets_standard_initial_version_number(version)
            newp_g4, newp_m4 = meets_repo_url_requirement(pkg;
                                                          registry_head = registry_head)
            newp_g5, newp_m5 = meets_compat_for_all_deps(registry_head,
                                                         pkg,
                                                         version)
            newp_g0through5 = [g0, newp_g1, newp_g2, newp_g3, newp_g4, newp_g5]
            @info("Only modifies the files that it's allowed to modify", meets_this_guideline = g0, message = m0)
            @info("Normal capitalization", meets_this_guideline = newp_g1, message = newp_m1)
            @info("Name not too short", meets_this_guideline = newp_g2, message = newp_m2)
            @info("Standard initial version number ", meets_this_guideline = newp_g3, message = newp_m3)
            @info("Repo URL ends with /name.jl.git", meets_this_guideline = newp_g4, message = newp_m4)
            @info("Compat (with upper bound) for all dependencies", meets_this_guideline = newp_g5, message = newp_m5)
            if all(newp_g0through5)
                newp_g6and7, newp_m6and7 = meets_version_can_be_loaded(registry_head,
                                                                         pkg,
                                                                         version)
                @info("Version can be installed and loaded", meets_this_guideline = newp_g6and7, message = newp_m6and7)
                if newp_g6and7
                    newp_commenttextpass = comment_text_pass(NewPackage(),
                                                             suggest_onepointzero,
                                                             version)
                    my_retry(() -> delete_all_of_my_reviews!(registry, pr; auth = auth, whoami = whoami))
                    my_retry(() -> approve!(registry, pr, current_pr_head_commit_sha; auth = auth, body = newp_commenttextpass, whoami=whoami))
                    my_retry(() -> GitHub.create_status(registry, current_pr_head_commit_sha; auth=auth, params=Dict("state" => "success", "context" => "automerge/decision", "description" => "New package. Approved. sha=\"$(current_pr_head_commit_sha)\"")))
                    return nothing
                else
                    newp_commenttext6and7 = comment_text_fail(NewPackage(),
                                                              [newp_m6and7],
                                                              suggest_onepointzero,
                                                              version)
                    my_retry(() -> post_comment!(registry, pr, newp_commenttext6and7; auth = auth))
                    my_retry(() -> GitHub.create_status(registry, current_pr_head_commit_sha; auth=auth, params=Dict("state" => "failure", "context" => "automerge/decision", "description" => "New package. Failed.")))
                    error("The automerge guidelines were not met.")
                    return nothing
                end
            else
                newp_allmessages0through5 = [m0, newp_m1, newp_m2, newp_m3, newp_m4, newp_m5]
                newp_failingmessages0through5 = newp_allmessages0through5[.!newp_g0through5]
                newp_commenttext0through5 = comment_text_fail(NewPackage(),
                                                              newp_failingmessages0through5,
                                                              suggest_onepointzero,
                                                              version)
                my_retry(() -> post_comment!(registry, pr, newp_commenttext0through5; auth = auth))
                my_retry(() -> GitHub.create_status(registry, current_pr_head_commit_sha; auth=auth, params=Dict("state" => "failure", "context" => "automerge/decision", "description" => "New package. Failed.")))
                error("The automerge guidelines were not met.")
                return nothing
            end
        else
            @info("Author $(pr_author_login) is not authorized to automerge. Exiting...")
            return nothing
        end
    else
        @info("The pull request is not open. Exiting...")
        return nothing
    end
    return nothing
end


text:
function FillCoefficients!(I::Vector{Int64}, J::Vector{Int64}, V::Vector{Float64}, istart::Int64, icoef::SubArray{Int64}, jcoef::SubArray{Int64}, vcoef::SubArray{Float64})
    @tturbo for i=1:length(icoef)
        I[istart+i] = icoef[i]
        J[istart+i] = jcoef[i]
        V[istart+i] = vcoef[i]
    end
end
export FillCoefficients!

###########


text:
# Naming convention:
# ------------------
# {type}[S]_{dims}[_d][_n]
# type: C: Continuous, D: Discrete
# scalar: add "S" if any input is scalar
# dims: "nxnuny"
# feedthrough: append "_d" if `D` is present
# names: append "_n" if some inputs/outputs/states are named

@testset "test_statespace" begin
# SCALARS
a_2 = [-5 -3; 2 -9]
CS_111 = ss(-5, 2, 3, [0])
CS_111_d = ss([3], 2, 1, 1)
CS_211 = ss(a_2, [1; 2], [1 0], 0)
CS_221 = ss(a_2, [1 0; 0 2], [1 0], 0)
CS_222 = ss(a_2, [1 0; 0 2], eye_(2), 0)

# CONTINUOUS
a_1 = [-5]
C_111 = ss(a_1, [2], [3], [0])
C_211 = ss(a_2, [1; 2], [1 0], [0])
C_212 = ss(a_2, [1; 2], eye_(2), [0; 0])
C_221 = ss(a_2, [1 0; 0 2], [1 0], [0 0])
C_222 = ss(a_2, [1 0; 0 2], eye_(2), zeros(2,2))
C_222_d = ss(a_2, [1 0; 0 2], eye_(2), eye_(2))
C_022 = ss(4.0*eye_(2))

# DISCRETE
da_1 = [-0.5]
da_2 = [0.2 -0.8; -0.8 0.07]
D_111 = ss(da_1, [2], [3], [0], 0.005)
D_211 = ss(da_2, [1; 2], [1 0], [0], 0.005)
D_221 = ss(da_2, [1 0; 0 2], [1 0], [0 0], 0.005)
D_222 = ss(da_2, [1 0; 0 2], eye_(2), zeros(2,2), 0.005)
D_222_d = ss(da_2, [1 0; 0 2], eye_(2), eye_(2), 0.005)
D_022 = ss(4.0*eye_(2), 0.005)

# Definition of input, output and state names
C_222_d_n = ss(a_2, [1 0; 0 2], eye_(2), eye_(2))

# TESTS
# Contstuct with scalars
@test CS_111 == C_111
@test CS_111_d == ss([3],[2],[1],[1])
@test CS_211 == C_211
@test CS_221 == C_221
@test CS_222 == C_222

# Addition
@test C_111 + C_111 == ss([-5 0; 0 -5],[2; 2],[3 3],[0])
@test C_222 + C_222 == ss([-5 -3 0 0; 2 -9 0 0; 0 0 -5 -3;
        0 0 2 -9],[1 0; 0 2; 1 0; 0 2], [1 0 1 0; 0 1 0 1],[0 0; 0 0])
@test C_222 + 1 == ss([-5 -3; 2 -9],[1 0; 0 2],[1 0; 0 1],[1 1; 1 1])
@test D_111 + D_111 == ss([-0.5 0; 0 -0.5],[2; 2],[3 3],[0], 0.005)

# Subtraction
@test C_111 - C_211 == ss([-5 0 0; 0 -5 -3; 0 2 -9],[2; 1; 2],[3 -1 -0],[0])
@test 1 - C_222 == ss([-5 -3; 2 -9],[1 0; 0 2],[-1 -0; -0 -1],[1 1; 1 1])
@test D_111 - D_211 == ss([-0.5 0 0; 0 0.2 -0.8; 0 -0.8 0.07],[2; 1; 2],
        [3 -1 -0],[0], 0.005)

# Multiplication
@test C_111 * C_221 == ss([-5 2 0; 0 -5 -3; 0 2 -9],
        [0 0; 1 0; 0 2],[3 0 0],[0 0])
@test C_212 * C_111 == ss([-5 -3 3; 2 -9 6; 0 0 -5],
        [0; 0; 2],[1 0 0; 0 1 0],[0; 0])
@test 4*C_222 == ss([-5 -3; 2 -9],[1 0; 0 2],[4 0; 0 4],[0 0; 0 0])
@test D_111 * D_221 == ss([-0.5 2 0; 0 0.2 -0.8; 0 -0.8 0.07],
        [0 0; 1 0; 0 2],[3 0 0],[0 0],0.005)

# Division
@test 1/C_222_d == ss([-6 -3; 2 -11],[1 0; 0 2],[-1 0; -0 -1],[1 -0; 0 1])
@test C_221/C_222_d == ss([-5 -3 -1 0; 2 -9 -0 -2; 0 0 -6 -3;
        0 0 2 -11],[1 0; 0 2; 1 0; 0 2],[1 0 0 0],[0 0])
@test 1/D_222_d == ss([-0.8 -0.8; -0.8 -1.93],[1 0; 0 2],[-1 0; -0 -1],
        [1 -0; 0 1],0.005)

# Indexing
@test size(C_222) == (2, 2)
@test size(C_212) == (2, 1)
@test C_222[1,1] == ss([-5 -3; 2 -9],[1; 0],[1 0],[0])
@test C_222[1:1,1] == ss([-5 -3; 2 -9],[1; 0],[1 0],[0])
@test C_222[1,1:2] == C_221
@test size(C_222[1,[]]) == (1,0)

A = [-1.0 -2.0; 0.0 -1.0]
B = [0.0; -2.0]
C = [1.0 1.0]
D = 1.0
sys = ss(A, B, C, D)

@test sys + 1.0 == ss(A, B, C, D + 1.0)
@test 2.0 + sys == ss(A, B, C, D + 2.0)

@test -sys == ss(A, B, -C, -D)

# Printing
res = ("StateSpace{Float64,Array{Float64,2}}\nA = \n -5.0  -3.0\n  2.0  -9.0\nB = \n 1.0  0.0\n 0.0  2.0\nC = \n 1.0  0.0\n 0.0  1.0\nD = \n 0.0  0.0\n 0.0  0.0\n\nContinuous-time state-space model")
@test sprint(show, C_222) == res
res = ("StateSpace{Float64,Array{Float64,2}}\nA = \n  0.2  -0.8 \n -0.8   0.07\nB = \n 1.0  0.0\n 0.0  2.0\nC = \n 1.0  0.0\n 0.0  1.0\nD = \n 0.0  0.0\n 0.0  0.0\n\nSample Time: 0.005 (seconds)\nDiscrete-time state-space model")
@test sprint(show, D_222) == res
res = ("StateSpace{Float64,Array{Float64,2}}\nD = \n 4.0  0.0\n 0.0  4.0\n\nStatic gain")
@test sprint(show, C_022) == res
res = "StateSpace{Float64,Array{Float64,2}}\nD = \n 4.0  0.0\n 0.0  4.0\n\nSample Time: 0.005 (seconds)\nStatic gain"
@test sprint(show, D_022) == res

# Errors
@test_throws ErrorException C_111 + C_222             # Dimension mismatch
@test_throws ErrorException C_111 - C_222             # Dimension mismatch
@test_throws ErrorException C_111 * C_222             # Dimension mismatch
@test_throws ErrorException D_111 + C_111             # Sampling time mismatch
@test_throws ErrorException D_111 - C_111             # Sampling time mismatch
@test_throws ErrorException D_111 * C_111             # Sampling time mismatch
D_diffTs = ss([1], [2], [3], [4], 0.1)
@test_throws ErrorException D_111 + D_diffTs            # Sampling time mismatch
@test_throws ErrorException D_111 - D_diffTs            # Sampling time mismatch
@test_throws ErrorException D_111 * D_diffTs            # Sampling time mismatch
@test_throws ErrorException 1/C_222                     # Not invertible
@test_throws ErrorException 1/C_212                     # Not invertible
@test_throws ErrorException ss([1 2], [1], [2], [3])      # Not square A
@test_throws ErrorException ss([1], [2 0], [1], [2])      # I/0 dim mismatch
@test_throws ErrorException ss([1], [2], [3 4], [1])      # I/0 dim mismatch
@test_throws ErrorException ss([1], [2], [3], [4], -0.1)  # Negative samping time
@test_throws ErrorException ss(eye_(2), eye_(2), eye_(2), [0]) # Dimension mismatch
end


text:
module ICD10Utilities

using CSV, TypedTables, Dates, CategoricalArrays, Missings
using FileIO

import Base: isless, show, (==)

export ICDOPTS

export AbstractICD10
export ICD10
export ICD10AM, ACHI
export ICD10CA, ICD10CM, ICD10GM
export ICD10CM
export ICD10AMAge
export icd3
export isvalidcode
export icd3, icd4
export seticdpunct
export ICDO3
export ICD10AMcodes, ACHIcodes

"""
    ICDOPTS

Dictionary of options for ICD10Utilities.

`:punct`: if `true` (default), ICD-10 codes are displayed with `.` after third character,
and ACHI codes are displayed with `-` after fifth character.
"""
const ICDOPTS = Dict{Symbol,Any}(:punct => true)
# easily set options
"""
    seticdpunct(v::Bool)

Set `ICDOPTS[:punct]` to new value.
"""
seticdpunct(v::Bool) = ICDOPTS[:punct] = v

"""
    AbstractICD10

Abstract type for ICD-10 codes
"""
abstract type AbstractICD10 end

Broadcast.broadcastable(icdcode::T) where {T<:AbstractICD10} = Ref(icdcode)

include("icd10.jl")
include("icd10fns.jl")
include("othericd10.jl")
include("icd10cm.jl")
include("achi.jl")
include("icd10amfns.jl")
include("icdO3.jl")

end # module


text:
"""
    Chain(layers...)

Chain multiple layers / functions together, so that they are called in sequence
on a given input.

`Chain` also supports indexing and slicing, e.g. `m[2]` or `m[1:end-1]`.
`m[1:3](x)` will calculate the output of the first three layers.

# Examples
```jldoctest
julia> m = Chain(x -> x^2, x -> x+1);

julia> m(5) == 26
true

julia> m = Chain(Dense(10, 5), Dense(5, 2));

julia> x = rand(10);

julia> m(x) == m[2](m[1](x))
true
```
"""
struct Chain{T<:Tuple}
  layers::T
  Chain(xs...) = new{typeof(xs)}(xs)
end

@forward Chain.layers Base.getindex, Base.length, Base.first, Base.last,
  Base.iterate, Base.lastindex

functor(::Type{<:Chain}, c) = c.layers, ls -> Chain(ls...)

applychain(::Tuple{}, x) = x
applychain(fs::Tuple, x) = applychain(tail(fs), first(fs)(x))

(c::Chain)(x) = applychain(c.layers, x)

Base.getindex(c::Chain, i::AbstractArray) = Chain(c.layers[i]...)

testmode!(m::Chain, mode = true) = (map(x -> testmode!(x, mode), m.layers); m)

function Base.show(io::IO, c::Chain)
  print(io, "Chain(")
  join(io, c.layers, ", ")
  print(io, ")")
end

"""
    outdims(c::Chain, isize)

Calculate the output dimensions given the input dimensions, `isize`.

```julia
m = Chain(Conv((3, 3), 3 => 16), Conv((3, 3), 16 => 32))
outdims(m, (10, 10)) == (6, 6)
```
"""
outdims(c::Chain, isize) = foldl(∘, map(l -> (x -> outdims(l, x)), c.layers))(isize)

# This is a temporary and naive implementation
# it might be replaced in the future for better performance
# see issue https://github.com/FluxML/Flux.jl/issues/702
# Johnny Chen -- @johnnychen94
# only slightly changed to better handle interaction with Zygote @dsweber2
"""
    activations(c::Chain, input)

Calculate the forward results of each layers in Chain `c` with `input` as model input.
"""
function activations(c::Chain, input)
    extraChain(c.layers, input)
end

function extraChain(fs::Tuple, x)
    res = first(fs)(x)
    return (res, extraChain(Base.tail(fs), res)...)
end

extraChain(::Tuple{}, x) = ()

"""
    Dense(in::Integer, out::Integer, σ = identity)

Create a traditional `Dense` layer with parameters `W` and `b`.

    y = σ.(W * x .+ b)

The input `x` must be a vector of length `in`, or a batch of vectors represented
as an `in × N` matrix. The out `y` will be a vector or batch of length `out`.

# Examples
```jldoctest; setup = :(using Random; Random.seed!(0))
julia> d = Dense(5, 2)
Dense(5, 2)

julia> d(rand(5))
2-element Array{Float32,1}:
  -0.16210233
   0.12311903```
"""
struct Dense{F,S,T}
  W::S
  b::T
  σ::F
end

Dense(W, b) = Dense(W, b, identity)

function Dense(in::Integer, out::Integer, σ = identity;
               initW = glorot_uniform, initb = zeros)
  return Dense(initW(out, in), initb(out), σ)
end

@functor Dense

function (a::Dense)(x::AbstractArray)
  W, b, σ = a.W, a.b, a.σ
  σ.(W*x .+ b)
end

function Base.show(io::IO, l::Dense)
  print(io, "Dense(", size(l.W, 2), ", ", size(l.W, 1))
  l.σ == identity || print(io, ", ", l.σ)
  print(io, ")")
end

# Try to avoid hitting generic matmul in some simple cases
# Base's matmul is so slow that it's worth the extra conversion to hit BLAS
(a::Dense{<:Any,W})(x::AbstractArray{T}) where {T <: Union{Float32,Float64}, W <: AbstractArray{T}} =
  invoke(a, Tuple{AbstractArray}, x)

(a::Dense{<:Any,W})(x::AbstractArray{<:AbstractFloat}) where {T <: Union{Float32,Float64}, W <: AbstractArray{T}} =
  a(T.(x))

"""
    outdims(l::Dense, isize)

Calculate the output dimensions given the input dimensions, `isize`.

```julia
m = Dense(10, 5)
outdims(m, (5, 2)) == (5,)
outdims(m, (10,)) == (5,)
```
"""
outdims(l::Dense, isize) = (size(l.W)[1],)

"""
    Diagonal(in::Integer)

Create an element-wise linear transformation layer with learnable
vectors `α` and `β`:

    y = α .* x .+ β

The input `x` must be a array where `size(x, 1) == in`.
"""
struct Diagonal{T}
  α::T
  β::T
end

Diagonal(in::Integer; initα = ones, initβ = zeros) =
  Diagonal(initα(in), initβ(in))

@functor Diagonal

function (a::Diagonal)(x)
  α, β = a.α, a.β
  α.*x .+ β
end

function Base.show(io::IO, l::Diagonal)
  print(io, "Diagonal(", length(l.α), ")")
end

outdims(l::Diagonal, isize) = (length(l.α),)

"""
    Maxout(over)

The [Maxout](https://arxiv.org/pdf/1302.4389.pdf) layer has a number of
internal layers which all receive the same input. It returns the elementwise
maximum of the internal layers' outputs.

Maxout over linear dense layers satisfies the univeral approximation theorem.
"""
struct Maxout{FS<:Tuple}
    over::FS
end

"""
    Maxout(f, n_alts)

Construct a Maxout layer over `n_alts` instances of the layer given by `f`.
The function takes no arguments and should return some callable layer.
Conventionally, this is a linear dense layer.

# Examples

This constructs a `Maxout` layer over 4 internal dense linear layers, each
identical in structure (784 inputs, 128 outputs):
```julia
insize = 784
outsize = 128
Maxout(()->Dense(insize, outsize), 4)
```
"""
function Maxout(f, n_alts)
  over = Tuple(f() for _ in 1:n_alts)
  return Maxout(over)
end

@functor Maxout

function (mo::Maxout)(input::AbstractArray)
    mapreduce(f -> f(input), (acc, out) -> max.(acc, out), mo.over)
end

outdims(l::Maxout, isize) = outdims(first(l.over), isize)

"""
    SkipConnection(layer, connection)

Create a skip connection which consists of a layer or `Chain` of consecutive
layers and a shortcut connection linking the block's input to the output
through a user-supplied 2-argument callable. The first argument to the callable
will be propagated through the given `layer` while the second is the unchanged,
"skipped" input.

The simplest "ResNet"-type connection is just `SkipConnection(layer, +)`,
and requires the output of the layers to be the same shape as the input.
Here is a more complicated example:
```julia
m = Conv((3,3), 4=>7, pad=(1,1))
x = ones(5,5,4,10);
size(m(x)) == (5, 5, 7, 10)

sm = SkipConnection(m, (mx, x) -> cat(mx, x, dims=3))
size(sm(x)) == (5, 5, 11, 10)
```
"""
struct SkipConnection
  layers
  connection  #user can pass arbitrary connections here, such as (a,b) -> a + b
end

@functor SkipConnection

function (skip::SkipConnection)(input)
  skip.connection(skip.layers(input), input)
end

function Base.show(io::IO, b::SkipConnection)
  print(io, "SkipConnection(", b.layers, ", ", b.connection, ")")
end


text:
#=
    NEMO: Next Energy Modeling system for Optimization.
    https://github.com/sei-international/NemoMod.jl

    Copyright © 2019: Stockholm Environment Institute U.S.

	File description: Tests of NemoMod package using Mosek solver.
=#

# Tests will be skipped if MosekTools package is not installed.
try
    using MosekTools
catch e
    @info "Error when initializing Mosek. Error message: " * sprint(showerror, e) * "."
    @info "Skipping Mosek tests."
end

if @isdefined MosekTools
    @info "Testing scenario solution with Mosek."

    @testset "Solving storage_test with Mosek" begin
        dbfile = joinpath(@__DIR__, "storage_test.sqlite")
        #dbfile = "c:/temp/storage_test.sqlite"
        chmod(dbfile, 0o777)  # Make dbfile read-write. Necessary because after Julia 1.0, Pkg.add makes all package files read-only

        # Test with default outputs
        NemoMod.calculatescenario(dbfile; jumpmodel = JuMP.Model(Mosek.Optimizer),
            numprocs=1, restrictvars=false, quiet = false)

        db = SQLite.DB(dbfile)
        testqry = SQLite.DBInterface.execute(db, "select * from vtotaldiscountedcost") |> DataFrame

        @test testqry[1,:y] == "2020"
        @test testqry[2,:y] == "2021"
        @test testqry[3,:y] == "2022"
        @test testqry[4,:y] == "2023"
        @test testqry[5,:y] == "2024"
        @test testqry[6,:y] == "2025"
        @test testqry[7,:y] == "2026"
        @test testqry[8,:y] == "2027"
        @test testqry[9,:y] == "2028"
        @test testqry[10,:y] == "2029"

        @test isapprox(testqry[1,:val], 3845.15702798338; atol=TOL)
        @test isapprox(testqry[2,:val], 146.552320420142; atol=TOL)
        @test isapprox(testqry[3,:val], 139.573628995129; atol=TOL)
        @test isapprox(testqry[4,:val], 132.927265707472; atol=TOL)
        @test isapprox(testqry[5,:val], 126.597395911861; atol=TOL)
        @test isapprox(testqry[6,:val], 120.568948035628; atol=TOL)
        @test isapprox(testqry[7,:val], 114.827569988096; atol=TOL)
        @test isapprox(testqry[8,:val], 109.359590464853; atol=TOL)
        @test isapprox(testqry[9,:val], 104.151990918907; atol=TOL)
        @test isapprox(testqry[10,:val], 99.1907267770023; atol=TOL)

        # Test with optional outputs
        NemoMod.calculatescenario(dbfile; jumpmodel = JuMP.Model(Mosek.Optimizer),
            varstosave =
                "vrateofproductionbytechnologybymode, vrateofusebytechnologybymode, vrateofdemand, vproductionbytechnology, vtotaltechnologyannualactivity, "
                * "vtotaltechnologymodelperiodactivity, vusebytechnology, vmodelperiodcostbyregion, vannualtechnologyemissionpenaltybyemission, "
                * "vtotaldiscountedcost",
            numprocs=1, restrictvars=false, quiet = false)

        testqry = SQLite.DBInterface.execute(db, "select * from vtotaldiscountedcost") |> DataFrame

        @test testqry[1,:y] == "2020"
        @test testqry[2,:y] == "2021"
        @test testqry[3,:y] == "2022"
        @test testqry[4,:y] == "2023"
        @test testqry[5,:y] == "2024"
        @test testqry[6,:y] == "2025"
        @test testqry[7,:y] == "2026"
        @test testqry[8,:y] == "2027"
        @test testqry[9,:y] == "2028"
        @test testqry[10,:y] == "2029"

        @test isapprox(testqry[1,:val], 3845.15702798338; atol=TOL)
        @test isapprox(testqry[2,:val], 146.552320420142; atol=TOL)
        @test isapprox(testqry[3,:val], 139.573628995129; atol=TOL)
        @test isapprox(testqry[4,:val], 132.927265707472; atol=TOL)
        @test isapprox(testqry[5,:val], 126.597395911861; atol=TOL)
        @test isapprox(testqry[6,:val], 120.568948035628; atol=TOL)
        @test isapprox(testqry[7,:val], 114.827569988096; atol=TOL)
        @test isapprox(testqry[8,:val], 109.359590464853; atol=TOL)
        @test isapprox(testqry[9,:val], 104.151990918907; atol=TOL)
        @test isapprox(testqry[10,:val], 99.1907267770023; atol=TOL)

        # Test with restrictvars
        NemoMod.calculatescenario(dbfile; jumpmodel = JuMP.Model(Mosek.Optimizer),
            varstosave = "vrateofproductionbytechnologybymode, vrateofusebytechnologybymode, vproductionbytechnology, vusebytechnology, "
                * "vtotaldiscountedcost",
            targetprocs=[1], restrictvars = true, quiet = false)

        testqry = SQLite.DBInterface.execute(db, "select * from vtotaldiscountedcost") |> DataFrame

        @test testqry[1,:y] == "2020"
        @test testqry[2,:y] == "2021"
        @test testqry[3,:y] == "2022"
        @test testqry[4,:y] == "2023"
        @test testqry[5,:y] == "2024"
        @test testqry[6,:y] == "2025"
        @test testqry[7,:y] == "2026"
        @test testqry[8,:y] == "2027"
        @test testqry[9,:y] == "2028"
        @test testqry[10,:y] == "2029"

        @test isapprox(testqry[1,:val], 3845.15702798338; atol=TOL)
        @test isapprox(testqry[2,:val], 146.552320420142; atol=TOL)
        @test isapprox(testqry[3,:val], 139.573628995129; atol=TOL)
        @test isapprox(testqry[4,:val], 132.927265707472; atol=TOL)
        @test isapprox(testqry[5,:val], 126.597395911861; atol=TOL)
        @test isapprox(testqry[6,:val], 120.568948035628; atol=TOL)
        @test isapprox(testqry[7,:val], 114.827569988096; atol=TOL)
        @test isapprox(testqry[8,:val], 109.359590464853; atol=TOL)
        @test isapprox(testqry[9,:val], 104.151990918907; atol=TOL)
        @test isapprox(testqry[10,:val], 99.1907267770023; atol=TOL)

        # Test with storage net zero constraints
        SQLite.DBInterface.execute(db, "update STORAGE set netzeroyear = 1")
        NemoMod.calculatescenario(dbfile; jumpmodel = JuMP.Model(Mosek.Optimizer), restrictvars=false, numprocs=1)
        testqry = SQLite.DBInterface.execute(db, "select * from vtotaldiscountedcost") |> DataFrame

        @test testqry[1,:y] == "2020"
        @test testqry[2,:y] == "2021"
        @test testqry[3,:y] == "2022"
        @test testqry[4,:y] == "2023"
        @test testqry[5,:y] == "2024"
        @test testqry[6,:y] == "2025"
        @test testqry[7,:y] == "2026"
        @test testqry[8,:y] == "2027"
        @test testqry[9,:y] == "2028"
        @test testqry[10,:y] == "2029"

        @test isapprox(testqry[1,:val], 3840.94023817774; atol=TOL)
        @test isapprox(testqry[2,:val], 459.292614204884; atol=TOL)
        @test isapprox(testqry[3,:val], 437.421537337638; atol=TOL)
        @test isapprox(testqry[4,:val], 416.59194032216; atol=TOL)
        @test isapprox(testqry[5,:val], 396.754228877952; atol=TOL)
        @test isapprox(testqry[6,:val], 377.861170359965; atol=TOL)
        @test isapprox(testqry[7,:val], 359.867781295205; atol=TOL)
        @test isapprox(testqry[8,:val], 342.731220281158; atol=TOL)
        @test isapprox(testqry[9,:val], 326.410685982046; atol=TOL)
        @test isapprox(testqry[10,:val], 310.8673199829; atol=TOL)

        SQLite.DBInterface.execute(db, "update STORAGE set netzeroyear = 0")

        # Test with calcyears
        NemoMod.calculatescenario(dbfile; jumpmodel = Model(Mosek.Optimizer), restrictvars=true,
            numprocs=1, calcyears=[2020,2029])
        testqry = SQLite.DBInterface.execute(db, "select * from vtotaldiscountedcost") |> DataFrame

        @test testqry[1,:y] == "2020"
        @test testqry[2,:y] == "2029"

        @test isapprox(testqry[1,:val], 3840.93779774689; atol=TOL)
        @test isapprox(testqry[2,:val], 3427.79849864389; atol=TOL)

        # Test MinimumUtilization
        SQLite.DBInterface.execute(db, "insert into MinimumUtilization select ROWID, '1', 'gas', val, 2025, 0.5 from TIMESLICE")
        NemoMod.calculatescenario(dbfile; jumpmodel = Model(Mosek.Optimizer), numprocs=1, varstosave="vproductionbytechnologyannual")
        testqry = SQLite.DBInterface.execute(db, "select * from vproductionbytechnologyannual where t = 'gas' and y = 2025") |> DataFrame

        @test isapprox(testqry[1,:val], 15.768; atol=TOL)

        SQLite.DBInterface.execute(db, "delete from MinimumUtilization")

        # Delete test results and re-compact test database
        NemoMod.dropresulttables(db)
        testqry = SQLite.DBInterface.execute(db, "VACUUM")
    end  # "Solving storage_test with Mosek"

    @testset "Solving storage_transmission_test with Mosek" begin
        dbfile = joinpath(@__DIR__, "storage_transmission_test.sqlite")
        #dbfile = "c:/temp/storage_transmission_test.sqlite"
        chmod(dbfile, 0o777)  # Make dbfile read-write. Necessary because after Julia 1.0, Pkg.add makes all package files read-only

        NemoMod.calculatescenario(dbfile; jumpmodel = JuMP.Model(Mosek.Optimizer),
            varstosave =
                "vdemandnn, vnewcapacity, vtotalcapacityannual, vproductionbytechnologyannual, vproductionnn, vusebytechnologyannual, vusenn, vtotaldiscountedcost, "
                * "vtransmissionbuilt, vtransmissionexists, vtransmissionbyline, vtransmissionannual",
            numprocs=1, restrictvars=false, quiet = false)

        db = SQLite.DB(dbfile)
        testqry = SQLite.DBInterface.execute(db, "select * from vtotaldiscountedcost") |> DataFrame

        @test testqry[1,:y] == "2020"
        @test testqry[2,:y] == "2021"
        @test testqry[3,:y] == "2022"
        @test testqry[4,:y] == "2023"
        @test testqry[5,:y] == "2024"
        @test testqry[6,:y] == "2025"
        @test testqry[7,:y] == "2026"
        @test testqry[8,:y] == "2027"
        @test testqry[9,:y] == "2028"
        @test testqry[10,:y] == "2029"

        @test isapprox(testqry[1,:val], 9786.56609769189; atol=TOL)
        @test isapprox(testqry[2,:val], 239.495174532817; atol=TOL)
        @test isapprox(testqry[3,:val], 228.090642412206; atol=TOL)
        @test isapprox(testqry[4,:val], 217.22918324972; atol=TOL)
        @test isapprox(testqry[5,:val], 206.88494534077; atol=TOL)
        @test isapprox(testqry[6,:val], 197.033272788862; atol=TOL)
        @test isapprox(testqry[7,:val], 187.650797384118; atol=TOL)
        @test isapprox(testqry[8,:val], 178.714986656564; atol=TOL)
        @test isapprox(testqry[9,:val], 170.204749196728; atol=TOL)
        @test isapprox(testqry[10,:val], 162.099761139741; atol=TOL)

        # Test MinimumUtilization
        SQLite.DBInterface.execute(db, "insert into MinimumUtilization select ROWID, '1', 'gas', val, 2025, 0.2 from TIMESLICE")
        NemoMod.calculatescenario(dbfile; jumpmodel = Model(Mosek.Optimizer), numprocs=1, varstosave="vproductionbytechnologyannual", calcyears=[2020,2025,2029])
        testqry = SQLite.DBInterface.execute(db, "select * from vproductionbytechnologyannual where t = 'gas' and y = 2025") |> DataFrame

        @test isapprox(testqry[1,:val], 16.3149963697108; atol=TOL)

        SQLite.DBInterface.execute(db, "delete from MinimumUtilization")

        # Test interest rates
        SQLite.DBInterface.execute(db, "insert into InterestRateStorage select rowid, 1, 'storage1', y.val, 0.05 from year y")
        SQLite.DBInterface.execute(db, "insert into InterestRateTechnology select rowid, 1, 'solar', y.val, 0.05 from year y")
        SQLite.DBInterface.execute(db, "update TransmissionLine set interestrate = 0.05 where id = 2")
        NemoMod.calculatescenario(dbfile; jumpmodel = Model(Mosek.Optimizer), numprocs=1, varstosave="vtotaldiscountedcost", calcyears=[2020,2025,2029])
        testqry = SQLite.DBInterface.execute(db, "select * from vtotaldiscountedcost") |> DataFrame

        @test testqry[1,:y] == "2020"
        @test testqry[2,:y] == "2025"
        @test testqry[3,:y] == "2029"

        @test isapprox(testqry[1,:val], 12672.8114730171; atol=TOL)
        @test isapprox(testqry[2,:val], 2510.44589185654; atol=TOL)
        @test isapprox(testqry[3,:val], 1611.02255483292; atol=TOL)

        SQLite.DBInterface.execute(db, "delete from InterestRateStorage")
        SQLite.DBInterface.execute(db, "delete from InterestRateTechnology")
        SQLite.DBInterface.execute(db, "update TransmissionLine set interestrate = null where id = 2")

        # Delete test results and re-compact test database
        NemoMod.dropresulttables(db)
        testqry = SQLite.DBInterface.execute(db, "VACUUM")
    end  # "Solving storage_transmission_test with Mosek"

    @testset "Solving ramp_test with Mosek" begin
        dbfile = joinpath(@__DIR__, "ramp_test.sqlite")
        #dbfile = "c:/temp/ramp_test.sqlite"
        chmod(dbfile, 0o777)  # Make dbfile read-write. Necessary because after Julia 1.0, Pkg.add makes all package files read-only

        NemoMod.calculatescenario(dbfile; jumpmodel = JuMP.Model(Mosek.Optimizer),
            numprocs=1, restrictvars=false, quiet = false)

        db = SQLite.DB(dbfile)
        testqry = SQLite.DBInterface.execute(db, "select * from vtotaldiscountedcost") |> DataFrame

        @test testqry[1,:y] == "2020"
        @test testqry[2,:y] == "2021"
        @test testqry[3,:y] == "2022"
        @test testqry[4,:y] == "2023"
        @test testqry[5,:y] == "2024"
        @test testqry[6,:y] == "2025"
        @test testqry[7,:y] == "2026"
        @test testqry[8,:y] == "2027"
        @test testqry[9,:y] == "2028"
        @test testqry[10,:y] == "2029"

        @test isapprox(testqry[1,:val], 4163.77771390785; atol=TOL)
        @test isapprox(testqry[2,:val], 1606.0180127748; atol=TOL)
        @test isapprox(testqry[3,:val], 1529.54096454743; atol=TOL)
        @test isapprox(testqry[4,:val], 1456.70568052136; atol=TOL)
        @test isapprox(testqry[5,:val], 1387.33874335368; atol=TOL)
        @test isapprox(testqry[6,:val], 1321.28388353478; atol=TOL)
        @test isapprox(testqry[7,:val], 1258.35996962756; atol=TOL)
        @test isapprox(testqry[8,:val], 1198.44067149609; atol=TOL)
        @test isapprox(testqry[9,:val], 1141.373808485; atol=TOL)
        @test isapprox(testqry[10,:val], 1084.83956451882; atol=TOL)

        # Delete test results and re-compact test database
        NemoMod.dropresulttables(db)
        testqry = SQLite.DBInterface.execute(db, "VACUUM")
    end  # "Solving ramp_test with Mosek"
end  # @isdefined MosekTools


text:
@testset "Ntfv.jl" begin
    filename = "data\\test.ntfv"
    dc = DDR2import.Ntfv.read(filename)
    @test dc["ABSON"].name == ""
    @test dc["ALASO"].category == "_"
    @test dc["AERODS1"].reflocname == "EDGGADS"
    @test dc["AERODS2"].refloctype == "AS"
    @test dc["AEROMN"].reflocrole == "G"
    #LAST
    @test dc["ALGERIA"].airblocks[1].name == "ELI>GC"
    @test dc["ALGERIA"].airblocks[2].upperFL == "IN"
    @test dc["ALGERIA"].airblocks[2].name == "GC>ELI"
end


text:
module stest
using LinearAlgebra
using BenchmarkTools
using ..sparse

function test_create_sparse_vec()
    data = [10, 20, 1e-11, -30, 1e-5, 50, 1e-50]
    index = [1, 10, 15, 20]
    # vec = sparse.SparseVector(data, index)

    vec2 = sparse.create_sparse(data)

end

function test_create_sparse_mat()

    n_row = 5
    n_col = 3

    data        = [10,  20, 30, 40, 50, 60, 70, 80, 90]
    row_index   = [1,   2,  3,  4,  5,  1,  1,  4,  5, n_row]
    col_index   = [1,   1,  1,  2,  2,  2,  3,  3,  3, n_col]

    spmat = sparse.SparseMatrixCOO(data, row_index, col_index)

end

function test_add()

    data_x  = [3.3, 4, 5, 6, 10, 4,  5,  2,  3]
    index_x = [3,   4, 6, 8, 12, 15, 16, 18, 20, 25]

    sp_x = sparse.SparseVector(data_x, index_x)

    data_y  = [5.0, 2, 3, 5, 2, 3,  5,  2,  3]
    index_y = [4,   5, 6, 7, 8, 14, 15, 16, 18, 25]

    sp_y = sparse.SparseVector(data_y, index_y)

    sp_z = sparse.add(sp_x, sp_y)

    show(sp_z)

end

function test_matrix_vector_mult()

    n_row = 5
    n_col = 3

    data        = [10,  20, 30, 40, 50, 60, 70, 80, 90]
    row_index   = [1,   2,  3,  4,  5,  1,  1,  4,  5, n_row]
    col_index   = [1,   1,  1,  2,  2,  2,  3,  3,  3, n_col]

    spmat = sparse.SparseMatrixCOO(data, row_index, col_index)

    fvec = [10.0, 5, 6]

    res = spmat * fvec # sparse.matrix_vector_prod(spmat, fvec)
    show(res)

    fmat = sparse.create_full(spmat)
    res_theo = fmat * fvec
    show(res_theo)

    show(res_theo - res)

end

function create_matrix(n_dim = 8)

    A = zeros(Float64, n_dim, n_dim)

    main_diag = [1; 4 * ones(n_dim-2); 4]
    side_diag = -ones(n_dim-1)
    sideside_diag = -ones(n_dim-3)

    main_diag_idx = diagind(A, 0)
    side_m1_diag_idx = diagind(A, 1)
    side_p1_diag_idx = diagind(A, -1)

    side_m3_diag_idx = diagind(A, 3)
    side_p3_diag_idx = diagind(A, -3)

    A[ main_diag_idx ] = main_diag
    A[ side_m1_diag_idx ] = side_diag
    A[ side_p1_diag_idx ] = side_diag
    A[ side_m3_diag_idx ] = sideside_diag
    A[ side_p3_diag_idx ] = sideside_diag

    return A
end

function test_csc()
    n_dim = 8
    A = create_matrix(n_dim)
    A[:,1] .= 0.0
    A[:,4] .= 0.0

    spmat = sparse.create_sparse_csc(A)

    Af = sparse.create_full(spmat)

    vec = rand(n_dim) # ones(n_dim)

    res = sparse.matrix_vector_prod(spmat, vec)

    res_theo = A * vec

    println( "$(norm(res - res_theo))" )

end

function test_matmult(spmat, vec)
    res = sparse.matrix_vector_prod(spmat, vec)
end

function benchmark_matmult()

    n_dim = 25000 # 25_000
    A = create_matrix(n_dim)
    spmat = sparse.create_sparse_csc(A)
    vec = rand(n_dim)

    b = @benchmark test_matmult($spmat, $vec)
    println("=== sparse matrix vec mult === ")
    show(
        stdout,
        MIME("text/plain"),
        b)
    println("")
    println("")

    b = @benchmark $A * $vec
    println("=== default full matrix vec mult === ")
    show(
        stdout,
        MIME("text/plain"),
        b)

    # res_sp = test_matmult(spmat, vec)
    # res_fu = A * vec;

    # println("$(norm(res_sp - res_fu))")

end

n_dim = 8
A = create_matrix(n_dim)
A[7, 6] = 3.14
A[8,7] = 2.6
spmat = sparse.create_sparse_csc(A)

(L, d, U) = sparse.matrix_split(spmat)

Lf = sparse.create_full(L)
Uf = sparse.create_full(U)
Df = diagm(d)

println(" $(norm( Lf + Uf + Df - A )) ")

# stest.L[8,6]

end # stest

max_stars_repo_name: HomoModelicus/julia

max_issues_repo_name: HomoModelicus/julia

max_forks_repo_name: HomoModelicus/julia


text:
@lexer TOMLLexer Dict(
    :name => "TOML",
    :description => "A lexer for TOML, a simple language for config files.",
    :aliases => ["TOML", "toml"],
    :filenames => ["*.toml"],
    :tokens => Dict(
        :root => [
            (r"\s+"sm, TEXT),
            (r"#.*?$"m, COMMENT),

            (r"\"\"\"", STRING, :triple_strings),
            (r"\"", STRING, :strings),
            (r"'''", STRING, :triple_literal_strings),
            (r"'", STRING, :literal_strings),

            (r"(true|false)$"m, KEYWORD_CONSTANT),
            (r"[a-zA-Z_][a-zA-Z0-9_\-]*", NAME),

            (r"([0-9]{4}-[0-9]{2}-[0-9]{2}T[0-9]{2}:[0-9]{2}:[0-9]{2}(?:\.[0-9]+)?(?:Z|[-+][0-9]{2}:[0-9]{2}|))", LITERAL_DATE),
            (r"([0-9]{4}-[0-9]{2}-[0-9]{2})", LITERAL_DATE),

            (r"([-+]?(?:[0-9_]+(?:\.[0-9_]+)?[eE][-+]?[0-9_]+|[0-9_]*\.[0-9_]+))", NUMBER_FLOAT),

            (r"([-+]?[0-9_]+)", NUMBER_INTEGER),

            (r"[\[\];.,:(){}]", PUNCTUATION),
            (r"\.", PUNCTUATION),

            (r"=", OPERATOR),
        ],
        :strings => [
            (r"\"", STRING, :__pop__),
            (r"\\([\\\"\'\$nrbtfav]|(x|u|U)[a-fA-F0-9]+|\d+)", STRING_ESCAPE),
            (r"."sm, STRING),
        ],
        :triple_strings => [
            (r"\"\"\"", STRING, :__pop__),
            (r"\\([\\\"\'\$nrbtfav]|(x|u|U)[a-fA-F0-9]+|\d+)", STRING_ESCAPE),
            (r"."sm, STRING),
        ],
        :literal_strings => [
            (r"'", STRING, :__pop__),
            (r"."sm, STRING),
        ],
        :triple_literal_strings => [
            (r"'''", STRING, :__pop__),
            (r"."sm, STRING),
        ],
    )
)


text:
# Copyright 2018 Erik Edin
#
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
#
#     http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.

module Experimental

import Base: (|), show

using Behavior.Gherkin: Given, When, Then, Scenario, ScenarioStep, AbstractScenario
using Behavior.Gherkin: Feature, FeatureHeader, Background, DataTableRow, DataTable
using Behavior.Gherkin: ScenarioOutline

struct GherkinSource
    lines::Vector{String}

    GherkinSource(source::String) = new(split(source, "\n"))
end

notempty(s) = !isempty(strip(s))
notcomment(s) = !startswith(strip(s), "#")
defaultlinecondition(s) = notempty(s) && notcomment(s)
anyline(s) = true

"""
    ParserInput

ParserInput encapsulates
- the Gherkin source
- the parser state (current line)
- the parser line operation (not implemented yet)
"""
struct ParserInput
    source::GherkinSource
    index::Int
    condition::Function

    ParserInput(source::String) = new(GherkinSource(source), 1, defaultlinecondition)
    ParserInput(input::ParserInput, index::Int) = new(input.source, index, defaultlinecondition)
    ParserInput(input::ParserInput, index::Int, condition::Function) = new(input.source, index, condition)
end

consume(input::ParserInput, n::Int) :: ParserInput = ParserInput(input, input.index + n)

function line(input::ParserInput) :: Tuple{Union{Nothing, String}, ParserInput}
    nextline = findfirst(input.condition, input.source.lines[input.index:end])
    if nextline === nothing
        return nothing, input
    end
    strip(input.source.lines[input.index + nextline - 1]), consume(input, nextline)
end

"""
    Parser{T}

A Parser{T} is a callable that takes a ParserInput and produces
a ParseResult{T}. The parameter T is the type of what it parses.
"""
abstract type Parser{T} end

"""
    ParseResult{T}

A ParseResult{T} is the either successful or failed result of trying to parse
a value of T from the input.
"""
abstract type ParseResult{T} end

"""
    OKParseResult{T}(value::T)

OKParseResult{T} is a successful parse result.
"""
struct OKParseResult{T} <: ParseResult{T}
    value::T
    newinput::ParserInput
end

"""
    BadParseResult{T}

BadParseResult{T} is an abstract supertype for all failed parse results.
"""
abstract type BadParseResult{T} <: ParseResult{T} end

"""
    isparseok(::ParseResult{T})

isparseok returns true for successful parse results, and false for unsuccessful ones.
"""

isparseok(::OKParseResult{T}) where {T} = true
isparseok(::BadParseResult{T}) where {T} = false

struct BadExpectationParseResult{T} <: BadParseResult{T}
    expected::String
    actual::String
    newinput::ParserInput
end

struct BadUnexpectedParseResult{T} <: BadParseResult{T}
    unexpected::String
    newinput::ParserInput
end

struct BadInnerParseResult{S, T} <: BadParseResult{T}
    inner::ParseResult{<:S}
    newinput::ParserInput
end

function Base.show(io::IO, result::BadInnerParseResult{S, T}) where {S, T}
    show(io, result.inner)
end
function Base.show(io::IO, mime::MIME"text/plain", result::BadInnerParseResult{S, T}) where {S, T}
    show(io, mime, result.inner)
end

struct BadCardinalityParseResult{S, T} <: BadParseResult{T}
    inner::BadParseResult{S}
    atleast::Int
    actual::Int
    newinput::ParserInput
end

struct BadUnexpectedEOFParseResult{T} <: BadParseResult{T}
    newinput::ParserInput
end

struct BadExpectedEOFParseResult{T} <: BadParseResult{T}
    newinput::ParserInput
end

function Base.show(io::IO, result::BadExpectedEOFParseResult{T}) where {T}
    s, _newinput = line(result.newinput)
    println(io, "Expected EOF but found at line $(result.newinput.index):")
    println(io)
    println(io, "  $s")
end

struct BadExceptionParseResult{T} <: BadParseResult{T}
    ex
end

function Base.show(io::IO, result::BadExceptionParseResult{T}) where {T}
    show(io, "Exception: $(result.ex)")
end

"""
    Line(expected::String)

Line is a parser that recognizes when a line is exactly an expected value.
It returns the expected value if it is found, or a bad parse result if not.
"""
struct Line <: Parser{String}
    expected::String
end

function (parser::Line)(input::ParserInput) :: ParseResult{String}
    s, newinput = line(input)

    if s === nothing
        BadUnexpectedEOFParseResult{String}(input)
    elseif s == parser.expected
        OKParseResult{String}(parser.expected, newinput)
    else
        BadExpectationParseResult{String}(parser.expected, s, input)
    end
end

"""
    Optionally{T}

Optionally{T} parses a value of type T if the inner parser succeeds, or nothing
if the inner parser fails. It always succeeds.
"""
struct Optionally{T} <: Parser{Union{Nothing, T}}
    inner::Parser{T}
end

function (parser::Optionally{T})(input::ParserInput) :: ParseResult{Union{Nothing, T}} where {T}
    result = parser.inner(input)
    if isparseok(result)
        OKParseResult{Union{Nothing, T}}(result.value, result.newinput)
    else
        OKParseResult{Union{Nothing, T}}(nothing, input)
    end
end

"""
    optionalordefault(value, default)

Return `value` if it is not nothing, default otherwise.
"""
function optionalordefault(value, default)
    if value === nothing
        default
    else
        value
    end
end

"""
    Or{T}(::Parser{T}, ::Parser{T})

Or matches either of the provided parsers. It short-circuits.
"""
struct Or{T} <: Parser{T}
    a::Parser{<:T}
    b::Parser{<:T}
end

function (parser::Or{T})(input::ParserInput) :: ParseResult{<:T} where {T}
    result = parser.a(input)
    if isparseok(result)
        result
    else
        parser.b(input)
    end
end

(|)(a::Parser{T}, b::Parser{T}) where {T} = Or{T}(a, b)

"""
    Transformer{S, T}

Transforms a parser of type S to a parser of type T, using a transform function.
"""
struct Transformer{S, T} <: Parser{T}
    inner::Parser{S}
    transform::Function
end

function (parser::Transformer{S, T})(input::ParserInput) where {S, T}
    result = parser.inner(input)
    if isparseok(result)
        newvalue = parser.transform(result.value)
        OKParseResult{T}(newvalue, result.newinput)
    else
        BadInnerParseResult{S, T}(result, input)
    end
end

"""
    Sequence{T}

Combine parsers into a sequence, that matches all of them in order.
"""
struct Sequence{T} <: Parser{Vector{T}}
    inner::Vector{Parser{<:T}}

    Sequence{T}(parsers...) where {T} = new(collect(parsers))
end

function (parser::Sequence{T})(input::ParserInput) :: ParseResult{Vector{T}} where {T}
    values = Vector{T}()
    currentinput = input

    for p in parser.inner
        result = p(currentinput)
        if isparseok(result)
            push!(values, result.value)
            currentinput = result.newinput
        else
            return BadInnerParseResult{T, Vector{T}}(result, input)
        end
    end

    OKParseResult{Vector{T}}(values, currentinput)
end

"""
    Joined

Joins a sequence of strings together into one string.
"""
Joined(inner::Parser{Vector{String}}) = Transformer{Vector{String}, String}(inner, x -> join(x, "\n"))

"""
    Repeating{T}

Repeats the provided parser until it no longer recognizes the input.
"""
struct Repeating{T} <: Parser{Vector{T}}
    inner::Parser{T}
    atleast::Int

    Repeating{T}(inner::Parser{T}; atleast::Int = 0) where {T} = new(inner, atleast)
end

function (parser::Repeating{T})(input::ParserInput) :: ParseResult{Vector{T}} where {T}
    values = Vector{T}()
    currentinput = input

    while true
        result = parser.inner(currentinput)
        if !isparseok(result)
            cardinality = length(values)
            if cardinality < parser.atleast
                return BadCardinalityParseResult{T, Vector{T}}(result, parser.atleast, cardinality, input)
            end
            break
        end
        push!(values, result.value)
        currentinput = result.newinput
    end

    OKParseResult{Vector{T}}(values, currentinput)
end

"""
    LineIfNot

Consumes a line if it does not match a given parser.
"""
struct LineIfNot <: Parser{String}
    inner::Parser{<:Any}
    linecondition::Function

    LineIfNot(inner::Parser{<:Any}) = new(inner, defaultlinecondition)
    LineIfNot(inner::Parser{<:Any}, condition::Function) = new(inner, condition)
end

function (parser::LineIfNot)(realinput::ParserInput) :: ParseResult{String}
    # This parser has support for providing the inner parser with a
    # non-default line condition.
    input = ParserInput(realinput, realinput.index, parser.linecondition)
    result = parser.inner(input)
    if isparseok(result)
        badline, _badinput = line(input)
        BadUnexpectedParseResult{String}(badline, input)
    else
        s, newinput = line(input)
        if s === nothing
            BadUnexpectedEOFParseResult{String}(input)
        else
            OKParseResult{String}(s, newinput)
        end
    end
end

"""
    StartsWith

Consumes a line if it starts with a given string.
"""
struct StartsWith <: Parser{String}
    prefix::String
end

function (parser::StartsWith)(input::ParserInput) :: ParseResult{String}
    s, newinput = line(input)
    if s === nothing
        BadUnexpectedEOFParseResult{String}(input)
    elseif startswith(s, parser.prefix)
        OKParseResult{String}(s, newinput)
    else
        BadExpectationParseResult{String}(parser.prefix, s, input)
    end
end

struct EOFParser <: Parser{Nothing} end

function (parser::EOFParser)(input::ParserInput) :: ParseResult{Nothing}
    s, newinput = line(input)
    if s === nothing
        OKParseResult{Nothing}(nothing, newinput)
    else
        BadExpectedEOFParseResult{Nothing}(input)
    end
end

##
## Gherkin-specific parser
##

takeelement(i::Int) = xs -> xs[i]

"""
    BlockText

Parses a Gherkin block text.
"""
BlockText() = Transformer{Vector{String}, String}(
    Sequence{String}(
        Line("\"\"\""),
        Joined(Repeating{String}(LineIfNot(Line("\"\"\""), anyline))),
        Line("\"\"\""),
    ),
    takeelement(2)
)

struct Keyword
    keyword::String
    rest::String
end

"""
    DataTableParser()

Consumes a data table in a scenario step.

## Example

    | Header 1 | Header 2 |
    | Foo      | Bar      |
    | Baz      | Quux     |
"""
struct DataTableRowParser <: Parser{DataTableRow} end

function (parser::DataTableRowParser)(input::ParserInput) :: ParseResult{DataTableRow}
    s, newinput = line(input)
    if s === nothing
        return BadUnexpectedEOFParseResult{DataTableRow}(input)
    end
    parts = split(s, "|")

    if length(parts) >= 3
        columns = collect([strip(p) for p in parts if strip(p) != ""])
        OKParseResult{DataTableRow}(columns, newinput)
    else
        BadExpectationParseResult{DataTableRow}("| column 1 | column 2 | ... | column n |", s, input)
    end
end

DataTableParser() = Transformer{Vector{DataTableRow}, DataTable}(
    Repeating{DataTableRow}(DataTableRowParser(), atleast=1),
    rows -> rows
)

struct AnyLine <: Parser{String} end

function (parser::AnyLine)(input::ParserInput) :: ParseResult{String}
    s, newinput = line(input)
    if s === nothing
        BadUnexpectedEOFParseResult{String}(input)
    else
        OKParseResult{String}(s, newinput)
    end
end

Splitter(inner::Parser{String}, delimiter) :: Parser{Vector{String}} = Transformer{String, Vector{String}}(
    inner,
    s -> split(s, delimiter, keepempty=false)
)

struct Validator{T} <: Parser{Vector{T}}
    inner::Parser{Vector{T}}
    f::Function
end

function (parser::Validator{T})(input::ParserInput) :: ParseResult{Vector{T}} where {T}
    result = parser.inner(input)
    if isparseok(result) && all(parser.f, result.value)
        OKParseResult{Vector{T}}(result.value, result.newinput)
    else
        BadInnerParseResult{Vector{T}, Vector{T}}(result, input)
    end
end

const TagList = Vector{String}

"""
    TagParser()

Consumes a single tag in the form `@tagname`.
"""
TagParser() :: Parser{TagList} = Validator{String}(Splitter(AnyLine(), isspace), x -> startswith(x, "@"))

"""
    TagLinesParser()

Read tags until they stop.
"""
TagLinesParser() :: Parser{TagList} = Transformer{Vector{TagList}, TagList}(
    Repeating{TagList}(TagParser()),
    taglines -> vcat(taglines...)
)

"""
    KeywordParser

Recognizes a keyword, and any following text on the same line.
"""
KeywordParser(word::String) = Transformer{String, Keyword}(
    StartsWith(word),
    s -> begin
        Keyword(word, strip(replace(s, word => "", count=1)))
    end
)

const TableOrBlockTextTypes = Union{String, DataTable}
const DataTableOrBlockText = (
    Sequence{TableOrBlockTextTypes}(DataTableParser(), BlockText()) |
    Sequence{TableOrBlockTextTypes}(BlockText(), DataTableParser()) |
    Sequence{TableOrBlockTextTypes}(BlockText()) |
    Sequence{TableOrBlockTextTypes}(DataTableParser())
)
const StepPieces = Union{Keyword, Union{Nothing, Vector{TableOrBlockTextTypes}}}

mutable struct StepBuilder{T}
    steptype::Type{T}
    keyword::Keyword
    blocktext::String
    datatable::DataTable

    StepBuilder{T}(steptype::Type{T}, keyword::Keyword) where {T} = new(steptype, keyword, "", DataTable())
end

accumulate!(sb::StepBuilder{T}, table::DataTable) where {T} = sb.datatable = table
accumulate!(sb::StepBuilder{T}, blocktext::String) where {T} = sb.blocktext = blocktext
accumulate!(sb::StepBuilder{T}, vs::Vector{TableOrBlockTextTypes}) where {T} = foreach(v -> accumulate!(sb, v), vs)
accumulate!(::StepBuilder{T}, ::Nothing) where {T} = nothing

buildstep(sb::StepBuilder{T}) where {T} = sb.steptype(sb.keyword.rest, block_text=sb.blocktext, datatable=sb.datatable)

function StepParser(steptype::Type{T}, keyword::String) :: Parser{T} where {T}
    Transformer{Vector{StepPieces}, T}(
        Sequence{StepPieces}(KeywordParser(keyword), Optionally(DataTableOrBlockText)),
        sequence -> begin
            keyword = sequence[1]
            stepbuilder = StepBuilder{T}(steptype, keyword)

            accumulate!(stepbuilder, sequence[2])

            buildstep(stepbuilder)
        end
    )
end

"""
    And

An `And` type scenario step
"""
struct And <: ScenarioStep
    text::String
    block_text::String
    datatable::DataTable

    And(text::AbstractString; block_text = "", datatable=DataTable()) = new(text, block_text, datatable)
end

"""
    GivenParser

Consumes a Given step.
"""
GivenParser() = StepParser(Given, "Given ")
WhenParser() = StepParser(When, "When ")
ThenParser() = StepParser(Then, "Then ")
AndParser() = StepParser(And, "And ")

# TODO Find a way to express this as
#      GivenParser() | WhenParser() | ThenParser()
const AnyStepParser = Or{ScenarioStep}(
    Or{ScenarioStep}(GivenParser(), WhenParser()),
    Or{ScenarioStep}(ThenParser(), AndParser())
)
"""
    StepsParser

Parses zero or more steps.
"""
StepsParser() = Repeating{ScenarioStep}(AnyStepParser)

const AnyKeyword = (
    KeywordParser("Given ") |
    KeywordParser("When ") |
    KeywordParser("Then ") |
    KeywordParser("Feature:") |
    KeywordParser("Scenario:") |
    KeywordParser("Scenario Outline:") |
    KeywordParser("Background:") |
    KeywordParser("Rule:")
)
const MaybeTags = Union{Nothing, Vector{String}}
const ScenarioBits = Union{Keyword, String, Vector{ScenarioStep}, MaybeTags}

const LongDescription = Joined(Repeating{String}(LineIfNot(AnyKeyword)))
"""
    ScenarioParser()

Consumes a Scenario.
"""
ScenarioParser() = Transformer{Vector{ScenarioBits}, Scenario}(
    Sequence{ScenarioBits}(
        Optionally(TagLinesParser()),
        KeywordParser("Scenario:"),
        Optionally(LongDescription),
        StepsParser()),
    sequence -> begin
        tags = optionalordefault(sequence[1], [])
        keyword = sequence[2]
        longdescription = optionalordefault(sequence[3], "")
        Scenario(keyword.rest, tags, Vector{ScenarioStep}(sequence[4]), long_description=longdescription)
    end
)

const ScenarioOutlineBits = Union{Keyword, String, Vector{ScenarioStep}, MaybeTags, DataTable}
"""
    ScenarioOutlineParser()

Consumes a Scenario Outline.
"""
ScenarioOutlineParser() = Transformer{Vector{ScenarioOutlineBits}, ScenarioOutline}(
    Sequence{ScenarioOutlineBits}(
        Optionally(TagLinesParser()),
        KeywordParser("Scenario Outline:"),
        Optionally(LongDescription),
        StepsParser(),
        Line("Examples:"),
        DataTableParser()
    ),
    sequence -> begin
        tags = optionalordefault(sequence[1], [])
        keyword = sequence[2]
        longdescription = optionalordefault(sequence[3], "")
        steps = sequence[4]
        examples = sequence[6]
        # The DataTableParser guarantees at least 1 row.
        placeholders = examples[1]
        ScenarioOutline(
            keyword.rest,
            tags,
            steps,
            placeholders,
            examples[2:end],
            long_description=longdescription
        )
    end
)

const BackgroundBits = ScenarioBits
"""
    BackgroundParser()

Consume a Background.
"""
BackgroundParser() = Transformer{Vector{BackgroundBits}, Background}(
    Sequence{BackgroundBits}(
        KeywordParser("Background:"),
        Optionally(LongDescription),
        StepsParser()
    ),
    sequence -> begin
        keyword = sequence[1]
        longdescription = optionalordefault(sequence[2], "")
        Background(keyword.rest, sequence[3], long_description=longdescription)
    end
)

struct Rule <: AbstractScenario
    description::String
    longdescription::String
    scenarios::Vector{AbstractScenario}
    tags::Vector{String}
end

const ScenarioList = Vector{Scenario}
const RuleBits = Union{Keyword, MaybeTags, Nothing, String, ScenarioList}

"""
    RuleParser

Consumes a Rule and its child scenarios.
"""
RuleParser() = Transformer{Vector{RuleBits}, Rule}(
    Sequence{RuleBits}(
        Optionally(TagLinesParser()),
        KeywordParser("Rule:"),
        Optionally(LongDescription),
        Repeating{Scenario}(ScenarioParser())),
    sequence -> begin
        tags = optionalordefault(sequence[1], String[])
        keyword = sequence[2]
        longdescription = optionalordefault(sequence[3], "")
        scenarios = sequence[4]
        Rule(keyword.rest, longdescription, scenarios, tags)
    end
)

const FeatureBits = Union{Keyword, String, Background, Nothing, Vector{AbstractScenario}, MaybeTags}
"""
    FeatureParser

Consumes a full feature file.
"""
const ScenarioOrRule = Or{AbstractScenario}(
    Or{AbstractScenario}(ScenarioParser(), ScenarioOutlineParser()),
    RuleParser()
)
FeatureParser() = Transformer{Vector{FeatureBits}, Feature}(
    Sequence{FeatureBits}(
        Optionally(TagLinesParser()),
        KeywordParser("Feature:"),
        Optionally(LongDescription),
        Optionally(BackgroundParser()),
        Repeating{AbstractScenario}(ScenarioOrRule)),
    sequence -> begin
        keyword = sequence[2]
        tags = optionalordefault(sequence[1], String[])
        longdescription = optionalordefault(sequence[3], "")
        background = optionalordefault(sequence[4], Background())
        Feature(FeatureHeader(keyword.rest, [longdescription], tags), background, sequence[5])
    end
)

const FeatureFileBits = Union{Feature, Nothing}
FeatureFileParser() = Transformer{Vector{FeatureFileBits}, Feature}(
    Sequence{FeatureFileBits}(FeatureParser(), EOFParser()),
    takeelement(1)
)

##
## Exports
##

export ParserInput, OKParseResult, BadParseResult, isparseok

# Basic combinators
export Line, Optionally, Or, Transformer, Sequence
export Joined, Repeating, LineIfNot, StartsWith, EOFParser

# Gherkin combinators
export BlockText, KeywordParser
export StepsParser, GivenParser, WhenParser, ThenParser
export ScenarioParser, RuleParser, FeatureParser, FeatureFileParser, BackgroundParser
export DataTableParser, TagParser, TagLinesParser, ScenarioOutlineParser

# Data carrier types
export Keyword, Rule

end


text:
using Images, TestImages

function do_edge_detection(sourceImage:: Array{GrayA{N0f8}, 2},
                           kernelSize::Int = 3) :: Array{GrayA{N0f8}, 2}
  return sourceImage
end

function build_kernels(kernelSize::Int) :: Array{UInt8, 2}
  if kernelSize & 1 == 0
    throw(ArgumentError("Kernel size must be odd"))
  end

  kernelBase = fill(-1, kernelSize, 1)
  kernelRows = [if (x == ceil(kernelSize/2)) 0 else -1 end for x in 1:kernelSize]

  kernelY = kernelBase * kernelRows

  (rotl90(kernelY), kernelY)
end

println(build_kernels(9))


text:
#######################################################################
# Macros for printing information depending on verbosity
#
# ASSUMPTION: the package defines variable `VERBOSE`
#######################################################################

const STATUS_BEGIN = "..."
const STATUS_END   = " completed"

# If verbose, prints information about the start and completion of
# of `action` named `name`
macro status(name :: String, action)
    quote
        @statusb $name
        $action
        @statuse $name
    end
end

# If verbose, prints information about the start of `name`
macro statusb(name :: String)
    :(if VERBOSE ; @info($name * $STATUS_BEGIN) end)
end

# If verbose, prints information about the completion of `name`
macro statuse(name :: String)
    :(if VERBOSE ; @info($name * $STATUS_END) end)
end

max_stars_repo_name: julbinb/julia-pkgsrc-downloader

max_issues_repo_name: julbinb/julia-pkgsrc-downloader

max_forks_repo_name: julbinb/julia-pkgsrc-downloader


text:
using BinaryBuilder

name = "Vulkan_Headers"
version = v"1.3.207"

source = "https://github.com/KhronosGroup/Vulkan-Headers.git"
commit = "2c45218b9008dd80f6004b419777758e6459ac16"

sources = [
    GitSource(source, commit)
]

script = raw"""
cd Vulkan-Headers
install_license LICENSE.txt

CMAKE_FLAGS=()

# Install things into $prefix
CMAKE_FLAGS+=(-DCMAKE_INSTALL_PREFIX=${prefix})

# Install things into $prefix
CMAKE_FLAGS+=(-DCMAKE_INSTALL_LIBDIR=${libdir})

cmake -B build -S . -GNinja ${CMAKE_FLAGS[@]}
ninja -C build -j ${nproc} install
"""

# platforms = supported_platforms()
platforms = [AnyPlatform()]

# The products that we will ensure are always built
products = Product[
    FileProduct("share/vulkan/registry/vk.xml", :vk_xml)
]

# Dependencies that must be installed before this package can be built
dependencies = Dependency[]

build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies)


text:
module HCIToolbox

export rotate,
       rotate!,
       derotate,
       derotate!,
       flatten,
       expand,
       combine,
       mask_circle!,
       mask_circle

# Utilities for dealing with cubes like derotating and median combining
include("morphology.jl")

# Utilities for mask design
include("masking.jl")

end


text:
export Vtree, vtree, variable, find_leaf,
    varsubset, varsubset_left, varsubset_right,
    depth,
    respects_vtree

import Base: parent, in # extend
import .Utils: find_leaf, find_inode, lca, depth # extend

#############
# Vtree
#############

"Root of the vtree node hiearchy"
abstract type Vtree <: Tree end

#############
# Methods
#############

"Get the variable in a vtree leaf"
function variable end

in(v::Var, vtree::Vtree) = v ∈ variables(vtree)

"Find the leaf in the vtree that represents the given variable"
function find_leaf(v::Var, n::Vtree)
    leaf = find_leaf(n, x -> v ∈ x)
    @assert variable(leaf) == v
    return leaf
end

"Compute the path length from vtree node `n` to leaf node for variable `var`"
depth(n::Vtree, v::Var)::Int = depth(n, x -> v ∈ x)

"""
Find an inner vtree node that has `left` in its left subtree and `right` in its right subtree.
Supports `nothing` as a catch-all for any node
"""
find_inode(left, right) = find_inode(left, right, varsubset)

# performance critical in SDD compilation:
"Are the variables in `n` contained in the variables in `m`?"
@inline varsubset(n::Vtree, m::Vtree) =
    (n===m || varsubset(n, m, NodeType(n), NodeType(m)))
@inline varsubset(n::Vtree, m::Vtree, ::Leaf, ::Leaf) =
    variable(n) == variable(m)
@inline varsubset(n::Vtree, m::Vtree, ::Inner, ::Leaf) =
    false
@inline varsubset(n::Vtree, m::Vtree, ::Leaf, ::Inner) =
    variable(n) ∈ m
@inline varsubset(n::Vtree, m::Vtree, ::Inner, ::Inner) =
    subseteq_fast(variables(n), variables(m)) # very slow if done with Base.BitSet

"Are the variables in `n` contained in the left branch of `m`?"
@inline varsubset_left(n::Vtree, m::Vtree)::Bool = varsubset(n, m.left)

"Are the variables in `n` contained in the right branch of `m`?"
@inline varsubset_right(n::Vtree, m::Vtree)::Bool = varsubset(n, m.right)

"""
Compute the lowest common ancestor of two vtree nodes
Warning: this method uses an incomplete `varsubset` check for `descends_from` and is only correct when `v` and `w` are part of the same larger vtree.
"""
lca(v::Union{Nothing,Vtree}, w::Union{Nothing,Vtree}) = lca(v, w, varsubset)

# "Compute the LCA of two vtree nodes, allowing either one to be `nothing`"
# lca_safe(x,y)::Union{Nothing,Vtree} =
#     (isnothing(x) ? y : (isnothing(y) ? x : lca(x,y)) )

"""
Does the circuit respect the given vtree?
This function allows for constants in conjunctions, but only when a vtree node can be found where the left and right conjunct can be assigned to the left and right vtree.
"""
function respects_vtree(circuit::LogicCircuit, vtree::Vtree)
    result::Bool = true
    f_con(n) = nothing
    f_lit(n) = if variable(n) ∉ vtree
        result = false
        return nothing
    else
        find_leaf(variable(n), vtree)
    end
    f_a(n, call) = if num_children(n) != 2
            result = false
            return nothing
        else
            vtree_left = call(children(n)[1])
            vtree_right = call(children(n)[2])
            isnothing(vtree_left) && isnothing(vtree_right) && return nothing
            vtree_and = find_inode(vtree_left, vtree_right)
            result &= issomething(vtree_and)
            vtree_and
        end
    f_o(n, call) = mapreduce(call, lca, children(n))
    foldup(circuit, f_con, f_lit, f_a, f_o, Union{Vtree,Nothing})
    return result
end

# print vtree nodes by showing their type and variable scope
Base.show(io::IO, c::Vtree) = print(io, "$(typeof(c))($(join(variables(c), ',')))")

#############
# Constructors
#############

# Syntactic sugar to compile circuits using a vtree
(vtree::Vtree)(arg) = compile(vtree, arg)

# construct vtrees from other vtrees
function (::Type{V})(vtree::Vtree)::V where V<:Vtree
    f_leaf(l) = V(variable(l))
    f_inner(i,call) = V(call(i.left), call(i.right))
    foldup(vtree,f_leaf, f_inner, V)
end

# construct a vtree for a given number of variables
function (::Type{V})(num_vars::Int, structure;
                    f = noop,
                    ordered_leafs = true)::V where V<:Vtree
    vars = Var.(ordered_leafs ? (1:num_vars) : randperm(num_vars))
    leaves = V.(vars)
    V(leaves, structure; f)
end

using Random: rand, randperm

# construct a vtree for a given set of leaves
function (::Type{V})(leafs::AbstractVector{<:V}, structure;
                    f = noop)::V where V<:Vtree
    length(leafs) == 1 && return leafs[1]
    if structure ∉ [:topdown, :bottomup]
        # check some predefined top-down structures
        if structure == :rightlinear
            g = x -> (x[1:1],x[2:end])
        elseif structure == :leftlinear
            g = x -> (x[1:length(x)-1],x[length(x):end])
        elseif structure == :balanced
            g = x -> (x[1:length(x)÷2],x[length(x)÷2+1:end])
        elseif structure == :random
            g = x -> begin
                split = rand(1:length(x)-1)
                (x[1:split],x[split+1:end])
            end
        else
            error("Vtree structure $(structure) not supported.")
        end
        return V(leafs, :topdown; f=g)
    end
    if structure == :topdown
        l, r = f(leafs)
        left = V(l, structure; f)
        right = V(r,structure; f)
        return V(left, right)
    else
        @assert structure == :bottomup
        pairs = f(leafs)
        leafs = map(x -> (x isa Tuple) ? V(x...) : x, pairs)
        return V(leafs, structure; f)
    end
end


text:
module test_julia_base

M = Base.require(Main, :Test)
using .M

@test true

@test basename(Base.active_project()) == "Project.toml"

end # module test_julia_base


text:
# Note that this test assumes TPF properly works
pm = PoolModel("ss1")
pm <= Setting(:data_vintage, "190822")
filepath = dirname(@__FILE__)
pm <= Setting(:dataroot, "$(filepath)/../reference/")
data = df_to_matrix(pm, load_data(pm))

# This commented code produces the saved output
tuning = Dict(:r_star => 2., :c_init => 0.3, :target_accept_rate => 0.4,
              :resampling_method => :systematic, :n_mh_steps => 1,
              :n_particles => 1000, :n_presample_periods => 0,
              :allout => true)
pm <= Setting(:tuning, tuning, "tuning parameters for TPF")
pm[:ρ].value = 0.5
pm[:μ].value = 0.
pm[:σ].value = 1.
Random.seed!(1793)
Φfilt, Ψfilt, F_ϵfilt, F_ufilt, F_λfilt = compute_system(pm)
s_init = reshape(rand(F_λfilt, tuning[:n_particles]), 1, 1000)
s_init = [s_init; 1 .- s_init]
tpf_out, ~, ~ = tempered_particle_filter(data, Φfilt, Ψfilt, F_ϵfilt, F_ufilt,
                                   s_init; tuning..., verbose = :none,
                                   fixed_sched = [1.], parallel = false, poolmodel = true)

Random.seed!(1793)
filt_tpf_out, ~, ~ = DSGE.filter(pm, data; tuning = get_setting(pm, :tuning))
Random.seed!(1793)
filt_lik_tpf_out = sum(DSGE.filter_likelihood(pm, data; tuning = get_setting(pm, :tuning)))

@testset "Check call to tempered particle filter without providing initial states" begin
    @test tpf_out == filt_tpf_out
    @test tpf_out == filt_lik_tpf_out
end

Random.seed!(1793)
s_init = reshape(rand(F_λfilt, tuning[:n_particles]), 1, 1000)
s_init = [s_init; 1 .- s_init] # this tpf output should be saved later
filt_tpf_out, ~, ~ = DSGE.filter(pm, data, s_init; tuning = get_setting(pm, :tuning))
Random.seed!(1793)
s_init = reshape(rand(F_λfilt, tuning[:n_particles]), 1, 1000)
s_init = [s_init; 1 .- s_init] # this tpf output should be saved later
filt_lik_tpf_out = sum(DSGE.filter_likelihood(pm, data, s_init; tuning = get_setting(pm, :tuning)))

@testset "Check call to tempered particle filter when providing initial states" begin
    @test tpf_out == filt_tpf_out
    @test tpf_out == filt_lik_tpf_out
end

pm = PoolModel("ss1", weight_type = :equal)
filepath = dirname(@__FILE__)
pm <= Setting(:dataroot, "$(filepath)/../reference/")
pm <= Setting(:data_vintage, "190822")
data = df_to_matrix(pm, load_data(pm))
filt_out, ~ = DSGE.filter(pm, data[:,1:2])
filtlik_out = sum(DSGE.filter_likelihood(pm, data[:,1:2]))
truefilt_out = sum(log.(sum(data[:,1:2] ./ 2, dims = 1)))
@testset "Check equal weights computes properly" begin
    @test filt_out == truefilt_out
    @test filtlik_out == truefilt_out
end

pm = PoolModel("ss1", weight_type = :static)
pm <= Setting(:data_vintage, "190822")
filepath = dirname(@__FILE__)
pm <= Setting(:dataroot, "$(filepath)/../reference/")
pm <= Setting(:data_vintage, "190822")
data = df_to_matrix(pm, load_data(pm))
s_init = [0.5, 0.5]
filt_out, ~  = DSGE.filter(pm, data[:,1:2], s_init)
filtlik_out = sum(DSGE.filter_likelihood(pm, data[:,1:2], s_init))
truefilt_out = sum(log.(sum(data[:,1:2] ./ 2, dims = 1)))
@testset "Check static weights computes properly" begin
    @test filt_out == truefilt_out
    @test filtlik_out == truefilt_out
end

@testset "Check warning and error conditions" begin
    pm = PoolModel("ss1")
    pm <= Setting(:data_vintage, "190822")
    filepath = dirname(@__FILE__)
    pm <= Setting(:dataroot, "$(filepath)/../reference/")
    tuning = get_setting(pm, :tuning)
    tuning[:n_particles] = 10
    pm <= Setting(:tuning, tuning)
    @test_throws ErrorException DSGE.filter(pm, data, vec([1. 1.; 1. 1.]))
    delete!(pm.settings, :tuning)
    @test_logs (:warn, "no tuning parameters provided; using default tempered particle filter values") DSGE.filter(pm, data)
    pm = PoolModel("ss1", weight_type = :bma)
    pm <= Setting(:data_vintage, "190822")
    filepath = dirname(@__FILE__)
    pm <= Setting(:dataroot, "$(filepath)/../reference/")
    @test_throws ErrorException DSGE.filter(pm, data)
end

nothing


text:
import ..UncertainValues.TheoreticalDistributionScalarValue
import ..UncertainValues.UncertainScalarBinomialDistributed
import ..UncertainValues.UncertainScalarBetaBinomialDistributed
import ..UncertainValues: AbstractUncertainValue

import Distributions.pdf

import ..SamplingConstraints:
    SamplingConstraint,
	constrain

function get_density(uv::AbstractUncertainValue)
	some_sample = resample(uv, 10000)
	xmin = minimum(some_sample) * 0.97
	xmax = maximum(some_sample) * 1.03
	step = (xmax-xmin)/300
	xvals = xmin:step:xmax+step
	density = pdf.(uv.distribution, xvals)

	xvals, density ./ sum(density)
end

function get_density(uv::UncertainScalarBinomialDistributed)
	some_sample = resample(uv, 10000)
	xmin = minimum(some_sample)
	xmax = maximum(some_sample)

	xvals = xmin:1:xmax
	density = pdf.(uv.distribution, xvals)

	xvals, density ./ sum(density)
end

function get_density(uv::UncertainScalarBetaBinomialDistributed)
	some_sample = resample(uv, 10000)
	xmin = minimum(some_sample)
	xmax = maximum(some_sample)

	xvals = xmin:1:xmax
	density = pdf.(uv.distribution, xvals)

	xvals, density ./ sum(density)
end

@recipe function plot_theoretical(uv::TheoreticalDistributionScalarValue,
		density = true, n_samples = 1000)

	if density
		@series begin
			get_density(uv)
		end
	else
		@series begin
			label --> ""
			resample(uv, n_samples)
		end
	end
end

@recipe function plot_theoretical(uv::TheoreticalDistributionScalarValue,
	constraint::SamplingConstraint, density = true, n_samples = 1000)
	cuv = constrain(uv, constraint)

	if density
		@series begin
			get_density(cuv)
		end
	else
		@series begin
			label --> ""
			resample(cuv, n_samples)
		end
	end
end


text:
@testset "CUDNN" begin

if !isdefined(CuArrays, :CUDNN)
@warn "Not testing CUDNN"
else
using CuArrays.CUDNN
@info "Testing CUDNN $(CUDNN.version())"

@testset "NNlib" begin
  using NNlib
  using NNlib: ∇conv_data, ∇conv_filter,
               maxpool, meanpool, ∇maxpool, ∇meanpool,
               softmax, ∇softmax, logsoftmax, ∇logsoftmax
  a, b, c = rand(Float64, 10, 10, 3, 1), rand(Float64, 2, 2, 3, 4), rand(Float64, 9, 9, 4, 1)
  da, db, dc = CuArray(a), CuArray(b), CuArray(c)
  cdims = DenseConvDims(a, b)
  @test NNlib.conv(a, b, cdims) ≈ collect(NNlib.conv(da, db, cdims))
  @test ∇conv_data(c, b, cdims) ≈ collect(∇conv_data(dc, db, cdims))
  @test ∇conv_filter(a, c, cdims) ≈ collect(∇conv_filter(da, dc, cdims))

  # Test for agreement between CPU NNlib and CuDNN versions, across a variety of kwargs
  for num_spatial_dims in (2, 3)
    # Initialize data we'll run our tests over
    C_in = 3
    C_out = 4
    batch_size = 1
    x = rand(Float64, repeat([8], num_spatial_dims)..., C_in, batch_size)
    w = rand(Float64, repeat([2], num_spatial_dims)..., C_in, C_out)
    b = rand(Float64, repeat([1], num_spatial_dims)..., C_in, C_out)
    options = (Dict(), Dict(:dilation => 2), Dict(:flipkernel => true), Dict(:stride => 2),)
    algos = (1, 0, 1, 1,)

    for (opts, algo) in zip(options, algos)
      cdims = DenseConvDims(x, w; opts...)
      y = NNlib.conv(x, w, cdims)

      # Test that basic convolution is equivalent across GPU/CPU
      @test testf((x, w) -> NNlib.conv(x, w, cdims), x, w)
      @test testf((y, w) -> ∇conv_data(y, w, cdims), y, w)
      @test testf((x, y) -> ∇conv_filter(x, y, cdims), x, y)
      # Test that we can use an alternative algorithm without dying
      @test_nowarn NNlib.conv!(cu(y), cu(x), cu(w), cdims; algo=algo)
      @test_nowarn NNlib.∇conv_data!(cu(x), cu(y), cu(w), cdims; algo=algo)
      @test_nowarn NNlib.∇conv_filter!(cu(w), cu(x), cu(y), cdims; algo=algo)
    end

    # Test that pooling is equivalent across GPU/CPU
    pdims = PoolDims(x, 2)
    y = maxpool(x, pdims)
    dy = ones(size(y))
    @test testf(x -> maxpool(x, pdims), x)
    @test testf((dy, y, x) -> ∇maxpool(dy, y, x, pdims), dy, y, x)
    @test testf(x -> maxpool(x, pdims), x)
    @test testf((dy, y, x) -> ∇maxpool(dy, y, x, pdims), dy, y, x)

    # CPU implementation of ∇conv_bias!
    db = zeros(Float64, 1, 1, 3, 1)
    function CuArrays.CUDNN.∇conv_bias!(db, y)
      db .= sum(y, dims=(1:(ndims(y)-2)))
      return db
    end
    #@test testf(CuArrays.CUDNN.∇conv_bias!, db, y)
  end

  for dims in [(5,5), (5,)]
    @test testf(softmax, rand(Float64, dims))
    @test testf(∇softmax, rand(Float64, dims), rand(Float64, dims))
    @test testf(logsoftmax, rand(Float64, dims))
    @test testf(∇logsoftmax, rand(Float64, dims), rand(Float64, dims))
  end
end

@testset "Activations and Other Ops" begin
  @test testf(CuArrays.CUDNN.cudnnAddTensor, cu(rand(Float64, 10, 10, 3, 1)), cu(rand(Float64, 10, 10, 3, 1)))
  @test testf(CuArrays.CUDNN.cudnnActivationForward, cu(rand(Float64, 10, 10, 3, 1)), cu(rand(Float64, 10, 10, 3, 1)))
  @test testf(CuArrays.CUDNN.cudnnActivationBackward, cu(rand(Float64, 10, 10, 3, 1)), cu(rand(Float64, 10, 10, 3, 1)), cu(rand(Float64, 10, 10, 3, 1)), cu(rand(Float64, 10, 10, 3, 1)))
end

end

end


text:
#-------------------------------------------------------------------
#* EMSO Model Library (EML) Copyright (C) 2004 - 2007 ALSOC.
#*
#* This LIBRARY is free software; you can distribute it and/or modify
#* it under the therms of the ALSOC FREE LICENSE as available at
#* http://www.enq.ufrgs.br/alsoc.
#*
#* EMSO Copyright (C) 2004 - 2007 ALSOC, original code
#* from http://www.rps.eng.br Copyright (C) 2002-2004.
#* All rights reserved.
#*
#* EMSO is distributed under the therms of the ALSOC LICENSE as
#* available at http://www.enq.ufrgs.br/alsoc.
#*
#*----------------------------------------------------------------------
#* File containg models of columns: distillation, stripping, absorbers
#* rectifier, ....
#*
#* The default nomenclature is:
#*		Type_Column_reboilertype_condensertyper
#*
#* where:
#*	Type = refluxed or reboiled or section
#*	Column = Stripping, Absorption, Rectifier, Distillation
#*	Reboiler type (if exists) = kettle or thermosyphon
#*	Condenser type (if exists) = with subccoling or without subcooling
#*
#*-----------------------------------------------------------------------
#* Author: Paula B. Staudt
#* $Id: column.mso 511 2008-05-12 17:25:33Z paula $
#*---------------------------------------------------------------------
type ControlSection
	ControlSection()=begin
		NComp=outers.NComp
		new(
			DanaInteger(Dict{Symbol,Any}(
				:Brief=>"Number of components"
			)),
			DanaInteger(Dict{Symbol,Any}(
				:Brief=>"Temperature Indicator in the tray"
			)),
			DanaInteger(Dict{Symbol,Any}(
				:Brief=>"Pressure Indicator in the tray"
			)),
			DanaInteger(Dict{Symbol,Any}(
				:Brief=>"Composition Indicator in the tray"
			)),
			DanaInteger(Dict{Symbol,Any}(
				:Brief=>"Component to indicate"
			)),
			[:NComp,:Tindicator_TrayNumber,:Pindicator_TrayNumber,:zindicator_TrayNumber,:zindicator_Component,],
		)
	end
	NComp::DanaInteger
	Tindicator_TrayNumber::DanaInteger
	Pindicator_TrayNumber::DanaInteger
	zindicator_TrayNumber::DanaInteger
	zindicator_Component::DanaInteger
	parameters::Array{Symbol,1}
	attributes::Dict{Symbol,Any}
end
export ControlSection
function atributes(in::ControlSection,_::Dict{Symbol,Any})
	fields::Dict{Symbol,Any}=Dict{Symbol,Any}()
	fields[:Pallete]=false
	fields[:Brief]="Control variables for column section."
	drive!(fields,_)
	return fields
end
ControlSection(_::Dict{Symbol,Any})=begin
	newModel=ControlSection()
	newModel.attributes=atributes(newModel,_)
	newModel
end


text:
# Use baremodule to shave off a few KB from the serialized `.ji` file
baremodule prrte_jll
using Base
using Base: UUID
import JLLWrappers

JLLWrappers.@generate_main_file_header("prrte")
JLLWrappers.@generate_main_file("prrte", UUID("eb928a42-fffd-568d-ab9c-3f5d54fc65b9"))
end  # module prrte_jll


text:
module MutualInformationImageRegistration

using FastHistograms
using OffsetArrays

include("mutual_information.jl")
include("register.jl")

export MutualInformationContainer
export mutual_information!, register!

"""
MutualInformationImageRegistration performs image registration using mutual information.

```jldoctest
using MutualInformationImageRegistration, MutualInformationImageRegistration.FastHistograms, Random

# Create the container used to hold intermediate variables for registration
mi = MutualInformationContainer(
    create_fast_histogram(
        FastHistograms.FixedWidth(),
        FastHistograms.Arithmetic(),
        FastHistograms.NoParallelization(),
        [(0x00, 0xff, 4), (0x00, 0xff, 4)],
    ),
)

# Create the full image that the smaller images to register will be pulled from
full_image = rand(UInt8, 500, 300)
view(full_image, 300:330, 200:220) .= 0xff

# The fixed image is the image that the other images are registered against
fixed = full_image[(300-10):(330+10), (200-10):(220+10)]

# The max shift is the maximum search range
MAX_SHIFT = 11
# Padding is used to grow the bbox for higher quality registration
padding = [-10, -10, 10, 10]

# A buffer is needed to hold intermediate data
buffer = Array{UInt8}(undef, (size(fixed) .+ (MAX_SHIFT * 2))...)

# Introduce a random shift to the moving bbox
expected_x = rand(-5:5)
expected_y = rand(-5:5)

# Register the image given by the bbox (called the moving bbox) against the fixed image
shift, mm, mms = register!(
    mi,
    full_image,
    fixed,
    [300, 200, 330, 220] .+ padding .+ [expected_x, expected_y, expected_x, expected_y],
    MAX_SHIFT,
    MAX_SHIFT,
    buffer
)

# The shift we get out should be equal and opposite of the shift we applied
shift == (-expected_x, -expected_y)

# output

true
```
"""
MutualInformationImageRegistration

end


text:
using NLSolvers, LinearAlgebra
n = 4
A = randn(n, n) + im * randn(n, n)
A = A'A + I
b = randn(n) + im * randn(n)
μ = 1.0

fcomplex(x) = real(dot(x, A * x) / 2 - dot(b, x)) + μ * sum(abs.(x) .^ 4)
gcomplex(x) = A * x - b + 4μ * (abs.(x) .^ 2) .* x
gcomplex!(stor, x) = copyto!(stor, gcomplex(x))
function fgcomplex!(stor, x)
    copyto!(stor, gcomplex(x))
    fcomplex(x), stor
end
obj = ScalarObjective(
    fcomplex,
    gcomplex!,
    fgcomplex!,
    nothing,
    nothing,
    nothing,
    nothing,
    nothing,
)
prob = OptimizationProblem(obj)
prob_oop = OptimizationProblem(obj; inplace = false)
x0 = randn(n) + im * randn(n)

####### WTF ! Inplace takes 1000 iters and oop doesn't
# solve(prob, x0, ParticleSwarm(), OptimizationOptions()) # fails
solve(prob, x0, SimulatedAnnealing(), OptimizationOptions())
solve(prob, copy(x0), LineSearch(LBFGS(), Backtracking()), OptimizationOptions()) # fails due to bounds
solve(
    prob,
    copy(x0),
    LineSearch(LBFGS(), Backtracking()),
    OptimizationOptions(g_abstol = 1e-6),
) # fails due to bounds
solve(prob, copy(x0), LineSearch(LBFGS()), OptimizationOptions(g_abstol = 1e-6)) # fails due to bounds
solve(prob, copy(x0), LineSearch(LBFGS(), HZAW()), OptimizationOptions(g_abstol = 1e-6)) # fails due to bounds
solve(
    prob,
    copy(x0),
    LineSearch(BFGS(), Backtracking()),
    OptimizationOptions(g_abstol = 1e-6),
) # fails due to bounds
solve(prob, copy(x0), LineSearch(BFGS()), OptimizationOptions(g_abstol = 1e-6)) # fails due to bounds
solve(prob, copy(x0), LineSearch(BFGS(), HZAW()), OptimizationOptions(g_abstol = 1e-6)) # fails due to bounds
solve(
    prob_oop,
    x0,
    LineSearch(BFGS(), Backtracking()),
    OptimizationOptions(g_abstol = 1e-6),
) # fails due to bounds
solve(prob_oop, x0, LineSearch(BFGS()), OptimizationOptions(g_abstol = 1e-6)) # fails due to bounds
solve(prob_oop, x0, LineSearch(BFGS(), HZAW()), OptimizationOptions(g_abstol = 1e-6)) # fails due to bounds

solve(
    prob,
    copy(x0),
    LineSearch(BFGS(), Backtracking()),
    OptimizationOptions(g_abstol = 1e-6),
) # fails due to bounds
solve(
    prob_oop,
    x0,
    LineSearch(BFGS(), Backtracking()),
    OptimizationOptions(g_abstol = 1e-6),
) # fails due to bounds
#=

@testset "Finite difference setup" begin
    oda1 = OnceDifferentiable(fcomplex, x0)
    fx, gx = NLSolversBase.value_gradient!(oda1, x0)
    @test fx == fcomplex(x0)
    @test gcomplex(x0) ≈ NLSolversBase.gradient(oda1)
end

@testset "Zeroth and second order methods" begin
    for method in (Optim.NelderMead,Optim.ParticleSwarm,Optim.Newton)
        @test_throws Any Optim.optimize(fcomplex, gcomplex!, x0, method())
    end
    #not supposed to converge, but it should at least go through without errors
    res = Optim.optimize(fcomplex, gcomplex!, x0, Optim.SimulatedAnnealing())
end

@testset "First order methods" begin
    options = Optim.Options(allow_f_increases=true)
    # TODO: AcceleratedGradientDescent fail to converge?
    for method in (Optim.GradientDescent(), Optim.ConjugateGradient(), Optim.LBFGS(), Optim.BFGS(),
                    Optim.NGMRES(), Optim.OACCEL(),Optim.MomentumGradientDescent(mu=0.1))

        debug_printing && printstyled("Solver: $(summary(method))\n", color=:green)
        res = Optim.optimize(fcomplex, gcomplex!, x0, method, options)
        debug_printing && printstyled("Iter\tf-calls\tg-calls\n", color=:green)
        debug_printing && printstyled("$(Optim.iterations(res))\t$(Optim.f_calls(res))\t$(Optim.g_calls(res))\n", color=:red)
        if !Optim.converged(res)
            @warn("$(summary(method)) failed.")
            display(res)
            println("########################")
        end
        ressum = summary(res) # Just check that no errors arise when doing display(res)
        @test typeof(fcomplex(x0)) == typeof(Optim.minimum(res))
        @test eltype(x0) == eltype(Optim.minimizer(res))
        @test Optim.converged(res)
        @test Optim.minimizer(res) ≈ xref rtol=1e-4

        res = Optim.optimize(fcomplex, x0, method, options)
        @test Optim.converged(res)
        @test Optim.minimizer(res) ≈ xref rtol=1e-4

        # # To compare with the equivalent real solvers
        # to_cplx(x) = x[1:n] + im*x[n+1:2n]
        # from_cplx(x) = [real(x);imag(x)]
        # freal(x) = fcomplex(to_cplx(x))
        # greal!(stor,x) = copyto!(stor, from_cplx(gcomplex(to_cplx(x))))
        # opt = Optim.Options(allow_f_increases=true,show_trace=true)
        # println("$(summary(method)) cplx")
        # res_cplx = Optim.optimize(fcomplex,gcomplex!,x0,method,opt)
        # println("$(summary(method)) real")
        # res_real = Optim.optimize(freal,greal!,from_cplx(x0),method,opt)
    end
end
=#

max_stars_repo_name: JuliaNLSolvers/NLSolvers

max_issues_repo_name: JuliaNLSolvers/NLSolvers

max_forks_repo_name: JuliaNLSolvers/NLSolvers