text:
using KitBase, FluxReconstruction, OrdinaryDiffEq, LinearAlgebra, Plots, OffsetArrays
using ProgressMeter: @showprogress
using Base.Threads: @threads

begin
    x0 = 0
    x1 = 1
    nx = 15
    y0 = 0
    y1 = 1
    ny = 15
    deg = 2
    nsp = deg + 1
    inK = 1
    γ = 5 / 3
    knudsen = 0.001
    muref = ref_vhs_vis(knudsen, 1.0, 0.5)
    cfl = 0.1
    dx = (x1 - x0) / nx
    dy = (y1 - y0) / ny
    dt = cfl * min(dx, dy) / (3.0)
    t = 0.0
    tmax = 0.15
    tspan = (0.0, tmax)
    nt = tmax ÷ dt |> Int
end

ps = FRPSpace2D(x0, x1, nx, y0, y1, ny, deg, 1, 1)

μᵣ = ref_vhs_vis(knudsen, 1.0, 0.5)
gas = Gas(knudsen, 0.0, 1.0, 1.0, γ, 0.81, 1.0, 0.5, μᵣ)

u0 = OffsetArray{Float64}(undef, 4, nsp, nsp, 0:ny+1, 0:nx+1)
#for i = 1:nsp, j = 1:nsp, k = 0:ny+1, l = 0:nx+1
for i = 0:nx+1, j = 0:ny+1, k = 1:nsp, l = 1:nsp
    u0[:, l, k, j, i] .= prim_conserve([1.0, 0.0, 0.0, 1.0], gas.γ)

    #ρ = max(exp(-50 * ((ps.xpg[l, k, j, i, 1] - 0.5)^2 + (ps.xpg[l, k, j, i, 2] - 0.5)^2)), 1e-2)
    #u0[:, i, j, k, l] .= prim_conserve([ρ, 0.0, 0.0, 1.0], gas.γ)

    #=if ps.x[i, j] < 0.5
        prim = [1.0, 0.0, 0.0, 0.5]
    else
        prim = [0.3, 0.0, 0.0, 0.625]
    end
    u0[:, l, k, j, i] .= prim_conserve(prim, gas.γ)=#
end

function KitBase.flux_gks!(
    fw::AbstractVector{T1},
    w::AbstractVector{T2},
    inK::Real,
    γ::Real,
    μᵣ::Real,
    ω::Real,
    sw = zero(w)::AbstractVector{T2},
) where {T1<:AbstractFloat,T2<:Real}
    prim = conserve_prim(w, γ)
    Mu, Mv, Mxi, MuL, MuR = gauss_moments(prim, inK)
    tau = vhs_collision_time(prim, μᵣ, ω)

    a = pdf_slope(prim, sw, inK)
    ∂ft = -prim[1] .* moments_conserve_slope(a, Mu, Mv, Mxi, 1, 0)
    A = pdf_slope(prim, ∂ft, inK)

    Muv = moments_conserve(Mu, Mv, Mxi, 1, 0, 0)
    Mau = moments_conserve_slope(a, Mu, Mv, Mxi, 2, 0)
    Mtu = moments_conserve_slope(A, Mu, Mv, Mxi, 1, 0)

    @. fw = prim[1] * (Muv - tau * Mau - tau * Mtu)

    return nothing
end

function KitBase.flux_gks!(
    fw::X,
    wL::Y,
    wR::Y,
    inK::Real,
    γ::Real,
    μᵣ::Real,
    ω::Real,
    dt::Real,
    swL = zero(wL)::Y,
    swR = zero(wR)::Y,
) where {X<:AbstractArray{<:AbstractFloat,1},Y<:AbstractArray{<:AbstractFloat,1}}
    primL = conserve_prim(wL, γ)
    primR = conserve_prim(wR, γ)

    Mu1, Mv1, Mxi1, MuL1, MuR1 = gauss_moments(primL, inK)
    Mu2, Mv2, Mxi2, MuL2, MuR2 = gauss_moments(primR, inK)

    w =
        primL[1] .* moments_conserve(MuL1, Mv1, Mxi1, 0, 0, 0) .+
        primR[1] .* moments_conserve(MuR2, Mv2, Mxi2, 0, 0, 0)
    prim = conserve_prim(w, γ)
    tau =
        vhs_collision_time(prim, μᵣ, ω) +
        2.0 * dt * abs(primL[1] / primL[end] - primR[1] / primR[end]) /
        (primL[1] / primL[end] + primR[1] / primR[end])

    if minimum(swL .* swR) < 0
        #swL .= 0.0
        #swR .= 0.0
    end

    faL = pdf_slope(primL, swL, inK)
    sw = -primL[1] .* moments_conserve_slope(faL, Mu1, Mv1, Mxi1, 1, 0)
    faTL = pdf_slope(primL, sw, inK)

    faR = pdf_slope(primR, swR, inK)
    sw = -primR[1] .* moments_conserve_slope(faR, Mu2, Mv1, Mxi2, 1, 0)
    faTR = pdf_slope(primR, sw, inK)

    Mu, Mv, Mxi, MuL, MuR = gauss_moments(prim, inK)

    # time-integration constants
    Mt = zeros(5)
    Mt[4] = dt#tau * (1.0 - exp(-dt / tau))
    Mt[5] = -tau * dt * exp(-dt / tau) + tau * Mt[4]
    Mt[1] = dt - Mt[4]
    Mt[2] = -tau * Mt[1] + Mt[5]
    Mt[3] = 0.5 * dt^2 - tau * Mt[1]

    # flux related to central distribution
    Muv = moments_conserve(Mu, Mv, Mxi, 1, 0, 0)
    fw .= Mt[1] .* prim[1] .* Muv

    # flux related to upwind distribution
    MuvL = moments_conserve(MuL1, Mv1, Mxi1, 1, 0, 0)
    MauL = moments_conserve_slope(faL, MuL1, Mv1, Mxi1, 2, 0)
    MauLT = moments_conserve_slope(faTL, MuL1, Mv1, Mxi1, 1, 0)

    MuvR = moments_conserve(MuR2, Mv2, Mxi2, 1, 0, 0)
    MauR = moments_conserve_slope(faR, MuR2, Mv2, Mxi2, 2, 0)
    MauRT = moments_conserve_slope(faTR, MuR2, Mv2, Mxi2, 1, 0)

    @. fw +=
        Mt[4] * primL[1] * MuvL - tau * Mt[4] * primL[1] * MauL -
        tau * Mt[4] * primL[1] * MauLT + Mt[4] * primR[1] * MuvR -
        tau * Mt[4] * primR[1] * MauR - tau * Mt[4] * primR[1] * MauRT
    fw ./= dt

    return nothing
end

function dudt!(du, u, p, t)
    boundary!(u, p, 1.0)

    fx,
    fy,
    ux_face,
    uy_face,
    fx_face,
    fy_face,
    fx_interaction,
    fy_interaction,
    rhs1,
    rhs2,
    ps,
    gas,
    dt = p

    nx = size(u, 5) - 2
    ny = size(u, 4) - 2
    nr = size(u, 3)
    ns = size(u, 2)

    @inbounds for i = 1:nx
        for j = 1:ny, k = 1:nr, l = 1:ns
            fw = @view fx[:, l, k, j, i]
            flux_gks!(fw, u[:, l, k, j, i], gas.K, gas.γ, gas.μᵣ, gas.ω, zeros(4))
            fw ./= ps.J[i, j][1]

            #fx[:, l, k, j, i] .= euler_flux(u[:, l, k, j, i], gas.γ)[1] ./ ps.J[i, j][1]
        end
    end
    @inbounds @threads for i = 1:nx
        for j = 1:ny, k = 1:nr, l = 1:ns
            fw = @view fy[:, l, k, j, i]
            ul = local_frame(u[:, l, k, j, i], 0.0, 1.0)
            flux_gks!(fw, ul, gas.K, gas.γ, gas.μᵣ, gas.ω, zeros(4))
            fy[:, l, k, j, i] .= global_frame(fw, 0.0, 1.0) ./ ps.J[i, j][2]

            #fy[:, l, k, j, i] .= euler_flux(u[:, l, k, j, i], gas.γ)[2] ./ ps.J[i, j][2]
        end
    end

    @inbounds for i = 0:nx+1
        for j = 1:ny, l = 1:ns, m = 1:4
            ux_face[m, 1, l, j, i] = dot(u[m, l, :, j, i], ps.ll)
            ux_face[m, 2, l, j, i] = dot(u[m, l, :, j, i], ps.lr)

            fx_face[m, 1, l, j, i] = dot(fx[m, l, :, j, i], ps.ll)
            fx_face[m, 2, l, j, i] = dot(fx[m, l, :, j, i], ps.lr)
        end
    end
    @inbounds for i = 1:nx
        for j = 0:ny+1, k = 1:nr, m = 1:4
            uy_face[m, 1, k, j, i] = dot(u[m, :, k, j, i], ps.ll)
            uy_face[m, 2, k, j, i] = dot(u[m, :, k, j, i], ps.lr)

            fy_face[m, 1, k, j, i] = dot(fy[m, :, k, j, i], ps.ll)
            fy_face[m, 2, k, j, i] = dot(fy[m, :, k, j, i], ps.lr)
        end
    end

    @inbounds for i = 1:nx+1
        for j = 1:ny, l = 1:ns
            swL = zeros(4)
            swR = zeros(4)
            for m in eachindex(swL)
                swL[m] = dot(u[m, l, :, j, i-1], ps.dll) / ps.J[i-1, j][1]
                swR[m] = dot(u[m, l, :, j, i], ps.dlr) / ps.J[i, j][1]
            end

            fw = @view fx_interaction[:, l, j, i]
            flux_gks!(
                fw,
                ux_face[:, 2, l, j, i-1],
                ux_face[:, 1, l, j, i],
                gas.K,
                gas.γ,
                gas.μᵣ,
                gas.ω,
                dt,
                swL,
                swR,
            )
            #=flux_hll!(
                fw,
                ux_face[:, 2, l, j, i-1],
                ux_face[:, 1, l, j, i],
                gas.γ,
                1.0,
            )=#
        end
    end
    @inbounds for i = 1:nx
        for j = 1:ny+1, k = 1:nr
            swL = zeros(4)
            swR = zeros(4)
            for m in eachindex(swL)
                swL[m] = dot(u[m, :, k, j-1, i], ps.dll) / ps.J[i, j-1][2]
                swR[m] = dot(u[m, :, k, j, i], ps.dlr) / ps.J[i, j][2]
            end

            fw = @view fy_interaction[:, k, j, i]
            uL = local_frame(uy_face[:, 2, k, j-1, i], 0.0, 1.0)
            uR = local_frame(uy_face[:, 1, k, j, i], 0.0, 1.0)

            flux_gks!(fw, uL, uR, gas.K, gas.γ, gas.μᵣ, gas.ω, dt, swL, swR)
            #=flux_hll!(
                fw,
                uL,
                uR,
                gas.γ,
                1.0,
            )=#
            fy_interaction[:, k, j, i] .= global_frame(fy_interaction[:, k, j, i], 0.0, 1.0)
        end
    end

    @inbounds for i = 1:nx, j = 1:ny, k = 1:nr, l = 1:ns, m = 1:4
        rhs1[m, l, k, j, i] = dot(fx[m, l, :, j, i], ps.dl[k, :])
        rhs2[m, l, k, j, i] = dot(fy[m, :, k, j, i], ps.dl[l, :])
    end

    @inbounds for i = 1:nx, j = 1:ny, k = 1:nr, l = 1:ns, m = 1:4
        du[m, l, k, j, i] =
            -(
                rhs1[m, l, k, j, i] +
                rhs2[m, l, k, j, i] +
                (fx_interaction[m, l, j, i] / ps.J[i, j][1] - fx_face[m, 1, l, j, i]) *
                ps.dhl[k] +
                (fx_interaction[m, l, j, i+1] / ps.J[i, j][1] - fx_face[m, 2, l, j, i]) *
                ps.dhr[k] +
                (fy_interaction[m, k, j, i] / ps.J[i, j][2] - fy_face[m, 1, k, j, i]) *
                ps.dhl[l] +
                (fy_interaction[m, k, j+1, i] / ps.J[i, j][2] - fy_face[m, 2, k, j, i]) *
                ps.dhr[l]
            )
    end
    du[:, :, :, :, 0] .= 0.0
    du[:, :, :, :, nx+1] .= 0.0
    du[:, :, :, 0, :] .= 0.0
    du[:, :, :, ny+1, :] .= 0.0

    return nothing
end

function boundary!(u, p, λ0)
    gas = p[end-1]

    nx = size(u, 5) - 2
    ny = size(u, 4) - 2
    nr = size(u, 3)
    ns = size(u, 2)

    pb = zeros(4)
    for j = 1:ny, k = 1:nr, l = 1:ns
        prim = conserve_prim(u[:, l, k, j, 1], gas.γ)

        pb[end] = 2 * λ0 - prim[end]
        tmp = (prim[end] - λ0) / λ0
        pb[1] = (1.0 - tmp) / (1.0 + tmp) * prim[1]
        pb[2] = -prim[2]
        pb[3] = -prim[3]

        u[:, l, nr+1-k, j, 0] .= prim_conserve(pb, gas.γ)
    end
    for j = 1:ny, k = 1:nr, l = 1:ns
        prim = conserve_prim(u[:, l, k, j, nx], gas.γ)

        pb[end] = 2 * λ0 - prim[end]
        tmp = (prim[end] - λ0) / λ0
        pb[1] = (1 - tmp) / (1 + tmp) * prim[1]
        pb[2] = -prim[2]
        pb[3] = -prim[3]

        u[:, l, nr+1-k, j, nx+1] .= prim_conserve(pb, gas.γ)
    end
    for i = 1:nx, k = 1:nr, l = 1:ns
        prim = conserve_prim(u[:, l, k, 1, i], gas.γ)

        pb[end] = 2 * λ0 - prim[end]
        tmp = (prim[end] - λ0) / λ0
        pb[1] = (1 - tmp) / (1 + tmp) * prim[1]
        pb[2] = -prim[2]
        pb[3] = -prim[3]

        u[:, ns+1-l, k, 0, i] .= prim_conserve(pb, gas.γ)
    end
    for i = 1:nx, k = 1:nr, l = 1:ns
        prim = conserve_prim(u[:, l, k, ny, i], gas.γ)

        pb[end] = 2 * λ0 - prim[end]
        tmp = (prim[end] - λ0) / λ0
        pb[1] = (1 - tmp) / (1 + tmp) * prim[1]
        pb[2] = 0.15#-prim[2] + 0.3
        pb[3] = -prim[3]

        u[:, ns+1-l, k, ny+1, i] .= prim_conserve(pb, gas.γ)
    end

    return nothing
end

begin
    du = zero(u0)
    fx = zero(u0)
    fy = zero(u0)
    ux_face = OffsetArray{Float64}(undef, 4, 2, ps.deg + 1, ps.ny, 0:ps.nx+1) |> zero
    uy_face = OffsetArray{Float64}(undef, 4, 2, ps.deg + 1, 0:ps.ny+1, ps.nx) |> zero
    fx_face = zero(ux_face)
    fy_face = zero(uy_face)
    fx_interaction = zeros(4, ps.deg + 1, ps.ny, ps.nx + 1)
    fy_interaction = zeros(4, ps.deg + 1, ps.ny + 1, ps.nx)
    rhs1 = zero(u0)
    rhs2 = zero(u0)
end

p = (
    fx,
    fy,
    ux_face,
    uy_face,
    fx_face,
    fy_face,
    fx_interaction,
    fy_interaction,
    rhs1,
    rhs2,
    ps,
    gas,
    dt,
)

#u = deepcopy(u0)
#dudt!(du, u, p, 0.0)

prob = ODEProblem(dudt!, u0, tspan, p)
itg = init(prob, Euler(), save_everystep = false, adaptive = false, dt = dt)

@showprogress for iter = 1:100#nt
    step!(itg)
end

contourf(ps.xpg[1:nx, 1, 1, 1, 1], ps.xpg[1, 1:ny, 1, 1, 2], itg.u[1, 1, 1, 1:ny, 1:nx])
contourf(ps.xpg[1:nx, 1, 1, 1, 1], ps.xpg[1, 1:ny, 1, 1, 2], itg.u[2, 1, 1, 1:ny, 1:nx])
contourf(ps.xpg[1:nx, 1, 1, 1, 1], ps.xpg[1, 1:ny, 1, 1, 2], itg.u[4, 1, 1, 1:ny, 1:nx])

plot(ps.xpg[1:nx, 1, 1, 1, 1], itg.u[1, 1, 1, ny÷2, 1:nx])
plot(ps.xpg[1:nx, 1, 1, 1, 1], itg.u[2, 1, 1, ny÷2, 1:nx])

plot(ps.xpg[1, 1:ny, 1, 1, 2], itg.u[1, 1, 1, 1:ny, nx÷2])
plot(ps.xpg[1, 1:ny, 1, 1, 2], itg.u[3, 1, 1, 1:ny, nx÷2])

begin
    coord = zeros(nx * nsp, ny * nsp, 2)
    prim = zeros(nx * nsp, ny * nsp, 4)
    for i = 1:nx, j = 1:ny
        idx0 = (i - 1) * nsp
        idy0 = (j - 1) * nsp

        for k = 1:nsp, l = 1:nsp
            idx = idx0 + k
            idy = idy0 + l
            coord[idx, idy, 1] = ps.xpg[i, j, k, l, 1]
            coord[idx, idy, 2] = ps.xpg[i, j, k, l, 2]

            _w = itg.u[:, l, k, j, i]
            prim[idx, idy, :] .= conserve_prim(_w, γ)
            prim[idx, idy, 4] = 1 / prim[idx, idy, 4]
        end
    end
end

contourf(coord[:, 1, 1], coord[1, :, 2], prim[:, :, 2]')

plot(coord[:, 1, 1], prim[:, 1, 1])
plot(coord[:, 1, 1], prim[:, 1, 4])

begin
    using PyCall
    itp = pyimport("scipy.interpolate")

    x_uni =
        coord[1, 1, 1]:(coord[end, 1, 1]-coord[1, 1, 1])/(nx*nsp-1):coord[end, 1, 1] |>
        collect
    y_uni =
        coord[1, 1, 2]:(coord[1, end, 2]-coord[1, 1, 2])/(ny*nsp-1):coord[1, end, 2] |>
        collect

    n_ref = itp.interp2d(coord[:, 1, 1], coord[1, :, 2], prim[:, :, 1], kind = "cubic")
    n_uni = n_ref(x_uni, y_uni)

    u_ref = itp.interp2d(coord[:, 1, 1], coord[1, :, 2], prim[:, :, 2], kind = "cubic")
    u_uni = u_ref(x_uni, y_uni)

    v_ref = itp.interp2d(coord[:, 1, 1], coord[1, :, 2], prim[:, :, 3], kind = "cubic")
    v_uni = v_ref(x_uni, y_uni)

    t_ref = itp.interp2d(coord[:, 1, 1], coord[1, :, 2], prim[:, :, 4], kind = "cubic")
    t_uni = t_ref(x_uni, y_uni)
end

contourf(x_uni, y_uni, n_uni')
contourf(x_uni, y_uni, u_uni')
contourf(x_uni, y_uni, t_uni')

plot(x_uni[:, 1], n_uni[:, 1])


text:
#---------------------------------------------------------------------#
#This function computes the Initial & Exact Solutions.
#Written by F.X. Giraldo on April 19, 2019
#           Department of Applied Mathematics
#           Naval Postgraduate School
#           Monterey; CA 93943-5216
#---------------------------------------------------------------------#
function exact_solution(coord,Npoin,time,case,DFloat)

    #Set some constants
    xmin=-1
    xmax=+1
    xc=0.5*(xmax+xmin)
    xl=xmax-xmin
    rc=0.125
    sigma=1/16; #best for right convergence rates
    # sigma=1/8
    # sigma=1/4; #fatter Gaussian to compare against 2D code
    u=DFloat(2.0)

    #Initialize
    qe=zeros(DFloat,Npoin)

    #timec=time - floor(time)
    timec=time
    #Generate Grid Points
    for I=1:Npoin
        x=coord[I]
        xbar=xc + u*timec
        if (xbar > xmax)
            xbar=xmin + (xbar-xmax)
        end
        r=x-xbar
        if (case == 1) #Gaussian
            qe[I]=exp( -(x-xbar)^2/(2*sigma)^2 )
        elseif (case == 2) #Square Wave
      	    if ( abs(r) <= rc )
                qe[I]=1
            end
        end
    end #I
    return (qe,u)
end

ext: jl


text: println("Hello from included.jl!")


text:
using Documenter, ExtensibleMCMC

makedocs(;
    modules=[ExtensibleMCMC],
    format=Documenter.HTML(
        mathengine = Documenter.MathJax(
            Dict(
                :TeX => Dict(
                    :equationNumbers => Dict(
                        :autoNumber => "AMS"
                    ),
                    :Macros => Dict(
                        :dd => "{\\textrm d}",
                        :RR => "\\mathbb{R}",
                        :wt => ["\\widetilde{#1}", 1]
                    ),
                )
            )
        ),
        collapselevel = 1,
    ),
    pages=[
        "Home" => "index.md",
        "Get started" => joinpath("get_started", "basic_use.md"),
        "User manual" => Any[
            "Internal structure" => joinpath("manual", "internal_structure.md"),
            "MCMC" => joinpath("manual", "mcmc.md"),
            "Updates & Decorators" => joinpath("manual", "updates_and_decorators.md"),
            "Callbacks" => joinpath("manual", "callbacks.md"),
            "Priors" => joinpath("manual", "priors.md"),
            "Workspaces" => joinpath("manual", "workspaces.md"),
            "MCMC Schedule" => joinpath("manual", "mcmc_schedule.md"),
        ],
        "How to..." => Any[
            "(TODO) First how to" => joinpath("how_to_guides", "first_how_to.md"),
            "(TODO) Definy my own workspace" => joinpath("how_to_guides", "define_my_own_workspace.md"),
        ],
        "Tutorials" => Any[
            "Estimate mean of a bivariate Gaussian" => joinpath("tutorials", "mean_of_bivariate_gaussian.md"),
        ],
        "Index" => "module_index.md"
    ],
    repo="https://github.com/JuliaDiffusionBayes/ExtensibleMCMC.jl/blob/{commit}{path}#L{line}",
    sitename="ExtensibleMCMC.jl",
    authors="Sebastiano Grazzi, Frank van der Meulen, Marcin Mider, Moritz Schauer",
    #assets=String[],
)

deploydocs(;
    repo="github.com/JuliaDiffusionBayes/ExtensibleMCMC.jl",
)


text:
using ChainRulesCore
using ChainRulesTestUtils

function two2three(x1::Float64, x2::Float64)
    return 1.0, 2.0 * x1, 3.0 * x2
end

function ChainRulesCore.frule((Δf, Δx1, Δx2), ::typeof(two2three), x1, x2)
    y = two2three(x1, x2)
    ∂y = Tangent{Tuple{Float64, Float64, Float64}}(
        ZeroTangent(),
        2.0 * Δx1,
        3.0 * Δx2,
    )
    return y, ∂y
end

function ChainRulesCore.rrule(::typeof(two2three), x1, x2)
    y = two2three(x1, x2)
    function two2three_pullback(Ȳ)
        return (NoTangent(), 2.0 * Ȳ[2], 3.0 * Ȳ[3])
    end
    return y, two2three_pullback
end

max_stars_repo_name: sandreza/GradientLearning

max_issues_repo_name: sandreza/GradientLearning

max_forks_repo_name: sandreza/GradientLearning


text:
is_constant_cache(cache::OrdinaryDiffEqConstantCache) = true
is_constant_cache(cache::OrdinaryDiffEqCache) = false
is_constant_cache(cache::CompositeCache) = is_constant_cache(cache.caches[1])


text:
using DataFrames

# Accuracy holds classification accuracy info
# which can be used to compute precision, recall, etc.
type Accuracy
   predictions::Int
   correct::Int
   incorrect::Int
   prediction::Array{Int}
   classes::Int

   Accuracy(classes) = new(0,0,0,round(Int,zeros(classes,classes)), classes)
end

# predict! adds a prediction to the accuracy
# type
function predict!(acc::Accuracy, predicted::Int, actual::Int)
    @assert predicted <= acc.classes
    @assert predicted > 0
    @assert actual <= acc.classes
    @assert actual > 0
    acc.predictions += 1
    if predicted == actual
        acc.correct += 1
    else
        acc.incorrect += 1
    end
    acc.prediction[actual,predicted] += 1
end

# accuracy prints and returns mean accuracy
function accuracy(acc::Accuracy)
    println("==> Mean Accuracy:  $(acc.correct/acc.predictions)")
    return acc.correct/acc.predictions
end

# precision returns and prints individual
# class precision and prints mean precision
function precision(acc::Accuracy)
    prec = diag(acc.prediction)./sum(acc.prediction,1)'
    println("==> Mean Precision:  $(mean(prec))")
    println("==> Class Precision: $(prec')")
    return prec
end

# recall returns and prints individual
# class recall and prints mean recall
function recall(acc::Accuracy)
    rec = diag(acc.prediction)./sum(acc.prediction,2)
    println("==> Mean Recall:  $(mean(rec))")
    println("==> Class Recall: $(rec')")
    return rec
end

max_stars_repo_name: cdipaolo/inflection

max_issues_repo_name: cdipaolo/inflection

max_forks_repo_name: cdipaolo/inflection


text:
"""
    f∂f(f,x)

returns f and ∂f/∂x evaluated in `x`, using `ForwardDiff.jl`, `DiffResults.jl` and `StaticArrays.jl` to calculate everything in one pass.
"""
function f∂f(f,x::T) where T
    _f(z) = f(only(z))
    x_vec =   SVector(x)
    ∂result = DiffResults.GradientResult(x_vec)
    _∂f =  ForwardDiff.gradient!(∂result, _f,x_vec)
    fx =  DiffResults.value(_∂f)
    ∂f∂x = only(DiffResults.gradient(_∂f))
    return fx,∂f∂x
end

"""
    f∂f∂2f(f,x)

returns f,∂f/∂x,and ∂²f/∂²x and evaluated in `x`, using `ForwardDiff.jl`, `DiffResults.jl` and `StaticArrays.jl` to calculate everything in one pass.
"""
function f∂f∂2f(f,x::T) where T
    _f(z) = f(only(z))
    x_vec =   SVector(x)
    ∂result = DiffResults.HessianResult(x_vec)
    _∂f =  ForwardDiff.hessian!(∂result, _f,x_vec)
    fx =  DiffResults.value(_∂f)
    ∂f∂x = only(DiffResults.gradient(_∂f))
    ∂²f∂²x =  only(DiffResults.hessian(_∂f))
    return fx,∂f∂x,∂²f∂²x
end

"""
    ad_newton(f,x0;kwargs...)

    Performs newton root finding (via Roots.jl), but the derivatives are calculated with ForwardDiff.jl
    kwargs are passed to `Roots.find_zero`
"""
function ad_newton(f::F,x0;xatol=nothing, xrtol=nothing, maxevals = 100) where F
    function f0(_x)
        fx,∂f∂x = f∂f(f,_x)
        return fx,fx/∂f∂x
    end
    x = float(x0)
    T = typeof(x)
    atol = xatol !== nothing ? xatol : oneunit(T) * (eps(one(T)))^(4/5)
    rtol = xrtol !== nothing ? xrtol : eps(one(T))^(4/5)
    xo = T(Inf)
    for i in 1:maxevals
        fx, Δx = f0(x)
        iszero(fx) && return x

        x -= Δx

        if isapprox(x, xo, atol=atol, rtol=rtol)
            return x
        end
        xo = x
    end
    return zero(x0)/zero(x0)
end


text:
# Distributions with compound variates
#
#	NormalGamma
#	NormalInverseGamma
#	NormalWishart
#	NormalInverseWishart
#
# and their applications for conjugates & posteriors
#

using Base.Test
using Distributions

### NormalGamma

n = 100

mu0 = 2.
nu0 = 3.
shape0 = 5.
rate0 = 2.
ng = NormalGamma(mu0, nu0, shape0, rate0)

# Random number generation
mu, tau2 = rand(ng)

# Did it generate something valid?
@test insupport(NormalGamma, mu, tau2)

mu = 2.5
tau2 = 3.

# pdf
npdf = pdf(Normal(mu0, 1./sqrt(nu0*tau2)), mu)
gpdf = pdf(Gamma(shape0, 1./rate0), tau2)
lnpdf = logpdf(Normal(mu0, 1./sqrt(nu0*tau2)), mu)
lgpdf = logpdf(Gamma(shape0, 1./rate0), tau2)

@test_approx_eq_eps pdf(ng, mu, tau2) (npdf*gpdf) 1e-8
@test_approx_eq_eps logpdf(ng, mu, tau2) (lnpdf+lgpdf) 1e-8

### NormalInverseGamma

n = 100

mu0 = 2.
v0 = 3.
shape0 = 5.
scale0 = 2.
nig = NormalInverseGamma(mu0, v0, shape0, scale0)

# Random number generation
mu, sig2 = rand(nig)

# Did it generate something valid?
@test insupport(NormalInverseGamma, mu, sig2)

mu = 2.5
sig2 = 3.

# pdf
npdf = pdf(Normal(mu0, sqrt(v0*sig2)), mu)
gpdf = pdf(InverseGamma(shape0, scale0), sig2)
lnpdf = logpdf(Normal(mu0, sqrt(v0*sig2)), mu)
lgpdf = logpdf(InverseGamma(shape0, scale0), sig2)

@test_approx_eq_eps pdf(nig, mu, sig2) (npdf*gpdf) 1e-8
@test_approx_eq_eps logpdf(nig, mu, sig2) (lnpdf+lgpdf) 1e-8

### NormalWishart

n = 100

mu0 = [2., 3.]
kappa0 = 3.
nu0 = 4.
T0 = eye(2)
T0[1,2] = T0[2,1] = .5
nw = NormalWishart(mu0, kappa0, T0, nu0)

# Random number generation
mu, Lam = rand(nw)

# Did it generate something valid?
@test insupport(NormalWishart, mu, Lam)

mu = [1.5, 2.5]
T = 0.75*eye(2)
T[1,2] = T[2,1] = 0.6

# pdf
npdf = pdf(MultivariateNormal(mu0, inv(kappa0*T)), mu)
wpdf = pdf(Wishart(nu0, T0), T)
lnpdf = logpdf(MultivariateNormal(mu0, inv(kappa0*T)), mu)
lwpdf = logpdf(Wishart(nu0, T0), T)

@test_approx_eq_eps pdf(nw, mu, T) (npdf*wpdf) 1e-8
@test_approx_eq_eps logpdf(nw, mu, T) (lnpdf+lwpdf) 1e-8

### NormalInverseWishart

n = 100

mu0 = [2., 3.]
kappa0 = 3.
nu0 = 4.
T0 = eye(2)
T0[1,2] = T0[2,1] = .5
niw = NormalInverseWishart(mu0, kappa0, T0, nu0)

# Random number generation
mu, Sig = rand(niw)

# Did it generate something valid?
@test insupport(NormalInverseWishart, mu, Sig)

mu = [1.5, 2.5]
T = 0.75*eye(2)
T[1,2] = T[2,1] = 0.6

# pdf
npdf = pdf(MultivariateNormal(mu0, 1/kappa0*T), mu)
wpdf = pdf(InverseWishart(nu0, T0), T)
lnpdf = logpdf(MultivariateNormal(mu0, 1/kappa0*T), mu)
lwpdf = logpdf(InverseWishart(nu0, T0), T)

@test_approx_eq_eps pdf(niw, mu, T) (npdf*wpdf) 1e-8
@test_approx_eq_eps logpdf(niw, mu, T) (lnpdf+lwpdf) 1e-8


text:
function odes!(dc, c, p, t)
    # unpack parameters
    imax, Ns, kpij, kmij, kp, kb, kub, kn, kg1, kd1, kg2, kd2, ka, ke, kc  = p

    # Define important indices for 2D array
    rows    = imax + 3
    columns = Ns[imax] + 3
    istart  = 3
    idimer  = 4
    iend    = rows - 1
    j0      = 2
    jmax    = columns - 1

    # array element (1,1) is assigned to M+ or c1,0+
    dc[1,1] = -kp*c[1,1]

    # array element (3,2) is assigned to M or c1,0
    c_kg2 = 0.0
    for j = j0:jmax
        for i = idimer:iend-1
            c_kg2 += kpij[i,j]*c[i,j]
        end
    end
    dc[3,2] = kp*c[1,1] - kb*c[3,2]*c[2,3] + kub*c[3,3] -
        2*kn*c[3,2]*c[3,2] - kn*c[3,2]*c[3,3] - kg2*c[3,2]*c_kg2

    # array element (3,3) is assigned to ML or c1,1
    c_kg1 = 0.0
    for j = j0:jmax
        for i = idimer:iend-1
            c_kg1 += kpij[i,j]*c[i,j]
        end
    end
    c_kd1 = 0.0
    for j = j0:jmax
        for i = idimer+1:iend
            c_kd1 += kmij[i,j]*c[i,j]
        end
    end
    dc[3,3] = kb*c[3,2]*c[2,3] - kub*c[3,3] - 2*kn*c[3,3]*c[3,3] - kn*c[3,3]*c[3,2] -
        kg1*c[3,3]*c_kg1 + kd1*c_kd1

    # array element (2,3) is assigned to L or c0,1
    c_ka = 0.0
    for j = j0:jmax
        for i = idimer:iend
            c_ka += kpij[i,j]*c[i,j]
        end
    end
    c_ke = 0.0
    for j = j0:jmax
        for i = idimer:iend
            c_ke += kmij[i,j]*c[i,j]
        end
    end
    dc[2,3] = -kb*c[3,2]*c[2,3] + kub*c[3,3] - ka*c[2,3]*c_ka + ke*c_ke

    # nested loop provides expression for ci,j starting at c2,0
    # ligands are the outer loop for memory efficiency in julia
    for j = j0:Ns[imax]+2
        for i = idimer:iend
            # # of monomers is index i minus 2
            # # of ligands on NC is the index j minus 2
            monomers = i - 2
            ligands  = j - 2
            # don't access array elements that correspond to non-existent clusters
            if ( ligands <= Ns[monomers] )

                # sum up clusters consumed by coalescence c_cns = kc * kpij *ci,j * sum_k,l kpkl * ck,l
                # k, l are dummy indices
                # conditions: i+k <= imax, j+l <= Ns[i+k]
                c_cns = 0.0
                if i < iend - 1
                for l = j0:jmax
                    for k = idimer:iend-monomers+2
                        if ( (i-2)+(k-2) <= imax)
                            if ( (j-2)+(l-2) <= Ns[(i-2)+(k-2)] )
                                c_cns += kpij[i,j]*c[i,j]*kpij[k,l]*c[k,l]
                            end
                        end
                    end
                end
                end

                # sum up clusters produced by coalescence c_prd = kc * sum_w,y sum_x,z kpwx * cw,x * kpyz * cy,z
                # w, x, y, z are dummy indices
                # conditions: i=w+y, j=x+z
                # sum over-counts by factor of 4, so divide by 4 later
                c_prd = 0.0
                if i > idimer + 1
                for x = j0:jmax
                    for z = j0:jmax
                        for w = idimer:i
                            for y = idimer:i
                                if ( (i-2) == (w-2)+(y-2) )
                                    if ( (j-2) == (x-2)+(z-2) )
                                        c_prd += kpij[w,x]*c[w,x]*kpij[y,z]*c[y,z]
                                    end
                                end
                            end
                        end
                    end
                end
                end

                # if clause provides special cases for ci,j
                if i == idimer
                    if j == j0
                        dc[i,j] = kn*c[3,2]*c[3,2] -
                            kg1*kpij[i,j]*c[3,3]*c[i,j] + kd1*kmij[i+1,j+1]*c[i+1,j+1] -
                            kg2*kpij[i,j]*c[3,2]*c[i,j] -
                            ka*kpij[i,j]*c[2,3]*c[i,j] +
                            ke*kmij[i,j+1]*c[i,j+1] -
                            kc*c_cns
                    elseif j == j0 + 1
                        dc[i,j] = kn*c[3,3]*c[3,2] -
                            kg1*kpij[i,j]*c[3,3]*c[i,j] + kd1*kmij[i+1,j+1]*c[i+1,j+1] -
                            kg2*kpij[i,j]*c[3,2]*c[i,j] +
                            ka*kpij[i,j-1]*c[2,3]*c[i,j-1] - ka*kpij[i,j]*c[2,3]*c[i,j] +
                            ke*kmij[i,j+1]*c[i,j+1] - ke*kmij[i,j]*c[i,j] -
                            kc*c_cns
                    elseif j == j0 + 2
                        dc[i,j] = kn*c[3,3]*c[3,3] -
                            kg1*kpij[i,j]*c[3,3]*c[i,j] + kd1*kmij[i+1,j+1]*c[i+1,j+1] -
                            kg2*kpij[i,j]*c[3,2]*c[i,j] +
                            ka*kpij[i,j-1]*c[2,3]*c[i,j-1] - ka*kpij[i,j]*c[2,3]*c[i,j] +
                            ke*kmij[i,j+1]*c[i,j+1] - ke*kmij[i,j]*c[i,j] -
                            kc*c_cns
                    elseif j == j0 + 3
                        dc[i,j] = -kg1*kpij[i,j]*c[3,3]*c[i,j] + kd1*kmij[i+1,j+1]*c[i+1,j+1] -
                            kg2*kpij[i,j]*c[3,2]*c[i,j] +
                            ka*kpij[i,j-1]*c[2,3]*c[i,j-1] - ka*kpij[i,j]*c[2,3]*c[i,j] +
                            ke*kmij[i,j+1]*c[i,j+1] - ke*kmij[i,j]*c[i,j] -
                            kc*c_cns
                    end
                elseif i == idimer + 1
                    dc[i,j] = kg1*kpij[i-1,j-1]*c[3,3]*c[i-1,j-1] - kg1*kpij[i,j]*c[3,3]*c[i,j] +
                        kd1*kmij[i+1,j+1]*c[i+1,j+1] - kd1*kmij[i,j]*c[i,j] +
                        kg2*kpij[i-1,j]*c[3,2]*c[i-1,j] - kg2*kpij[i,j]*c[3,2]*c[i,j] +
                        ka*kpij[i,j-1]*c[2,3]*c[i,j-1] - ka*kpij[i,j]*c[2,3]*c[i,j] +
                        ke*kmij[i,j+1]*c[i,j+1] - ke*kmij[i,j]*c[i,j] -
                        kc*c_cns
                elseif i > idimer + 1 && i < iend - 1
                    dc[i,j] = kg1*kpij[i-1,j-1]*c[3,3]*c[i-1,j-1] - kg1*kpij[i,j]*c[3,3]*c[i,j] +
                        kd1*kmij[i+1,j+1]*c[i+1,j+1] - kd1*kmij[i,j]*c[i,j] +
                        kg2*kpij[i-1,j]*c[3,2]*c[i-1,j] - kg2*kpij[i,j]*c[3,2]*c[i,j] +
                        ka*kpij[i,j-1]*c[2,3]*c[i,j-1] - ka*kpij[i,j]*c[2,3]*c[i,j] +
                        ke*kmij[i,j+1]*c[i,j+1] - ke*kmij[i,j]*c[i,j] -
                        kc*c_cns +
                        kc*c_prd/4.0
                elseif i == iend - 1
                    dc[i,j] = kg1*kpij[i-1,j-1]*c[3,3]*c[i-1,j-1] - kg1*kpij[i,j]*c[3,3]*c[i,j] +
                        kd1*kmij[i+1,j+1]*c[i+1,j+1] - kd1*kmij[i,j]*c[i,j] +
                        kg2*kpij[i-1,j]*c[3,2]*c[i-1,j] - kg2*kpij[i,j]*c[3,2]*c[i,j] +
                        ka*kpij[i,j-1]*c[2,3]*c[i,j-1] - ka*kpij[i,j]*c[2,3]*c[i,j] +
                        ke*kmij[i,j+1]*c[i,j+1] - ke*kmij[i,j]*c[i,j] +
                        kc*c_prd/4.0
                elseif i == iend
                    dc[i,j] = kg1*kpij[i-1,j-1]*c[3,3]*c[i-1,j-1] +
                        kd1*kmij[i+1,j+1]*c[i+1,j+1] - kd1*kmij[i,j]*c[i,j] +
                        kg2*kpij[i-1,j]*c[3,2]*c[i-1,j] +
                        ka*kpij[i,j-1]*c[2,3]*c[i,j-1] - ka*kpij[i,j]*c[2,3]*c[i,j] +
                        ke*kmij[i,j+1]*c[i,j+1] - ke*kmij[i,j]*c[i,j] +
                        kc*c_prd/4.0
                end
            end
        end
    end
end


text:
using Complementarity
using CSV
using DataFrames
using JuMP

cge = MCPModel();

smallval = 0.0

@variable(cge,X>=smallval,start=1)
@variable(cge,Y>=smallval,start=1)
@variable(cge,W>=smallval,start=1)
@variable(cge,PX>=smallval,start=1)
@variable(cge,PY>=smallval,start=1)
@variable(cge,PW>=smallval,start=1)
@variable(cge,PL>=smallval,start=1)
@variable(cge,PK>=smallval,start=1)
@variable(cge,CONS>=smallval,start=1)

#parameters for counterfactuals
TX = 0.0
LENDOW = 1.0

#need to be indexed by s
@mapping(cge,PRF_X,100 * PL^0.25 * PK^0.75 * (1+TX) - 100*PX)
@mapping(cge,PRF_Y,100 * PL^0.75* PK^0.25 - 100 * PY)
@mapping(cge,PRF_W,200 * PX^0.5 * PY^0.5 - 200 * PW)

#need to be indexed by m
@mapping(cge,MKT_X, 100 * X - (100 * W * PX^0.5 * PY^0.5) / PX)
@mapping(cge,MKT_Y, 100 * Y - (100 * W * PX^0.5 * PY^0.5) / PY)
@mapping(cge,MKT_W, 200 * W - CONS / PW)
@mapping(cge,MKT_L, 100 * LENDOW - ((25 * X * (PL^0.25) * (PK^0.75)) / PL + (75 * Y * PL^0.75 * PK^0.25) / PL))
@mapping(cge,MKT_K, 100 - (75 * X * PL^0.25 * PK^0.75 / PK + 25 * Y * PL^0.75 * PK^0.25 / PK))

#INCOME Constraint
@mapping(cge,I_CONS, CONS == (100*LENDOW*PL + TX * 100 * X * (PL^0.25) * (PK^0.75) ) )

#@complementarity(d, capcon, pk);

@complementarity(cge,PRF_X,X)
@complementarity(cge,PRF_Y,Y)
@complementarity(cge,PRF_W,W)
@complementarity(cge,MKT_X,PX)
@complementarity(cge,MKT_Y,PY)
@complementarity(cge,MKT_L,PL)
@complementarity(cge,MKT_K,PK)
@complementarity(cge,MKT_W,PW)
@complementarity(cge,I_CONS,CONS)

set_start_value(X,1)
set_start_value(Y,1)
set_start_value(W,1)
set_start_value(PX,1)
set_start_value(PY,1)
set_start_value(PW,1)
set_start_value(PL,1)
set_start_value(PK,1)
set_start_value(CONS,200)

PATHSolver.options(convergence_tolerance=1e-8, output=:yes, time_limit=3600)
status = solveMCP(cge)

#retrieve the values of x
result_value(X)
result_value(Y)
result_value(W)
result_value(P)
result_value(PY)
result_value(PW)
result_value(PL)
result_value(PK)
result_value(CONS)

max_stars_repo_name: matbesancon/SLiDE

max_issues_repo_name: matbesancon/SLiDE

max_forks_repo_name: matbesancon/SLiDE


text:
"Holds the tableau of an explicit Runge-Kutta method."
struct TableauERK{T} <: AbstractTableauERK{T}
    @HeaderTableau

    q::CoefficientsRK{T}

    function TableauERK{T}(q) where {T}
        @assert q.c[1] == 0
        @assert istrilstrict(q.a)
        @assert !(q.s==1 && q.a[1,1] ≠ 0)
        new(q.name, q.o, q.s, q)
    end
end

function TableauERK(q::CoefficientsRK{T}) where {T}
    TableauERK{T}(q)
end

function TableauERK(name::Symbol, order::Int, a::Matrix{T}, b::Vector{T}, c::Vector{T}) where {T}
    TableauERK{T}(CoefficientsRK(name, order, a, b, c))
end

"Read explicit Runge-Kutta tableau from file."
function readTableauERKFromFile(dir::AbstractString, name::AbstractString)
    file = string(dir, "/", name, ".tsv")

    o, s, T = readTableauRKHeaderFromFile(file)

    # TODO Read data in original format (e.g., Rational).
    #      For this we need to save tableaus as jld or hdf5.
#    tab_array = readdlm(file, T)
    tab_array = readdlm(file, comments=true)

    if s == 0
        s = size(tab_array, 1)
    end

    @assert s == size(tab_array, 1) == size(tab_array, 2)-2

    c = tab_array[1:s, 1]
    b = tab_array[1:s, 2]
    a = tab_array[1:s, 3:s+2]

    get_config(:verbosity) > 1 ? @info("Reading explicit Runge-Kutta tableau $(name) with $(s) stages and order $(o) from file\n$(file)") : nothing

    TableauERK(Symbol(name), o, a, b, c)
end

"Parameters for right-hand side function of explicit Runge-Kutta methods."
struct ParametersERK{DT, TT, ET <: ODE{DT,TT}, D, S} <: Parameters{DT,TT}
    equ::ET
    tab::TableauERK{TT}
    Δt::TT
end

function ParametersERK(equ::ET, tab::TableauERK{TT}, Δt::TT) where {DT, TT, ET <: ODE{DT,TT}}
    ParametersERK{DT, TT, ET, equ.d, tab.s}(equ, tab, Δt)
end

"Explicit Runge-Kutta integrator cache."
struct IntegratorCacheERK{DT,D,S} <: ODEIntegratorCache{DT,D}
    Q::Vector{Vector{DT}}
    V::Vector{Vector{DT}}

    function IntegratorCacheERK{DT,D,S}() where {DT,D,S}
        Q = create_internal_stage_vector(DT, D, S)
        V = create_internal_stage_vector(DT, D, S)
        new(Q, V)
    end
end

"Explicit Runge-Kutta integrator."
struct IntegratorERK{DT, TT, PT <: ParametersERK{DT,TT}, D, S} <: IntegratorRK{DT,TT}
    params::PT
    cache::IntegratorCacheERK{DT,D,S}

    function IntegratorERK(equation::ODE{DT,TT,FT}, tableau::TableauERK{TT}, Δt::TT) where {DT,TT,FT}
        D = equation.d
        S = tableau.s

        # create params
        params = ParametersERK(equation, tableau, Δt)

        # create cache
        cache = IntegratorCacheERK{DT,D,S}()

        # create integrators
        new{DT, TT, typeof(params), D, S}(params, cache)
    end
end

"Integrate ODE with explicit Runge-Kutta integrator."
function integrate_step!(int::IntegratorERK{DT,TT}, sol::AtomicSolutionODE{DT,TT}) where {DT,TT}
    local tᵢ::TT
    local yᵢ::DT

    # reset cache
    reset!(sol, timestep(int))

    # compute internal stages
    for i in eachstage(int)
        @inbounds for k in eachindex(int.cache.Q[i], int.cache.V[i])
            yᵢ = 0
            for j in 1:i-1
                yᵢ += tableau(int).q.a[i,j] * int.cache.V[j][k]
            end
            int.cache.Q[i][k] = sol.q̅[k] + timestep(int) * yᵢ
        end
        tᵢ = sol.t̅ + timestep(int) * tableau(int).q.c[i]
        equation(int).v(tᵢ, int.cache.Q[i], int.cache.V[i])
    end

    # compute final update
    update_solution!(sol.q, sol.q̃, int.cache.V, tableau(int).q.b, tableau(int).q.b̂, timestep(int))
end


text:
const TASK_NUM = 8
const CHUNK_CHANNEL_SIZE = 2

include("multithreads.jl")
#include("multiprocesses.jl")
include("sequential.jl")
#include("sharedarray.jl")
include("taskthreads.jl")


text:
using Test
using JustSayIt
using JustSayIt.API
using PyCall
import JustSayIt: DEFAULT_MODEL_NAME, TYPE_MODEL_NAME, MODELDIR_PREFIX, DEFAULT_NOISES, COMMAND_RECOGNIZER_ID
import JustSayIt: init_jsi, finalize_jsi, recognizer, noises, reader, start_reading, stop_reading, read_wav, set_default_streamer, reset_all, _are_next

# Test setup
const SAMPLEDIR_TYPE    = joinpath("samples", "type")
const SAMPLEDIR_SILENCE = joinpath("samples", "silence")

commands = Dict("help"      => Help.help,
                "type"      => Keyboard.type)
modeldirs = Dict(DEFAULT_MODEL_NAME => joinpath(MODELDIR_PREFIX, "vosk-model-small-en-us-0.15"),
                 TYPE_MODEL_NAME    => joinpath(MODELDIR_PREFIX, "vosk-model-small-en-us-0.15"))
init_jsi(commands, modeldirs, DEFAULT_NOISES)

and           = read_wav(joinpath(SAMPLEDIR_TYPE, "and.wav"))
a             = read_wav(joinpath(SAMPLEDIR_TYPE, "a.wav"))
colon         = read_wav(joinpath(SAMPLEDIR_TYPE, "colon.wav"))
comma         = read_wav(joinpath(SAMPLEDIR_TYPE, "comma.wav"))
c             = read_wav(joinpath(SAMPLEDIR_TYPE, "c.wav"))
digits        = read_wav(joinpath(SAMPLEDIR_TYPE, "digits.wav"))
dot           = read_wav(joinpath(SAMPLEDIR_TYPE, "dot.wav"))
eight         = read_wav(joinpath(SAMPLEDIR_TYPE, "eight.wav"))
exclamation   = read_wav(joinpath(SAMPLEDIR_TYPE, "exclamation.wav"))
five          = read_wav(joinpath(SAMPLEDIR_TYPE, "five.wav"))
four          = read_wav(joinpath(SAMPLEDIR_TYPE, "four.wav"))
interrogation = read_wav(joinpath(SAMPLEDIR_TYPE, "interrogation.wav"))
i             = read_wav(joinpath(SAMPLEDIR_TYPE, "i.wav"))
julia         = read_wav(joinpath(SAMPLEDIR_TYPE, "julia.wav"))
j             = read_wav(joinpath(SAMPLEDIR_TYPE, "j.wav"))
language      = read_wav(joinpath(SAMPLEDIR_TYPE, "language.wav"))
let_s         = read_wav(joinpath(SAMPLEDIR_TYPE, "let_s.wav"))
letters       = read_wav(joinpath(SAMPLEDIR_TYPE, "letters.wav"))
lowercase     = read_wav(joinpath(SAMPLEDIR_TYPE, "lowercase.wav"))
l             = read_wav(joinpath(SAMPLEDIR_TYPE, "l.wav"))
nine          = read_wav(joinpath(SAMPLEDIR_TYPE, "nine.wav"))
one           = read_wav(joinpath(SAMPLEDIR_TYPE, "one.wav"))
paragraph     = read_wav(joinpath(SAMPLEDIR_TYPE, "paragraph.wav"))
point         = read_wav(joinpath(SAMPLEDIR_TYPE, "point.wav"))
programming   = read_wav(joinpath(SAMPLEDIR_TYPE, "programming.wav"))
redo          = read_wav(joinpath(SAMPLEDIR_TYPE, "redo.wav"))
semicolon     = read_wav(joinpath(SAMPLEDIR_TYPE, "semicolon.wav"))
seven         = read_wav(joinpath(SAMPLEDIR_TYPE, "seven.wav"))
six           = read_wav(joinpath(SAMPLEDIR_TYPE, "six.wav"))
some          = read_wav(joinpath(SAMPLEDIR_TYPE, "some.wav"))
space         = read_wav(joinpath(SAMPLEDIR_TYPE, "space.wav"))
s             = read_wav(joinpath(SAMPLEDIR_TYPE, "s.wav"))
terminus      = read_wav(joinpath(SAMPLEDIR_TYPE, "terminus.wav"))
text          = read_wav(joinpath(SAMPLEDIR_TYPE, "text.wav"))
the           = read_wav(joinpath(SAMPLEDIR_TYPE, "the.wav"))
three         = read_wav(joinpath(SAMPLEDIR_TYPE, "three.wav"))
_try          = read_wav(joinpath(SAMPLEDIR_TYPE, "try.wav"))
two           = read_wav(joinpath(SAMPLEDIR_TYPE, "two.wav"))
type          = read_wav(joinpath(SAMPLEDIR_TYPE, "type.wav"))
undo          = read_wav(joinpath(SAMPLEDIR_TYPE, "undo.wav"))
uppercase     = read_wav(joinpath(SAMPLEDIR_TYPE, "uppercase.wav"))
u             = read_wav(joinpath(SAMPLEDIR_TYPE, "u.wav"))
words         = read_wav(joinpath(SAMPLEDIR_TYPE, "words.wav"))
_zero         = read_wav(joinpath(SAMPLEDIR_TYPE, "zero.wav"))

_05 = read_wav(joinpath(SAMPLEDIR_SILENCE, "silence_501ms.wav"))
_1  = read_wav(joinpath(SAMPLEDIR_SILENCE, "silence_1001ms.wav"))
_2  = read_wav(joinpath(SAMPLEDIR_SILENCE, "silence_2001ms.wav"))
_5  = read_wav(joinpath(SAMPLEDIR_SILENCE, "silence_5001ms.wav"))

tests_words  = Dict("the julia programming language" => [words; _5;   the; julia; programming; language;                                                    _5; terminus; _5],
                    "the julia programming"          => [words; _5;   the; julia; programming; _5; language; _5; undo;                                      _5; terminus; _5],
                    "the julia"                      => [words; _5;   the; julia; _5; undo; _5; undo; _5; redo; _5; redo; _5; redo;                         _5; terminus; _5],
                    "the Julia programming language" => [words; _5;   the; _5; uppercase; _5; julia; _5; uppercase; lowercase; _5; programming; language;   _5; terminus; _5],
                   );

tests_letters = Dict("julia"  => [letters; _5;   j; u; l; i; _05; a;                                              _5; terminus; _5],
                     "cscs"   => [letters; _5;   c; s; c; s;                                                      _5; terminus; _5],
                     "juli a" => [letters; _5;   j; u; l; i; _5; space; _5;  a;                                   _5; terminus; _5],
                     "cs"     => [letters; _5;   c; s; _5; j; u; _5; undo;                                        _5; terminus; _5],
                     "csc"    => [letters; _5;   c; s; _5; c; _5; undo; _5; undo; _5; redo; _5; redo; _5; redo;   _5; terminus; _5],
                    );

tests_digits  = Dict("05261.78 394" => [digits; _5;   _zero; five; two; six; one; _5; dot; _5; seven; eight; _5; space; _5; three; nine; four;   _5; terminus; _5],
                     "05"           => [digits; _5;   _zero; five; _5; six; _5; undo;                                                            _5; terminus; _5],
                     "056"          => [digits; _5;   _zero; five; _5; six; _5; undo; _5; undo; _5; redo; _5; redo; _5; redo;                    _5; terminus; _5],
                    );

tests_text    = Dict("Let's type some words: the Julia programming language."                => [text; _5;   let_s; type; some; words;   _5; colon; _5;               the; _5; uppercase; _5;  julia; _5; uppercase; lowercase; _5; programming; language; _5; point;                _5; terminus; _5],
                     "Let's type some letters: cscs."                                        => [text; _5;   let_s; type; some; letters; _5; colon; _5; letters; _5;  space; c; s; c; s; _5; point;                                                                                _5; terminus; _5],
                     "Let's type some digits: 07.6."                                         => [text; _5;   let_s; type; some; digits;  _5; colon; _5; digits;  _5;  space; _zero; seven; dot; six; _5; point;                                                               _5; terminus; _5],
                     "Let's try some interrogation exclamation point and paragraph:.,:;!?\n" => [text; _5;   let_s; _try; some; interrogation; exclamation; point; and; paragraph; _5; colon;  point; comma; colon; semicolon; exclamation; interrogation; _5; paragraph;            _5; terminus; _5],
                     "Let's try some interrogation exclamation point and paragraph:.,:;!?"   => [text; _5;   let_s; _try; some; interrogation; exclamation; point; and; paragraph; _5; colon;  point; comma; colon; semicolon; exclamation; interrogation; _5; paragraph; _5; undo;  _5; terminus; _5],
                     "Let's try some interrogation exclamation point and paragraph."         => [text; _5;   let_s; _try; some; interrogation; exclamation; point; and; paragraph; _5;         undo; _5; undo; _5; redo; _5; redo; _5; redo; _5; point;                              _5; terminus; _5],
                    );

@testset "$(basename(@__FILE__))" begin
    @testset "1. type words (\"$test\")" for test in keys(tests_words)
        id = test
        start_reading(tests_words[test]; id=id)
        set_default_streamer(reader, id)
        @test test == Keyboard.type(; do_keystrokes=false)
        stop_reading(id=id)
    end;
    @testset "2. type letters (\"$test\")" for test in keys(tests_letters)
        id = test
        start_reading(tests_letters[test]; id=id)
        set_default_streamer(reader, id)
        @test test == Keyboard.type(; do_keystrokes=false)
        stop_reading(id=id)
    end;
    @testset "3. type digits (\"$test\")" for test in keys(tests_digits)
        id = test
        start_reading(tests_digits[test]; id=id)
        set_default_streamer(reader, id)
        @test test == Keyboard.type(; do_keystrokes=false)
        stop_reading(id=id)
    end;
    @testset "4. type text (\"$test\")" for test in keys(tests_text)
        id = test
        start_reading(tests_text[test]; id=id)
        set_default_streamer(reader, id)
        @test test == Keyboard.type(; do_keystrokes=false)
        stop_reading(id=id)
    end;
end;

# Test tear down
finalize_jsi()


text:
"""
    map_expression(e, f)

If `e` is an expression, map `f` over the arguments in `e`.
"""
map_expression(e, f) = e
map_expression(e::Expr, f) = Expr(f(e.head), map(f, e.args)...)

export with
"""
    @with(e, df = \_)

Extract any symbols `e`, such as `:colA`, from `df`, e.g. `df[:colA]`.

Anything wrapped in `^`, such as `^(expr)`, gets passed through untouched.

### Examples
```julia
a = 1
_ = DataFrames.DataFrame(a = 2)

@test DataFrames.DataFrame(b = _[:a] + a, c = :d) ==
   @with DataFrames.DataFrame(b = :a + a, c = ^(:d))
```
"""
with(e, df = :_) = MacroTools.@match e begin
    ^(e_) => e
    :(e_) => Expr(:ref, df, Meta.quot(e) )
    a_.b_ => Expr(:., a, Meta.quot(b) )
    e_ => map_expression(e, e -> with(e, df) )
end
export @with
@nonstandard with

"""
    double_match(e, first, second)

Test whether `e` calls `first`, and the first argument calls `second`

#Examples
```julia
e = Expr(:parameters, Expr(:..., :a) )
first = :parameters
second = :...

@test ChainMap.double_match(e, first, second)
@test !ChainMap.double_match(:b, :parameters, :...)
```
"""
double_match(e, first, second) =
    MacroTools.isexpr(e, first) &&
    length(e.args) > 0 &&
    MacroTools.isexpr(e.args[1], second)

"""
    replace_key(e, symbol = gensym())

Replace `e` with `gensym()`.

`e` can be wrapped in `...` or `...` and `;`.

#Examples
```julia
e = Expr(:parameters, Expr(:..., :a) )
@test ChainMap.replace_key(e, :b) == Expr(:parameters, Expr(:..., :b) )

e = Expr(:..., :a)
@test ChainMap.replace_key(e, :b) == Expr(:..., :b)

@test ChainMap.replace_key(:a, :b) == :b
```
"""
function replace_key(e, symbol = gensym())
    if double_match(e, :parameters, :...)
        @chain symbol Expr(:..., _) Expr(:parameters, _)
    elseif MacroTools.isexpr(e, :...)
        Expr(:..., symbol)
    else
        symbol
    end
end

"""
    unparameterize(e)

Remove `;` from `e`.

#Examples
```julia
a = Expr(:..., :a)

@test ChainMap.unparameterize(Expr(:parameters, a) ) == a

@test ChainMap.unparameterize(:b) == :b
```
"""
unparameterize(e) =
     if MacroTools.isexpr(e, :parameters) && length(e.args) == 1
         e.args[1]
     else
         e
     end

"""
    add_key!(d, e, symbol = gensym() )

'replace_key' of `e` and add `e` mapped to `symbol` to `d`. Return `e`
unparameterized.

# Examples
```julia
e = :a
symbol = :z
d = Dict(:a => :b, :(a + 1) => :(b + 1))

ChainMap.add_key!(d, e, symbol)
@test d[:a] == :b

ChainMap.add_key!(d, :c, symbol)
@test d[:c] == :z

e = Expr(:parameters, Expr(:..., :d) )
@test ChainMap.add_key!(d, e, symbol) == Expr(:..., :z)
@test d[e] == Expr(:parameters, Expr(:..., :z) )
```
"""
function add_key!(d, e, symbol = gensym() )
    if @chain d haskey(_, e) !
        d[e] = replace_key(e, symbol)
    end
    unparameterize( d[e] )
end

"""
# Examples
```julia
d = Dict()
e = :( 1 + ~(a))
ChainMap.replace_record!(e, d)
@test :a in keys(d)
```
"""
replace_record!(e, d) =
    MacroTools.@match e begin
        ~(e_) => add_key!(d, e)
        e_ => map_expression(e, e -> replace_record!(e, d) )
    end


text:
__precompile__(true)

module Quantum
using FEDVR

include("potentials.jl")
include("hamiltonians.jl")
include("eigenstates.jl")

end # module


text:
# Tabix Overlap Iterator
# ======================

struct TabixOverlapIterator{T}
    reader::T
    interval::Interval
end

function Base.eltype(::Type{TabixOverlapIterator{T}}) where T
    return eltype(T)
end

function Base.IteratorSize(::Type{TabixOverlapIterator{T}}) where T
    return Base.SizeUnknown()
end

mutable struct TabixOverlapIteratorState{T}
    chunks::Vector{Indexes.Chunk}
    chunkid::Int
    done::Bool
    record::T
end

function Base.iterate(iter::TabixOverlapIterator)
    @assert iter.reader.index !== nothing
    # TODO: Use a method that resets the reading position.
    buffer = BioGenerics.IO.stream(iter.reader)
    iter.reader.state = BioGenerics.Automa.State(TranscodingStreams.NoopStream(buffer.stream), 1, 1, false)
    state = TabixOverlapIteratorState(Indexes.overlapchunks(iter.reader.index, iter.interval), 0, false, eltype(iter)())

    return iterate(iter, state)
end

function done(iter::TabixOverlapIterator, state)
    buffer = BioGenerics.IO.stream(iter.reader)
    source = buffer.stream
    if state.chunkid == 0
        if isempty(state.chunks)
            return true
        end
        state.chunkid += 1
        seek(source, state.chunks[state.chunkid].start)
    end
    while state.chunkid ≤ lastindex(state.chunks)
        chunk = state.chunks[state.chunkid]
        #=
        The `virtualoffset(source)` is not synchronized with the current reading position because data are buffered in `buffer` for parsing text.
        So we need to check not only `virtualoffset` but also `nb_available`, which returns the current buffered data size.
        =#
        while bytesavailable(buffer) > 0 || CodecBGZF.virtualoffset(source) < chunk.stop
            read!(iter.reader, state.record)
            c = icmp(state.record, iter.interval)
            if c == 0  # overlapping
                return false
            end
            if c > 0
                # no more overlapping records in this chunk
                break
            end
        end
        state.chunkid += 1
        if state.chunkid ≤ lastindex(state.chunks)
            seek(source, state.chunks[state.chunkid].start)
        end
    end
    # no more overlapping records
    return true
end

function Base.iterate(iter::TabixOverlapIterator, state)
    if done(iter, state)
        return nothing
    end

    return copy(state.record), state
end

function icmp(record, interval)
    c = cmp(BioGenerics.seqname(record), BioGenerics.seqname(interval))

    if c < 0 || (c == 0 && BioGenerics.rightposition(record) < BioGenerics.leftposition(interval))
        return -1
    end

    if c > 0 || (c == 0 && BioGenerics.leftposition(record) > BioGenerics.rightposition(interval))
        return +1
    end

    return 0
end


text:
using PyPlot
# Tests different ways of keeping track of time.
# Looks like multiplication is just as accurate as
# compensated summation, yet is much faster.  Not
# that this is a limiting factor.
function test_time(nt)
#nt = 10000000
T = Float64
s2 = zero(T)
h = 0.04
hbig = big(h)
t = zeros(T,nt)
tb = zeros(BigFloat,nt)
tc = zeros(T,nt)
tm = zeros(T,nt)
for i=2:nt
  tc[i],s2 = comp_sum(tc[i-1],s2,h)
  t[i] = t[i-1] + h
  tm[i] = (i-1)*h
  tb[i] = tb[i-1] + hbig
end
# Add compensated error at the end:
println("Compensated: ",tc[nt]," big: ",tb[nt]," diff: ",tc[nt]-tb[nt]," diff+err: ",big(tc[nt])+big(s2)-tb[nt])

#tref = convert(Array{Float64,1},tb)
#semilogy(abs.(t-tref),label="Added")
#println("Maximum difference added: ",maximum(abs.(t-tref)))
#plot(abs.(tc-tref),":",label="Compensated")
#println("Maximum difference compensated: ",maximum(abs.(tc-tref)))
#plot(abs.(tm-tref),"--",label="Multiplied")
#println("Maximum difference multiplied: ",maximum(abs.(tm-tref)))
#legend()
return
end


text:
module AbstractActuators

using AbstractDAQs

export devname
export AbstractActuator, AbstractRobot, AbstractCartesianRobot, numaxes
export moveto
export move, moveX, moveY, moveZ
export rmove, rmoveX, rmoveY, rmoveZ
export axesnames

export absposition, abspositionX, abspositionY, abspositionZ
export devposition, positionX, positionY, positionZ

export setreference, setreferenceX, setreferenceY, setreferenceZ
export setabsreference, setabsreferenceX, setabsreferenceY, setabsreferenceZ

export home, homeX, homeY, homeZ
export stopmotion, waituntildone

export TestRobot1d, TestRobot

"Abstract base type for actuators"
abstract type AbstractActuator <: AbstractDAQs.AbstractDevice end

"Abstract base type for Robots"
abstract type AbstractRobot <: AbstractActuator end

"Abstract base type for cartesian robots"
abstract type AbstractCartesianRobot <: AbstractRobot end

include("move_interface.jl")
include("actuatorset.jl")
include("experimentmatrix.jl")
include("manualactuator.jl")
include("testdevices.jl")
include("hdf5io.jl")
end


text:
abstract type SLMixingRule end

function init_model(model::SLMixingRule,components,userlocations,verbose)
    return model
end

function init_model(model::Type{<:SLMixingRule},components,userlocations,verbose)
    verbose && @info("""Now creating Sanchez-Lacombe mixing rule:
    $model""")
    return model(components;userlocations,verbose)
end


text:
# A Collection of static analysis functions.
#
# These functions implement some static type-checking of Julia programs.
# There are also other type-related functions to allow human-readable output
# of things like the types of function-local variables.
module TypeCheck
export checkreturntypes, checklooptypes, checkmethodcalls,
  methodswithdescendants

## Modifying functions from Base

# Base.code_typed requires a function and tuple of argument types.
# It can be more convenient to pass in a Function or Method.
# A method specifies the same information as a Function + argument types,
# so the implementation is obvious.
# For a Function, code_typed should return a list of Expr, on per method.

# Return the type-inferred AST for each method of a generic Function.
function Base.code_typed(f::Function)
  Expr[code_typed(m) for m in f.env]
end

# Return the type-inferred AST for one Method (of a generic Function)
function Base.code_typed(m::Method)
 linfo = m.func.code
 (tree,ty) = Base.typeinf(linfo,m.sig,())
 if !isa(tree,Expr)
     ccall(:jl_uncompress_ast, Any, (Any,Any), linfo, tree)
 else
    tree
 end
end

# The whos function takes a namespace (a Module or the current one) and
# prints out the bound names and types (variables, types, functions, consts).
# I'm adding a version to print out function-local variables.
# Specify a method via a Method, Function, or Function+argument-type-tuple.
function Base.whos(f,args...)
  for e in code_typed(f,args...)
    display(MethodSignature(e))
    for x in sort(e.args[2][2];by=x->x[1])
      println("\t",x[1],"\t",x[2])
    end
    println("")
  end
end

## Basic Operations on Function Exprs

# Return the inferred return of type of an Expr representing a Method
returntype(e::Expr) =  e.args[3].typ

# Return the body of a Method.
# Takes an Expr representing a Method, returns Vector{Expr}.
body(e::Expr) = e.args[3].args

# Return all of the return statements in the body of a Method
# Takes an Expr representing a Method, return Vector{Expr}
returns(e::Expr) = filter(x-> typeof(x) == Expr && x.head==:return,body(e))

# Returns all function calls contained in return statements in the body of a Method.
# Takes an Expr representing a Method, returns Vector{Expr}
function extractcallsfromreturns(e::Expr)
  rs = returns(e)
  rs_with_calls = filter(x->typeof(x.args[1]) == Expr && x.args[1].head == :call,rs)
  Expr[expr.args[1] for expr in rs_with_calls]
end

# A type that covers both Types and TypeVars
AType = Union(Type,TypeVar)

# Returns the type of each argument in a Method's signature
# Takes an Expr representing a Method, returns Vector{AType}
function argumenttypes(e::Expr)
  argnames = Symbol[typeof(x) == Symbol ? x : x.args[1] for x in e.args[1]]
  argtuples = filter(x->x[1] in argnames, e.args[2][2]) #only arguments, no local vars
  AType[t[2] for t in argtuples]
end

# Returns true if the Expr is a :call to a TopNode
# Takes an Expr (or Any) and returns a Bool
istop(e) = Base.is_expr(e,:call) && typeof(e.args[1]) == TopNode

# Returns the return type of an Expr (:call, :call1, :new)
# Takes an Expr(call) and the Expr representing the method it occurs in
# Returns an AType
function returntype(e::Expr,context::Expr) #must be :call,:new,:call1
  if Base.is_expr(e,:new); return e.typ; end
  if Base.is_expr(e,:call1) && isa(e.args[1], TopNode); return e.typ; end
  if !Base.is_expr(e,:call); error("Expected :call Expr"); end

  if istop(e)
    return e.typ
  end

  callee = e.args[1]
  if istop(callee)
    return returntype(callee,context)
  elseif isa(callee,SymbolNode) # only seen (func::F), so non-generic function
    return Any
  elseif is(callee,Symbol)
    if e.typ != Any || any([isa(x,LambdaStaticData) for x in e.args[2:end]])
      return e.typ
    end

    if isdefined(Base,callee)
      f = eval(Base,callee)
      if !isa(f,Function) || !isgeneric(f)
        return e.typ
      end
      fargtypes = tuple([argumenttype(ea,context) for ea in e.args[2:end]])
      return Union([returntype(ef) for ef in code_typed(f,fargtypes)]...)
    else
      return @show e.typ
    end
  end

  return e.typ
end

# Returns the type of an Expr
# Takes an Expr and the Expr representing a Method it occurs in
function argumenttype(e::Expr,context::Expr)
 if Base.is_expr(e,:call) || Base.is_expr(e,:new) || Base.is_expr(e,:call1)
   return returntype(e,context)
 end

 @show e
 return Any
end

# Returns the type of a Symbol used as an argument to a function call
# Takes a Symbol and an Expr representing a Method it occurs in
function argumenttype(s::Symbol,e::Expr)
  vartypes = [x[1] => x[2] for x in e.args[2][2]]
  s in vartypes ? (vartypes[@show s]) : Any
end

# as above, but for different call argument types
argumenttype(s::SymbolNode,e::Expr) = s.typ
argumenttype(t::TopNode,e::Expr) = Any
argumenttype(l::LambdaStaticData,e::Expr) = Function
argumenttype(q::QuoteNode,e::Expr) = argumenttype(q.value,e)

# as above, but for various literal values
argumenttype(n::Number,e::Expr) = typeof(n)
argumenttype(c::Char,e::Expr) = typeof(c)
argumenttype(s::String,e::Expr) = typeof(s)
argumenttype(i,e::Expr) = typeof(i) #catch all, hopefully for more literals

# Make Types Iterable
# start, next, and done are the functions for-loops use to iterate.
# This defines iterating over a DataType to mean iterating over its descendants
# (breadth-first-search ordering)
Base.start(t::DataType) = [t]
function Base.next(t::DataType,arr::Vector{DataType})
  c = pop!(arr)
  append!(arr,[x for x in subtypes(c)])
  (c,arr)
end
Base.done(t::DataType,arr::Vector{DataType}) = length(arr) == 0

# Returns a list of Function names to implementation percentage.
# By considering all the descendents of the provided type,
# this function collects a list of functions implemented by at least one
# considered type and calculates what percentage of the considered types
# implements each function.
#
# Takes a DataType and returns a Vector{(Symbol,Float64)}
# Keyword args:
#   onlyleaves: If this is true, only concrete types are considered;
#     if this is false, all descendent types are considered.
#     (defaults to false)
#   lim: The maximum length of the returned list. The lim functions with
#     the highest implementation rate will be returned.
function methodswithdescendants(t::DataType;onlyleaves::Bool=false,lim::Int=10)
  d = Dict{Symbol,Int}()
  count = 0
  for s in t
    if !onlyleaves || (onlyleaves && isleaftype(s))
      count += 1
      fs = Set{Symbol}()
      for m in methodswith(s)
        push!(fs,m.func.code.name)
      end
      for sym in fs
        d[sym] = get(d,sym,0) + 1
      end
    end
  end
  l = [(k,v/count) for (k,v) in d]
  sort!(l,by=(x->x[2]),rev=true)
  l[1:min(lim,end)]
end

# Run a function on each generic function of a Module.
# Takes a Module, returns nothing.
# The output of each run of the test function will be displayed.
function checkallmodule(m::Module;test=checkreturntypes,kwargs...)
  score = 0
  for n in names(m)
    f = eval(m,n)
    if isgeneric(f) && typeof(f) == Function
      fm = test(f;mod=m,kwargs...)
      score += length(fm.methods)
      display(fm)
    end
  end
  println("The total number of failed methods in $m is $score")
end

# Use checkallmodule to implement the Module version of other checks
checkreturntypes(m::Module;kwargs...) = checkallmodule(m;test=checkreturntypes,kwargs...)
checklooptypes(m::Module) = checkallmodule(m;test=checklooptypes)
checkmethodcalls(m::Module) = checkallmodule(m;test=checkmethodcalls)

## Checking that return values are base only on input *types*, not values.

type MethodSignature
  typs::Vector{AType}
  returntype::Union(Type,TypeVar) # v0.2 has TypeVars as returntypes; v0.3 does not
end
MethodSignature(e::Expr) = MethodSignature(argumenttypes(e),returntype(e))
function Base.writemime(io, ::MIME"text/plain", x::MethodSignature)
  println(io,"(",join([string(t) for t in x.typs],","),")::",x.returntype)
end

type FunctionSignature
  methods::Vector{MethodSignature}
  name::Symbol
end

function Base.writemime(io, ::MIME"text/plain", x::FunctionSignature)
  for m in x.methods
    print(io,string(x.name))
    display(m)
  end
end

# given a function, run checkreturntypes on each method
function checkreturntypes(f::Function;kwargs...)
  results = MethodSignature[]
  for e in code_typed(f)
    (ms,b) = checkreturntype(e;kwargs...)
    if b push!(results,ms) end
  end
  FunctionSignature(results,f.env.name)
end

# given an Expr representing a Method,
# determine whether its return type is based
# only on the arugment types or whether it is
# also influenced by argument values
# (the Method fails the check if the return type depends on values)
function checkreturntype(e::Expr;kwargs...)
  (typ,b) = isreturnbasedonvalues(e;kwargs...)
  (MethodSignature(argumenttypes(e),typ),b)
end

# Determine whether this method's return type might change based on input values rather than input types
function isreturnbasedonvalues(e::Expr;mod=Base)
  rt = returntype(e)
  ts = argumenttypes(e)
  if isleaftype(rt) || rt == None return (rt,false) end

  for t in ts
    if !isleaftype(t)
      return (rt,false)
    end
  end

  cs = [returntype(c,e) for c in extractcallsfromreturns(e)]
  for c in cs
    if rt == c
       return (rt,false)
    end
  end

  return (rt,true) # return is not concrete type; all args are concrete types
end

## Checking that variables in loops have concrete types

type LoopResult
  msig::MethodSignature
  lines::Vector{(Symbol,Type)} #TODO should this be a specialized type? SymbolNode?
  LoopResult(ms::MethodSignature,ls::Vector{(Symbol,Type)}) = new(ms,unique(ls))
end

function Base.writemime(io, ::MIME"text/plain", x::LoopResult)
  display(x.msig)
  for (s,t) in x.lines
    println(io,"\t",string(s),"::",string(t))
  end
end

type LoopResults
  name::Symbol
  methods::Vector{LoopResult}
end

function Base.writemime(io, ::MIME"text/plain", x::LoopResults)
  for lr in x.methods
    print(io,string(x.name))
    display(lr)
  end
end

# for a given Function, run checklooptypes on each Method
function checklooptypes(f::Function;kwargs...)
  lrs = LoopResult[]
  for e in code_typed(f)
    lr = checklooptypes(e)
    if length(lr.lines) > 0 push!(lrs,lr) end
  end
  LoopResults(f.env.name,lrs)
end

# for an Expr representing a Method,
# check that the type of each variable used in a loop
# has a concrete type
checklooptypes(e::Expr;kwargs...) = LoopResult(MethodSignature(e),loosetypes(loopcontents(e)))

# This is a function for trying to detect loops in the body of a Method
# Returns lines that are inside one or more loops
function loopcontents(e::Expr)
  b = body(e)
  loops = Int[]
  nesting = 0
  lines = {}
  for i in 1:length(b)
    if typeof(b[i]) == LabelNode
      l = b[i].label
      jumpback = findnext(
        x-> (typeof(x) == GotoNode && x.label == l) || (Base.is_expr(x,:gotoifnot) && x.args[end] == l),
        b, i)
      if jumpback != 0
        push!(loops,jumpback)
        nesting += 1
      end
    end
    if nesting > 0
      push!(lines,(i,b[i]))
    end

    if typeof(b[i]) == GotoNode && in(i,loops)
      splice!(loops,findfirst(loops,i))
      nesting -= 1
    end
  end
  lines
end

# given `lr`, a Vector of expressions (Expr + literals, etc)
# try to find all occurances of a variables in `lr`
# and determine their types
# `method` is only used as part of constructing the return type
function loosetypes(lr::Vector)
  lines = (Symbol,Type)[]
  for (i,e) in lr
    if typeof(e) == Expr
      es = copy(e.args)
      while !isempty(es)
        e1 = pop!(es)
        if typeof(e1) == Expr
          append!(es,e1.args)
        elseif typeof(e1) == SymbolNode && !isleaftype(e1.typ) && typeof(e1.typ) == UnionType
          push!(lines,(e1.name,e1.typ))
        end
      end
    end
  end
  return lines
end

## Check method calls

type CallSignature
  name::Symbol
  argumenttypes::Vector{AType}
end
function Base.writemime(io, ::MIME"text/plain", x::CallSignature)
  println(io,string(x.name),"(",join([string(t) for t in x.argumenttypes],","),")")
end

type MethodCalls
  m::MethodSignature
  calls::Vector{CallSignature}
end

function Base.writemime(io, ::MIME"text/plain", x::MethodCalls)
  display(x.m)
  for c in x.calls
    print(io,"\t")
    display(c)
  end
end

type FunctionCalls
  name::Symbol
  methods::Vector{MethodCalls}
end

function Base.writemime(io, ::MIME"text/plain", x::FunctionCalls)
  for mc in x.methods
    print(io,string(x.name))
    display(mc)
  end
end

# given a Function, run `checkmethodcalls` on each Method
# and collect the results into a FunctionCalls
function checkmethodcalls(f::Function;kwargs...)
  calls = MethodCalls[]
  for m in f.env
    e = code_typed(m)
    mc = checkmethodcalls(e,m;kwargs...)
    if !isempty(mc.calls)
      push!(calls, mc)
    end
  end
  FunctionCalls(f.env.name,calls)
end

# given an Expr representing a Method,
# and the Method it represents,
# check the Method body for calls to non-existant Methods
function checkmethodcalls(e::Expr,m::Method;kwargs...)
  if Base.arg_decl_parts(m)[3] == symbol("deprecated.jl")
    CallSignature[]
  end
  nomethoderrors(e,methodcalls(e);kwargs...)
end

# Find any methods that match the given CallSignature
function hasmatches(mod::Module,cs::CallSignature)
  if isdefined(mod,cs.name)
    f = eval(mod,cs.name)
    if isgeneric(f)
      opts = methods(f,tuple(cs.argumenttypes...))
      if isempty(opts)
        return false
      end
    end
  else
    #println("$mod.$(cs.name) is undefined")
  end
  return true
end

# Find any CallSignatures that indicate potential NoMethodErrors
function nomethoderrors(e::Expr,cs::Vector{CallSignature};mod=Base)
  output = CallSignature[]
  for callsig in cs
    if !hasmatches(mod,callsig)
      push!(output,callsig)
    end
  end
  MethodCalls(MethodSignature(e),output)
end

# Look through the body of the function for `:call`s
function methodcalls(e::Expr)
  b = body(e)
  lines = CallSignature[]
  for s in b
    if typeof(s) == Expr
      if s.head == :return
        append!(b, s.args)
      elseif s.head == :call
        if typeof(s.args[1]) == Symbol
          push!(lines,CallSignature(s.args[1], [argumenttype(e1,e) for e1 in s.args[2:end]]))
        end
      end
    end
  end
  lines
end

# Return a list of all variables used on the left-hand-side of assignment (=)
#
# Arguments:
#   e: an Expr representing a Method, as from code_typed
#
# Returns:
#   a Set{Symbol}, where each element appears on the lhs of an assignment in e.
#
function find_lhs_variables(e::Expr)
  output = Set{Symbol}()
  for ex in body(e)
   isa(ex,Expr) && ex.head == symbol("=") && push!(output,ex.args[1])
  end
  return output
end

# Given an Expression, finds variables used in it (on right-hand-side)
#
# Arguments: e: an Expr
#
# Returns: a Set{Symbol}, where each e is used in a rhs expression in e
#
function find_rhs_variables(e::Expr)
  output = Set{Symbol}()

  if e.head == :lambda
    for ex in body(e)
      union!(output,find_rhs_variables(ex))
    end
  elseif e.head == :(=)
    for ex in e.args[2:end]  # skip lhs
      union!(output,find_rhs_variables(ex))
    end
  elseif e.head == :return
    output = find_rhs_variables(e.args[1])
  elseif e.head == :call
    start = 2  # skip function name
    e.args[1] == TopNode(:box) && (start = 3)  # skip type name
    for ex in e.args[start:end]
      union!(output,find_rhs_variables(ex))
    end
  elseif e.head == :if
   for ex in e.args # want to check condition, too
     union!(output,find_rhs_variables(ex))
   end
  elseif e.head == :(::)
    output = find_rhs_variables(e.args[1])
  end

  return output
end

# Recursive Base Cases, to simplify control flow in the Expr version
find_rhs_variables(a) = Set{Symbol}()  # unhandled, should be an immediate value, like an Int.
find_rhs_variables(s::Symbol) = Set{Symbol}([s])
find_rhs_variables(s::SymbolNode) = Set{Symbol}([s.name])

function unused_locals(e::Expr)
  lhs = find_lhs_variables(e)
  rhs = find_rhs_variables(e)
  setdiff(lhs,rhs)
end

check_locals(f::Function) = all([check_locals(e) for e in code_typed(f)])
check_locals(e::Expr) = isempty(unused_locals(e))

end  #end module


text:
function plot_vf(velocities)
    """
        velocities is supposed to be of shape nxnx2
    """
    n = size(velocities,1)
    u = velocities[:,:,1]
    v = velocities[:,:,2]
    x,y = meshgrid(0:n-1, 0:n-1)
    quiver(x,y,u,v,color="b")
    xlim(0,n)
    ylim(0,n)
    axis("off")
    show()
end


text:
"""
Blending is blah
"""

module Blending

# package code goes here

"""
sdcsdcsdc
"""
include("Direct.jl")
include("Indirect.jl")

end # module


text:
type abstractMap
	range #the numercal value of the subState
	sub #array of max subState index
	nRange #size of the range
	nSize #how many total states
end

function makeMap(range, nSeg)
	n=length(range)
	return abstractMap(range, int64(ones(nSeg)*n), n, n^nSeg)
end

function randEval(space::abstractMap) #returns some random evaluated point in the space (useful for getting random state)
	return ind2eval(rand(1:space.nSize), space)
end

function ind2eval(ind, space::abstractMap) #maps space index to value
	#function body
	sub=ind2sub(tuple(space.sub...), ind)
	return [space.range[sub[i]]	for i=1:length(sub)]
end

function eval2ind(value, space::abstractMap) #maps space value to index
	sub=[find(space.range .== value[i])[] for i=1:length(value)]
	#return int64(sub)
	return sub2ind(space.sub, int64(sub))
end

function detTrans(stateInd, actionInd, stateMap::abstractMap, actionMap::abstractMap)
	state=ind2eval(stateInd, stateMap)
	action=ind2eval(actionInd, actionMap)

	##purly determinmistic transision
	state_=state+action

	#check bound
	sMax=maximum(stateMap.range)
	sMin=minimum(stateMap.range)
	for i=1:length(state_)
		if state_[i]>sMax
			state_[i]=sMax
		elseif state_[i]<sMin
			state_[i]=sMin
		end
	end
	#return new index
	return eval2ind(state_, stateMap)
end

println("Import successful: floopMap.jl")


text:
"""
Function that solves the first unit commitment optimization
"""
function bucr1_model(day::Int64, hour::Int64, df_gens::DataFrame,
    df_peakers::DataFrame, df_storage::DataFrame, fuelprice::Matrix{Float64},
    fuelprice_peakers::Matrix{Float64}, map_gens::Matrix{Int64},
	map_peakers::Matrix{Int64}, map_storage::Matrix{Int64},
	bucr1_results::UC_Results, sucr_results::UC_Results,
	bucr2_results::UC_Results, fucr_to_bucr1::UC_Results, bucr1_prepoc::DemandHrPreprocGens,
	bucr_commit_LB::Vector{Int64}, bucr_commit_UB::Vector{Int64},
	bucr_commit_peaker_LB::Vector{Int64}, bucr_commit_peaker_UB::Vector{Int64})

	t1_BUCR1model = time_ns()
    BUCR1model = direct_model(CPLEX.Optimizer())
    set_optimizer_attribute(BUCR1model, "CPX_PARAM_EPGAP", SOLVER_EPGAP)

    # Declaring the decision variables for conventional generators
    @variable(BUCR1model, bucrm_genOnOff[1:GENS], Bin) #Bin
    @variable(BUCR1model, bucrm_genStartUp[1:GENS], Bin) # startup variable
    @variable(BUCR1model, bucrm_genShutDown[1:GENS], Bin) # shutdown variable
    @variable(BUCR1model, bucrm_genOut[1:GENS]>=0) # Generator's output schedule
    @variable(BUCR1model, bucrm_genOut_Block[1:GENS, 1:BLOCKS]>=0) # Generator's output schedule from each block (block rfers to IHR curve blocks/segments)
    @variable(BUCR1model, bucrm_TotGenVioP[g=1:GENS]>=0)
    @variable(BUCR1model, bucrm_TotGenVioN[g=1:GENS]>=0)
    @variable(BUCR1model, bucrm_MaxGenVioP[g=1:GENS]>=0)
    @variable(BUCR1model, bucrm_MinGenVioP[g=1:GENS]>=0)
    #@variable(BUCR1model, bucrm_genResUp[1:GENS]>=0) # Generators' up reserve schedule
    #@variable(BUCR1model, bucrm_genResDn[1:GENS]>=0) # Generator's down rserve schedule

    # Declaring the decision variables for peakers
    @variable(BUCR1model, 0<=bucrm_peakerOnOff[1:PEAKERS]<=1) #Bin, relaxed
    @variable(BUCR1model, 0<=bucrm_peakerStartUp[1:PEAKERS]<=1) # startup variable
    @variable(BUCR1model, 0<=bucrm_peakerShutDown[1:PEAKERS]<=1) # shutdown variable
    @variable(BUCR1model, bucrm_peakerOut[1:PEAKERS]>=0) # Generator's output schedule
    @variable(BUCR1model, bucrm_peakerOut_Block[1:PEAKERS, 1:BLOCKS]>=0) # Generator's output schedule from each block (block rfers to IHR curve blocks/segments)

    # declaring decision variables for storage Units
    @variable(BUCR1model, bucrm_storgChrg[1:STORG_UNITS], Bin) #Bin variable equal to 1 if unit runs in the charging mode
    @variable(BUCR1model, bucrm_storgDisc[1:STORG_UNITS], Bin) #Bin variable equal to 1 if unit runs in the discharging mode
    @variable(BUCR1model, bucrm_storgIdle[1:STORG_UNITS], Bin) ##Bin variable equal to 1 if unit runs in the idle mode
    @variable(BUCR1model, bucrm_storgChrgPwr[1:STORG_UNITS]>=0) #Chargung power
    @variable(BUCR1model, bucrm_storgDiscPwr[1:STORG_UNITS]>=0) # Discharging Power
    @variable(BUCR1model, bucrm_storgSOC[1:STORG_UNITS]>=0) # state of charge (stored energy level for storage unit at time t)
    #@variable(BUCR1model, bucrm_storgResUp[1:STORG_UNITS]>=0) # Scheduled up reserve
    #@variable(BUCR1model, bucrm_storgResDn[1:STORG_UNITS]>=0) # Scheduled down reserve

    # declaring decision variables for renewable generation
    @variable(BUCR1model, bucrm_solarG[1:N_ZONES]>=0) # solar energy schedules
    @variable(BUCR1model, bucrm_windG[1:N_ZONES]>=0) # wind energy schedules
    @variable(BUCR1model, bucrm_hydroG[1:N_ZONES]>=0) # hydro energy schedules
    @variable(BUCR1model, bucrm_solarGSpil[1:N_ZONES]>=0) # solar energy schedules
    @variable(BUCR1model, bucrm_windGSpil[1:N_ZONES]>=0) # wind energy schedules
    @variable(BUCR1model, bucrm_hydroGSpil[1:N_ZONES]>=0) # hydro energy schedules

    # Declaring decision variables for hourly dispatched and curtailed demand
    @variable(BUCR1model, bucrm_Demand[1:N_ZONES]>=0) # Hourly scheduled demand
    @variable(BUCR1model, bucrm_Demand_Curt[1:N_ZONES]>=0) # Hourly schedule demand
    @variable(BUCR1model, bucrm_OverGen[1:N_ZONES]>=0) #

    # declaring variables for transmission system
    @variable(BUCR1model, bucrm_voltAngle[1:N_ZONES]) #voltage angle at zone/bus n in t//
    @variable(BUCR1model, bucrm_powerFlow[1:N_ZONES, 1:M_ZONES]) #transmission Flow from zone n to zone m//

    # Defining the objective function that minimizes the overal cost of supplying electricity (=total variable, no-load, startup, and shutdown costs)

    #@objective(BUCR1model, Min, sum(sum(df_gens.VariableCost[g]*bucrm_genOut[g]+df_gens.NoLoadCost[g]*bucrm_genOnOff[g] +df_gens.StartUpCost[g]*bucrm_genStartUp[g] + df_gens.ShutdownCost[g]*bucrm_genShutDown[g] for g in 1:GENS)))
    @objective(BUCR1model, Min,
		sum(df_gens.IHRC_B1_HR[g]*fuelprice[g,day]*(bucrm_genOut_Block[g,1])
	   +df_gens.IHRC_B2_HR[g]*fuelprice[g,day]*(bucrm_genOut_Block[g,2])
	   +df_gens.IHRC_B3_HR[g]*fuelprice[g,day]*(bucrm_genOut_Block[g,3])
	   +df_gens.IHRC_B4_HR[g]*fuelprice[g,day]*(bucrm_genOut_Block[g,4])
	   +df_gens.IHRC_B5_HR[g]*fuelprice[g,day]*(bucrm_genOut_Block[g,5])
	   +df_gens.IHRC_B6_HR[g]*fuelprice[g,day]*(bucrm_genOut_Block[g,6])
	   +df_gens.IHRC_B7_HR[g]*fuelprice[g,day]*(bucrm_genOut_Block[g,7])
	   +df_gens.NoLoadHR[g]*fuelprice[g,day]*(bucrm_genOnOff[g])
	   +((df_gens.HotStartU_FixedCost[g]+(df_gens.HotStartU_HeatRate[g]*fuelprice[g,day]))*(bucrm_genStartUp[g]))
	   +df_gens.ShutdownCost[g]*(bucrm_genShutDown[g])
	   +(bucrm_TotGenVioP[g]*VIOLATION_PENALTY)+(bucrm_TotGenVioN[g]*VIOLATION_PENALTY)+(bucrm_MaxGenVioP[g]*VIOLATION_PENALTY)+(bucrm_MinGenVioP[g]*VIOLATION_PENALTY) for g in 1:GENS)
	   +sum(df_peakers.IHRC_B1_HR[k]*fuelprice_peakers[k,day]*(bucrm_peakerOut_Block[k,1])
	   +df_peakers.IHRC_B2_HR[k]*fuelprice_peakers[k,day]*(bucrm_peakerOut_Block[k,2])
	   +df_peakers.IHRC_B3_HR[k]*fuelprice_peakers[k,day]*(bucrm_peakerOut_Block[k,3])
	   +df_peakers.IHRC_B4_HR[k]*fuelprice_peakers[k,day]*(bucrm_peakerOut_Block[k,4])
	   +df_peakers.IHRC_B5_HR[k]*fuelprice_peakers[k,day]*(bucrm_peakerOut_Block[k,5])
	   +df_peakers.IHRC_B6_HR[k]*fuelprice_peakers[k,day]*(bucrm_peakerOut_Block[k,6])
	   +df_peakers.IHRC_B7_HR[k]*fuelprice_peakers[k,day]*(bucrm_peakerOut_Block[k,7])
	   +df_peakers.NoLoadHR[k]*fuelprice_peakers[k,day]*(bucrm_peakerOnOff[k])
	   +((df_peakers.HotStartU_FixedCost[k]+(df_peakers.HotStartU_HeatRate[k]*fuelprice_peakers[k,day]))*(bucrm_peakerStartUp[k]))
	   +df_peakers.ShutdownCost[k]*(bucrm_peakerShutDown[k]) for k in 1:PEAKERS)
	   +sum((bucrm_Demand_Curt[n]*LOAD_SHED_PENALTY)+(bucrm_OverGen[n]*OVERGEN_PENALTY) for n=1:N_ZONES) )

    # Baseload Operation of nuclear units
    @constraint(BUCR1model, conNuckBaseLoad[g=1:GENS], bucrm_genOnOff[g]>=df_gens.Nuclear[g]) #
    @constraint(BUCR1model, conNuclearTotGenZone[n=1:N_ZONES],
    	sum((bucrm_genOut[g]*map_gens[g,n]*df_gens.Nuclear[g]) for g=1:GENS) - bucr1_prepoc.nuclear[n] ==0)

    #Limits on generation of cogen units
    @constraint(BUCR1model, conCoGenBaseLoad[g=1:GENS], bucrm_genOnOff[g]>=df_gens.Cogen[g]) #
    @constraint(BUCR1model, conCoGenTotGenZone[n=1:N_ZONES], sum((bucrm_genOut[g]*map_gens[g,n]*df_gens.Cogen[g]) for g=1:GENS) -bucr1_prepoc.cogen[n] ==0)

    # Constraints representing technical limits of conventional generators
    #Status transition trajectory of
    @constraint(BUCR1model, conStartUpAndDn[g=1:GENS], (bucrm_genOnOff[g] - bucr1_results.gens.onoff_init[g] - bucrm_genStartUp[g] + bucrm_genShutDown[g])==0)
    # Max Power generation limit in Block 1
    @constraint(BUCR1model, conMaxPowBlock1[g=1:GENS],  bucrm_genOut_Block[g,1] <= df_gens.IHRC_B1_Q[g]*bucrm_genOnOff[g] )
    # Max Power generation limit in Block 2
    @constraint(BUCR1model, conMaxPowBlock2[g=1:GENS],  bucrm_genOut_Block[g,2] <= df_gens.IHRC_B2_Q[g]*bucrm_genOnOff[g] )
    # Max Power generation limit in Block 3
    @constraint(BUCR1model, conMaxPowBlock3[g=1:GENS],  bucrm_genOut_Block[g,3] <= df_gens.IHRC_B3_Q[g]*bucrm_genOnOff[g] )
    # Max Power generation limit in Block 4
    @constraint(BUCR1model, conMaxPowBlock4[g=1:GENS],  bucrm_genOut_Block[g,4] <= df_gens.IHRC_B4_Q[g]*bucrm_genOnOff[g] )
    # Max Power generation limit in Block 5
    @constraint(BUCR1model, conMaxPowBlock5[g=1:GENS],  bucrm_genOut_Block[g,5] <= df_gens.IHRC_B5_Q[g]*bucrm_genOnOff[g] )
    # Max Power generation limit in Block 6
    @constraint(BUCR1model, conMaxPowBlock6[g=1:GENS],  bucrm_genOut_Block[g,6] <= df_gens.IHRC_B6_Q[g]*bucrm_genOnOff[g] )
    # Max Power generation limit in Block 7
    @constraint(BUCR1model, conMaxPowBlock7[g=1:GENS],  bucrm_genOut_Block[g,7] <= df_gens.IHRC_B7_Q[g]*bucrm_genOnOff[g] )
    # Total Production of each generation equals the sum of generation from its all blocks
    @constraint(BUCR1model, conTotalGen[g=1:GENS],  sum(bucrm_genOut_Block[g,b] for b=1:BLOCKS)+ bucrm_TotGenVioP[g]-bucrm_TotGenVioN[g]==bucrm_genOut[g])
    #Max power generation limit
    @constraint(BUCR1model, conMaxPow[g=1:GENS],  bucrm_genOut[g] - bucrm_MaxGenVioP[g] <= df_gens.MaxPowerOut[g]*bucrm_genOnOff[g])
    # Min power generation limit
    @constraint(BUCR1model, conMinPow[g=1:GENS],  bucrm_genOut[g] + bucrm_MinGenVioP[g] >= df_gens.MinPowerOut[g]*bucrm_genOnOff[g])
    #Up ramp rate limit
    @constraint(BUCR1model, conRampRateUp[g=1:GENS], (bucrm_genOut[g] - bucr1_results.gens.power_out_init[g] <=(df_gens.RampUpLimit[g]*bucr1_results.gens.onoff_init[g]) + (df_gens.RampStartUpLimit[g]*bucrm_genStartUp[g])))
    # Down ramp rate limit
    @constraint(BUCR1model, conRampRateDown[g=1:GENS], (bucr1_results.gens.power_out_init[g] - bucrm_genOut[g] <=(df_gens.RampDownLimit[g]*bucrm_genOnOff[g]) + (df_gens.RampShutDownLimit[g]*bucrm_genShutDown[g])))
    # Min Up Time limit with alternative formulation
    #The next two constraints enforce limits on binary commitment variables of slow and fast generators
    # scheduled slow units are forced to remain on, offline slow units remain off, and fast start units
    # could change their commitment dependent on their MUT and MDT
    @constraint(BUCR1model, conCommitmentUB[g=1:GENS], (bucrm_genOnOff[g] <= bucr_commit_UB[g]))
    # if the generator is slow start and scheduled "on" in the FUCR,  is fixed by the following constraint
    @constraint(BUCR1model, conCommitmentLB[g=1:GENS], (bucrm_genOnOff[g] >= bucr_commit_LB[g]))

    # Constraints representing technical limits of peakers
    #Status transition trajectory of
    @constraint(BUCR1model, conStartUpAndDn_Peaker[k=1:PEAKERS], (bucrm_peakerOnOff[k] - bucr1_results.peakers.onoff_init[k] - bucrm_peakerStartUp[k] + bucrm_peakerShutDown[k])==0)
    # Max Power generation limit in Block 1
    @constraint(BUCR1model, conMaxPowBlock1_Peaker[k=1:PEAKERS],  bucrm_peakerOut_Block[k,1] <= df_peakers.IHRC_B1_Q[k]*bucrm_peakerOnOff[k] )
    # Max Power generation limit in Block 2
    @constraint(BUCR1model, conMaxPowBlock2_Peaker[k=1:PEAKERS],  bucrm_peakerOut_Block[k,2] <= df_peakers.IHRC_B2_Q[k]*bucrm_peakerOnOff[k] )
    # Max Power generation limit in Block 3
    @constraint(BUCR1model, conMaxPowBlock3_Peaker[k=1:PEAKERS],  bucrm_peakerOut_Block[k,3] <= df_peakers.IHRC_B3_Q[k]*bucrm_peakerOnOff[k] )
    # Max Power generation limit in Block 4
    @constraint(BUCR1model, conMaxPowBlock4_Peaker[k=1:PEAKERS],  bucrm_peakerOut_Block[k,4] <= df_peakers.IHRC_B4_Q[k]*bucrm_peakerOnOff[k] )
    # Max Power generation limit in Block 5
    @constraint(BUCR1model, conMaxPowBlock5_Peaker[k=1:PEAKERS],  bucrm_peakerOut_Block[k,5] <= df_peakers.IHRC_B5_Q[k]*bucrm_peakerOnOff[k] )
    # Max Power generation limit in Block 6
    @constraint(BUCR1model, conMaxPowBlock6_Peaker[k=1:PEAKERS],  bucrm_peakerOut_Block[k,6] <= df_peakers.IHRC_B6_Q[k]*bucrm_peakerOnOff[k] )
    # Max Power generation limit in Block 7
    @constraint(BUCR1model, conMaxPowBlock7_Peaker[k=1:PEAKERS],  bucrm_peakerOut_Block[k,7] <= df_peakers.IHRC_B7_Q[k]*bucrm_peakerOnOff[k] )
    # Total Production of each generation equals the sum of generation from its all blocks
    @constraint(BUCR1model, conTotalGen_Peaker[k=1:PEAKERS], sum(bucrm_peakerOut_Block[k,b] for b=1:BLOCKS)>= bucrm_peakerOut[k])
    #Max power generation limit
    @constraint(BUCR1model, conMaxPow_Peaker[k=1:PEAKERS],  bucrm_peakerOut[k] <= df_peakers.MaxPowerOut[k]*bucrm_peakerOnOff[k])
    # Min power generation limit
    @constraint(BUCR1model, conMinPow_Peaker[k=1:PEAKERS],  bucrm_peakerOut[k] >= df_peakers.MinPowerOut[k]*bucrm_peakerOnOff[k])
    #Up ramp rate limit
    @constraint(BUCR1model, conRampRateUp_Peaker[k=1:PEAKERS], (bucrm_peakerOut[k] - bucr1_results.peakers.power_out_init[k] <=(df_peakers.RampUpLimit[k]*bucr1_results.peakers.onoff_init[k]) + (df_peakers.RampStartUpLimit[k]*bucrm_peakerStartUp[k])))
    # Down ramp rate limit
    @constraint(BUCR1model, conRampRateDown_Peaker[k=1:PEAKERS], (bucr1_results.peakers.power_out_init[k] - bucrm_peakerOut[k] <=(df_peakers.RampDownLimit[k]*bucrm_peakerOnOff[k]) + (df_peakers.RampShutDownLimit[k]*bucrm_peakerShutDown[k])))
    # Min Up Time limit with alternative formulation
    #The next two constraints enforce limits on binary commitment variables of slow and fast generators
    # scheduled slow units are forced to remain on, offline slow units remain off, and fast start units
    # could change their commitment dependent on their MUT and MDT
    @constraint(BUCR1model, conCommitmentUB_Peaker[k=1:PEAKERS], (bucrm_peakerOnOff[k] <= bucr_commit_peaker_UB[k]))
    # if the generator is slow start and scheduled "on" in the FUCR,  is fixed by the following constraint
    @constraint(BUCR1model, conCommitmentLB_Peaker[k=1:PEAKERS], (bucrm_peakerOnOff[k] >= bucr_commit_peaker_LB[k]))

    # Renewable generation constraints
    @constraint(BUCR1model, conSolarLimit[n=1:N_ZONES], bucrm_solarG[n] + bucrm_solarGSpil[n] <= bucr1_prepoc.solar[n])
    @constraint(BUCR1model, conWindLimit[n=1:N_ZONES], bucrm_windG[n] + bucrm_windGSpil[n] <= bucr1_prepoc.wind[n])
    @constraint(BUCR1model, conHydroLimit[n=1:N_ZONES], bucrm_hydroG[n] + bucrm_hydroGSpil[n] <= bucr1_prepoc.hydro[n])

    # Constraints representing technical characteristics of storage units
    # the next three constraints fix the balancing charging/discharging/Idle status to their optimal outcomes as determined by FUCR
    @constraint(BUCR1model, conStorgChrgStatusFixed[p=1:STORG_UNITS], (bucrm_storgChrg[p]==fucr_to_bucr1.storg.chrg[p,hour]))
    @constraint(BUCR1model, conStorgDisChrgStatusFixed[p=1:STORG_UNITS], (bucrm_storgDisc[p]==fucr_to_bucr1.storg.disc[p,hour]))
    @constraint(BUCR1model, conStorgIdleStatusFixed[p=1:STORG_UNITS], (bucrm_storgIdle[p]==fucr_to_bucr1.storg.idle[p,hour]))

    # charging power limit
    @constraint(BUCR1model, conStrgChargPowerLimit[p=1:STORG_UNITS], (bucrm_storgChrgPwr[p])<=df_storage.Power[p]*bucrm_storgChrg[p])
    # Discharging power limit
    @constraint(BUCR1model, conStrgDisChgPowerLimit[p=1:STORG_UNITS], (bucrm_storgDiscPwr[p])<=df_storage.Power[p]*bucrm_storgDisc[p])
    # State of charge at t
    @constraint(BUCR1model, conStorgSOC[p=1:STORG_UNITS], bucrm_storgSOC[p]==bucr1_results.storg.soc_init[p]-(bucrm_storgDiscPwr[p]/df_storage.TripEfficDown[p])+(bucrm_storgChrgPwr[p]*df_storage.TripEfficUp[p])-(bucrm_storgSOC[p]*df_storage.SelfDischarge[p]))
    # minimum energy limit
    @constraint(BUCR1model, conMinEnrgStorgLimi[p=1:STORG_UNITS], bucrm_storgSOC[p]>=0)
    # Maximum energy limit
    @constraint(BUCR1model, conMaxEnrgStorgLimi[p=1:STORG_UNITS], bucrm_storgSOC[p]<=(df_storage.Power[p]/df_storage.PowerToEnergRatio[p]))

    # Constraints representing transmission grid capacity constraints
    # DC Power Flow Calculation
    #@constraint(BUCR1model, conDCPowerFlowPos[n=1:N_ZONES, m=1:N_ZONES], bucrm_powerFlow[n,m]-(TranS[n,m]*(bucrm_voltAngle[n]-bucrm_voltAngle[m])) ==0)
    @constraint(BUCR1model, conDCPowerFlowNeg[n=1:N_ZONES, m=1:N_ZONES], bucrm_powerFlow[n,m]+bucrm_powerFlow[m,n]==0)
    # Tranmission flow bounds (from n to m and from m to n)
    @constraint(BUCR1model, conPosFlowLimit[n=1:N_ZONES, m=1:N_ZONES], bucrm_powerFlow[n,m]<=TranC[n,m])
    @constraint(BUCR1model, conNegFlowLimit[n=1:N_ZONES, m=1:N_ZONES], bucrm_powerFlow[n,m]>=-TranC[n,m])
    # Voltage Angle bounds and reference point
    #@constraint(BUCR1model, conVoltAnglUB[n=1:N_ZONES], bucrm_voltAngle[n]<=π)
    #@constraint(BUCR1model, conVoltAnglLB[n=1:N_ZONES], bucrm_voltAngle[n]>=-π)
    #@constraint(BUCR1model, conVoltAngRef, bucrm_voltAngle[1]==0)

    # Demand-side Constraints
    @constraint(BUCR1model, conDemandLimit[n=1:N_ZONES], bucrm_Demand[n]+ bucrm_Demand_Curt[n] == bucr1_prepoc.demand[n])

    # Demand Curtailment and wind generation limits
    @constraint(BUCR1model, conDemandCurtLimit[n=1:N_ZONES], bucrm_Demand_Curt[n] <= LOAD_SHED_MAX);
    @constraint(BUCR1model, conOverGenLimit[n=1:N_ZONES], bucrm_OverGen[n] <= OVERGEN_MAX);

    # System-wide Constraints
    #nodal balance constraint
    @constraint(BUCR1model, conNodBalanc[n=1:N_ZONES],
        sum((bucrm_genOut[g]*map_gens[g,n]) for g=1:GENS) +
        sum((bucrm_peakerOut[k]*map_peakers[k,n]) for k=1:PEAKERS) +
        sum((bucrm_storgDiscPwr[p]*map_storage[p,n]) for p=1:STORG_UNITS) -
        sum((bucrm_storgChrgPwr[p]*map_storage[p,n]) for p=1:STORG_UNITS) + bucrm_solarG[n] + bucrm_windG[n] +bucrm_hydroG[n] - bucrm_OverGen[n]- bucrm_Demand[n] == sum(bucrm_powerFlow[n,m] for m=1:M_ZONES))
    # Minimum up reserve requirement
    #    @constraint(BUCR1model, conMinUpReserveReq[t=1:N_Hrs_BUCR], sum(genResUp[g,t] for g=1:GENS) + sum(storgResUp[p,t] for p=1:STORG_UNITS) >= Reserve_Req_Up[t] )

    # Minimum down reserve requirement
    #    @constraint(BUCR1model, conMinDnReserveReq[t=1:N_Hrs_BUCR], sum(genResDn[g,t] for g=1:GENS) + sum(storgResDn[p,t] for p=1:STORG_UNITS) >= Reserve_Req_Dn[t] )

    t2_BUCR1model = time_ns()
    time_BUCR1model = (t2_BUCR1model -t1_BUCR1model)/1.0e9;
    @info "BUCR1model for day: $day, hour $(hour+INIT_HR_FUCR) setup executed in (s):  $time_BUCR1model";

    open(".//outputs//time_performance.csv", FILE_ACCESS_APPEND) do io
    		writedlm(io, hcat("BUCR1model", time_BUCR1model, "day: $day",
    				"hour $(hour+INIT_HR_FUCR)", "", "Model Setup"), ',')
    end; # closes file

    # solve the First WAUC model (BUCR)
    JuMP.optimize!(BUCR1model)

    # Pricing general results in the terminal window
    println("Objective value: ", JuMP.objective_value(BUCR1model))
    open(".//outputs//objective_values_v76.csv", FILE_ACCESS_APPEND) do io
    		writedlm(io, hcat("BUCR1model", "day: $day",
    				"hour $(hour+INIT_HR_FUCR)", "", JuMP.objective_value(BUCR1model)), ',')
    end;

    println("------------------------------------")
    println("------- BAUC1 OBJECTIVE VALUE -------")
    println("Objective value for day ", day, " and hour ", hour+INIT_HR_FUCR," is: ", JuMP.objective_value(BUCR1model))
    println("------------------------------------")
    println("------- BAUC1 PRIMAL STATUS -------")
    println(primal_status(BUCR1model))
    println("------------------------------------")
    println("------- BAUC1 DUAL STATUS -------")
    println(JuMP.dual_status(BUCR1model))
    println("Day: ", day, " and hour ", hour+INIT_HR_FUCR, ": solved")
    println("---------------------------")
    println("BUCR1model Number of variables: ", JuMP.num_variables(BUCR1model))
    @info "BUCR1model Number of variables: " JuMP.num_variables(BUCR1model)
    open(".//outputs//time_performance.csv", FILE_ACCESS_APPEND) do io
    		writedlm(io, hcat("BUCR1model", JuMP.num_variables(BUCR1model), "day: $day",
    				"hour $(hour+INIT_HR_FUCR)", "", "Variables"), ',')
    end;

    @debug "BUCR1model for day: $day, hour $(hour+INIT_HR_FUCR) optimized executed in (s): $(solve_time(BUCR1model))";

    t1_write_BUCR1model_results = time_ns()
    open(".//outputs//time_performance.csv", FILE_ACCESS_APPEND) do io
    		writedlm(io, hcat("BUCR1model", solve_time(BUCR1model), "day: $day",
    				"hour $(hour+INIT_HR_FUCR)", "", "Model Optimization"), ',')
    end; # closes file

    ## Write the optimal outcomes into spreadsheets
    #TODO: Later we need to include a variable for day so the cell number in
    # which the results are printed is updated accordingly
    # Write the conventional generators' schedules
    open(".//outputs//BUCR_GenOutputs.csv", FILE_ACCESS_APPEND) do io
    	for g=1:GENS
    		writedlm(io, hcat(day, hour+INIT_HR_FUCR, g, df_gens.UNIT_ID[g],
    			df_gens.MinPowerOut[g], df_gens.MaxPowerOut[g],
    			JuMP.value.(bucrm_genOut[g]), JuMP.value.(bucrm_genOnOff[g]),
    			JuMP.value.(bucrm_genShutDown[g]),
    			JuMP.value.(bucrm_genStartUp[g]),
    			JuMP.value.(bucrm_TotGenVioP[g]),
    			JuMP.value.(bucrm_TotGenVioN[g]), JuMP.value.(bucrm_MaxGenVioP[g]),
    			JuMP.value.(bucrm_MinGenVioP[g]) ), ',')
    	end # ends the loop
    end; # closes file

    # Write the conventional peakers' schedules
    open(".//outputs//BUCR_PeakerOutputs.csv", FILE_ACCESS_APPEND) do io
    	for k=1:PEAKERS
    		writedlm(io, hcat(day, hour+INIT_HR_FUCR, k, df_peakers.UNIT_ID[k],
    			df_peakers.MinPowerOut[k], df_peakers.MaxPowerOut[k],
    			JuMP.value.(bucrm_peakerOut[k]),
    			JuMP.value.(bucrm_peakerOnOff[k]),
    			JuMP.value.(bucrm_peakerShutDown[k]),
    			JuMP.value.(bucrm_peakerStartUp[k]) ), ',')
    	end # ends the loop
    end; # closes file

    # Writing storage units' optimal schedules in CSV file
    open(".//outputs//BUCR_StorageOutputs.csv", FILE_ACCESS_APPEND) do io
    	for p=1:STORG_UNITS
    		writedlm(io, hcat(day, hour+INIT_HR_FUCR, p, df_storage.Name[p],
    			df_storage.Power[p],
    			df_storage.Power[p]/df_storage.PowerToEnergRatio[p],
    			JuMP.value.(bucrm_storgChrg[p]),
    			JuMP.value.(bucrm_storgDisc[p]),
    			JuMP.value.(bucrm_storgIdle[p]),
    			JuMP.value.(bucrm_storgChrgPwr[p]),
    			JuMP.value.(bucrm_storgDiscPwr[p]),
    			JuMP.value.(bucrm_storgSOC[p]) ), ',')
    	end # ends the loop
    end; # closes file

    # Writing the transmission flow schedules in CSV file
    open(".//outputs//BUCR_TranFlowOutputs.csv", FILE_ACCESS_APPEND) do io
    	for n=1:N_ZONES, m=1:M_ZONES
    		writedlm(io, hcat(day, hour+INIT_HR_FUCR, n, m,
    			JuMP.value.(bucrm_powerFlow[n,m]), TranC[n,m] ), ',')
    	end # ends the loop
    end; # closes file

    # Writing the curtilment, overgeneration, and spillage outcomes in CSV file
    open(".//outputs//BUCR_Curtail.csv", FILE_ACCESS_APPEND) do io
    	for n=1:N_ZONES
    	   writedlm(io, hcat(day, hour+INIT_HR_FUCR, n,
    		   JuMP.value.(bucrm_OverGen[n]), JuMP.value.(bucrm_Demand_Curt[n]),
    		   JuMP.value.(bucrm_windGSpil[n]), JuMP.value.(bucrm_solarGSpil[n]),
    		   JuMP.value.(bucrm_hydroGSpil[n]) ), ',')
    	end # ends the loop
    end; # closes file

    t2_write_BUCR1model_results = time_ns()
    time_write_BUCR1model_results = (t2_write_BUCR1model_results -t1_write_BUCR1model_results)/1.0e9;
    @info "Write BUCR1model results for day $day and hour $(hour+INIT_HR_FUCR) executed in (s): $time_write_BUCR1model_results";

    open(".//outputs//time_performance.csv", FILE_ACCESS_APPEND) do io
    		writedlm(io, hcat("BUCR1model", time_write_BUCR1model_results, "day: $day",
    				"hour $(hour+INIT_HR_FUCR)", "", "Write CSV files"), ',')
    end; #closes file
    ## Initilization of the next UC Run
    # This must be updated later when we run two WAUCs every day and then RTUCs
    # Setting initial values for bucr1_results (next hour), SUCR1, and bucr2_results
    t1_bucrm_init_next_UC = time_ns();

    for g=1:GENS
    	global bucr1_results.gens.onoff_init[g] = JuMP.value.(bucrm_genOnOff[g]);
    	global bucr1_results.gens.power_out_init[g] = JuMP.value.(bucrm_genOut[g]);
    	if hour==(INIT_HR_SUCR-INIT_HR_FUCR)
    		bucr2_results.gens.onoff_init[g] = JuMP.value.(bucrm_genOnOff[g]);
    		bucr2_results.gens.power_out_init[g] = JuMP.value.(bucrm_genOut[g]);
    		sucr_results.gens.onoff_init[g] = JuMP.value.(bucrm_genOnOff[g]);
    		sucr_results.gens.power_out_init[g] = JuMP.value.(bucrm_genOut[g]);
    	end
    end

    for k=1:PEAKERS
    	bucr1_results.peakers.onoff_init[k] = round(JuMP.value.(bucrm_peakerOnOff[k]));
    	bucr1_results.peakers.power_out_init[k] = JuMP.value.(bucrm_peakerOut[k]);
    	if hour==(INIT_HR_SUCR-INIT_HR_FUCR)
    		bucr2_results.peakers.onoff_init[k] = round(JuMP.value.(bucrm_peakerOnOff[k]));
    		bucr2_results.peakers.power_out_init[k] = JuMP.value.(bucrm_peakerOut[k]);
    		sucr_results.peakers.onoff_init[k] = round(JuMP.value.(bucrm_peakerOnOff[k]));
    		sucr_results.peakers.power_out_init[k] = JuMP.value.(bucrm_peakerOut[k]);
    	end
    end

    for p=1:STORG_UNITS
    	bucr1_results.storg.soc_init[p]=JuMP.value.(bucrm_storgSOC[p]);
    	if hour==(INIT_HR_SUCR-INIT_HR_FUCR)
    		bucr2_results.storg.soc_init[p] = JuMP.value.(bucrm_storgSOC[p]);
    		sucr_results.storg.soc_init[p] = JuMP.value.(bucrm_storgSOC[p]);
    	end
    end

	#NOTE: No need to use round function (round(JuMP.value.(BUCR1_genStartUp[g])))==1
    for g=1:GENS
    	if JuMP.value.(bucrm_genStartUp[g]) > 0.5
    		bucr1_results.gens.uptime_init[g]= 1;
    		bucr1_results.gens.dntime_init[g] = 0;
    	elseif JuMP.value.(bucrm_genShutDown[g]) > 0.5
    		bucr1_results.gens.uptime_init[g]= 0;
    		bucr1_results.gens.dntime_init[g]= 1;
    	else
    		if JuMP.value.(bucrm_genOnOff[g]) > 0.5
    			bucr1_results.gens.uptime_init[g]= bucr1_results.gens.uptime_init[g]+1;
    			bucr1_results.gens.dntime_init[g]= 0;
    		else
    			bucr1_results.gens.uptime_init[g]= 0;
    			bucr1_results.gens.dntime_init[g]= bucr1_results.gens.dntime_init[g]+1;
    		end
    	end
    	if hour == (INIT_HR_SUCR-INIT_HR_FUCR)
    		bucr2_results.gens.uptime_init[g]= bucr1_results.gens.uptime_init[g];
    		bucr2_results.gens.dntime_init[g]= bucr1_results.gens.dntime_init[g];
    		sucr_results.gens.uptime_init[g]= bucr1_results.gens.uptime_init[g];
    		sucr_results.gens.dntime_init[g]= bucr1_results.gens.dntime_init[g];
    	end
    end
    #TODO:No need to use the round function.  Since 0<= bucrm_peakerStartUp<=1,
    # bucrm_genStartUp>0 in the If Statement is enough
    for k=1:PEAKERS
    	if (JuMP.value.(bucrm_peakerStartUp[k]))>0
    		bucr1_results.peakers.uptime_init[k]= 1;
    		bucr1_results.peakers.dntime_init[k] = 0;
    	elseif (JuMP.value.(bucrm_peakerShutDown[k]))>0
    		bucr1_results.peakers.uptime_init[k]= 0;
    		bucr1_results.peakers.dntime_init[k]= 1;
    	else
    		if (JuMP.value.(bucrm_peakerOnOff[k]))>0
    			bucr1_results.peakers.uptime_init[k]= bucr1_results.peakers.uptime_init[k]+1;
    			bucr1_results.peakers.dntime_init[k]= 0;
    		else
    			bucr1_results.peakers.uptime_init[k]= 0;
    			bucr1_results.peakers.dntime_init[k]= bucr1_results.peakers.dntime_init[k]+1;
    		end
    	end
    	if hour==(INIT_HR_SUCR-INIT_HR_FUCR)
    		bucr2_results.peakers.uptime_init[k]= bucr1_results.peakers.uptime_init[k];
    		bucr2_results.peakers.dntime_init[k]= bucr1_results.peakers.dntime_init[k];
    		sucr_results.peakers.uptime_init[k]= bucr1_results.peakers.uptime_init[k];
    		sucr_results.peakers.dntime_init[k]= bucr1_results.peakers.dntime_init[k];
    	end
    end

    t2_bucrm_init_next_UC = time_ns();

    time_bucrm_init_next_UC = (t2_bucrm_init_next_UC -t1_bucrm_init_next_UC)/1.0e9;
    @info "bucrm_init_next_UC data handling for day $day and hour $(hour+INIT_HR_FUCR) executed in (s): $time_bucrm_init_next_UC";

    open(".//outputs//time_performance.csv", FILE_ACCESS_APPEND) do io
    		writedlm(io, hcat("BUCR1model", time_bucrm_init_next_UC, "day: $day",
    				"hour $(hour+INIT_HR_FUCR)", "bucrm_init_next_UC", "Data Manipulation"), ',')
    end; #closes file

    return bucr1_results, sucr_results, bucr2_results
end; # function

max_stars_repo_name: power-systems-op/grace-arpae

max_issues_repo_name: power-systems-op/grace-arpae

max_forks_repo_name: power-systems-op/grace-arpae